Actual source code: bdf.c

  1: /*
  2:   Code for timestepping with BDF methods
  3: */
  4: #include <petsc/private/tsimpl.h>
  5: #include <petscdm.h>

  7: static PetscBool  cited = PETSC_FALSE;
  8: static const char citation[] =
  9:   "@book{Brenan1995,\n"
 10:   "  title     = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n"
 11:   "  author    = {Brenan, K. and Campbell, S. and Petzold, L.},\n"
 12:   "  publisher = {Society for Industrial and Applied Mathematics},\n"
 13:   "  year      = {1995},\n"
 14:   "  doi       = {10.1137/1.9781611971224},\n}\n";

 16: typedef struct {
 17:   PetscInt  k,n;
 18:   PetscReal time[6+2];
 19:   Vec       work[6+2];
 20:   Vec       tvwork[6+2];
 21:   PetscReal shift;
 22:   Vec       vec_dot;            /* Xdot when !transientvar, else Cdot where C(X) is the transient variable. */
 23:   Vec       vec_wrk;
 24:   Vec       vec_lte;

 26:   PetscBool    transientvar;
 27:   PetscInt     order;
 28:   TSStepStatus status;
 29: } TS_BDF;

 31: /* Compute Lagrange polynomials on T[:n] evaluated at t.
 32:  * If one has data (T[i], Y[i]), then the interpolation/extrapolation is f(t) = \sum_i L[i]*Y[i].
 33:  */
 34: PETSC_STATIC_INLINE void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])
 35: {
 36:   PetscInt k,j;
 37:   for (k=0; k<n; k++)
 38:     for (L[k]=1, j=0; j<n; j++)
 39:       if (j != k)
 40:         L[k] *= (t - T[j])/(T[k] - T[j]);
 41: }

 43: PETSC_STATIC_INLINE void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[])
 44: {
 45:   PetscInt  k,j,i;
 46:   for (k=0; k<n; k++)
 47:     for (dL[k]=0, j=0; j<n; j++)
 48:       if (j != k) {
 49:         PetscReal L = 1/(T[k] - T[j]);
 50:         for (i=0; i<n; i++)
 51:           if (i != j && i != k)
 52:             L *= (t - T[i])/(T[k] - T[i]);
 53:         dL[k] += L;
 54:       }
 55: }

 57: static PetscErrorCode TSBDF_GetVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
 58: {
 59:   TS_BDF         *bdf = (TS_BDF*)ts->data;

 63:   if (dm && dm != ts->dm) {
 64:     DMGetNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot);
 65:     DMGetNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot);
 66:   } else {
 67:     *Xdot = bdf->vec_dot;
 68:     *Ydot = bdf->vec_wrk;
 69:   }
 70:   return(0);
 71: }

 73: static PetscErrorCode TSBDF_RestoreVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
 74: {
 75:   TS_BDF         *bdf = (TS_BDF*)ts->data;

 79:   if (dm && dm != ts->dm) {
 80:     DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot);
 81:     DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot);
 82:   } else {
 83:     if (*Xdot != bdf->vec_dot) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Vec does not match the cache");
 84:     if (*Ydot != bdf->vec_wrk) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Vec does not match the cache");
 85:     *Xdot = NULL;
 86:     *Ydot = NULL;
 87:   }
 88:   return(0);
 89: }

 91: static PetscErrorCode DMCoarsenHook_TSBDF(DM fine,DM coarse,void *ctx)
 92: {
 94:   return(0);
 95: }

 97: static PetscErrorCode DMRestrictHook_TSBDF(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
 98: {
 99:   TS             ts = (TS)ctx;
100:   Vec            Ydot,Ydot_c;
101:   Vec            Xdot,Xdot_c;

105:   TSBDF_GetVecs(ts,fine,&Xdot,&Ydot);
106:   TSBDF_GetVecs(ts,coarse,&Xdot_c,&Ydot_c);

108:   MatRestrict(restrct,Ydot,Ydot_c);
109:   VecPointwiseMult(Ydot_c,rscale,Ydot_c);

111:   TSBDF_RestoreVecs(ts,fine,&Xdot,&Ydot);
112:   TSBDF_RestoreVecs(ts,coarse,&Xdot_c,&Ydot_c);
113:   return(0);
114: }

116: static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X)
117: {
118:   TS_BDF         *bdf = (TS_BDF*)ts->data;
119:   PetscInt       i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec));
120:   Vec            tail = bdf->work[n-1],tvtail = bdf->tvwork[n-1];

