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;
41: } TS_Alpha;
43: static PetscErrorCode TSAlpha_StageTime(TS ts)
44: {
45: TS_Alpha *th = (TS_Alpha *)ts->data;
46: PetscReal t = ts->ptime;
47: PetscReal dt = ts->time_step;
48: PetscReal Alpha_m = th->Alpha_m;
49: PetscReal Alpha_f = th->Alpha_f;
50: PetscReal Gamma = th->Gamma;
51: PetscReal Beta = th->Beta;
53: PetscFunctionBegin;
54: th->stage_time = t + Alpha_f * dt;
55: th->shift_V = Gamma / (dt * Beta);
56: th->shift_A = Alpha_m / (Alpha_f * dt * dt * Beta);
57: th->scale_F = 1 / Alpha_f;
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode TSAlpha_StageVecs(TS ts, Vec X)
62: {
63: TS_Alpha *th = (TS_Alpha *)ts->data;
64: Vec X1 = X, V1 = th->V1, A1 = th->A1;
65: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
66: Vec X0 = th->X0, V0 = th->V0, A0 = th->A0;
67: PetscReal dt = ts->time_step;
68: PetscReal Alpha_m = th->Alpha_m;
69: PetscReal Alpha_f = th->Alpha_f;
70: PetscReal Gamma = th->Gamma;
71: PetscReal Beta = th->Beta;
73: PetscFunctionBegin;
74: /* A1 = ... */
75: PetscCall(VecWAXPY(A1, -1.0, X0, X1));
76: PetscCall(VecAXPY(A1, -dt, V0));
77: PetscCall(VecAXPBY(A1, -(1 - 2 * Beta) / (2 * Beta), 1 / (dt * dt * Beta), A0));
78: /* V1 = ... */
79: PetscCall(VecWAXPY(V1, (1.0 - Gamma) / Gamma, A0, A1));
80: PetscCall(VecAYPX(V1, dt * Gamma, V0));
81: /* Xa = X0 + Alpha_f*(X1-X0) */
82: PetscCall(VecWAXPY(Xa, -1.0, X0, X1));
83: PetscCall(VecAYPX(Xa, Alpha_f, X0));
84: /* Va = V0 + Alpha_f*(V1-V0) */
85: PetscCall(VecWAXPY(Va, -1.0, V0, V1));
86: PetscCall(VecAYPX(Va, Alpha_f, V0));
87: /* Aa = A0 + Alpha_m*(A1-A0) */
88: PetscCall(VecWAXPY(Aa, -1.0, A0, A1));
89: PetscCall(VecAYPX(Aa, Alpha_m, A0));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
94: {
95: PetscInt nits, lits;
97: PetscFunctionBegin;
98: PetscCall(SNESSolve(ts->snes, b, x));
99: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
100: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
101: ts->snes_its += nits;
102: ts->ksp_its += lits;
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: /*
107: Compute a consistent initial state for the generalized-alpha method.
108: - Solve two successive backward Euler steps with halved time step.
109: - Compute the initial second time derivative using backward differences.
110: - If using adaptivity, estimate the LTE of the initial step.
111: */
112: static PetscErrorCode TSAlpha_Restart(TS ts, PetscBool *initok)
113: {
114: TS_Alpha *th = (TS_Alpha *)ts->data;
115: PetscReal time_step;
116: PetscReal alpha_m, alpha_f, gamma, beta;
117: Vec X0 = ts->vec_sol, X1, X2 = th->X1;
118: Vec V0 = ts->vec_dot, V1, V2 = th->V1;
119: PetscBool stageok;
121: PetscFunctionBegin;
122: PetscCall(VecDuplicate(X0, &X1));
123: PetscCall(VecDuplicate(V0, &V1));
125: /* Setup backward Euler with halved time step */
126: PetscCall(TSAlpha2GetParams(ts, &alpha_m, &alpha_f, &gamma, &beta));
127: PetscCall(TSAlpha2SetParams(ts, 1, 1, 1, 0.5));
128: PetscCall(TSGetTimeStep(ts, &time_step));
129: ts->time_step = time_step / 2;
130: PetscCall(TSAlpha_StageTime(ts));
131: th->stage_time = ts->ptime;
132: PetscCall(VecZeroEntries(th->A0));
134: /* First BE step, (t0,X0,V0) -> (t1,X1,V1) */
135: th->stage_time += ts->time_step;
136: PetscCall(VecCopy(X0, th->X0));
137: PetscCall(VecCopy(V0, th->V0));
138: PetscCall(TSPreStage(ts, th->stage_time));
139: PetscCall(VecCopy(th->X0, X1));
140: PetscCall(TSAlpha_SNESSolve(ts, NULL, X1));
141: PetscCall(VecCopy(th->V1, V1));
142: PetscCall(TSPostStage(ts, th->stage_time, 0, &X1));
143: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
144: if (!stageok) goto finally;
146: /* Second BE step, (t1,X1,V1) -> (t2,X2,V2) */
147: th->stage_time += ts->time_step;
148: PetscCall(VecCopy(X1, th->X0));
149: PetscCall(VecCopy(V1, th->V0));
150: PetscCall(TSPreStage(ts, th->stage_time));
151: PetscCall(VecCopy(th->X0, X2));
152: PetscCall(TSAlpha_SNESSolve(ts, NULL, X2));
153: PetscCall(VecCopy(th->V1, V2));
154: PetscCall(TSPostStage(ts, th->stage_time, 0, &X2));
155: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
156: if (!stageok) goto finally;
158: /* Compute A0 ~ dV/dt at t0 with backward differences */
159: PetscCall(VecZeroEntries(th->A0));
160: PetscCall(VecAXPY(th->A0, -3 / ts->time_step, V0));
161: PetscCall(VecAXPY(th->A0, +4 / ts->time_step, V1));
162: PetscCall(VecAXPY(th->A0, -1 / ts->time_step, V2));
164: /* Rough, lower-order estimate LTE of the initial step */
165: if (th->vec_lte_work[0]) {
166: PetscCall(VecZeroEntries(th->vec_lte_work[0]));
167: PetscCall(VecAXPY(th->vec_lte_work[0], +2, X2));
168: PetscCall(VecAXPY(th->vec_lte_work[0], -4, X1));
169: PetscCall(VecAXPY(th->vec_lte_work[0], +2, X0));
170: }
171: if (th->vec_lte_work[1]) {
172: PetscCall(VecZeroEntries(th->vec_lte_work[1]));
173: PetscCall(VecAXPY(th->vec_lte_work[1], +2, V2));
174: PetscCall(VecAXPY(th->vec_lte_work[1], -4, V1));
175: PetscCall(VecAXPY(th->vec_lte_work[1], +2, V0));
176: }
178: finally:
179: /* Revert TSAlpha to the initial state (t0,X0,V0) */
180: if (initok) *initok = stageok;
181: PetscCall(TSSetTimeStep(ts, time_step));
182: PetscCall(TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta));
183: PetscCall(VecCopy(ts->vec_sol, th->X0));
184: PetscCall(VecCopy(ts->vec_dot, th->V0));
186: PetscCall(VecDestroy(&X1));
187: PetscCall(VecDestroy(&V1));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: static PetscErrorCode TSStep_Alpha(TS ts)
192: {
193: TS_Alpha *th = (TS_Alpha *)ts->data;
194: PetscInt rejections = 0;
195: PetscBool stageok, accept = PETSC_TRUE;
196: PetscReal next_time_step = ts->time_step;
198: PetscFunctionBegin;
199: PetscCall(PetscCitationsRegister(citation, &cited));
201: if (!ts->steprollback) {
202: if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
203: if (th->vec_dot_prev) PetscCall(VecCopy(th->V0, th->vec_dot_prev));
204: PetscCall(VecCopy(ts->vec_sol, th->X0));
205: PetscCall(VecCopy(ts->vec_dot, th->V0));
206: PetscCall(VecCopy(th->A1, th->A0));
207: }
209: th->status = TS_STEP_INCOMPLETE;
210: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
211: if (ts->steprestart) {
212: PetscCall(TSAlpha_Restart(ts, &stageok));
213: if (!stageok) goto reject_step;
214: }
216: PetscCall(TSAlpha_StageTime(ts));
217: PetscCall(VecCopy(th->X0, th->X1));
218: PetscCall(TSPreStage(ts, th->stage_time));
219: PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
220: PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
221: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
222: if (!stageok) goto reject_step;
224: th->status = TS_STEP_PENDING;
225: PetscCall(VecCopy(th->X1, ts->vec_sol));
226: PetscCall(VecCopy(th->V1, ts->vec_dot));
227: PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
228: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
229: if (!accept) {
230: PetscCall(VecCopy(th->X0, ts->vec_sol));
231: PetscCall(VecCopy(th->V0, ts->vec_dot));
232: ts->time_step = next_time_step;
233: goto reject_step;
234: }
236: ts->ptime += ts->time_step;
237: ts->time_step = next_time_step;
238: break;
240: reject_step:
241: ts->reject++;
242: accept = PETSC_FALSE;
243: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
244: ts->reason = TS_DIVERGED_STEP_REJECTED;
245: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
246: }
247: }
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
252: {
253: TS_Alpha *th = (TS_Alpha *)ts->data;
254: Vec X = th->X1; /* X = solution */
255: Vec V = th->V1; /* V = solution */
256: Vec Y = th->vec_lte_work[0]; /* Y = X + LTE */
257: Vec Z = th->vec_lte_work[1]; /* Z = V + LTE */
258: PetscReal enormX, enormV, enormXa, enormVa, enormXr, enormVr;
260: PetscFunctionBegin;
261: if (!th->vec_sol_prev) {
262: *wlte = -1;
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
265: if (!th->vec_dot_prev) {
266: *wlte = -1;
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
269: if (!th->vec_lte_work[0]) {
270: *wlte = -1;
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
273: if (!th->vec_lte_work[1]) {
274: *wlte = -1;
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
277: if (ts->steprestart) {
278: /* th->vec_lte_prev is set to the LTE in TSAlpha_Restart() */
279: PetscCall(VecAXPY(Y, 1, X));
280: PetscCall(VecAXPY(Z, 1, V));
281: } else {
282: /* Compute LTE using backward differences with non-constant time step */
283: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
284: PetscReal a = 1 + h_prev / h;
285: PetscScalar scal[3];
286: Vec vecX[3], vecV[3];
287: scal[0] = +1 / a;
288: scal[1] = -1 / (a - 1);
289: scal[2] = +1 / (a * (a - 1));
290: vecX[0] = th->X1;
291: vecX[1] = th->X0;
292: vecX[2] = th->vec_sol_prev;
293: vecV[0] = th->V1;
294: vecV[1] = th->V0;
295: vecV[2] = th->vec_dot_prev;
296: PetscCall(VecCopy(X, Y));
297: PetscCall(VecMAXPY(Y, 3, scal, vecX));
298: PetscCall(VecCopy(V, Z));
299: PetscCall(VecMAXPY(Z, 3, scal, vecV));
300: }
301: /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
302: PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, &enormX, &enormXa, &enormXr));
303: PetscCall(TSErrorWeightedNorm(ts, V, Z, wnormtype, &enormV, &enormVa, &enormVr));
304: if (wnormtype == NORM_2) *wlte = PetscSqrtReal(PetscSqr(enormX) / 2 + PetscSqr(enormV) / 2);
305: else *wlte = PetscMax(enormX, enormV);
306: if (order) *order = 2;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode TSRollBack_Alpha(TS ts)
311: {
312: TS_Alpha *th = (TS_Alpha *)ts->data;
314: PetscFunctionBegin;
315: PetscCall(VecCopy(th->X0, ts->vec_sol));
316: PetscCall(VecCopy(th->V0, ts->vec_dot));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*
321: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X,Vec V)
322: {
323: TS_Alpha *th = (TS_Alpha*)ts->data;
324: PetscReal dt = t - ts->ptime;
326: PetscFunctionBegin;
327: PetscCall(VecCopy(ts->vec_dot,V));
328: PetscCall(VecAXPY(V,dt*(1-th->Gamma),th->A0));
329: PetscCall(VecAXPY(V,dt*th->Gamma,th->A1));
330: PetscCall(VecCopy(ts->vec_sol,X));
331: PetscCall(VecAXPY(X,dt,V));
332: PetscCall(VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0));
333: PetscCall(VecAXPY(X,dt*dt*th->Beta,th->A1));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
336: */
338: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
339: {
340: TS_Alpha *th = (TS_Alpha *)ts->data;
341: PetscReal ta = th->stage_time;
342: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
344: PetscFunctionBegin;
345: PetscCall(TSAlpha_StageVecs(ts, X));
346: /* F = Function(ta,Xa,Va,Aa) */
347: PetscCall(TSComputeI2Function(ts, ta, Xa, Va, Aa, F));
348: PetscCall(VecScale(F, th->scale_F));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
353: {
354: TS_Alpha *th = (TS_Alpha *)ts->data;
355: PetscReal ta = th->stage_time;
356: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
357: PetscReal dVdX = th->shift_V, dAdX = th->shift_A;
359: PetscFunctionBegin;
360: /* J,P = Jacobian(ta,Xa,Va,Aa) */
361: PetscCall(TSComputeI2Jacobian(ts, ta, Xa, Va, Aa, dVdX, dAdX, J, P));
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: static PetscErrorCode TSReset_Alpha(TS ts)
366: {
367: TS_Alpha *th = (TS_Alpha *)ts->data;
369: PetscFunctionBegin;
370: PetscCall(VecDestroy(&th->X0));
371: PetscCall(VecDestroy(&th->Xa));
372: PetscCall(VecDestroy(&th->X1));
373: PetscCall(VecDestroy(&th->V0));
374: PetscCall(VecDestroy(&th->Va));
375: PetscCall(VecDestroy(&th->V1));
376: PetscCall(VecDestroy(&th->A0));
377: PetscCall(VecDestroy(&th->Aa));
378: PetscCall(VecDestroy(&th->A1));
379: PetscCall(VecDestroy(&th->vec_sol_prev));
380: PetscCall(VecDestroy(&th->vec_dot_prev));
381: PetscCall(VecDestroy(&th->vec_lte_work[0]));
382: PetscCall(VecDestroy(&th->vec_lte_work[1]));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: static PetscErrorCode TSDestroy_Alpha(TS ts)
387: {
388: PetscFunctionBegin;
389: PetscCall(TSReset_Alpha(ts));
390: PetscCall(PetscFree(ts->data));
392: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", NULL));
393: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", NULL));
394: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", NULL));
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: static PetscErrorCode TSSetUp_Alpha(TS ts)
399: {
400: TS_Alpha *th = (TS_Alpha *)ts->data;
401: PetscBool match;
403: PetscFunctionBegin;
404: PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
405: PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
406: PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
407: PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
408: PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
409: PetscCall(VecDuplicate(ts->vec_sol, &th->V1));
410: PetscCall(VecDuplicate(ts->vec_sol, &th->A0));
411: PetscCall(VecDuplicate(ts->vec_sol, &th->Aa));
412: PetscCall(VecDuplicate(ts->vec_sol, &th->A1));
414: PetscCall(TSGetAdapt(ts, &ts->adapt));
415: PetscCall(TSAdaptCandidatesClear(ts->adapt));
416: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
417: if (!match) {
418: PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
419: PetscCall(VecDuplicate(ts->vec_sol, &th->vec_dot_prev));
420: PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work[0]));
421: PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work[1]));
422: }
424: PetscCall(TSGetSNES(ts, &ts->snes));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
429: {
430: TS_Alpha *th = (TS_Alpha *)ts->data;
432: PetscFunctionBegin;
433: PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
434: {
435: PetscBool flg;
436: PetscReal radius = 1;
437: PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlpha2SetRadius", radius, &radius, &flg));
438: if (flg) PetscCall(TSAlpha2SetRadius(ts, radius));
439: PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlpha2SetParams", th->Alpha_m, &th->Alpha_m, NULL));
440: PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlpha2SetParams", th->Alpha_f, &th->Alpha_f, NULL));
441: PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlpha2SetParams", th->Gamma, &th->Gamma, NULL));
442: PetscCall(PetscOptionsReal("-ts_alpha_beta", "Algorithmic parameter beta", "TSAlpha2SetParams", th->Beta, &th->Beta, NULL));
443: PetscCall(TSAlpha2SetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma, th->Beta));
444: }
445: PetscOptionsHeadEnd();
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
450: {
451: TS_Alpha *th = (TS_Alpha *)ts->data;
452: PetscBool iascii;
454: PetscFunctionBegin;
455: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
456: 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));
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts, PetscReal radius)
461: {
462: PetscReal alpha_m, alpha_f, gamma, beta;
464: PetscFunctionBegin;
465: PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
466: alpha_m = (2 - radius) / (1 + radius);
467: alpha_f = 1 / (1 + radius);
468: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
469: beta = (PetscReal)0.5 * (1 + alpha_m - alpha_f);
470: beta *= beta;
471: PetscCall(TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
476: {
477: TS_Alpha *th = (TS_Alpha *)ts->data;
478: PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
479: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
481: PetscFunctionBegin;
482: th->Alpha_m = alpha_m;
483: th->Alpha_f = alpha_f;
484: th->Gamma = gamma;
485: th->Beta = beta;
486: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
491: {
492: TS_Alpha *th = (TS_Alpha *)ts->data;
494: PetscFunctionBegin;
495: if (alpha_m) *alpha_m = th->Alpha_m;
496: if (alpha_f) *alpha_f = th->Alpha_f;
497: if (gamma) *gamma = th->Gamma;
498: if (beta) *beta = th->Beta;
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: /*MC
503: TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method for second-order systems {cite}`chung1993`
505: Level: beginner
507: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
508: M*/
509: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
510: {
511: TS_Alpha *th;
513: PetscFunctionBegin;
514: ts->ops->reset = TSReset_Alpha;
515: ts->ops->destroy = TSDestroy_Alpha;
516: ts->ops->view = TSView_Alpha;
517: ts->ops->setup = TSSetUp_Alpha;
518: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
519: ts->ops->step = TSStep_Alpha;
520: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
521: ts->ops->rollback = TSRollBack_Alpha;
522: /*ts->ops->interpolate = TSInterpolate_Alpha;*/
523: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
524: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
525: ts->default_adapt_type = TSADAPTNONE;
527: ts->usessnes = PETSC_TRUE;
529: PetscCall(PetscNew(&th));
530: ts->data = (void *)th;
532: th->Alpha_m = 0.5;
533: th->Alpha_f = 0.5;
534: th->Gamma = 0.5;
535: th->Beta = 0.25;
536: th->order = 2;
538: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", TSAlpha2SetRadius_Alpha));
539: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", TSAlpha2SetParams_Alpha));
540: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", TSAlpha2GetParams_Alpha));
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /*@
545: TSAlpha2SetRadius - sets the desired spectral radius of the method for `TSALPHA2`
546: (i.e. high-frequency numerical damping)
548: Logically Collective
550: Input Parameters:
551: + ts - timestepping context
552: - radius - the desired spectral radius
554: Options Database Key:
555: . -ts_alpha_radius <radius> - set the desired spectral radius
557: Level: intermediate
559: Notes:
561: The algorithmic parameters $\alpha_m$ and $\alpha_f$ of the generalized-$\alpha$ method can
562: be computed in terms of a specified spectral radius $\rho$ in `[0, 1]` for infinite time step
563: in order to control high-frequency numerical damping\:
565: $$
566: \begin{align*}
567: \alpha_m = (2-\rho)/(1+\rho) \\
568: \alpha_f = 1/(1+\rho)
569: \end{align*}
570: $$
572: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetParams()`, `TSAlpha2GetParams()`
573: @*/
574: PetscErrorCode TSAlpha2SetRadius(TS ts, PetscReal radius)
575: {
576: PetscFunctionBegin;
579: PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
580: PetscTryMethod(ts, "TSAlpha2SetRadius_C", (TS, PetscReal), (ts, radius));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: /*@
585: TSAlpha2SetParams - sets the algorithmic parameters for `TSALPHA2`
587: Logically Collective
589: Input Parameters:
590: + ts - timestepping context
591: . alpha_m - algorithmic parameter
592: . alpha_f - algorithmic parameter
593: . gamma - algorithmic parameter
594: - beta - algorithmic parameter
596: Options Database Keys:
597: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
598: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
599: . -ts_alpha_gamma <gamma> - set gamma
600: - -ts_alpha_beta <beta> - set beta
602: Level: advanced
604: Notes:
605: Second-order accuracy can be obtained so long as\:
607: $$
608: \begin{align*}
609: \gamma = 1/2 + \alpha_m - \alpha_f \\
610: \beta = 1/4 (1 + \alpha_m - \alpha_f)^2.
611: \end{align*}
612: $$
614: Unconditional stability requires\:
615: $$
616: \alpha_m >= \alpha_f >= 1/2.
617: $$
619: Use of this function is normally only required to hack `TSALPHA2` to use a modified
620: integration scheme. Users should call `TSAlpha2SetRadius()` to set the desired spectral
621: radius of the methods (i.e. high-frequency damping) in order so select optimal values for
622: these parameters.
624: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2GetParams()`
625: @*/
626: PetscErrorCode TSAlpha2SetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
627: {
628: PetscFunctionBegin;
634: PetscTryMethod(ts, "TSAlpha2SetParams_C", (TS, PetscReal, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma, beta));
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: /*@
639: TSAlpha2GetParams - gets the algorithmic parameters for `TSALPHA2`
641: Not Collective
643: Input Parameter:
644: . ts - timestepping context
646: Output Parameters:
647: + alpha_m - algorithmic parameter
648: . alpha_f - algorithmic parameter
649: . gamma - algorithmic parameter
650: - beta - algorithmic parameter
652: Level: advanced
654: Note:
655: Use of this function is normally only required to hack `TSALPHA2` to use a modified
656: integration scheme. Users should call `TSAlpha2SetRadius()` to set the high-frequency damping
657: (i.e. spectral radius of the method) in order so select optimal values for these parameters.
659: .seealso: [](ch_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
660: @*/
661: PetscErrorCode TSAlpha2GetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
662: {
663: PetscFunctionBegin;
665: if (alpha_m) PetscAssertPointer(alpha_m, 2);
666: if (alpha_f) PetscAssertPointer(alpha_f, 3);
667: if (gamma) PetscAssertPointer(gamma, 4);
668: if (beta) PetscAssertPointer(beta, 5);
669: PetscUseMethod(ts, "TSAlpha2GetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma, beta));
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }