Actual source code: bdf.c

petsc-3.7.7 2017-09-25
Report Typos and Errors
  1: /*
  2:   Code for timestepping with BDF methods
  3: */
  4: #include <petsc/private/tsimpl.h>  /*I "petscts.h" I*/

  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";

 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;

 23:   PetscInt     order;
 24:   PetscBool    adapt;
 25:   TSStepStatus status;
 26: } TS_BDF;


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

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

 54: static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X)
 55: {
 56:   TS_BDF         *bdf = (TS_BDF*)ts->data;
 57:   PetscInt       i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec));
 58:   Vec            tail = bdf->work[n-1];

 62:   for (i=n-1; i>=2; i--) {
 63:     bdf->time[i] = bdf->time[i-1];
 64:     bdf->work[i] = bdf->work[i-1];
 65:   }
 66:   bdf->n       = PetscMin(bdf->n+1,n-1);
 67:   bdf->time[1] = t;
 68:   bdf->work[1] = tail;
 69:   VecCopy(X,tail);
 70:   return(0);
 71: }

 75: static PetscErrorCode TSBDF_VecDot(TS ts,PetscInt order,PetscReal t,Vec X,Vec Xdot,PetscReal *shift)
 76: {
 77:   TS_BDF         *bdf = (TS_BDF*)ts->data;
 78:   PetscInt       i,n = order+1;
 79:   PetscReal      time[7];
 80:   Vec            vecs[7];
 81:   PetscScalar    alpha[7];

 85:   n = PetscMax(n,2);
 86:   time[0] = t; for (i=1; i<n; i++) time[i] = bdf->time[i];
 87:   vecs[0] = X; for (i=1; i<n; i++) vecs[i] = bdf->work[i];
 88:   LagrangeBasisDers(n,t,time,alpha);
 89:   VecZeroEntries(Xdot);
 90:   VecMAXPY(Xdot,n,alpha,vecs);
 91:   if (shift) *shift = PetscRealPart(alpha[0]);
 92:   return(0);
 93: }

 97: static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)
 98: {
 99:   TS_BDF         *bdf = (TS_BDF*)ts->data;
100:   PetscInt       i,n = order+1;
101:   PetscReal      *time = bdf->time;
102:   Vec            *vecs = bdf->work;
103:   PetscScalar    a[8],b[8],alpha[8];

107:   LagrangeBasisDers(n+0,time[0],time,a); a[n] =0;
108:   LagrangeBasisDers(n+1,time[0],time,b);
109:   for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0];
110:   VecZeroEntries(lte);
111:   VecMAXPY(lte,n+1,alpha,vecs);
112:   return(0);
113: }

117: static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)
118: {
119:   TS_BDF         *bdf = (TS_BDF*)ts->data;
120:   PetscInt       n = order+1;
121:   PetscReal      *time = bdf->time+1;
122:   Vec            *vecs = bdf->work+1;
123:   PetscScalar    alpha[7];

127:   n = PetscMin(n,bdf->n);
128:   LagrangeBasisVals(n,t,time,alpha);
129:   VecZeroEntries(X);
130:   VecMAXPY(X,n,alpha,vecs);
131:   return(0);
132: }

136: static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)
137: {
138:   TS_BDF         *bdf = (TS_BDF*)ts->data;
139:   PetscInt       n = order+1;
140:   PetscReal      *time = bdf->time;
141:   Vec            *vecs = bdf->work;
142:   PetscScalar    alpha[7];

146:   LagrangeBasisVals(n,t,time,alpha);
147:   VecZeroEntries(X);
148:   VecMAXPY(X,n,alpha,vecs);
149:   return(0);
150: }

154: static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x)
155: {
156:   PetscInt       nits,lits;

160:   SNESSolve(ts->snes,b,x);
161:   SNESGetIterationNumber(ts->snes,&nits);
162:   SNESGetLinearSolveIterations(ts->snes,&lits);
163:   ts->snes_its += nits; ts->ksp_its += lits;
164:   return(0);
165: }

169: static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
170: {
171:   TS_BDF         *bdf = (TS_BDF*)ts->data;

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

178:   bdf->time[0] = ts->ptime + ts->time_step/2;
179:   VecCopy(bdf->work[1],bdf->work[0]);
180:   TSPreStage(ts,bdf->time[0]);
181:   TS_SNESSolve(ts,NULL,bdf->work[0]);
182:   TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
183:   TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);
184:   if (!*accept) return(0);

186:   bdf->k = PetscMin(2,bdf->order); bdf->n++;
187:   VecCopy(bdf->work[0],bdf->work[2]);
188:   bdf->time[2] = bdf->time[0];
189:   return(0);
190: }

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

