Actual source code: bdf.c

petsc-master 2019-10-15
Report Typos and Errors
  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:   PetscReal shift;
 21:   Vec       vec_dot;
 22:   Vec       vec_wrk;
 23:   Vec       vec_lte;

 25:   PetscInt     order;
 26:   TSStepStatus status;
 27: } TS_BDF;


 30: PETSC_STATIC_INLINE void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])
 31: {
 32:   PetscInt k,j;
 33:   for (k=0; k<n; k++)
 34:     for (L[k]=1, j=0; j<n; j++)
 35:       if (j != k)
 36:         L[k] *= (t - T[j])/(T[k] - T[j]);
 37: }

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

 53: static PetscErrorCode TSBDF_GetVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
 54: {
 55:   TS_BDF         *bdf = (TS_BDF*)ts->data;

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

 69: static PetscErrorCode TSBDF_RestoreVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
 70: {
 71:   TS_BDF         *bdf = (TS_BDF*)ts->data;

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

 87: static PetscErrorCode DMCoarsenHook_TSBDF(DM fine,DM coarse,void *ctx)
 88: {
 90:   return(0);
 91: }

 93: static PetscErrorCode DMRestrictHook_TSBDF(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
 94: {
 95:   TS             ts = (TS)ctx;
 96:   Vec            Ydot,Ydot_c;
 97:   Vec            Xdot,Xdot_c;

101:   TSBDF_GetVecs(ts,fine,&Xdot,&Ydot);
102:   TSBDF_GetVecs(ts,coarse,&Xdot_c,&Ydot_c);

104:   MatRestrict(restrct,Ydot,Ydot_c);
105:   VecPointwiseMult(Ydot_c,rscale,Ydot_c);

107:   TSBDF_RestoreVecs(ts,fine,&Xdot,&Ydot);
108:   TSBDF_RestoreVecs(ts,coarse,&Xdot_c,&Ydot_c);
109:   return(0);
110: }

112: static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X)
113: {
114:   TS_BDF         *bdf = (TS_BDF*)ts->data;
115:   PetscInt       i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec));
116:   Vec            tail = bdf->work[n-1];

120:   for (i=n-1; i>=2; i--) {
121:     bdf->time[i] = bdf->time[i-1];
122:     bdf->work[i] = bdf->work[i-1];
123:   }
124:   bdf->n       = PetscMin(bdf->n+1,n-1);
125:   bdf->time[1] = t;
126:   bdf->work[1] = tail;
127:   VecCopy(X,tail);
128:   return(0);
129: }

131: static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)
132: {
133:   TS_BDF         *bdf = (TS_BDF*)ts->data;
134:   PetscInt       i,n = order+1;
135:   PetscReal      *time = bdf->time;
136:   Vec            *vecs = bdf->work;
137:   PetscScalar    a[8],b[8],alpha[8];

141:   LagrangeBasisDers(n+0,time[0],time,a); a[n] =0;
142:   LagrangeBasisDers(n+1,time[0],time,b);
143:   for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0];
144:   VecZeroEntries(lte);
145:   VecMAXPY(lte,n+1,alpha,vecs);
146:   return(0);
147: }

149: static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)
150: {
151:   TS_BDF         *bdf = (TS_BDF*)ts->data;
152:   PetscInt       n = order+1;
153:   PetscReal      *time = bdf->time+1;
154:   Vec            *vecs = bdf->work+1;
155:   PetscScalar    alpha[7];

159:   n = PetscMin(n,bdf->n);
160:   LagrangeBasisVals(n,t,time,alpha);
161:   VecZeroEntries(X);
162:   VecMAXPY(X,n,alpha,vecs);
163:   return(0);
164: }

166: static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)
167: {
168:   TS_BDF         *bdf = (TS_BDF*)ts->data;
169:   PetscInt       n = order+1;
170:   PetscReal      *time = bdf->time;
171:   Vec            *vecs = bdf->work;
172:   PetscScalar    alpha[7];

176:   LagrangeBasisVals(n,t,time,alpha);
177:   VecZeroEntries(X);
178:   VecMAXPY(X,n,alpha,vecs);
179:   return(0);
180: }

182: static PetscErrorCode TSBDF_PreSolve(TS ts)
183: {
184:   TS_BDF         *bdf = (TS_BDF*)ts->data;
185:   PetscInt       i,n = PetscMax(bdf->k,1) + 1;
186:   Vec            V,V0;
187:   Vec            vecs[7];
188:   PetscScalar    alpha[7];

192:   TSBDF_GetVecs(ts,NULL,&V,&V0);
193:   LagrangeBasisDers(n,bdf->time[0],bdf->time,alpha);
194:   for (i=1; i<n; i++) vecs[i] = bdf->work[i];
195:   VecZeroEntries(V0);
196:   VecMAXPY(V0,n-1,alpha+1,vecs+1);
197:   bdf->shift = PetscRealPart(alpha[0]);
198:   TSBDF_RestoreVecs(ts,NULL,&V,&V0);
199:   return(0);
200: }

