Actual source code: alpha2.c
1: /*
2: Code for timestepping with implicit generalized-\alpha method
3: for second order systems.
4: */
5: #include <petsc/private/tsimpl.h>
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";
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;
28: Vec vec_dot;
30: PetscReal Alpha_m;
31: PetscReal Alpha_f;
32: PetscReal Gamma;
33: PetscReal Beta;
34: PetscInt order;
36: Vec vec_sol_prev;
37: Vec vec_dot_prev;
38: Vec vec_lte_work[2];
40: TSStepStatus status;
42: TSAlpha2PredictorFn *predictor;
43: void *predictor_ctx;
44: } TS_Alpha;
46: /*@C
47: TSAlpha2SetPredictor - sets the callback for computing a predictor (i.e., initial guess
48: for the nonlinear solver).
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
55: Level: intermediate
57: Notes:
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.
63: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2PredictorFn`
64: @*/
65: PetscErrorCode TSAlpha2SetPredictor(TS ts, TSAlpha2PredictorFn *predictor, void *ctx)
66: {
67: TS_Alpha *th = (TS_Alpha *)ts->data;
69: PetscFunctionBegin;
70: th->predictor = predictor;
71: th->predictor_ctx = ctx;
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
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;
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: }
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;
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: }
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;
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: }
136: static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
137: {
138: PetscInt nits, lits;
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: }
149: /*
150: Compute a consistent initial state for the generalized-alpha method.
151: - Solve two successive backward Euler 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: */
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;
164: PetscFunctionBegin;
165: PetscCall(VecDuplicate(X0, &X1));
166: PetscCall(VecDuplicate(V0, &V1));
168: /* Setup backward Euler 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));
177: /* First BE 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;
189: /* Second BE 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, X1, &stageok));
199: if (!stageok) goto finally;
201: /* Compute A0 ~ dV/dt at t0 with backward differences */
202: PetscCall(VecZeroEntries(th->A0));
203: PetscCall(VecAXPY(th->A0, -3 / ts->time_step, V0));
204: PetscCall(VecAXPY(th->A0, +4 / ts->time_step, V1));
205: PetscCall(VecAXPY(th->A0, -1 / ts->time_step, V2));
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: }
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));
230: PetscCall(VecDestroy(&X1));
231: PetscCall(VecDestroy(&V1));
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
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;
242: PetscFunctionBegin;
243: PetscCall(PetscCitationsRegister(citation, &cited));
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: }
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: }
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;
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: }
280: ts->ptime += ts->time_step;
281: ts->time_step = next_time_step;
282: break;
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: }
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;
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: }
354: static PetscErrorCode TSRollBack_Alpha(TS ts)
355: {
356: TS_Alpha *th = (TS_Alpha *)ts->data;
358: PetscFunctionBegin;
359: PetscCall(VecCopy(th->X0, ts->vec_sol));
360: PetscCall(VecCopy(th->V0, ts->vec_dot));
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
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;
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: */
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;
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: }
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;
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: }
409: static PetscErrorCode TSReset_Alpha(TS ts)
410: {
411: TS_Alpha *th = (TS_Alpha *)ts->data;
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: }
430: static PetscErrorCode TSDestroy_Alpha(TS ts)
431: {
432: PetscFunctionBegin;
433: PetscCall(TSReset_Alpha(ts));
434: PetscCall(PetscFree(ts->data));
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: }
442: static PetscErrorCode TSSetUp_Alpha(TS ts)
443: {
444: TS_Alpha *th = (TS_Alpha *)ts->data;
445: PetscBool match;
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));
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: }
468: PetscCall(TSGetSNES(ts, &ts->snes));
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
473: {
474: TS_Alpha *th = (TS_Alpha *)ts->data;
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: }
493: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
494: {
495: TS_Alpha *th = (TS_Alpha *)ts->data;
496: PetscBool iascii;
498: PetscFunctionBegin;
499: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
500: if (iascii) 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: }
504: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts, PetscReal radius)
505: {
506: PetscReal alpha_m, alpha_f, gamma, beta;
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: }
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;
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: }
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;
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: }
546: /*MC
547: TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method for second-order systems {cite}`chung1993`
549: Level: beginner
551: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
552: M*/
553: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
554: {
555: TS_Alpha *th;
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;
571: ts->usessnes = PETSC_TRUE;
573: PetscCall(PetscNew(&th));
574: ts->data = (void *)th;
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;
582: th->predictor = NULL;
583: th->predictor_ctx = NULL;
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: }
591: /*@
592: TSAlpha2SetRadius - sets the desired spectral radius of the method for `TSALPHA2`
593: (i.e. high-frequency numerical damping)
595: Logically Collective
597: Input Parameters:
598: + ts - timestepping context
599: - radius - the desired spectral radius
601: Options Database Key:
602: . -ts_alpha_radius <radius> - set the desired spectral radius
604: Level: intermediate
606: Notes:
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\:
612: $$
613: \begin{align*}
614: \alpha_m = (2-\rho)/(1+\rho) \\
615: \alpha_f = 1/(1+\rho)
616: \end{align*}
617: $$
619: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetParams()`, `TSAlpha2GetParams()`
620: @*/
621: PetscErrorCode TSAlpha2SetRadius(TS ts, PetscReal radius)
622: {
623: PetscFunctionBegin;
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: }
631: /*@
632: TSAlpha2SetParams - sets the algorithmic parameters for `TSALPHA2`
634: Logically Collective
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
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
649: Level: advanced
651: Notes:
652: Second-order accuracy can be obtained so long as\:
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: $$
661: Unconditional stability requires\:
662: $$
663: \alpha_m >= \alpha_f >= 1/2.
664: $$
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.
671: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2GetParams()`
672: @*/
673: PetscErrorCode TSAlpha2SetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
674: {
675: PetscFunctionBegin;
681: PetscTryMethod(ts, "TSAlpha2SetParams_C", (TS, PetscReal, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma, beta));
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: /*@
686: TSAlpha2GetParams - gets the algorithmic parameters for `TSALPHA2`
688: Not Collective
690: Input Parameter:
691: . ts - timestepping context
693: Output Parameters:
694: + alpha_m - algorithmic parameter
695: . alpha_f - algorithmic parameter
696: . gamma - algorithmic parameter
697: - beta - algorithmic parameter
699: Level: advanced
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.
706: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
707: @*/
708: PetscErrorCode TSAlpha2GetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
709: {
710: PetscFunctionBegin;
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: }