124:   for (i=n-1; i>=2; i--) {
125:     bdf->time[i] = bdf->time[i-1];
126:     bdf->work[i] = bdf->work[i-1];
127:     bdf->tvwork[i] = bdf->tvwork[i-1];
128:   }
129:   bdf->n       = PetscMin(bdf->n+1,n-1);
130:   bdf->time[1] = t;
131:   bdf->work[1] = tail;
132:   bdf->tvwork[1] = tvtail;
133:   VecCopy(X,tail);
134:   TSComputeTransientVariable(ts,tail,tvtail);
135:   return(0);
136: }

138: static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)
139: {
140:   TS_BDF         *bdf = (TS_BDF*)ts->data;
141:   PetscInt       i,n = order+1;
142:   PetscReal      *time = bdf->time;
143:   Vec            *vecs = bdf->work;
144:   PetscScalar    a[8],b[8],alpha[8];

148:   LagrangeBasisDers(n+0,time[0],time,a); a[n] =0;
149:   LagrangeBasisDers(n+1,time[0],time,b);
150:   for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0];
151:   VecZeroEntries(lte);
152:   VecMAXPY(lte,n+1,alpha,vecs);
153:   return(0);
154: }

156: static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)
157: {
158:   TS_BDF         *bdf = (TS_BDF*)ts->data;
159:   PetscInt       n = order+1;
160:   PetscReal      *time = bdf->time+1;
161:   Vec            *vecs = bdf->work+1;
162:   PetscScalar    alpha[7];

166:   n = PetscMin(n,bdf->n);
167:   LagrangeBasisVals(n,t,time,alpha);
168:   VecZeroEntries(X);
169:   VecMAXPY(X,n,alpha,vecs);
170:   return(0);
171: }

173: static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)
174: {
175:   TS_BDF         *bdf = (TS_BDF*)ts->data;
176:   PetscInt       n = order+1;
177:   PetscReal      *time = bdf->time;
178:   Vec            *vecs = bdf->work;
179:   PetscScalar    alpha[7];

183:   LagrangeBasisVals(n,t,time,alpha);
184:   VecZeroEntries(X);
185:   VecMAXPY(X,n,alpha,vecs);
186:   return(0);
187: }

189: /* Compute the affine term V0 such that Xdot = shift*X + V0.
190:  *
191:  * When using transient variables, we're computing Cdot = shift*C(X) + V0, and thus choose a linear combination of tvwork.
192:  */
193: static PetscErrorCode TSBDF_PreSolve(TS ts)
194: {
195:   TS_BDF         *bdf = (TS_BDF*)ts->data;
196:   PetscInt       i,n = PetscMax(bdf->k,1) + 1;
197:   Vec            V,V0;
198:   Vec            vecs[7];
199:   PetscScalar    alpha[7];

203:   TSBDF_GetVecs(ts,NULL,&V,&V0);
204:   LagrangeBasisDers(n,bdf->time[0],bdf->time,alpha);
205:   for (i=1; i<n; i++) {
206:     vecs[i] = bdf->transientvar ? bdf->tvwork[i] : bdf->work[i];
207:   }
208:   VecZeroEntries(V0);
209:   VecMAXPY(V0,n-1,alpha+1,vecs+1);
210:   bdf->shift = PetscRealPart(alpha[0]);
211:   TSBDF_RestoreVecs(ts,NULL,&V,&V0);
212:   return(0);
213: }

215: static PetscErrorCode TSBDF_SNESSolve(TS ts,Vec b,Vec x)
216: {
217:   PetscInt       nits,lits;

221:   TSBDF_PreSolve(ts);
222:   SNESSolve(ts->snes,b,x);
223:   SNESGetIterationNumber(ts->snes,&nits);
224:   SNESGetLinearSolveIterations(ts->snes,&lits);
225:   ts->snes_its += nits; ts->ksp_its += lits;
226:   return(0);
227: }

229: static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
230: {
231:   TS_BDF         *bdf = (TS_BDF*)ts->data;

235:   bdf->k = 1; bdf->n = 0;
236:   TSBDF_Advance(ts,ts->ptime,ts->vec_sol);

238:   bdf->time[0] = ts->ptime + ts->time_step/2;
239:   VecCopy(bdf->work[1],bdf->work[0]);
240:   TSPreStage(ts,bdf->time[0]);
241:   TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
242:   TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
243:   TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);
244:   if (!*accept) return(0);

246:   bdf->k = PetscMin(2,bdf->order); bdf->n++;
247:   VecCopy(bdf->work[0],bdf->work[2]);
248:   bdf->time[2] = bdf->time[0];
249:   TSComputeTransientVariable(ts,bdf->work[2],bdf->tvwork[2]);
250:   return(0);
251: }

253: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};