202: static PetscErrorCode TSBDF_SNESSolve(TS ts,Vec b,Vec x)
203: {
204:   PetscInt       nits,lits;

208:   TSBDF_PreSolve(ts);
209:   SNESSolve(ts->snes,b,x);
210:   SNESGetIterationNumber(ts->snes,&nits);
211:   SNESGetLinearSolveIterations(ts->snes,&lits);
212:   ts->snes_its += nits; ts->ksp_its += lits;
213:   return(0);
214: }

216: static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
217: {
218:   TS_BDF         *bdf = (TS_BDF*)ts->data;

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

225:   bdf->time[0] = ts->ptime + ts->time_step/2;
226:   VecCopy(bdf->work[1],bdf->work[0]);
227:   TSPreStage(ts,bdf->time[0]);
228:   TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
229:   TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
230:   TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);
231:   if (!*accept) return(0);

233:   bdf->k = PetscMin(2,bdf->order); bdf->n++;
234:   VecCopy(bdf->work[0],bdf->work[2]);
235:   bdf->time[2] = bdf->time[0];
236:   return(0);
237: }

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

241: static PetscErrorCode TSStep_BDF(TS ts)
242: {
243:   TS_BDF         *bdf = (TS_BDF*)ts->data;
244:   PetscInt       rejections = 0;
245:   PetscBool      stageok,accept = PETSC_TRUE;
246:   PetscReal      next_time_step = ts->time_step;

250:   PetscCitationsRegister(citation,&cited);

252:   if (!ts->steprollback && !ts->steprestart) {
253:     bdf->k = PetscMin(bdf->k+1,bdf->order);
254:     TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
255:   }

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

260:     if (ts->steprestart) {
261:       TSBDF_Restart(ts,&stageok);
262:       if (!stageok) goto reject_step;
263:     }

265:     bdf->time[0] = ts->ptime + ts->time_step;
266:     TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);
267:     TSPreStage(ts,bdf->time[0]);
268:     TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
269:     TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
270:     TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);
271:     if (!stageok) goto reject_step;

273:     bdf->status = TS_STEP_PENDING;
274:     TSAdaptCandidatesClear(ts->adapt);
275:     TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);
276:     TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
277:     bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
278:     if (!accept) { ts->time_step = next_time_step; goto reject_step; }

280:     VecCopy(bdf->work[0],ts->vec_sol);
281:     ts->ptime += ts->time_step;
282:     ts->time_step = next_time_step;
283:     break;

285:   reject_step:
286:     ts->reject++; accept = PETSC_FALSE;
287:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
288:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
289:       ts->reason = TS_DIVERGED_STEP_REJECTED;
290:     }
291:   }
292:   return(0);
293: }

295: static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
296: {
297:   TS_BDF         *bdf = (TS_BDF*)ts->data;

301:   TSBDF_Interpolate(ts,bdf->k,t,X);
302:   return(0);
303: }

305: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
306: {
307:   TS_BDF         *bdf = (TS_BDF*)ts->data;
308:   PetscInt       k = bdf->k;
309:   PetscReal      wltea,wlter;
310:   Vec            X = bdf->work[0], Y = bdf->vec_lte;

314:   k = PetscMin(k,bdf->n-1);
315:   TSBDF_VecLTE(ts,k,Y);
316:   VecAXPY(Y,1,X);
317:   TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
318:   if (order) *order = k + 1;
319:   return(0);
320: }

322: static PetscErrorCode TSRollBack_BDF(TS ts)
323: {
324:   TS_BDF         *bdf = (TS_BDF*)ts->data;

328:   VecCopy(bdf->work[1],ts->vec_sol);
329:   return(0);
330: }

332: static PetscErrorCode SNESTSFormFunction_BDF(SNES snes,Vec X,Vec F,TS ts)
333: {
334:   TS_BDF         *bdf = (TS_BDF*)ts->data;
335:   DM             dm, dmsave = ts->dm;
336:   PetscReal      t = bdf->time[0];
337:   PetscReal      shift = bdf->shift;
338:   Vec            V,V0;

342:   SNESGetDM(snes,&dm);
343:   TSBDF_GetVecs(ts,dm,&V,&V0);
344:   VecWAXPY(V,shift,X,V0);

346:   /* F = Function(t,X,V) */
347:   ts->dm = dm;
348:   TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);
349:   ts->dm = dmsave;

351:   TSBDF_RestoreVecs(ts,dm,&V,&V0);
352:   return(0);
353: }

355: static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes,Vec X,Mat J,Mat P,TS ts)
356: {
357:   TS_BDF         *bdf = (TS_BDF*)ts->data;
358:   DM             dm, dmsave = ts->dm;
359:   PetscReal      t = bdf->time[0];
360:   PetscReal      shift = bdf->shift;
361:   Vec            V,V0;

365:   SNESGetDM(snes,&dm);
366:   TSBDF_GetVecs(ts,dm,&V,&V0);

368:   /* J,P = Jacobian(t,X,V) */
369:   ts->dm = dm;
370:   TSComputeIJacobian(ts,t,X,V,shift,J,P,PETSC_FALSE);
371:   ts->dm = dmsave;

373:   TSBDF_RestoreVecs(ts,dm,&V,&V0);
374:   return(0);
375: }

