xref: /petsc/src/ts/impls/implicit/alpha/alpha2.c (revision 4e8208cbcbc709572b8abe32f33c78b69c819375) 
1 /*
2   Code for timestepping with implicit generalized-\alpha method
3   for second order systems.
4 */
5 #include <petsc/private/tsimpl.h> /*I   "petscts.h"   I*/
6 
7 static PetscBool  cited      = PETSC_FALSE;
8 static const char citation[] = "@article{Chung1993,\n"
9                                "  title   = {A Time Integration Algorithm for Structural Dynamics with Improved Numerical Dissipation: The Generalized-$\\alpha$ Method},\n"
10                                "  author  = {J. Chung, G. M. Hubert},\n"
11                                "  journal = {ASME Journal of Applied Mechanics},\n"
12                                "  volume  = {60},\n"
13                                "  number  = {2},\n"
14                                "  pages   = {371--375},\n"
15                                "  year    = {1993},\n"
16                                "  issn    = {0021-8936},\n"
17                                "  doi     = {http://dx.doi.org/10.1115/1.2900803}\n}\n";
18 
19 typedef struct {
20   PetscReal stage_time;
21   PetscReal shift_V;
22   PetscReal shift_A;
23   PetscReal scale_F;
24   Vec       X0, Xa, X1;
25   Vec       V0, Va, V1;
26   Vec       A0, Aa, A1;
27 
28   Vec vec_dot;
29 
30   PetscReal Alpha_m;
31   PetscReal Alpha_f;
32   PetscReal Gamma;
33   PetscReal Beta;
34   PetscInt  order;
35 
36   Vec vec_sol_prev;
37   Vec vec_dot_prev;
38   Vec vec_lte_work[2];
39 
40   TSStepStatus status;
41 
42   TSAlpha2PredictorFn *predictor;
43   void                *predictor_ctx;
44 } TS_Alpha;
45 
46 /*@C
47   TSAlpha2SetPredictor - sets the callback for computing a predictor (i.e., initial guess
48   for the nonlinear solver).
49 
50   Input Parameters:
51 + ts        - timestepping context
52 . predictor - callback to set the predictor in each step
53 - ctx       - the application context, which may be set to `NULL` if not used
54 
55   Level: intermediate
56 
57   Notes:
58 
59   If this function is never called, a same-state-vector predictor will be used, i.e.,
60   the initial guess will be the converged solution from the previous time step, without regard
61   for the previous velocity or acceleration.
62 
63 .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2PredictorFn`
64 @*/
TSAlpha2SetPredictor(TS ts,TSAlpha2PredictorFn * predictor,PetscCtx ctx)65 PetscErrorCode TSAlpha2SetPredictor(TS ts, TSAlpha2PredictorFn *predictor, PetscCtx ctx)
66 {
67   TS_Alpha *th = (TS_Alpha *)ts->data;
68 
69   PetscFunctionBegin;
70   th->predictor     = predictor;
71   th->predictor_ctx = ctx;
72   PetscFunctionReturn(PETSC_SUCCESS);
73 }
74 
TSAlpha_ApplyPredictor(TS ts,Vec X1)75 static PetscErrorCode TSAlpha_ApplyPredictor(TS ts, Vec X1)
76 {
77   /* Apply a custom predictor if set, or default to same-displacement. */
78   TS_Alpha *th = (TS_Alpha *)ts->data;
79 
80   PetscFunctionBegin;
81   if (th->predictor) PetscCall(th->predictor(ts, th->X0, th->V0, th->A0, X1, th->predictor_ctx));
82   else PetscCall(VecCopy(th->X0, X1));
83   PetscFunctionReturn(PETSC_SUCCESS);
84 }
85 
TSAlpha_StageTime(TS ts)86 static PetscErrorCode TSAlpha_StageTime(TS ts)
87 {
88   TS_Alpha *th      = (TS_Alpha *)ts->data;
89   PetscReal t       = ts->ptime;
90   PetscReal dt      = ts->time_step;
91   PetscReal Alpha_m = th->Alpha_m;
92   PetscReal Alpha_f = th->Alpha_f;
93   PetscReal Gamma   = th->Gamma;
94   PetscReal Beta    = th->Beta;
95 
96   PetscFunctionBegin;
97   th->stage_time = t + Alpha_f * dt;
98   th->shift_V    = Gamma / (dt * Beta);
99   th->shift_A    = Alpha_m / (Alpha_f * dt * dt * Beta);
100   th->scale_F    = 1 / Alpha_f;
101   PetscFunctionReturn(PETSC_SUCCESS);
102 }
103 
TSAlpha_StageVecs(TS ts,Vec X)104 static PetscErrorCode TSAlpha_StageVecs(TS ts, Vec X)
105 {
106   TS_Alpha *th = (TS_Alpha *)ts->data;
107   Vec       X1 = X, V1 = th->V1, A1 = th->A1;
108   Vec       Xa = th->Xa, Va = th->Va, Aa = th->Aa;
109   Vec       X0 = th->X0, V0 = th->V0, A0 = th->A0;
110   PetscReal dt      = ts->time_step;
111   PetscReal Alpha_m = th->Alpha_m;
112   PetscReal Alpha_f = th->Alpha_f;
113   PetscReal Gamma   = th->Gamma;
114   PetscReal Beta    = th->Beta;
115 
116   PetscFunctionBegin;
117   /* A1 = ... */
118   PetscCall(VecWAXPY(A1, -1.0, X0, X1));
119   PetscCall(VecAXPY(A1, -dt, V0));
120   PetscCall(VecAXPBY(A1, -(1 - 2 * Beta) / (2 * Beta), 1 / (dt * dt * Beta), A0));
121   /* V1 = ... */
122   PetscCall(VecWAXPY(V1, (1.0 - Gamma) / Gamma, A0, A1));
123   PetscCall(VecAYPX(V1, dt * Gamma, V0));
124   /* Xa = X0 + Alpha_f*(X1-X0) */
125   PetscCall(VecWAXPY(Xa, -1.0, X0, X1));
126   PetscCall(VecAYPX(Xa, Alpha_f, X0));
127   /* Va = V0 + Alpha_f*(V1-V0) */
128   PetscCall(VecWAXPY(Va, -1.0, V0, V1));
129   PetscCall(VecAYPX(Va, Alpha_f, V0));
130   /* Aa = A0 + Alpha_m*(A1-A0) */
131   PetscCall(VecWAXPY(Aa, -1.0, A0, A1));
132   PetscCall(VecAYPX(Aa, Alpha_m, A0));
133   PetscFunctionReturn(PETSC_SUCCESS);
134 }
135 
TSAlpha_SNESSolve(TS ts,Vec b,Vec x)136 static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
137 {
138   PetscInt nits, lits;
139 
140   PetscFunctionBegin;
141   PetscCall(SNESSolve(ts->snes, b, x));
142   PetscCall(SNESGetIterationNumber(ts->snes, &nits));
143   PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
144   ts->snes_its += nits;
145   ts->ksp_its += lits;
146   PetscFunctionReturn(PETSC_SUCCESS);
147 }
148 
149 /*
150   Compute a consistent initial state for the generalized-alpha method.
151   - Solve two successive first-order-accurate steps with halved time step.
152   - Compute the initial second time derivative using backward differences.
153   - If using adaptivity, estimate the LTE of the initial step.
154 */
TSAlpha_Restart(TS ts,PetscBool * initok)155 static PetscErrorCode TSAlpha_Restart(TS ts, PetscBool *initok)
156 {
157   TS_Alpha *th = (TS_Alpha *)ts->data;
158   PetscReal time_step;
159   PetscReal alpha_m, alpha_f, gamma, beta;
160   Vec       X0 = ts->vec_sol, X1, X2 = th->X1;
161   Vec       V0 = ts->vec_dot, V1, V2 = th->V1;
162   PetscBool stageok;
163 
164   PetscFunctionBegin;
165   PetscCall(VecDuplicate(X0, &X1));
166   PetscCall(VecDuplicate(V0, &V1));
167 
168   /* Setup first-order-accurate method with halved time step */
169   PetscCall(TSAlpha2GetParams(ts, &alpha_m, &alpha_f, &gamma, &beta));
170   PetscCall(TSAlpha2SetParams(ts, 1, 1, 1, 0.5));
171   PetscCall(TSGetTimeStep(ts, &time_step));
172   ts->time_step = time_step / 2;
173   PetscCall(TSAlpha_StageTime(ts));
174   th->stage_time = ts->ptime;
175   PetscCall(VecZeroEntries(th->A0));
176 
177   /* First half step, (t0,X0,V0) -> (t1,X1,V1) */
178   th->stage_time += ts->time_step;
179   PetscCall(VecCopy(X0, th->X0));
180   PetscCall(VecCopy(V0, th->V0));
181   PetscCall(TSPreStage(ts, th->stage_time));
182   PetscCall(TSAlpha_ApplyPredictor(ts, X1));
183   PetscCall(TSAlpha_SNESSolve(ts, NULL, X1));
184   PetscCall(VecCopy(th->V1, V1));
185   PetscCall(TSPostStage(ts, th->stage_time, 0, &X1));
186   PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
187   if (!stageok) goto finally;
188 
189   /* Second half step, (t1,X1,V1) -> (t2,X2,V2) */
190   th->stage_time += ts->time_step;
191   PetscCall(VecCopy(X1, th->X0));
192   PetscCall(VecCopy(V1, th->V0));
193   PetscCall(TSPreStage(ts, th->stage_time));
194   PetscCall(TSAlpha_ApplyPredictor(ts, X2));
195   PetscCall(TSAlpha_SNESSolve(ts, NULL, X2));
196   PetscCall(VecCopy(th->V1, V2));
197   PetscCall(TSPostStage(ts, th->stage_time, 0, &X2));
198   PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X2, &stageok));
199   if (!stageok) goto finally;
200 
201   /* Compute A0 ~ dV/dt at t0 with backward differences */
202   PetscCall(VecZeroEntries(th->A0));
203   PetscCall(VecAXPY(th->A0, -3 / time_step, V0));
204   PetscCall(VecAXPY(th->A0, +4 / time_step, V1));
205   PetscCall(VecAXPY(th->A0, -1 / time_step, V2));
206 
207   /* Rough, lower-order estimate LTE of the initial step */
208   if (th->vec_lte_work[0]) {
209     PetscCall(VecZeroEntries(th->vec_lte_work[0]));
210     PetscCall(VecAXPY(th->vec_lte_work[0], +2, X2));
211     PetscCall(VecAXPY(th->vec_lte_work[0], -4, X1));
212     PetscCall(VecAXPY(th->vec_lte_work[0], +2, X0));
213   }
214   if (th->vec_lte_work[1]) {
215     PetscCall(VecZeroEntries(th->vec_lte_work[1]));
216     PetscCall(VecAXPY(th->vec_lte_work[1], +2, V2));
217     PetscCall(VecAXPY(th->vec_lte_work[1], -4, V1));
218     PetscCall(VecAXPY(th->vec_lte_work[1], +2, V0));
219   }
220 
221 finally:
222   /* Revert TSAlpha to the initial state (t0,X0,V0), but retain
223      potential time step reduction by factor after failed solve. */
224   if (initok) *initok = stageok;
225   PetscCall(TSSetTimeStep(ts, 2 * ts->time_step));
226   PetscCall(TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta));
227   PetscCall(VecCopy(ts->vec_sol, th->X0));
228   PetscCall(VecCopy(ts->vec_dot, th->V0));
229 
230   PetscCall(VecDestroy(&X1));
231   PetscCall(VecDestroy(&V1));
232   PetscFunctionReturn(PETSC_SUCCESS);
233 }
234 
TSStep_Alpha(TS ts)235 static PetscErrorCode TSStep_Alpha(TS ts)
236 {
237   TS_Alpha *th         = (TS_Alpha *)ts->data;
238   PetscInt  rejections = 0;
239   PetscBool stageok, accept = PETSC_TRUE;
240   PetscReal next_time_step = ts->time_step;
241 
242   PetscFunctionBegin;
243   PetscCall(PetscCitationsRegister(citation, &cited));
244 
245   if (!ts->steprollback) {
246     if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
247     if (th->vec_dot_prev) PetscCall(VecCopy(th->V0, th->vec_dot_prev));
248     PetscCall(VecCopy(ts->vec_sol, th->X0));
249     PetscCall(VecCopy(ts->vec_dot, th->V0));
250     PetscCall(VecCopy(th->A1, th->A0));
251   }
252 
253   th->status = TS_STEP_INCOMPLETE;
254   while (!ts->reason && th->status != TS_STEP_COMPLETE) {
255     if (ts->steprestart) {
256       PetscCall(TSAlpha_Restart(ts, &stageok));
257       if (!stageok) goto reject_step;
258     }
259 
260     PetscCall(TSAlpha_StageTime(ts));
261     PetscCall(TSAlpha_ApplyPredictor(ts, th->X1));
262     PetscCall(TSPreStage(ts, th->stage_time));
263     PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
264     PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
265     PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
266     if (!stageok) goto reject_step;
267 
268     th->status = TS_STEP_PENDING;
269     PetscCall(VecCopy(th->X1, ts->vec_sol));
270     PetscCall(VecCopy(th->V1, ts->vec_dot));
271     PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
272     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
273     if (!accept) {
274       PetscCall(VecCopy(th->X0, ts->vec_sol));
275       PetscCall(VecCopy(th->V0, ts->vec_dot));
276       ts->time_step = next_time_step;
277       goto reject_step;
278     }
279 
280     ts->ptime += ts->time_step;
281     ts->time_step = next_time_step;
282     break;
283 
284   reject_step:
285     ts->reject++;
286     accept = PETSC_FALSE;
287     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
288       ts->reason = TS_DIVERGED_STEP_REJECTED;
289       PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
290     }
291   }
292   PetscFunctionReturn(PETSC_SUCCESS);
293 }
294 
TSEvaluateWLTE_Alpha(TS ts,NormType wnormtype,PetscInt * order,PetscReal * wlte)295 static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
296 {
297   TS_Alpha *th = (TS_Alpha *)ts->data;
298   Vec       X  = th->X1;              /* X = solution */
299   Vec       V  = th->V1;              /* V = solution */
300   Vec       Y  = th->vec_lte_work[0]; /* Y = X + LTE  */
301   Vec       Z  = th->vec_lte_work[1]; /* Z = V + LTE  */
302   PetscReal enormX, enormV, enormXa, enormVa, enormXr, enormVr;
303 
304   PetscFunctionBegin;
305   if (!th->vec_sol_prev) {
306     *wlte = -1;
307     PetscFunctionReturn(PETSC_SUCCESS);
308   }
309   if (!th->vec_dot_prev) {
310     *wlte = -1;
311     PetscFunctionReturn(PETSC_SUCCESS);
312   }
313   if (!th->vec_lte_work[0]) {
314     *wlte = -1;
315     PetscFunctionReturn(PETSC_SUCCESS);
316   }
317   if (!th->vec_lte_work[1]) {
318     *wlte = -1;
319     PetscFunctionReturn(PETSC_SUCCESS);
320   }
321   if (ts->steprestart) {
322     /* th->vec_lte_prev is set to the LTE in TSAlpha_Restart() */
323     PetscCall(VecAXPY(Y, 1, X));
324     PetscCall(VecAXPY(Z, 1, V));
325   } else {
326     /* Compute LTE using backward differences with non-constant time step */
327     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
328     PetscReal   a = 1 + h_prev / h;
329     PetscScalar scal[3];
330     Vec         vecX[3], vecV[3];
331     scal[0] = +1 / a;
332     scal[1] = -1 / (a - 1);
333     scal[2] = +1 / (a * (a - 1));
334     vecX[0] = th->X1;
335     vecX[1] = th->X0;
336     vecX[2] = th->vec_sol_prev;
337     vecV[0] = th->V1;
338     vecV[1] = th->V0;
339     vecV[2] = th->vec_dot_prev;
340     PetscCall(VecCopy(X, Y));
341     PetscCall(VecMAXPY(Y, 3, scal, vecX));
342     PetscCall(VecCopy(V, Z));
343     PetscCall(VecMAXPY(Z, 3, scal, vecV));
344   }
345   /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
346   PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, &enormX, &enormXa, &enormXr));
347   PetscCall(TSErrorWeightedNorm(ts, V, Z, wnormtype, &enormV, &enormVa, &enormVr));
348   if (wnormtype == NORM_2) *wlte = PetscSqrtReal(PetscSqr(enormX) / 2 + PetscSqr(enormV) / 2);
349   else *wlte = PetscMax(enormX, enormV);
350   if (order) *order = 2;
351   PetscFunctionReturn(PETSC_SUCCESS);
352 }
353 
TSRollBack_Alpha(TS ts)354 static PetscErrorCode TSRollBack_Alpha(TS ts)
355 {
356   TS_Alpha *th = (TS_Alpha *)ts->data;
357 
358   PetscFunctionBegin;
359   PetscCall(VecCopy(th->X0, ts->vec_sol));
360   PetscCall(VecCopy(th->V0, ts->vec_dot));
361   PetscFunctionReturn(PETSC_SUCCESS);
362 }
363 
364 /*
365 static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X,Vec V)
366 {
367   TS_Alpha       *th = (TS_Alpha*)ts->data;
368   PetscReal      dt  = t - ts->ptime;
369 
370   PetscFunctionBegin;
371   PetscCall(VecCopy(ts->vec_dot,V));
372   PetscCall(VecAXPY(V,dt*(1-th->Gamma),th->A0));
373   PetscCall(VecAXPY(V,dt*th->Gamma,th->A1));
374   PetscCall(VecCopy(ts->vec_sol,X));
375   PetscCall(VecAXPY(X,dt,V));
376   PetscCall(VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0));
377   PetscCall(VecAXPY(X,dt*dt*th->Beta,th->A1));
378   PetscFunctionReturn(PETSC_SUCCESS);
379 }
380 */
381 
SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)382 static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
383 {
384   TS_Alpha *th = (TS_Alpha *)ts->data;
385   PetscReal ta = th->stage_time;
386   Vec       Xa = th->Xa, Va = th->Va, Aa = th->Aa;
387 
388   PetscFunctionBegin;
389   PetscCall(TSAlpha_StageVecs(ts, X));
390   /* F = Function(ta,Xa,Va,Aa) */
391   PetscCall(TSComputeI2Function(ts, ta, Xa, Va, Aa, F));
392   PetscCall(VecScale(F, th->scale_F));
393   PetscFunctionReturn(PETSC_SUCCESS);
394 }
395 
SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,TS ts)396 static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
397 {
398   TS_Alpha *th = (TS_Alpha *)ts->data;
399   PetscReal ta = th->stage_time;
400   Vec       Xa = th->Xa, Va = th->Va, Aa = th->Aa;
401   PetscReal dVdX = th->shift_V, dAdX = th->shift_A;
402 
403   PetscFunctionBegin;
404   /* J,P = Jacobian(ta,Xa,Va,Aa) */
405   PetscCall(TSComputeI2Jacobian(ts, ta, Xa, Va, Aa, dVdX, dAdX, J, P));
406   PetscFunctionReturn(PETSC_SUCCESS);
407 }
408 
TSReset_Alpha(TS ts)409 static PetscErrorCode TSReset_Alpha(TS ts)
410 {
411   TS_Alpha *th = (TS_Alpha *)ts->data;
412 
413   PetscFunctionBegin;
414   PetscCall(VecDestroy(&th->X0));
415   PetscCall(VecDestroy(&th->Xa));
416   PetscCall(VecDestroy(&th->X1));
417   PetscCall(VecDestroy(&th->V0));
418   PetscCall(VecDestroy(&th->Va));
419   PetscCall(VecDestroy(&th->V1));
420   PetscCall(VecDestroy(&th->A0));
421   PetscCall(VecDestroy(&th->Aa));
422   PetscCall(VecDestroy(&th->A1));
423   PetscCall(VecDestroy(&th->vec_sol_prev));
424   PetscCall(VecDestroy(&th->vec_dot_prev));
425   PetscCall(VecDestroy(&th->vec_lte_work[0]));
426   PetscCall(VecDestroy(&th->vec_lte_work[1]));
427   PetscFunctionReturn(PETSC_SUCCESS);
428 }
429 
TSDestroy_Alpha(TS ts)430 static PetscErrorCode TSDestroy_Alpha(TS ts)
431 {
432   PetscFunctionBegin;
433   PetscCall(TSReset_Alpha(ts));
434   PetscCall(PetscFree(ts->data));
435 
436   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", NULL));
437   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", NULL));
438   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", NULL));
439   PetscFunctionReturn(PETSC_SUCCESS);
440 }
441 
TSSetUp_Alpha(TS ts)442 static PetscErrorCode TSSetUp_Alpha(TS ts)
443 {
444   TS_Alpha *th = (TS_Alpha *)ts->data;
445   PetscBool match;
446 
447   PetscFunctionBegin;
448   PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
449   PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
450   PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
451   PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
452   PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
453   PetscCall(VecDuplicate(ts->vec_sol, &th->V1));
454   PetscCall(VecDuplicate(ts->vec_sol, &th->A0));
455   PetscCall(VecDuplicate(ts->vec_sol, &th->Aa));
456   PetscCall(VecDuplicate(ts->vec_sol, &th->A1));
457 
458   PetscCall(TSGetAdapt(ts, &ts->adapt));
459   PetscCall(TSAdaptCandidatesClear(ts->adapt));
460   PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
461   if (!match) {
462     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
463     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_dot_prev));
464     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work[0]));
465     PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work[1]));
466   }
467 
468   PetscCall(TSGetSNES(ts, &ts->snes));
469   PetscFunctionReturn(PETSC_SUCCESS);
470 }
471 
TSSetFromOptions_Alpha(TS ts,PetscOptionItems PetscOptionsObject)472 static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems PetscOptionsObject)
473 {
474   TS_Alpha *th = (TS_Alpha *)ts->data;
475 
476   PetscFunctionBegin;
477   PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
478   {
479     PetscBool flg;
480     PetscReal radius = 1;
481     PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlpha2SetRadius", radius, &radius, &flg));
482     if (flg) PetscCall(TSAlpha2SetRadius(ts, radius));
483     PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlpha2SetParams", th->Alpha_m, &th->Alpha_m, NULL));
484     PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlpha2SetParams", th->Alpha_f, &th->Alpha_f, NULL));
485     PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlpha2SetParams", th->Gamma, &th->Gamma, NULL));
486     PetscCall(PetscOptionsReal("-ts_alpha_beta", "Algorithmic parameter beta", "TSAlpha2SetParams", th->Beta, &th->Beta, NULL));
487     PetscCall(TSAlpha2SetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma, th->Beta));
488   }
489   PetscOptionsHeadEnd();
490   PetscFunctionReturn(PETSC_SUCCESS);
491 }
492 
TSView_Alpha(TS ts,PetscViewer viewer)493 static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
494 {
495   TS_Alpha *th = (TS_Alpha *)ts->data;
496   PetscBool isascii;
497 
498   PetscFunctionBegin;
499   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
500   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "  Alpha_m=%g, Alpha_f=%g, Gamma=%g, Beta=%g\n", (double)th->Alpha_m, (double)th->Alpha_f, (double)th->Gamma, (double)th->Beta));
501   PetscFunctionReturn(PETSC_SUCCESS);
502 }
503 
TSAlpha2SetRadius_Alpha(TS ts,PetscReal radius)504 static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts, PetscReal radius)
505 {
506   PetscReal alpha_m, alpha_f, gamma, beta;
507 
508   PetscFunctionBegin;
509   PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
510   alpha_m = (2 - radius) / (1 + radius);
511   alpha_f = 1 / (1 + radius);
512   gamma   = (PetscReal)0.5 + alpha_m - alpha_f;
513   beta    = (PetscReal)0.5 * (1 + alpha_m - alpha_f);
514   beta *= beta;
515   PetscCall(TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta));
516   PetscFunctionReturn(PETSC_SUCCESS);
517 }
518 
TSAlpha2SetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)519 static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
520 {
521   TS_Alpha *th  = (TS_Alpha *)ts->data;
522   PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
523   PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
524 
525   PetscFunctionBegin;
526   th->Alpha_m = alpha_m;
527   th->Alpha_f = alpha_f;
528   th->Gamma   = gamma;
529   th->Beta    = beta;
530   th->order   = (PetscAbsReal(res) < tol) ? 2 : 1;
531   PetscFunctionReturn(PETSC_SUCCESS);
532 }
533 
TSAlpha2GetParams_Alpha(TS ts,PetscReal * alpha_m,PetscReal * alpha_f,PetscReal * gamma,PetscReal * beta)534 static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
535 {
536   TS_Alpha *th = (TS_Alpha *)ts->data;
537 
538   PetscFunctionBegin;
539   if (alpha_m) *alpha_m = th->Alpha_m;
540   if (alpha_f) *alpha_f = th->Alpha_f;
541   if (gamma) *gamma = th->Gamma;
542   if (beta) *beta = th->Beta;
543   PetscFunctionReturn(PETSC_SUCCESS);
544 }
545 
546 /*MC
547   TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method for second-order systems {cite}`chung1993`
548 
549   Level: beginner
550 
551 .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
552 M*/
TSCreate_Alpha2(TS ts)553 PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
554 {
555   TS_Alpha *th;
556 
557   PetscFunctionBegin;
558   ts->ops->reset          = TSReset_Alpha;
559   ts->ops->destroy        = TSDestroy_Alpha;
560   ts->ops->view           = TSView_Alpha;
561   ts->ops->setup          = TSSetUp_Alpha;
562   ts->ops->setfromoptions = TSSetFromOptions_Alpha;
563   ts->ops->step           = TSStep_Alpha;
564   ts->ops->evaluatewlte   = TSEvaluateWLTE_Alpha;
565   ts->ops->rollback       = TSRollBack_Alpha;
566   /*ts->ops->interpolate  = TSInterpolate_Alpha;*/
567   ts->ops->snesfunction  = SNESTSFormFunction_Alpha;
568   ts->ops->snesjacobian  = SNESTSFormJacobian_Alpha;
569   ts->default_adapt_type = TSADAPTNONE;
570 
571   ts->usessnes = PETSC_TRUE;
572 
573   PetscCall(PetscNew(&th));
574   ts->data = (void *)th;
575 
576   th->Alpha_m = 0.5;
577   th->Alpha_f = 0.5;
578   th->Gamma   = 0.5;
579   th->Beta    = 0.25;
580   th->order   = 2;
581 
582   th->predictor     = NULL;
583   th->predictor_ctx = NULL;
584 
585   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", TSAlpha2SetRadius_Alpha));
586   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", TSAlpha2SetParams_Alpha));
587   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", TSAlpha2GetParams_Alpha));
588   PetscFunctionReturn(PETSC_SUCCESS);
589 }
590 
591 /*@
592   TSAlpha2SetRadius - sets the desired spectral radius of the method for `TSALPHA2`
593   (i.e. high-frequency numerical damping)
594 
595   Logically Collective
596 
597   Input Parameters:
598 + ts     - timestepping context
599 - radius - the desired spectral radius
600 
601   Options Database Key:
602 . -ts_alpha_radius <radius> - set the desired spectral radius
603 
604   Level: intermediate
605 
606   Notes:
607 
608   The algorithmic parameters $\alpha_m$ and $\alpha_f$ of the generalized-$\alpha$ method can
609   be computed in terms of a specified spectral radius $\rho$ in `[0, 1]` for infinite time step
610   in order to control high-frequency numerical damping\:
611 
612   $$
613   \begin{align*}
614   \alpha_m = (2-\rho)/(1+\rho) \\
615   \alpha_f = 1/(1+\rho)
616   \end{align*}
617   $$
618 
619 .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetParams()`, `TSAlpha2GetParams()`
620 @*/
TSAlpha2SetRadius(TS ts,PetscReal radius)621 PetscErrorCode TSAlpha2SetRadius(TS ts, PetscReal radius)
622 {
623   PetscFunctionBegin;
624   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
625   PetscValidLogicalCollectiveReal(ts, radius, 2);
626   PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
627   PetscTryMethod(ts, "TSAlpha2SetRadius_C", (TS, PetscReal), (ts, radius));
628   PetscFunctionReturn(PETSC_SUCCESS);
629 }
630 
631 /*@
632   TSAlpha2SetParams - sets the algorithmic parameters for `TSALPHA2`
633 
634   Logically Collective
635 
636   Input Parameters:
637 + ts      - timestepping context
638 . alpha_m - algorithmic parameter
639 . alpha_f - algorithmic parameter
640 . gamma   - algorithmic parameter
641 - beta    - algorithmic parameter
642 
643   Options Database Keys:
644 + -ts_alpha_alpha_m <alpha_m> - set alpha_m
645 . -ts_alpha_alpha_f <alpha_f> - set alpha_f
646 . -ts_alpha_gamma   <gamma>   - set gamma
647 - -ts_alpha_beta    <beta>    - set beta
648 
649   Level: advanced
650 
651   Notes:
652   Second-order accuracy can be obtained so long as\:
653 
654   $$
655   \begin{align*}
656   \gamma = 1/2 + \alpha_m - \alpha_f \\
657   \beta  = 1/4 (1 + \alpha_m - \alpha_f)^2.
658   \end{align*}
659   $$
660 
661   Unconditional stability requires\:
662   $$
663   \alpha_m >= \alpha_f >= 1/2.
664   $$
665 
666   Use of this function is normally only required to hack `TSALPHA2` to use a modified
667   integration scheme. Users should call `TSAlpha2SetRadius()` to set the desired spectral
668   radius of the methods (i.e. high-frequency damping) in order so select optimal values for
669   these parameters.
670 
671 .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2GetParams()`
672 @*/
TSAlpha2SetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma,PetscReal beta)673 PetscErrorCode TSAlpha2SetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
674 {
675   PetscFunctionBegin;
676   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
677   PetscValidLogicalCollectiveReal(ts, alpha_m, 2);
678   PetscValidLogicalCollectiveReal(ts, alpha_f, 3);
679   PetscValidLogicalCollectiveReal(ts, gamma, 4);
680   PetscValidLogicalCollectiveReal(ts, beta, 5);
681   PetscTryMethod(ts, "TSAlpha2SetParams_C", (TS, PetscReal, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma, beta));
682   PetscFunctionReturn(PETSC_SUCCESS);
683 }
684 
685 /*@
686   TSAlpha2GetParams - gets the algorithmic parameters for `TSALPHA2`
687 
688   Not Collective
689 
690   Input Parameter:
691 . ts - timestepping context
692 
693   Output Parameters:
694 + alpha_m - algorithmic parameter
695 . alpha_f - algorithmic parameter
696 . gamma   - algorithmic parameter
697 - beta    - algorithmic parameter
698 
699   Level: advanced
700 
701   Note:
702   Use of this function is normally only required to hack `TSALPHA2` to use a modified
703   integration scheme. Users should call `TSAlpha2SetRadius()` to set the high-frequency damping
704   (i.e. spectral radius of the method) in order so select optimal values for these parameters.
705 
706 .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
707 @*/
TSAlpha2GetParams(TS ts,PetscReal * alpha_m,PetscReal * alpha_f,PetscReal * gamma,PetscReal * beta)708 PetscErrorCode TSAlpha2GetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
709 {
710   PetscFunctionBegin;
711   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
712   if (alpha_m) PetscAssertPointer(alpha_m, 2);
713   if (alpha_f) PetscAssertPointer(alpha_f, 3);
714   if (gamma) PetscAssertPointer(gamma, 4);
715   if (beta) PetscAssertPointer(beta, 5);
716   PetscUseMethod(ts, "TSAlpha2GetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma, beta));
717   PetscFunctionReturn(PETSC_SUCCESS);
718 }
719