Actual source code: bdf.c
petsc-3.14.6 2021-03-30
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;
32: /* Compute Lagrange polynomials on T[:n] evaluated at t.
33: * If one has data (T[i], Y[i]), then the interpolation/extrapolation is f(t) = \sum_i L[i]*Y[i].
34: */
35: PETSC_STATIC_INLINE void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])
36: {
37: PetscInt k,j;
38: for (k=0; k<n; k++)
39: for (L[k]=1, j=0; j<n; j++)
40: if (j != k)
41: L[k] *= (t - T[j])/(T[k] - T[j]);
42: }
44: PETSC_STATIC_INLINE void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[])
45: {
46: PetscInt k,j,i;
47: for (k=0; k<n; k++)
48: for (dL[k]=0, j=0; j<n; j++)
49: if (j != k) {
50: PetscReal L = 1/(T[k] - T[j]);
51: for (i=0; i<n; i++)
52: if (i != j && i != k)
53: L *= (t - T[i])/(T[k] - T[i]);
54: dL[k] += L;
55: }
56: }
58: static PetscErrorCode TSBDF_GetVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
59: {
60: TS_BDF *bdf = (TS_BDF*)ts->data;
64: if (dm && dm != ts->dm) {
65: DMGetNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot);
66: DMGetNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot);
67: } else {
68: *Xdot = bdf->vec_dot;
69: *Ydot = bdf->vec_wrk;
70: }
71: return(0);
72: }
74: static PetscErrorCode TSBDF_RestoreVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot)
75: {
76: TS_BDF *bdf = (TS_BDF*)ts->data;
80: if (dm && dm != ts->dm) {
81: DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot);
82: DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot);
83: } else {
84: if (*Xdot != bdf->vec_dot) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Vec does not match the cache");
85: if (*Ydot != bdf->vec_wrk) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Vec does not match the cache");
86: *Xdot = NULL;
87: *Ydot = NULL;
88: }
89: return(0);
90: }
92: static PetscErrorCode DMCoarsenHook_TSBDF(DM fine,DM coarse,void *ctx)
93: {
95: return(0);
96: }
98: static PetscErrorCode DMRestrictHook_TSBDF(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
99: {
100: TS ts = (TS)ctx;
101: Vec Ydot,Ydot_c;
102: Vec Xdot,Xdot_c;
106: TSBDF_GetVecs(ts,fine,&Xdot,&Ydot);
107: TSBDF_GetVecs(ts,coarse,&Xdot_c,&Ydot_c);
109: MatRestrict(restrct,Ydot,Ydot_c);
110: VecPointwiseMult(Ydot_c,rscale,Ydot_c);
112: TSBDF_RestoreVecs(ts,fine,&Xdot,&Ydot);
113: TSBDF_RestoreVecs(ts,coarse,&Xdot_c,&Ydot_c);
114: return(0);
115: }
117: static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X)
118: {
119: TS_BDF *bdf = (TS_BDF*)ts->data;
120: PetscInt i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec));
121: Vec tail = bdf->work[n-1],tvtail = bdf->tvwork[n-1];
125: for (i=n-1; i>=2; i--) {
126: bdf->time[i] = bdf->time[i-1];
127: bdf->work[i] = bdf->work[i-1];
128: bdf->tvwork[i] = bdf->tvwork[i-1];
129: }
130: bdf->n = PetscMin(bdf->n+1,n-1);
131: bdf->time[1] = t;
132: bdf->work[1] = tail;
133: bdf->tvwork[1] = tvtail;
134: VecCopy(X,tail);
135: TSComputeTransientVariable(ts,tail,tvtail);
136: return(0);
137: }
139: static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte)
140: {
141: TS_BDF *bdf = (TS_BDF*)ts->data;
142: PetscInt i,n = order+1;
143: PetscReal *time = bdf->time;
144: Vec *vecs = bdf->work;
145: PetscScalar a[8],b[8],alpha[8];
149: LagrangeBasisDers(n+0,time[0],time,a); a[n] =0;
150: LagrangeBasisDers(n+1,time[0],time,b);
151: for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0];
152: VecZeroEntries(lte);
153: VecMAXPY(lte,n+1,alpha,vecs);
154: return(0);
155: }
157: static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X)
158: {
159: TS_BDF *bdf = (TS_BDF*)ts->data;
160: PetscInt n = order+1;
161: PetscReal *time = bdf->time+1;
162: Vec *vecs = bdf->work+1;
163: PetscScalar alpha[7];
167: n = PetscMin(n,bdf->n);
168: LagrangeBasisVals(n,t,time,alpha);
169: VecZeroEntries(X);
170: VecMAXPY(X,n,alpha,vecs);
171: return(0);
172: }
174: static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X)
175: {
176: TS_BDF *bdf = (TS_BDF*)ts->data;
177: PetscInt n = order+1;
178: PetscReal *time = bdf->time;
179: Vec *vecs = bdf->work;
180: PetscScalar alpha[7];
184: LagrangeBasisVals(n,t,time,alpha);
185: VecZeroEntries(X);
186: VecMAXPY(X,n,alpha,vecs);
187: return(0);
188: }
190: /* Compute the affine term V0 such that Xdot = shift*X + V0.
191: *
192: * When using transient variables, we're computing Cdot = shift*C(X) + V0, and thus choose a linear combination of tvwork.
193: */
194: static PetscErrorCode TSBDF_PreSolve(TS ts)
195: {
196: TS_BDF *bdf = (TS_BDF*)ts->data;
197: PetscInt i,n = PetscMax(bdf->k,1) + 1;
198: Vec V,V0;
199: Vec vecs[7];
200: PetscScalar alpha[7];
204: TSBDF_GetVecs(ts,NULL,&V,&V0);
205: LagrangeBasisDers(n,bdf->time[0],bdf->time,alpha);
206: for (i=1; i<n; i++) {
207: vecs[i] = bdf->transientvar ? bdf->tvwork[i] : bdf->work[i];
208: }
209: VecZeroEntries(V0);
210: VecMAXPY(V0,n-1,alpha+1,vecs+1);
211: bdf->shift = PetscRealPart(alpha[0]);
212: TSBDF_RestoreVecs(ts,NULL,&V,&V0);
213: return(0);
214: }
216: static PetscErrorCode TSBDF_SNESSolve(TS ts,Vec b,Vec x)
217: {
218: PetscInt nits,lits;
222: TSBDF_PreSolve(ts);
223: SNESSolve(ts->snes,b,x);
224: SNESGetIterationNumber(ts->snes,&nits);
225: SNESGetLinearSolveIterations(ts->snes,&lits);
226: ts->snes_its += nits; ts->ksp_its += lits;
227: return(0);
228: }
230: static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept)
231: {
232: TS_BDF *bdf = (TS_BDF*)ts->data;
236: bdf->k = 1; bdf->n = 0;
237: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
239: bdf->time[0] = ts->ptime + ts->time_step/2;
240: VecCopy(bdf->work[1],bdf->work[0]);
241: TSPreStage(ts,bdf->time[0]);
242: TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
243: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
244: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);
245: if (!*accept) return(0);
247: bdf->k = PetscMin(2,bdf->order); bdf->n++;
248: VecCopy(bdf->work[0],bdf->work[2]);
249: bdf->time[2] = bdf->time[0];
250: TSComputeTransientVariable(ts,bdf->work[2],bdf->tvwork[2]);
251: return(0);
252: }
254: static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"};
256: static PetscErrorCode TSStep_BDF(TS ts)
257: {
258: TS_BDF *bdf = (TS_BDF*)ts->data;
259: PetscInt rejections = 0;
260: PetscBool stageok,accept = PETSC_TRUE;
261: PetscReal next_time_step = ts->time_step;
265: PetscCitationsRegister(citation,&cited);
267: if (!ts->steprollback && !ts->steprestart) {
268: bdf->k = PetscMin(bdf->k+1,bdf->order);
269: TSBDF_Advance(ts,ts->ptime,ts->vec_sol);
270: }
272: bdf->status = TS_STEP_INCOMPLETE;
273: while (!ts->reason && bdf->status != TS_STEP_COMPLETE) {
275: if (ts->steprestart) {
276: TSBDF_Restart(ts,&stageok);
277: if (!stageok) goto reject_step;
278: }
280: bdf->time[0] = ts->ptime + ts->time_step;
281: TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);
282: TSPreStage(ts,bdf->time[0]);
283: TSBDF_SNESSolve(ts,NULL,bdf->work[0]);
284: TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);
285: TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);
286: if (!stageok) goto reject_step;
288: bdf->status = TS_STEP_PENDING;
289: TSAdaptCandidatesClear(ts->adapt);
290: TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);
291: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
292: bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
293: if (!accept) { ts->time_step = next_time_step; goto reject_step; }
295: VecCopy(bdf->work[0],ts->vec_sol);
296: ts->ptime += ts->time_step;
297: ts->time_step = next_time_step;
298: break;
300: reject_step:
301: ts->reject++; accept = PETSC_FALSE;
302: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
303: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
304: ts->reason = TS_DIVERGED_STEP_REJECTED;
305: }
306: }
307: return(0);
308: }
310: static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X)
311: {
312: TS_BDF *bdf = (TS_BDF*)ts->data;
316: TSBDF_Interpolate(ts,bdf->k,t,X);
317: return(0);
318: }
320: static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
321: {
322: TS_BDF *bdf = (TS_BDF*)ts->data;
323: PetscInt k = bdf->k;
324: PetscReal wltea,wlter;
325: Vec X = bdf->work[0], Y = bdf->vec_lte;
329: k = PetscMin(k,bdf->n-1);
330: TSBDF_VecLTE(ts,k,Y);
331: VecAXPY(Y,1,X);
332: TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
333: if (order) *order = k + 1;
334: return(0);
335: }
337: static PetscErrorCode TSRollBack_BDF(TS ts)
338: {
339: TS_BDF *bdf = (TS_BDF*)ts->data;
343: VecCopy(bdf->work[1],ts->vec_sol);
344: return(0);
345: }
347: static PetscErrorCode SNESTSFormFunction_BDF(SNES snes,Vec X,Vec F,TS ts)
348: {
349: TS_BDF *bdf = (TS_BDF*)ts->data;
350: DM dm, dmsave = ts->dm;
351: PetscReal t = bdf->time[0];
352: PetscReal shift = bdf->shift;
353: Vec V,V0;
357: SNESGetDM(snes,&dm);
358: TSBDF_GetVecs(ts,dm,&V,&V0);
359: if (bdf->transientvar) { /* shift*C(X) + V0 */
360: TSComputeTransientVariable(ts,X,V);
361: VecAYPX(V,shift,V0);
362: } else { /* shift*X + V0 */
363: VecWAXPY(V,shift,X,V0);
364: }
366: /* F = Function(t,X,V) */
367: ts->dm = dm;
368: TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);
369: ts->dm = dmsave;
371: TSBDF_RestoreVecs(ts,dm,&V,&V0);
372: return(0);
373: }
375: static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes,Vec X,Mat J,Mat P,TS ts)
376: {
377: TS_BDF *bdf = (TS_BDF*)ts->data;
378: DM dm, dmsave = ts->dm;
379: PetscReal t = bdf->time[0];
380: PetscReal shift = bdf->shift;
381: Vec V,V0;
385: SNESGetDM(snes,&dm);
386: TSBDF_GetVecs(ts,dm,&V,&V0);
388: /* J,P = Jacobian(t,X,V) */
389: ts->dm = dm;
390: TSComputeIJacobian(ts,t,X,V,shift,J,P,PETSC_FALSE);
391: ts->dm = dmsave;
393: TSBDF_RestoreVecs(ts,dm,&V,&V0);
394: return(0);
395: }
397: static PetscErrorCode TSReset_BDF(TS ts)
398: {
399: TS_BDF *bdf = (TS_BDF*)ts->data;
400: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
404: bdf->k = bdf->n = 0;
405: for (i=0; i<n; i++) {
406: VecDestroy(&bdf->work[i]);
407: VecDestroy(&bdf->tvwork[i]);
408: }
409: VecDestroy(&bdf->vec_dot);
410: VecDestroy(&bdf->vec_wrk);
411: VecDestroy(&bdf->vec_lte);
412: if (ts->dm) {DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);}
413: return(0);
414: }
416: static PetscErrorCode TSDestroy_BDF(TS ts)
417: {
421: TSReset_BDF(ts);
422: PetscFree(ts->data);
423: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);
424: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);
425: return(0);
426: }
428: static PetscErrorCode TSSetUp_BDF(TS ts)
429: {
430: TS_BDF *bdf = (TS_BDF*)ts->data;
431: size_t i,n = sizeof(bdf->work)/sizeof(Vec);
432: PetscReal low,high,two = 2;
436: TSHasTransientVariable(ts,&bdf->transientvar);
437: bdf->k = bdf->n = 0;
438: for (i=0; i<n; i++) {
439: VecDuplicate(ts->vec_sol,&bdf->work[i]);
440: if (i && bdf->transientvar) {
441: VecDuplicate(ts->vec_sol,&bdf->tvwork[i]);
442: }
443: }
444: VecDuplicate(ts->vec_sol,&bdf->vec_dot);
445: VecDuplicate(ts->vec_sol,&bdf->vec_wrk);
446: VecDuplicate(ts->vec_sol,&bdf->vec_lte);
447: TSGetDM(ts,&ts->dm);
448: DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts);
450: TSGetAdapt(ts,&ts->adapt);
451: TSAdaptCandidatesClear(ts->adapt);
452: TSAdaptGetClip(ts->adapt,&low,&high);
453: TSAdaptSetClip(ts->adapt,low,PetscMin(high,two));
455: TSGetSNES(ts,&ts->snes);
456: return(0);
457: }
459: static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts)
460: {
464: PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");
465: {
466: PetscBool flg;
467: PetscInt order;
468: TSBDFGetOrder(ts,&order);
469: PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);
470: if (flg) {TSBDFSetOrder(ts,order);}
471: }
472: PetscOptionsTail();
473: return(0);
474: }
476: static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer)
477: {
478: TS_BDF *bdf = (TS_BDF*)ts->data;
479: PetscBool iascii;
483: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
484: if (iascii) {
485: PetscViewerASCIIPrintf(viewer," Order=%D\n",bdf->order);
486: }
487: return(0);
488: }
490: /* ------------------------------------------------------------ */
492: static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order)
493: {
494: TS_BDF *bdf = (TS_BDF*)ts->data;
497: if (order == bdf->order) return(0);
498: if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order);
499: bdf->order = order;
500: return(0);
501: }
503: static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order)
504: {
505: TS_BDF *bdf = (TS_BDF*)ts->data;
508: *order = bdf->order;
509: return(0);
510: }
512: /* ------------------------------------------------------------ */
514: /*MC
515: TSBDF - DAE solver using BDF methods
517: Level: beginner
519: .seealso: TS, TSCreate(), TSSetType()
520: M*/
521: PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts)
522: {
523: TS_BDF *bdf;
527: ts->ops->reset = TSReset_BDF;
528: ts->ops->destroy = TSDestroy_BDF;
529: ts->ops->view = TSView_BDF;
530: ts->ops->setup = TSSetUp_BDF;
531: ts->ops->setfromoptions = TSSetFromOptions_BDF;
532: ts->ops->step = TSStep_BDF;
533: ts->ops->evaluatewlte = TSEvaluateWLTE_BDF;
534: ts->ops->rollback = TSRollBack_BDF;
535: ts->ops->interpolate = TSInterpolate_BDF;
536: ts->ops->snesfunction = SNESTSFormFunction_BDF;
537: ts->ops->snesjacobian = SNESTSFormJacobian_BDF;
538: ts->default_adapt_type = TSADAPTBASIC;
540: ts->usessnes = PETSC_TRUE;
542: PetscNewLog(ts,&bdf);
543: ts->data = (void*)bdf;
545: bdf->status = TS_STEP_COMPLETE;
547: PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);
548: PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);
549: TSBDFSetOrder(ts,2);
550: return(0);
551: }
553: /* ------------------------------------------------------------ */
555: /*@
556: TSBDFSetOrder - Set the order of the BDF method
558: Logically Collective on TS
560: Input Parameter:
561: + ts - timestepping context
562: - order - order of the method
564: Options Database:
565: . -ts_bdf_order <order>
567: Level: intermediate
569: @*/
570: PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order)
571: {
577: PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));
578: return(0);
579: }
581: /*@
582: TSBDFGetOrder - Get the order of the BDF method
584: Not Collective
586: Input Parameter:
587: . ts - timestepping context
589: Output Parameter:
590: . order - order of the method
592: Level: intermediate
594: @*/
595: PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order)
596: {
602: PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));
603: return(0);
604: }