255: static PetscErrorCode TSStep_BDF(TS ts)
256: {
257:   TS_BDF         *bdf = (TS_BDF*)ts->data;
258:   PetscInt       rejections = 0;
259:   PetscBool      stageok,accept = PETSC_TRUE;
260:   PetscReal      next_time_step = ts->time_step;

264:   PetscCitationsRegister(citation,&cited);

266:   if (!ts->steprollback && !ts->steprestart) {
267:     bdf->k = PetscMin(bdf->k+1,bdf->order);
268:     TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
269:   }

271:   bdf->status = TS_STEP_INCOMPLETE;
272:   while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {

274:     if (ts->steprestart) {
275:       TSBDF_Restart(ts,&stageok);
276:       if (!stageok) goto reject_step;
277:     }

279:     bdf->time[0] = ts->ptime + ts->time_step;
280:     TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);
281:     TSPreStage(ts,bdf->time[0]);
282:     TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
283:     TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
284:     TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);
285:     if (!stageok) goto reject_step;

287:     bdf->status = TS_STEP_PENDING;
288:     TSAdaptCandidatesClear(ts->adapt);
289:     TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);
290:     TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
291:     bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
292:     if (!accept) { ts->time_step = next_time_step; goto reject_step; }

294:     VecCopy(bdf->work[0],ts->vec_sol);
295:     ts->ptime += ts->time_step;
296:     ts->time_step = next_time_step;
297:     break;

299:   reject_step:
300:     ts->reject++; accept = PETSC_FALSE;
301:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
302:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
303:       ts->reason = TS_DIVERGED_STEP_REJECTED;
304:     }
305:   }
306:   return(0);
307: }

309: static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
310: {
311:   TS_BDF         *bdf = (TS_BDF*)ts->data;

315:   TSBDF_Interpolate(ts,bdf->k,t,X);
316:   return(0);
317: }

319: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
320: {
321:   TS_BDF         *bdf = (TS_BDF*)ts->data;
322:   PetscInt       k = bdf->k;
323:   PetscReal      wltea,wlter;
324:   Vec            X = bdf->work[0], Y = bdf->vec_lte;

328:   k = PetscMin(k,bdf->n-1);
329:   TSBDF_VecLTE(ts,k,Y);
330:   VecAXPY(Y,1,X);
331:   TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
332:   if (order) *order = k + 1;
333:   return(0);
334: }

336: static PetscErrorCode TSRollBack_BDF(TS ts)
337: {
338:   TS_BDF         *bdf = (TS_BDF*)ts->data;

342:   VecCopy(bdf->work[1],ts->vec_sol);
343:   return(0);
344: }

346: static PetscErrorCode SNESTSFormFunction_BDF(SNES snes,Vec X,Vec F,TS ts)
347: {
348:   TS_BDF         *bdf = (TS_BDF*)ts->data;
349:   DM             dm, dmsave = ts->dm;
350:   PetscReal      t = bdf->time[0];
351:   PetscReal      shift = bdf->shift;
352:   Vec            V,V0;

356:   SNESGetDM(snes,&dm);
357:   TSBDF_GetVecs(ts,dm,&V,&V0);
358:   if (bdf->transientvar) {      /* shift*C(X) + V0 */
359:     TSComputeTransientVariable(ts,X,V);
360:     VecAYPX(V,shift,V0);
361:   } else {                      /* shift*X + V0 */
362:     VecWAXPY(V,shift,X,V0);
363:   }

365:   /* F = Function(t,X,V) */
366:   ts->dm = dm;
367:   TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);
368:   ts->dm = dmsave;

370:   TSBDF_RestoreVecs(ts,dm,&V,&V0);
371:   return(0);
372: }

374: static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes,Vec X,Mat J,Mat P,TS ts)
375: {
376:   TS_BDF         *bdf = (TS_BDF*)ts->data;
377:   DM             dm, dmsave = ts->dm;
378:   PetscReal      t = bdf->time[0];
379:   PetscReal      shift = bdf->shift;
380:   Vec            V,V0;

384:   SNESGetDM(snes,&dm);
385:   TSBDF_GetVecs(ts,dm,&V,&V0);

387:   /* J,P = Jacobian(t,X,V) */
388:   ts->dm = dm;
389:   TSComputeIJacobian(ts,t,X,V,shift,J,P,PETSC_FALSE);
390:   ts->dm = dmsave;

392:   TSBDF_RestoreVecs(ts,dm,&V,&V0);
393:   return(0);
394: }

396: static PetscErrorCode TSReset_BDF(TS ts)
397: {
398:   TS_BDF         *bdf = (TS_BDF*)ts->data;
399:   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);

403:   bdf->k = bdf->n = 0;
404:   for (i=0; i<n; i++) {
405:     VecDestroy(&bdf->work[i]);
406:     VecDestroy(&bdf->tvwork[i]);
407:   }
408:   VecDestroy(&bdf->vec_dot);
409:   VecDestroy(&bdf->vec_wrk);
410:   VecDestroy(&bdf->vec_lte);
411:   if (ts->dm) {DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);}
412:   return(0);
413: }

