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 TSBDF_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 = TSBDF_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 = TSBDF_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 PetscReal low,high,two = 2; 332 PetscErrorCode ierr; 333 334 PetscFunctionBegin; 335 bdf->k = bdf->n = 0; 336 for (i=0; i<n; i++) {ierr = VecDuplicate(ts->vec_sol,&bdf->work[i]);CHKERRQ(ierr);} 337 ierr = VecDuplicate(ts->vec_sol,&bdf->vec_dot);CHKERRQ(ierr); 338 ierr = VecDuplicate(ts->vec_sol,&bdf->vec_lte);CHKERRQ(ierr); 339 340 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 341 ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr); 342 ierr = TSAdaptGetClip(ts->adapt,&low,&high);CHKERRQ(ierr); 343 ierr = TSAdaptSetClip(ts->adapt,low,PetscMin(high,two));CHKERRQ(ierr); 344 345 ierr = TSGetSNES(ts,&ts->snes);CHKERRQ(ierr); 346 PetscFunctionReturn(0); 347 } 348 349 static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts) 350 { 351 TS_BDF *bdf = (TS_BDF*)ts->data; 352 PetscErrorCode ierr; 353 354 PetscFunctionBegin; 355 ierr = PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");CHKERRQ(ierr); 356 { 357 PetscBool flg; 358 PetscInt order = bdf->order; 359 ierr = PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);CHKERRQ(ierr); 360 if (flg) {ierr = TSBDFSetOrder(ts,order);CHKERRQ(ierr);} 361 } 362 ierr = PetscOptionsTail();CHKERRQ(ierr); 363 PetscFunctionReturn(0); 364 } 365 366 static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer) 367 { 368 TS_BDF *bdf = (TS_BDF*)ts->data; 369 PetscBool iascii; 370 PetscErrorCode ierr; 371 372 PetscFunctionBegin; 373 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 374 if (iascii) { 375 ierr = PetscViewerASCIIPrintf(viewer," Order=%D\n",bdf->order);CHKERRQ(ierr); 376 } 377 PetscFunctionReturn(0); 378 } 379 380 /* ------------------------------------------------------------ */ 381 382 static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order) 383 { 384 TS_BDF *bdf = (TS_BDF*)ts->data; 385 386 PetscFunctionBegin; 387 if (order == bdf->order) PetscFunctionReturn(0); 388 if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order); 389 bdf->order = order; 390 PetscFunctionReturn(0); 391 } 392 393 static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order) 394 { 395 TS_BDF *bdf = (TS_BDF*)ts->data; 396 397 PetscFunctionBegin; 398 *order = bdf->order; 399 PetscFunctionReturn(0); 400 } 401 402 /* ------------------------------------------------------------ */ 403 404 /*MC 405 TSBDF - DAE solver using BDF methods 406 407 Level: beginner 408 409 .seealso: TS, TSCreate(), TSSetType() 410 M*/ 411 PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts) 412 { 413 TS_BDF *bdf; 414 PetscErrorCode ierr; 415 416 PetscFunctionBegin; 417 ts->ops->reset = TSReset_BDF; 418 ts->ops->destroy = TSDestroy_BDF; 419 ts->ops->view = TSView_BDF; 420 ts->ops->setup = TSSetUp_BDF; 421 ts->ops->setfromoptions = TSSetFromOptions_BDF; 422 ts->ops->step = TSStep_BDF; 423 ts->ops->evaluatewlte = TSEvaluateWLTE_BDF; 424 ts->ops->rollback = TSRollBack_BDF; 425 ts->ops->interpolate = TSInterpolate_BDF; 426 ts->ops->snesfunction = SNESTSFormFunction_BDF; 427 ts->ops->snesjacobian = SNESTSFormJacobian_BDF; 428 ts->default_adapt_type = TSADAPTBASIC; 429 430 ts->usessnes = PETSC_TRUE; 431 432 ierr = PetscNewLog(ts,&bdf);CHKERRQ(ierr); 433 ts->data = (void*)bdf; 434 435 bdf->order = 2; 436 bdf->status = TS_STEP_COMPLETE; 437 438 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);CHKERRQ(ierr); 439 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);CHKERRQ(ierr); 440 PetscFunctionReturn(0); 441 } 442 443 /* ------------------------------------------------------------ */ 444 445 /*@ 446 TSBDFSetOrder - Set the order of the BDF method 447 448 Logically Collective on TS 449 450 Input Parameter: 451 + ts - timestepping context 452 - order - order of the method 453 454 Options Database: 455 . -ts_bdf_order <order> 456 457 Level: intermediate 458 459 @*/ 460 PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order) 461 { 462 PetscErrorCode ierr; 463 464 PetscFunctionBegin; 465 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 466 PetscValidLogicalCollectiveInt(ts,order,2); 467 ierr = PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));CHKERRQ(ierr); 468 PetscFunctionReturn(0); 469 } 470 471 /*@ 472 TSBDFGetOrder - Get the order of the BDF method 473 474 Not Collective 475 476 Input Parameter: 477 . ts - timestepping context 478 479 Output Parameter: 480 . order - order of the method 481 482 Level: intermediate 483 484 @*/ 485 PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order) 486 { 487 PetscErrorCode ierr; 488 489 PetscFunctionBegin; 490 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 491 PetscValidIntPointer(order,2); 492 ierr = PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));CHKERRQ(ierr); 493 PetscFunctionReturn(0); 494 } 495