377: static PetscErrorCode TSReset_BDF(TS ts)
378: {
379:   TS_BDF         *bdf = (TS_BDF*)ts->data;
380:   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);

384:   bdf->k = bdf->n = 0;
385:   for (i=0; i<n; i++) {VecDestroy(&bdf->work[i]);}
386:   VecDestroy(&bdf->vec_dot);
387:   VecDestroy(&bdf->vec_wrk);
388:   VecDestroy(&bdf->vec_lte);
389:   if (ts->dm) {DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);}
390:   return(0);
391: }

393: static PetscErrorCode TSDestroy_BDF(TS ts)
394: {

398:   TSReset_BDF(ts);
399:   PetscFree(ts->data);
400:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);
401:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);
402:   return(0);
403: }

405: static PetscErrorCode TSSetUp_BDF(TS ts)
406: {
407:   TS_BDF         *bdf = (TS_BDF*)ts->data;
408:   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);
409:   PetscReal      low,high,two = 2;

413:   bdf->k = bdf->n = 0;
414:   for (i=0; i<n; i++) {VecDuplicate(ts->vec_sol,&bdf->work[i]);}
415:   VecDuplicate(ts->vec_sol,&bdf->vec_dot);
416:   VecDuplicate(ts->vec_sol,&bdf->vec_wrk);
417:   VecDuplicate(ts->vec_sol,&bdf->vec_lte);
418:   TSGetDM(ts,&ts->dm);
419:   DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);

421:   TSGetAdapt(ts,&ts->adapt);
422:   TSAdaptCandidatesClear(ts->adapt);
423:   TSAdaptGetClip(ts->adapt,&low,&high);
424:   TSAdaptSetClip(ts->adapt,low,PetscMin(high,two));

426:   TSGetSNES(ts,&ts->snes);
427:   return(0);
428: }

430: static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
431: {
432:   TS_BDF         *bdf = (TS_BDF*)ts->data;

436:   PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");
437:   {
438:     PetscBool flg;
439:     PetscInt  order = bdf->order;
440:     PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);
441:     if (flg) {TSBDFSetOrder(ts,order);}
442:   }
443:   PetscOptionsTail();
444:   return(0);
445: }

447: static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
448: {
449:   TS_BDF         *bdf = (TS_BDF*)ts->data;
450:   PetscBool      iascii;

454:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
455:   if (iascii) {
456:     PetscViewerASCIIPrintf(viewer,"  Order=%D\n",bdf->order);
457:   }
458:   return(0);
459: }

461: /* ------------------------------------------------------------ */

463: static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
464: {
465:   TS_BDF *bdf = (TS_BDF*)ts->data;

468:   if (order == bdf->order) return(0);
469:   if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order);
470:   bdf->order = order;
471:   return(0);
472: }

474: static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
475: {
476:   TS_BDF *bdf = (TS_BDF*)ts->data;

479:   *order = bdf->order;
480:   return(0);
481: }

483: /* ------------------------------------------------------------ */

485: /*MC
486:       TSBDF - DAE solver using BDF methods

488:   Level: beginner

490: .seealso:  TS, TSCreate(), TSSetType()
491: M*/
492: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
493: {
494:   TS_BDF         *bdf;

498:   ts->ops->reset          = TSReset_BDF;
499:   ts->ops->destroy        = TSDestroy_BDF;
500:   ts->ops->view           = TSView_BDF;
501:   ts->ops->setup          = TSSetUp_BDF;
502:   ts->ops->setfromoptions = TSSetFromOptions_BDF;
503:   ts->ops->step           = TSStep_BDF;
504:   ts->ops->evaluatewlte   = TSEvaluateWLTE_BDF;
505:   ts->ops->rollback       = TSRollBack_BDF;
506:   ts->ops->interpolate    = TSInterpolate_BDF;
507:   ts->ops->snesfunction   = SNESTSFormFunction_BDF;
508:   ts->ops->snesjacobian   = SNESTSFormJacobian_BDF;
509:   ts->default_adapt_type  = TSADAPTBASIC;

511:   ts->usessnes = PETSC_TRUE;

513:   PetscNewLog(ts,&bdf);
514:   ts->data = (void*)bdf;

516:   bdf->order  = 2;
517:   bdf->status = TS_STEP_COMPLETE;

519:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);
520:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);
521:   return(0);
522: }

524: /* ------------------------------------------------------------ */

526: /*@
527:   TSBDFSetOrder - Set the order of the BDF method

529:   Logically Collective on TS

531:   Input Parameter:
532: +  ts - timestepping context
533: -  order - order of the method

535:   Options Database:
536: .  -ts_bdf_order <order>

538:   Level: intermediate

540: @*/
541: PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
542: {

548:   PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));
549:   return(0);
550: }

552: /*@
553:   TSBDFGetOrder - Get the order of the BDF method

555:   Not Collective

557:   Input Parameter:
558: .  ts - timestepping context

560:   Output Parameter:
561: .  order - order of the method

563:   Level: intermediate

565: @*/
566: PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
567: {

573:   PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));
574:   return(0);
575: }