415: static PetscErrorCode TSDestroy_BDF(TS ts)
416: {

420:   TSReset_BDF(ts);
421:   PetscFree(ts->data);
422:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);
423:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);
424:   return(0);
425: }

427: static PetscErrorCode TSSetUp_BDF(TS ts)
428: {
429:   TS_BDF         *bdf = (TS_BDF*)ts->data;
430:   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);
431:   PetscReal      low,high,two = 2;

435:   TSHasTransientVariable(ts,&bdf->transientvar);
436:   bdf->k = bdf->n = 0;
437:   for (i=0; i<n; i++) {
438:     VecDuplicate(ts->vec_sol,&bdf->work[i]);
439:     if (i && bdf->transientvar) {
440:       VecDuplicate(ts->vec_sol,&bdf->tvwork[i]);
441:     }
442:   }
443:   VecDuplicate(ts->vec_sol,&bdf->vec_dot);
444:   VecDuplicate(ts->vec_sol,&bdf->vec_wrk);
445:   VecDuplicate(ts->vec_sol,&bdf->vec_lte);
446:   TSGetDM(ts,&ts->dm);
447:   DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);

449:   TSGetAdapt(ts,&ts->adapt);
450:   TSAdaptCandidatesClear(ts->adapt);
451:   TSAdaptGetClip(ts->adapt,&low,&high);
452:   TSAdaptSetClip(ts->adapt,low,PetscMin(high,two));

454:   TSGetSNES(ts,&ts->snes);
455:   return(0);
456: }

458: static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
459: {

463:   PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");
464:   {
465:     PetscBool flg;
466:     PetscInt  order;
467:     TSBDFGetOrder(ts,&order);
468:     PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);
469:     if (flg) {TSBDFSetOrder(ts,order);}
470:   }
471:   PetscOptionsTail();
472:   return(0);
473: }

475: static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
476: {
477:   TS_BDF         *bdf = (TS_BDF*)ts->data;
478:   PetscBool      iascii;

482:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
483:   if (iascii) {
484:     PetscViewerASCIIPrintf(viewer,"  Order=%D\n",bdf->order);
485:   }
486:   return(0);
487: }

489: /* ------------------------------------------------------------ */

491: static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
492: {
493:   TS_BDF *bdf = (TS_BDF*)ts->data;

496:   if (order == bdf->order) return(0);
497:   if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order);
498:   bdf->order = order;
499:   return(0);
500: }

502: static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
503: {
504:   TS_BDF *bdf = (TS_BDF*)ts->data;

507:   *order = bdf->order;
508:   return(0);
509: }

511: /* ------------------------------------------------------------ */

513: /*MC
514:       TSBDF - DAE solver using BDF methods

516:   Level: beginner

518: .seealso:  TS, TSCreate(), TSSetType()
519: M*/
520: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
521: {
522:   TS_BDF         *bdf;

526:   ts->ops->reset          = TSReset_BDF;
527:   ts->ops->destroy        = TSDestroy_BDF;
528:   ts->ops->view           = TSView_BDF;
529:   ts->ops->setup          = TSSetUp_BDF;
530:   ts->ops->setfromoptions = TSSetFromOptions_BDF;
531:   ts->ops->step           = TSStep_BDF;
532:   ts->ops->evaluatewlte   = TSEvaluateWLTE_BDF;
533:   ts->ops->rollback       = TSRollBack_BDF;
534:   ts->ops->interpolate    = TSInterpolate_BDF;
535:   ts->ops->snesfunction   = SNESTSFormFunction_BDF;
536:   ts->ops->snesjacobian   = SNESTSFormJacobian_BDF;
537:   ts->default_adapt_type  = TSADAPTBASIC;

539:   ts->usessnes = PETSC_TRUE;

541:   PetscNewLog(ts,&bdf);
542:   ts->data = (void*)bdf;

544:   bdf->status = TS_STEP_COMPLETE;

546:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);
547:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);
548:   TSBDFSetOrder(ts,2);
549:   return(0);
550: }

552: /* ------------------------------------------------------------ */

554: /*@
555:   TSBDFSetOrder - Set the order of the BDF method

557:   Logically Collective on TS

559:   Input Parameter:
560: +  ts - timestepping context
561: -  order - order of the method

563:   Options Database:
564: .  -ts_bdf_order <order>

566:   Level: intermediate

568: @*/
569: PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
570: {

576:   PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));
577:   return(0);
578: }

580: /*@
581:   TSBDFGetOrder - Get the order of the BDF method

583:   Not Collective

585:   Input Parameter:
586: .  ts - timestepping context

588:   Output Parameter:
589: .  order - order of the method

591:   Level: intermediate

593: @*/
594: PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
595: {

601:   PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));
602:   return(0);
603: }