1 /* 2 Code for timestepping with BDF methods 3 */ 4 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 5 6 static PetscBool cited = PETSC_FALSE; 7 static const char citation[] = 8 "@book{Brenan1995,\n" 9 " title = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n" 10 " author = {Brenan, K. and Campbell, S. and Petzold, L.},\n" 11 " publisher = {Society for Industrial and Applied Mathematics},\n" 12 " year = {1995},\n" 13 " doi = {10.1137/1.9781611971224},\n}\n"; 14 15 typedef struct { 16 PetscInt k,n; 17 PetscReal time[6+2]; 18 Vec work[6+2]; 19 PetscReal shift; 20 Vec vec_dot; 21 Vec vec_lte; 22 23 PetscInt order; 24 TSStepStatus status; 25 } TS_BDF; 26 27 28 PETSC_STATIC_INLINE void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[]) 29 { 30 PetscInt k,j; 31 for (k=0; k<n; k++) 32 for (L[k]=1, j=0; j<n; j++) 33 if (j != k) 34 L[k] *= (t - T[j])/(T[k] - T[j]); 35 } 36 37 PETSC_STATIC_INLINE void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[]) 38 { 39 PetscInt k,j,i; 40 for (k=0; k<n; k++) 41 for (dL[k]=0, j=0; j<n; j++) 42 if (j != k) { 43 PetscReal L = 1/(T[k] - T[j]); 44 for (i=0; i<n; i++) 45 if (i != j && i != k) 46 L *= (t - T[i])/(T[k] - T[i]); 47 dL[k] += L; 48 } 49 } 50 51 static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X) 52 { 53 TS_BDF *bdf = (TS_BDF*)ts->data; 54 PetscInt i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec)); 55 Vec tail = bdf->work[n-1]; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 for (i=n-1; i>=2; i--) { 60 bdf->time[i] = bdf->time[i-1]; 61 bdf->work[i] = bdf->work[i-1]; 62 } 63 bdf->n = PetscMin(bdf->n+1,n-1); 64 bdf->time[1] = t; 65 bdf->work[1] = tail; 66 ierr = VecCopy(X,tail);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 static PetscErrorCode TSBDF_VecDot(TS ts,PetscInt order,PetscReal t,Vec X,Vec Xdot,PetscReal *shift) 71 { 72 TS_BDF *bdf = (TS_BDF*)ts->data; 73 PetscInt i,n = order+1; 74 PetscReal time[7]; 75 Vec vecs[7]; 76 PetscScalar alpha[7]; 77 PetscErrorCode ierr; 78 79 PetscFunctionBegin; 80 n = PetscMax(n,2); 81 time[0] = t; for (i=1; i<n; i++) time[i] = bdf->time[i]; 82 vecs[0] = X; for (i=1; i<n; i++) vecs[i] = bdf->work[i]; 83 LagrangeBasisDers(n,t,time,alpha); 84 ierr = VecZeroEntries(Xdot);CHKERRQ(ierr); 85 ierr = VecMAXPY(Xdot,n,alpha,vecs);CHKERRQ(ierr); 86 if (shift) *shift = PetscRealPart(alpha[0]); 87 PetscFunctionReturn(0); 88 } 89 90 static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte) 91 { 92 TS_BDF *bdf = (TS_BDF*)ts->data; 93 PetscInt i,n = order+1; 94 PetscReal *time = bdf->time; 95 Vec *vecs = bdf->work; 96 PetscScalar a[8],b[8],alpha[8]; 97 PetscErrorCode ierr; 98 99 PetscFunctionBegin; 100 LagrangeBasisDers(n+0,time[0],time,a); a[n] =0; 101 LagrangeBasisDers(n+1,time[0],time,b); 102 for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0]; 103 ierr = VecZeroEntries(lte);CHKERRQ(ierr); 104 ierr = VecMAXPY(lte,n+1,alpha,vecs);CHKERRQ(ierr); 105 PetscFunctionReturn(0); 106 } 107 108 static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X) 109 { 110 TS_BDF *bdf = (TS_BDF*)ts->data; 111 PetscInt n = order+1; 112 PetscReal *time = bdf->time+1; 113 Vec *vecs = bdf->work+1; 114 PetscScalar alpha[7]; 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 n = PetscMin(n,bdf->n); 119 LagrangeBasisVals(n,t,time,alpha); 120 ierr = VecZeroEntries(X);CHKERRQ(ierr); 121 ierr = VecMAXPY(X,n,alpha,vecs);CHKERRQ(ierr); 122 PetscFunctionReturn(0); 123 } 124 125 static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X) 126 { 127 TS_BDF *bdf = (TS_BDF*)ts->data; 128 PetscInt n = order+1; 129 PetscReal *time = bdf->time; 130 Vec *vecs = bdf->work; 131 PetscScalar alpha[7]; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 LagrangeBasisVals(n,t,time,alpha); 136 ierr = VecZeroEntries(X);CHKERRQ(ierr); 137 ierr = VecMAXPY(X,n,alpha,vecs);CHKERRQ(ierr); 138 PetscFunctionReturn(0); 139 } 140 141 static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x) 142 { 143 PetscInt nits,lits; 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 ierr = SNESSolve(ts->snes,b,x);CHKERRQ(ierr); 148 ierr = SNESGetIterationNumber(ts->snes,&nits);CHKERRQ(ierr); 149 ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 150 ts->snes_its += nits; ts->ksp_its += lits; 151 PetscFunctionReturn(0); 152 } 153 154 static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept) 155 { 156 TS_BDF *bdf = (TS_BDF*)ts->data; 157 PetscErrorCode ierr; 158 159 PetscFunctionBegin; 160 bdf->k = 1; bdf->n = 0; 161 ierr = TSBDF_Advance(ts,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 162 163 bdf->time[0] = ts->ptime + ts->time_step/2; 164 ierr = VecCopy(bdf->work[1],bdf->work[0]);CHKERRQ(ierr); 165 ierr = TSPreStage(ts,bdf->time[0]);CHKERRQ(ierr); 166 ierr = TS_SNESSolve(ts,NULL,bdf->work[0]);CHKERRQ(ierr); 167 ierr = TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);CHKERRQ(ierr); 168 ierr = TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);CHKERRQ(ierr); 169 if (!*accept) PetscFunctionReturn(0); 170 171 bdf->k = PetscMin(2,bdf->order); bdf->n++; 172 ierr = VecCopy(bdf->work[0],bdf->work[2]);CHKERRQ(ierr); 173 bdf->time[2] = bdf->time[0]; 174 PetscFunctionReturn(0); 175 } 176 177 static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"}; 178 179 static PetscErrorCode TSStep_BDF(TS ts) 180 { 181 TS_BDF *bdf = (TS_BDF*)ts->data; 182 PetscInt rejections = 0; 183 PetscBool stageok,accept = PETSC_TRUE; 184 PetscReal next_time_step = ts->time_step; 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 189 190 if (!ts->steprollback && !ts->steprestart) { 191 bdf->k = PetscMin(bdf->k+1,bdf->order); 192 ierr = TSBDF_Advance(ts,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 193 } 194 195 bdf->status = TS_STEP_INCOMPLETE; 196 while (!ts->reason && bdf->status != TS_STEP_COMPLETE) { 197 198 if (ts->steprestart) { 199 ierr = TSBDF_Restart(ts,&stageok);CHKERRQ(ierr); 200 if (!stageok) goto reject_step; 201 } 202 203 bdf->time[0] = ts->ptime + ts->time_step; 204 ierr = TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);CHKERRQ(ierr); 205 ierr = TSPreStage(ts,bdf->time[0]);CHKERRQ(ierr); 206 ierr = TS_SNESSolve(ts,NULL,bdf->work[0]);CHKERRQ(ierr); 207 ierr = TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);CHKERRQ(ierr); 208 ierr = TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);CHKERRQ(ierr); 209 if (!stageok) goto reject_step; 210 211 bdf->status = TS_STEP_PENDING; 212 ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr); 213 ierr = TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);CHKERRQ(ierr); 214 ierr = TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 215 bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 216 if (!accept) { ts->time_step = next_time_step; goto reject_step; } 217 218 ierr = VecCopy(bdf->work[0],ts->vec_sol);CHKERRQ(ierr); 219 ts->ptime += ts->time_step; 220 ts->time_step = next_time_step; 221 break; 222 223 reject_step: 224 ts->reject++; accept = PETSC_FALSE; 225 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 226 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 227 ts->reason = TS_DIVERGED_STEP_REJECTED; 228 } 229 } 230 PetscFunctionReturn(0); 231 } 232 233 static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X) 234 { 235 TS_BDF *bdf = (TS_BDF*)ts->data; 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 ierr = TSBDF_Interpolate(ts,bdf->k,t,X);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 243 static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte) 244 { 245 TS_BDF *bdf = (TS_BDF*)ts->data; 246 PetscInt k = bdf->k; 247 PetscReal wltea,wlter; 248 Vec X = bdf->work[0], Y = bdf->vec_lte; 249 PetscErrorCode ierr; 250 251 PetscFunctionBegin; 252 k = PetscMin(k,bdf->n-1); 253 ierr = TSBDF_VecLTE(ts,k,Y);CHKERRQ(ierr); 254 ierr = VecAXPY(Y,1,X);CHKERRQ(ierr); 255 ierr = TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);CHKERRQ(ierr); 256 if (order) *order = k + 1; 257 PetscFunctionReturn(0); 258 } 259 260 static PetscErrorCode TSRollBack_BDF(TS ts) 261 { 262 TS_BDF *bdf = (TS_BDF*)ts->data; 263 PetscErrorCode ierr; 264 265 PetscFunctionBegin; 266 ierr = VecCopy(bdf->work[1],ts->vec_sol);CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 static PetscErrorCode SNESTSFormFunction_BDF(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts) 271 { 272 TS_BDF *bdf = (TS_BDF*)ts->data; 273 PetscInt k = bdf->k; 274 PetscReal t = bdf->time[0]; 275 Vec V = bdf->vec_dot; 276 PetscErrorCode ierr; 277 278 PetscFunctionBegin; 279 ierr = TSBDF_VecDot(ts,k,t,X,V,&bdf->shift);CHKERRQ(ierr); 280 /* F = Function(t,X,V) */ 281 ierr = TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);CHKERRQ(ierr); 282 PetscFunctionReturn(0); 283 } 284 285 static PetscErrorCode SNESTSFormJacobian_BDF(PETSC_UNUSED SNES snes, 286 PETSC_UNUSED Vec X, 287 Mat J,Mat P, 288 TS ts) 289 { 290 TS_BDF *bdf = (TS_BDF*)ts->data; 291 PetscReal t = bdf->time[0]; 292 Vec V = bdf->vec_dot; 293 PetscReal dVdX = bdf->shift; 294 PetscErrorCode ierr; 295 296 PetscFunctionBegin; 297 /* J,P = Jacobian(t,X,V) */ 298 ierr = TSComputeIJacobian(ts,t,X,V,dVdX,J,P,PETSC_FALSE);CHKERRQ(ierr); 299 PetscFunctionReturn(0); 300 } 301 302 static PetscErrorCode TSReset_BDF(TS ts) 303 { 304 TS_BDF *bdf = (TS_BDF*)ts->data; 305 size_t i,n = sizeof(bdf->work)/sizeof(Vec); 306 PetscErrorCode ierr; 307 308 PetscFunctionBegin; 309 for (i=0; i<n; i++) {ierr = VecDestroy(&bdf->work[i]);CHKERRQ(ierr);} 310 ierr = VecDestroy(&bdf->vec_dot);CHKERRQ(ierr); 311 ierr = VecDestroy(&bdf->vec_lte);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 static PetscErrorCode TSDestroy_BDF(TS ts) 316 { 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 ierr = TSReset_BDF(ts);CHKERRQ(ierr); 321 ierr = PetscFree(ts->data);CHKERRQ(ierr); 322 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);CHKERRQ(ierr); 323 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);CHKERRQ(ierr); 324 PetscFunctionReturn(0); 325 } 326 327 static PetscErrorCode TSSetUp_BDF(TS ts) 328 { 329 TS_BDF *bdf = (TS_BDF*)ts->data; 330 size_t i,n = sizeof(bdf->work)/sizeof(Vec); 331 PetscErrorCode ierr; 332 333 PetscFunctionBegin; 334 bdf->k = bdf->n = 0; 335 for (i=0; i<n; i++) {ierr = VecDuplicate(ts->vec_sol,&bdf->work[i]);CHKERRQ(ierr);} 336 ierr = VecDuplicate(ts->vec_sol,&bdf->vec_dot);CHKERRQ(ierr); 337 ierr = VecDuplicate(ts->vec_sol,&bdf->vec_lte);CHKERRQ(ierr); 338 339 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 340 ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr); 341 { 342 PetscBool isnone; 343 PetscReal low,high; 344 ierr = TSAdaptBasicGetClip(ts->adapt,&low,&high);CHKERRQ(ierr); 345 high = PetscMin(high,2.0); 346 ierr = TSAdaptBasicSetClip(ts->adapt,low,high);CHKERRQ(ierr); 347 ierr = PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&isnone);CHKERRQ(ierr); 348 if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) 349 ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP; 350 } 351 352 ierr = TSGetSNES(ts,&ts->snes);CHKERRQ(ierr); 353 PetscFunctionReturn(0); 354 } 355 356 static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts) 357 { 358 TS_BDF *bdf = (TS_BDF*)ts->data; 359 PetscErrorCode ierr; 360 361 PetscFunctionBegin; 362 ierr = PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");CHKERRQ(ierr); 363 { 364 PetscBool flg; 365 PetscInt order = bdf->order; 366 ierr = PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);CHKERRQ(ierr); 367 if (flg) {ierr = TSBDFSetOrder(ts,order);CHKERRQ(ierr);} 368 } 369 ierr = PetscOptionsTail();CHKERRQ(ierr); 370 PetscFunctionReturn(0); 371 } 372 373 static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer) 374 { 375 TS_BDF *bdf = (TS_BDF*)ts->data; 376 PetscBool iascii; 377 PetscErrorCode ierr; 378 379 PetscFunctionBegin; 380 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 381 if (iascii) { 382 ierr = PetscViewerASCIIPrintf(viewer," Order=%D\n",bdf->order);CHKERRQ(ierr); 383 } 384 if (ts->adapt) {ierr = TSAdaptView(ts->adapt,viewer);CHKERRQ(ierr);} 385 if (ts->snes) {ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr);} 386 PetscFunctionReturn(0); 387 } 388 389 /* ------------------------------------------------------------ */ 390 391 static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order) 392 { 393 TS_BDF *bdf = (TS_BDF*)ts->data; 394 395 PetscFunctionBegin; 396 if (order == bdf->order) PetscFunctionReturn(0); 397 if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order); 398 bdf->order = order; 399 PetscFunctionReturn(0); 400 } 401 402 static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order) 403 { 404 TS_BDF *bdf = (TS_BDF*)ts->data; 405 406 PetscFunctionBegin; 407 *order = bdf->order; 408 PetscFunctionReturn(0); 409 } 410 411 /* ------------------------------------------------------------ */ 412 413 /*MC 414 TSBDF - DAE solver using BDF methods 415 416 Level: beginner 417 418 .seealso: TS, TSCreate(), TSSetType() 419 M*/ 420 PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts) 421 { 422 TS_BDF *bdf; 423 PetscErrorCode ierr; 424 425 PetscFunctionBegin; 426 ts->ops->reset = TSReset_BDF; 427 ts->ops->destroy = TSDestroy_BDF; 428 ts->ops->view = TSView_BDF; 429 ts->ops->setup = TSSetUp_BDF; 430 ts->ops->setfromoptions = TSSetFromOptions_BDF; 431 ts->ops->step = TSStep_BDF; 432 ts->ops->evaluatewlte = TSEvaluateWLTE_BDF; 433 ts->ops->rollback = TSRollBack_BDF; 434 ts->ops->interpolate = TSInterpolate_BDF; 435 ts->ops->snesfunction = SNESTSFormFunction_BDF; 436 ts->ops->snesjacobian = SNESTSFormJacobian_BDF; 437 438 ierr = PetscNewLog(ts,&bdf);CHKERRQ(ierr); 439 ts->data = (void*)bdf; 440 441 bdf->order = 2; 442 bdf->status = TS_STEP_COMPLETE; 443 444 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);CHKERRQ(ierr); 445 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);CHKERRQ(ierr); 446 PetscFunctionReturn(0); 447 } 448 449 /* ------------------------------------------------------------ */ 450 451 /*@ 452 TSBDFSetOrder - Set the order of the BDF method 453 454 Logically Collective on TS 455 456 Input Parameter: 457 + ts - timestepping context 458 - order - order of the method 459 460 Options Database: 461 . -ts_bdf_order <order> 462 463 Level: intermediate 464 465 @*/ 466 PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order) 467 { 468 PetscErrorCode ierr; 469 470 PetscFunctionBegin; 471 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 472 PetscValidLogicalCollectiveInt(ts,order,2); 473 ierr = PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));CHKERRQ(ierr); 474 PetscFunctionReturn(0); 475 } 476 477 /*@ 478 TSBDFGetOrder - Get the order of the BDF method 479 480 Not Collective 481 482 Input Parameter: 483 . ts - timestepping context 484 485 Output Parameter: 486 . order - order of the method 487 488 Level: intermediate 489 490 @*/ 491 PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order) 492 { 493 PetscErrorCode ierr; 494 495 PetscFunctionBegin; 496 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 497 PetscValidIntPointer(order,2); 498 ierr = PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));CHKERRQ(ierr); 499 PetscFunctionReturn(0); 500 } 501