196: static PetscErrorCode TSStep_BDF(TS ts)
197: {
198:   TS_BDF         *bdf = (TS_BDF*)ts->data;
199:   PetscInt       rejections = 0;
200:   PetscBool      stageok,accept = PETSC_TRUE;
201:   PetscReal      next_time_step = ts->time_step;

205:   PetscCitationsRegister(citation,&cited);

207:   if (!ts->steprollback && !ts->steprestart) {
208:     bdf->k = PetscMin(bdf->k+1,bdf->order);
209:     TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
210:   }

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

215:     if (ts->steprestart) {
216:       TSBDF_Restart(ts,&stageok);
217:       if (!stageok) goto reject_step;
218:     }

220:     bdf->time[0] = ts->ptime + ts->time_step;
221:     TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);
222:     TSPreStage(ts,bdf->time[0]);
223:     TS_SNESSolve(ts,NULL,bdf->work[0]);
224:     TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
225:     TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);
226:     if (!stageok) goto reject_step;

228:     bdf->status = TS_STEP_PENDING;
229:     TSAdaptCandidatesClear(ts->adapt);
230:     TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);
231:     TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
232:     bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
233:     if (!accept) { ts->time_step = next_time_step; goto reject_step; }

235:     VecCopy(bdf->work[0],ts->vec_sol);
236:     ts->ptime += ts->time_step;
237:     ts->time_step = next_time_step;
238:     break;

240:   reject_step:
241:     ts->reject++; accept = PETSC_FALSE;
242:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
243:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
244:       ts->reason = TS_DIVERGED_STEP_REJECTED;
245:     }
246:   }
247:   return(0);
248: }

252: static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
253: {
254:   TS_BDF         *bdf = (TS_BDF*)ts->data;

258:   TSBDF_Interpolate(ts,bdf->k,t,X);
259:   return(0);
260: }

264: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
265: {
266:   TS_BDF         *bdf = (TS_BDF*)ts->data;
267:   PetscInt       k = bdf->k;
268:   Vec            X = bdf->work[0], Y = bdf->vec_lte;

272:   k = PetscMin(k,bdf->n-1);
273:   TSBDF_VecLTE(ts,k,Y);
274:   VecAXPY(Y,1,X);
275:   TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte);
276:   if (order) *order = k + 1;
277:   return(0);
278: }

282: static PetscErrorCode TSRollBack_BDF(TS ts)
283: {
284:   TS_BDF         *bdf = (TS_BDF*)ts->data;

288:   VecCopy(bdf->work[1],ts->vec_sol);
289:   return(0);
290: }

294: static PetscErrorCode SNESTSFormFunction_BDF(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
295: {
296:   TS_BDF         *bdf = (TS_BDF*)ts->data;
297:   PetscInt       k = bdf->k;
298:   PetscReal      t = bdf->time[0];
299:   Vec            V = bdf->vec_dot;

303:   TSBDF_VecDot(ts,k,t,X,V,&bdf->shift);
304:   /* F = Function(t,X,V) */
305:   TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);
306:   return(0);
307: }

311: static PetscErrorCode SNESTSFormJacobian_BDF(PETSC_UNUSED SNES snes,
312:                                              PETSC_UNUSED Vec X,
313:                                              Mat J,Mat P,
314:                                              TS ts)
315: {
316:   TS_BDF         *bdf = (TS_BDF*)ts->data;
317:   PetscReal      t = bdf->time[0];
318:   Vec            V = bdf->vec_dot;
319:   PetscReal      dVdX = bdf->shift;

323:   /* J,P = Jacobian(t,X,V) */
324:   TSComputeIJacobian(ts,t,X,V,dVdX,J,P,PETSC_FALSE);
325:   return(0);
326: }

330: static PetscErrorCode TSReset_BDF(TS ts)
331: {
332:   TS_BDF         *bdf = (TS_BDF*)ts->data;
333:   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);

337:   for (i=0; i<n; i++) {VecDestroy(&bdf->work[i]);}
338:   VecDestroy(&bdf->vec_dot);
339:   VecDestroy(&bdf->vec_lte);
340:   return(0);
341: }

345: static PetscErrorCode TSDestroy_BDF(TS ts)
346: {

350:   TSReset_BDF(ts);
351:   PetscFree(ts->data);
352:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);
353:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);
354:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",NULL);
355:   return(0);
356: }

360: static PetscErrorCode TSSetUp_BDF(TS ts)
361: {
362:   TS_BDF         *bdf = (TS_BDF*)ts->data;
363:   size_t         i,n = sizeof(bdf->work)/sizeof(Vec);

367:   bdf->k = bdf->n = 0;
368:   for (i=0; i<n; i++) {VecDuplicate(ts->vec_sol,&bdf->work[i]);}
369:   VecDuplicate(ts->vec_sol,&bdf->vec_dot);
370:   VecDuplicate(ts->vec_sol,&bdf->vec_lte);

372:   TSGetAdapt(ts,&ts->adapt);
373:   TSAdaptCandidatesClear(ts->adapt);
374:   if (!bdf->adapt) {
375:     TSAdaptSetType(ts->adapt,TSADAPTNONE);
376:   } else {
377:     PetscReal low,high;
378:     TSAdaptBasicGetClip(ts->adapt,&low,&high);
379:     high = PetscMin(high,2.0);
380:     TSAdaptBasicSetClip(ts->adapt,low,high);
381:     if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED)
382:       ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
383:   }

