Actual source code: bdf.c
petsc-3.12.5 2020-03-29
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: }