385:   TSGetSNES(ts,&ts->snes);
386:   return(0);
387: }

391: static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
392: {
393:   TS_BDF         *bdf = (TS_BDF*)ts->data;

397:   PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");
398:   {
399:     PetscBool flg;
400:     PetscInt  order = bdf->order;
401:     PetscBool adapt = bdf->adapt;
402:     PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);
403:     if (flg) {TSBDFSetOrder(ts,order);}
404:     PetscOptionsBool("-ts_bdf_adapt","Use time-step adaptivity with the BDF method","TSBDFUseAdapt",adapt,&adapt,&flg);
405:     if (flg) {TSBDFUseAdapt(ts,adapt);}
406:   }
407:   PetscOptionsTail();
408:   return(0);
409: }

413: static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
414: {
415:   TS_BDF         *bdf = (TS_BDF*)ts->data;
416:   PetscBool      iascii;

420:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
421:   if (iascii)    {
422:     PetscViewerASCIIPrintf(viewer,"  Order=%D\n",bdf->order);
423:   }
424:   if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
425:   if (ts->snes)  {SNESView(ts->snes,viewer);}
426:   return(0);
427: }

429: /* ------------------------------------------------------------ */

433: static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
434: {
435:   TS_BDF *bdf = (TS_BDF*)ts->data;

438:   if (order == bdf->order) return(0);
439:   if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order);
440:   bdf->order = order;
441:   return(0);
442: }

446: static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
447: {
448:   TS_BDF *bdf = (TS_BDF*)ts->data;

451:   *order = bdf->order;
452:   return(0);
453: }

457: static PetscErrorCode TSBDFUseAdapt_BDF(TS ts,PetscBool use)
458: {
459:   TS_BDF *bdf = (TS_BDF*)ts->data;

462:   if (use == bdf->adapt) return(0);
463:   if (ts->setupcalled) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ORDER,"Cannot change adaptivity after TSSetUp()");
464:   bdf->adapt = use;
465:   return(0);
466: }

468: /* ------------------------------------------------------------ */

470: /*MC
471:       TSBDF - DAE solver using BDF methods

473:   Level: beginner

475: .seealso:  TS, TSCreate(), TSSetType()
476: M*/
479: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
480: {
481:   TS_BDF         *bdf;

485:   ts->ops->reset          = TSReset_BDF;
486:   ts->ops->destroy        = TSDestroy_BDF;
487:   ts->ops->view           = TSView_BDF;
488:   ts->ops->setup          = TSSetUp_BDF;
489:   ts->ops->setfromoptions = TSSetFromOptions_BDF;
490:   ts->ops->step           = TSStep_BDF;
491:   ts->ops->evaluatewlte   = TSEvaluateWLTE_BDF;
492:   ts->ops->rollback       = TSRollBack_BDF;
493:   ts->ops->interpolate    = TSInterpolate_BDF;
494:   ts->ops->snesfunction   = SNESTSFormFunction_BDF;
495:   ts->ops->snesjacobian   = SNESTSFormJacobian_BDF;

497:   PetscNewLog(ts,&bdf);
498:   ts->data = (void*)bdf;

500:   bdf->order  = 2;
501:   bdf->adapt  = PETSC_FALSE;
502:   bdf->status = TS_STEP_COMPLETE;

504:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);
505:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);
506:   PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",TSBDFUseAdapt_BDF);
507:   return(0);
508: }

510: /* ------------------------------------------------------------ */

514: /*@
515:   TSBDFSetOrder - Set the order of the BDF method

517:   Logically Collective on TS

519:   Input Parameter:
520: +  ts - timestepping context
521: -  order - order of the method

523:   Options Database:
524: .  -ts_bdf_order <order>

526:   Level: intermediate

528: @*/
529: PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
530: {

536:   PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));
537:   return(0);
538: }

542: /*@
543:   TSBDFGetOrder - Get the order of the BDF method

545:   Not Collective

547:   Input Parameter:
548: .  ts - timestepping context

550:   Output Parameter:
551: .  order - order of the method

553:   Level: intermediate

555: @*/
556: PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
557: {

563:   PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));
564:   return(0);
565: }

569: /*@
570:   TSBDFUseAdapt - Use time-step adaptivity with the BDF method

572:   Logically Collective on TS

574:   Input Parameter:
575: +  ts - timestepping context
576: -  use - flag to use adaptivity

578:   Options Database:
579: .  -ts_bdf_adapt

581:   Level: intermediate

583: .seealso: TSAdapt
584: @*/
585: PetscErrorCode TSBDFUseAdapt(TS ts,PetscBool use)
586: {

592:   PetscTryMethod(ts,"TSBDFUseAdapt_C",(TS,PetscBool),(ts,use));
593:   return(0);
594: }