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: th->stage_time = t + Alpha_f * dt;
54: th->shift_V = Gamma / (dt * Beta);
55: th->shift_A = Alpha_m / (Alpha_f * dt * dt * Beta);
56: th->scale_F = 1 / Alpha_f;
57: return 0;
58: }
60: static PetscErrorCode TSAlpha_StageVecs(TS ts, Vec X)
61: {
62: TS_Alpha *th = (TS_Alpha *)ts->data;
63: Vec X1 = X, V1 = th->V1, A1 = th->A1;
64: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
65: Vec X0 = th->X0, V0 = th->V0, A0 = th->A0;
66: PetscReal dt = ts->time_step;
67: PetscReal Alpha_m = th->Alpha_m;
68: PetscReal Alpha_f = th->Alpha_f;
69: PetscReal Gamma = th->Gamma;
70: PetscReal Beta = th->Beta;
72: /* A1 = ... */
73: VecWAXPY(A1, -1.0, X0, X1);
74: VecAXPY(A1, -dt, V0);
75: VecAXPBY(A1, -(1 - 2 * Beta) / (2 * Beta), 1 / (dt * dt * Beta), A0);
76: /* V1 = ... */
77: VecWAXPY(V1, (1.0 - Gamma) / Gamma, A0, A1);
78: VecAYPX(V1, dt * Gamma, V0);
79: /* Xa = X0 + Alpha_f*(X1-X0) */
80: VecWAXPY(Xa, -1.0, X0, X1);
81: VecAYPX(Xa, Alpha_f, X0);
82: /* Va = V0 + Alpha_f*(V1-V0) */
83: VecWAXPY(Va, -1.0, V0, V1);
84: VecAYPX(Va, Alpha_f, V0);
85: /* Aa = A0 + Alpha_m*(A1-A0) */
86: VecWAXPY(Aa, -1.0, A0, A1);
87: VecAYPX(Aa, Alpha_m, A0);
88: return 0;
89: }
91: static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
92: {
93: PetscInt nits, lits;
95: SNESSolve(ts->snes, b, x);
96: SNESGetIterationNumber(ts->snes, &nits);
97: SNESGetLinearSolveIterations(ts->snes, &lits);
98: ts->snes_its += nits;
99: ts->ksp_its += lits;
100: return 0;
101: }
103: /*
104: Compute a consistent initial state for the generalized-alpha method.
105: - Solve two successive backward Euler steps with halved time step.
106: - Compute the initial second time derivative using backward differences.
107: - If using adaptivity, estimate the LTE of the initial step.
108: */
109: static PetscErrorCode TSAlpha_Restart(TS ts, PetscBool *initok)
110: {
111: TS_Alpha *th = (TS_Alpha *)ts->data;
112: PetscReal time_step;
113: PetscReal alpha_m, alpha_f, gamma, beta;
114: Vec X0 = ts->vec_sol, X1, X2 = th->X1;
115: Vec V0 = ts->vec_dot, V1, V2 = th->V1;
116: PetscBool stageok;
118: VecDuplicate(X0, &X1);
119: VecDuplicate(V0, &V1);
121: /* Setup backward Euler with halved time step */
122: TSAlpha2GetParams(ts, &alpha_m, &alpha_f, &gamma, &beta);
123: TSAlpha2SetParams(ts, 1, 1, 1, 0.5);
124: TSGetTimeStep(ts, &time_step);
125: ts->time_step = time_step / 2;
126: TSAlpha_StageTime(ts);
127: th->stage_time = ts->ptime;
128: VecZeroEntries(th->A0);
130: /* First BE step, (t0,X0,V0) -> (t1,X1,V1) */
131: th->stage_time += ts->time_step;
132: VecCopy(X0, th->X0);
133: VecCopy(V0, th->V0);
134: TSPreStage(ts, th->stage_time);
135: VecCopy(th->X0, X1);
136: TSAlpha_SNESSolve(ts, NULL, X1);
137: VecCopy(th->V1, V1);
138: TSPostStage(ts, th->stage_time, 0, &X1);
139: TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok);
140: if (!stageok) goto finally;
142: /* Second BE step, (t1,X1,V1) -> (t2,X2,V2) */
143: th->stage_time += ts->time_step;
144: VecCopy(X1, th->X0);
145: VecCopy(V1, th->V0);
146: TSPreStage(ts, th->stage_time);
147: VecCopy(th->X0, X2);
148: TSAlpha_SNESSolve(ts, NULL, X2);
149: VecCopy(th->V1, V2);
150: TSPostStage(ts, th->stage_time, 0, &X2);
151: TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok);
152: if (!stageok) goto finally;
154: /* Compute A0 ~ dV/dt at t0 with backward differences */
155: VecZeroEntries(th->A0);
156: VecAXPY(th->A0, -3 / ts->time_step, V0);
157: VecAXPY(th->A0, +4 / ts->time_step, V1);
158: VecAXPY(th->A0, -1 / ts->time_step, V2);
160: /* Rough, lower-order estimate LTE of the initial step */
161: if (th->vec_lte_work[0]) {
162: VecZeroEntries(th->vec_lte_work[0]);
163: VecAXPY(th->vec_lte_work[0], +2, X2);
164: VecAXPY(th->vec_lte_work[0], -4, X1);
165: VecAXPY(th->vec_lte_work[0], +2, X0);
166: }
167: if (th->vec_lte_work[1]) {
168: VecZeroEntries(th->vec_lte_work[1]);
169: VecAXPY(th->vec_lte_work[1], +2, V2);
170: VecAXPY(th->vec_lte_work[1], -4, V1);
171: VecAXPY(th->vec_lte_work[1], +2, V0);
172: }
174: finally:
175: /* Revert TSAlpha to the initial state (t0,X0,V0) */
176: if (initok) *initok = stageok;
177: TSSetTimeStep(ts, time_step);
178: TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta);
179: VecCopy(ts->vec_sol, th->X0);
180: VecCopy(ts->vec_dot, th->V0);
182: VecDestroy(&X1);
183: VecDestroy(&V1);
184: return 0;
185: }
187: static PetscErrorCode TSStep_Alpha(TS ts)
188: {
189: TS_Alpha *th = (TS_Alpha *)ts->data;
190: PetscInt rejections = 0;
191: PetscBool stageok, accept = PETSC_TRUE;
192: PetscReal next_time_step = ts->time_step;
194: PetscCitationsRegister(citation, &cited);
196: if (!ts->steprollback) {
197: if (th->vec_sol_prev) VecCopy(th->X0, th->vec_sol_prev);
198: if (th->vec_dot_prev) VecCopy(th->V0, th->vec_dot_prev);
199: VecCopy(ts->vec_sol, th->X0);
200: VecCopy(ts->vec_dot, th->V0);
201: VecCopy(th->A1, th->A0);
202: }
204: th->status = TS_STEP_INCOMPLETE;
205: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
206: if (ts->steprestart) {
207: TSAlpha_Restart(ts, &stageok);
208: if (!stageok) goto reject_step;
209: }
211: TSAlpha_StageTime(ts);
212: VecCopy(th->X0, th->X1);
213: TSPreStage(ts, th->stage_time);
214: TSAlpha_SNESSolve(ts, NULL, th->X1);
215: TSPostStage(ts, th->stage_time, 0, &th->Xa);
216: TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok);
217: if (!stageok) goto reject_step;
219: th->status = TS_STEP_PENDING;
220: VecCopy(th->X1, ts->vec_sol);
221: VecCopy(th->V1, ts->vec_dot);
222: TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
223: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
224: if (!accept) {
225: VecCopy(th->X0, ts->vec_sol);
226: VecCopy(th->V0, ts->vec_dot);
227: ts->time_step = next_time_step;
228: goto reject_step;
229: }
231: ts->ptime += ts->time_step;
232: ts->time_step = next_time_step;
233: break;
235: reject_step:
236: ts->reject++;
237: accept = PETSC_FALSE;
238: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
239: ts->reason = TS_DIVERGED_STEP_REJECTED;
240: PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections);
241: }
242: }
243: return 0;
244: }
246: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
247: {
248: TS_Alpha *th = (TS_Alpha *)ts->data;
249: Vec X = th->X1; /* X = solution */
250: Vec V = th->V1; /* V = solution */
251: Vec Y = th->vec_lte_work[0]; /* Y = X + LTE */
252: Vec Z = th->vec_lte_work[1]; /* Z = V + LTE */
253: PetscReal enormX, enormV, enormXa, enormVa, enormXr, enormVr;
255: if (!th->vec_sol_prev) {
256: *wlte = -1;
257: return 0;
258: }
259: if (!th->vec_dot_prev) {
260: *wlte = -1;
261: return 0;
262: }
263: if (!th->vec_lte_work[0]) {
264: *wlte = -1;
265: return 0;
266: }
267: if (!th->vec_lte_work[1]) {
268: *wlte = -1;
269: return 0;
270: }
271: if (ts->steprestart) {
272: /* th->vec_lte_prev is set to the LTE in TSAlpha_Restart() */
273: VecAXPY(Y, 1, X);
274: VecAXPY(Z, 1, V);
275: } else {
276: /* Compute LTE using backward differences with non-constant time step */
277: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
278: PetscReal a = 1 + h_prev / h;
279: PetscScalar scal[3];
280: Vec vecX[3], vecV[3];
281: scal[0] = +1 / a;
282: scal[1] = -1 / (a - 1);
283: scal[2] = +1 / (a * (a - 1));
284: vecX[0] = th->X1;
285: vecX[1] = th->X0;
286: vecX[2] = th->vec_sol_prev;
287: vecV[0] = th->V1;
288: vecV[1] = th->V0;
289: vecV[2] = th->vec_dot_prev;
290: VecCopy(X, Y);
291: VecMAXPY(Y, 3, scal, vecX);
292: VecCopy(V, Z);
293: VecMAXPY(Z, 3, scal, vecV);
294: }
295: /* XXX ts->atol and ts->vatol are not appropriate for computing enormV */
296: TSErrorWeightedNorm(ts, X, Y, wnormtype, &enormX, &enormXa, &enormXr);
297: TSErrorWeightedNorm(ts, V, Z, wnormtype, &enormV, &enormVa, &enormVr);
298: if (wnormtype == NORM_2) *wlte = PetscSqrtReal(PetscSqr(enormX) / 2 + PetscSqr(enormV) / 2);
299: else *wlte = PetscMax(enormX, enormV);
300: if (order) *order = 2;
301: return 0;
302: }
304: static PetscErrorCode TSRollBack_Alpha(TS ts)
305: {
306: TS_Alpha *th = (TS_Alpha *)ts->data;
308: VecCopy(th->X0, ts->vec_sol);
309: VecCopy(th->V0, ts->vec_dot);
310: return 0;
311: }
313: /*
314: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X,Vec V)
315: {
316: TS_Alpha *th = (TS_Alpha*)ts->data;
317: PetscReal dt = t - ts->ptime;
319: VecCopy(ts->vec_dot,V);
320: VecAXPY(V,dt*(1-th->Gamma),th->A0);
321: VecAXPY(V,dt*th->Gamma,th->A1);
322: VecCopy(ts->vec_sol,X);
323: VecAXPY(X,dt,V);
324: VecAXPY(X,dt*dt*((PetscReal)0.5-th->Beta),th->A0);
325: VecAXPY(X,dt*dt*th->Beta,th->A1);
326: return 0;
327: }
328: */
330: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
331: {
332: TS_Alpha *th = (TS_Alpha *)ts->data;
333: PetscReal ta = th->stage_time;
334: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
336: TSAlpha_StageVecs(ts, X);
337: /* F = Function(ta,Xa,Va,Aa) */
338: TSComputeI2Function(ts, ta, Xa, Va, Aa, F);
339: VecScale(F, th->scale_F);
340: return 0;
341: }
343: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
344: {
345: TS_Alpha *th = (TS_Alpha *)ts->data;
346: PetscReal ta = th->stage_time;
347: Vec Xa = th->Xa, Va = th->Va, Aa = th->Aa;
348: PetscReal dVdX = th->shift_V, dAdX = th->shift_A;
350: /* J,P = Jacobian(ta,Xa,Va,Aa) */
351: TSComputeI2Jacobian(ts, ta, Xa, Va, Aa, dVdX, dAdX, J, P);
352: return 0;
353: }
355: static PetscErrorCode TSReset_Alpha(TS ts)
356: {
357: TS_Alpha *th = (TS_Alpha *)ts->data;
359: VecDestroy(&th->X0);
360: VecDestroy(&th->Xa);
361: VecDestroy(&th->X1);
362: VecDestroy(&th->V0);
363: VecDestroy(&th->Va);
364: VecDestroy(&th->V1);
365: VecDestroy(&th->A0);
366: VecDestroy(&th->Aa);
367: VecDestroy(&th->A1);
368: VecDestroy(&th->vec_sol_prev);
369: VecDestroy(&th->vec_dot_prev);
370: VecDestroy(&th->vec_lte_work[0]);
371: VecDestroy(&th->vec_lte_work[1]);
372: return 0;
373: }
375: static PetscErrorCode TSDestroy_Alpha(TS ts)
376: {
377: TSReset_Alpha(ts);
378: PetscFree(ts->data);
380: PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", NULL);
381: PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", NULL);
382: PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", NULL);
383: return 0;
384: }
386: static PetscErrorCode TSSetUp_Alpha(TS ts)
387: {
388: TS_Alpha *th = (TS_Alpha *)ts->data;
389: PetscBool match;
391: VecDuplicate(ts->vec_sol, &th->X0);
392: VecDuplicate(ts->vec_sol, &th->Xa);
393: VecDuplicate(ts->vec_sol, &th->X1);
394: VecDuplicate(ts->vec_sol, &th->V0);
395: VecDuplicate(ts->vec_sol, &th->Va);
396: VecDuplicate(ts->vec_sol, &th->V1);
397: VecDuplicate(ts->vec_sol, &th->A0);
398: VecDuplicate(ts->vec_sol, &th->Aa);
399: VecDuplicate(ts->vec_sol, &th->A1);
401: TSGetAdapt(ts, &ts->adapt);
402: TSAdaptCandidatesClear(ts->adapt);
403: PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match);
404: if (!match) {
405: VecDuplicate(ts->vec_sol, &th->vec_sol_prev);
406: VecDuplicate(ts->vec_sol, &th->vec_dot_prev);
407: VecDuplicate(ts->vec_sol, &th->vec_lte_work[0]);
408: VecDuplicate(ts->vec_sol, &th->vec_lte_work[1]);
409: }
411: TSGetSNES(ts, &ts->snes);
412: return 0;
413: }
415: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
416: {
417: TS_Alpha *th = (TS_Alpha *)ts->data;
419: PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
420: {
421: PetscBool flg;
422: PetscReal radius = 1;
423: PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlpha2SetRadius", radius, &radius, &flg);
424: if (flg) TSAlpha2SetRadius(ts, radius);
425: PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlpha2SetParams", th->Alpha_m, &th->Alpha_m, NULL);
426: PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlpha2SetParams", th->Alpha_f, &th->Alpha_f, NULL);
427: PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlpha2SetParams", th->Gamma, &th->Gamma, NULL);
428: PetscOptionsReal("-ts_alpha_beta", "Algorithmic parameter beta", "TSAlpha2SetParams", th->Beta, &th->Beta, NULL);
429: TSAlpha2SetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma, th->Beta);
430: }
431: PetscOptionsHeadEnd();
432: return 0;
433: }
435: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
436: {
437: TS_Alpha *th = (TS_Alpha *)ts->data;
438: PetscBool iascii;
440: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
441: if (iascii) 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);
442: return 0;
443: }
445: static PetscErrorCode TSAlpha2SetRadius_Alpha(TS ts, PetscReal radius)
446: {
447: PetscReal alpha_m, alpha_f, gamma, beta;
450: alpha_m = (2 - radius) / (1 + radius);
451: alpha_f = 1 / (1 + radius);
452: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
453: beta = (PetscReal)0.5 * (1 + alpha_m - alpha_f);
454: beta *= beta;
455: TSAlpha2SetParams(ts, alpha_m, alpha_f, gamma, beta);
456: return 0;
457: }
459: static PetscErrorCode TSAlpha2SetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
460: {
461: TS_Alpha *th = (TS_Alpha *)ts->data;
462: PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
463: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
465: th->Alpha_m = alpha_m;
466: th->Alpha_f = alpha_f;
467: th->Gamma = gamma;
468: th->Beta = beta;
469: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
470: return 0;
471: }
473: static PetscErrorCode TSAlpha2GetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
474: {
475: TS_Alpha *th = (TS_Alpha *)ts->data;
477: if (alpha_m) *alpha_m = th->Alpha_m;
478: if (alpha_f) *alpha_f = th->Alpha_f;
479: if (gamma) *gamma = th->Gamma;
480: if (beta) *beta = th->Beta;
481: return 0;
482: }
484: /*MC
485: TSALPHA2 - ODE/DAE solver using the implicit Generalized-Alpha method
486: for second-order systems
488: Level: beginner
490: References:
491: . * - J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
492: Dynamics with Improved Numerical Dissipation: The Generalized-alpha
493: Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.
495: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
496: M*/
497: PETSC_EXTERN PetscErrorCode TSCreate_Alpha2(TS ts)
498: {
499: TS_Alpha *th;
501: ts->ops->reset = TSReset_Alpha;
502: ts->ops->destroy = TSDestroy_Alpha;
503: ts->ops->view = TSView_Alpha;
504: ts->ops->setup = TSSetUp_Alpha;
505: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
506: ts->ops->step = TSStep_Alpha;
507: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
508: ts->ops->rollback = TSRollBack_Alpha;
509: /*ts->ops->interpolate = TSInterpolate_Alpha;*/
510: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
511: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
512: ts->default_adapt_type = TSADAPTNONE;
514: ts->usessnes = PETSC_TRUE;
516: PetscNew(&th);
517: ts->data = (void *)th;
519: th->Alpha_m = 0.5;
520: th->Alpha_f = 0.5;
521: th->Gamma = 0.5;
522: th->Beta = 0.25;
523: th->order = 2;
525: PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetRadius_C", TSAlpha2SetRadius_Alpha);
526: PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2SetParams_C", TSAlpha2SetParams_Alpha);
527: PetscObjectComposeFunction((PetscObject)ts, "TSAlpha2GetParams_C", TSAlpha2GetParams_Alpha);
528: return 0;
529: }
531: /*@
532: TSAlpha2SetRadius - sets the desired spectral radius of the method for `TSALPHA2`
533: (i.e. high-frequency numerical damping)
535: Logically Collective
537: The algorithmic parameters \alpha_m and \alpha_f of the
538: generalized-\alpha method can be computed in terms of a specified
539: spectral radius \rho in [0,1] for infinite time step in order to
540: control high-frequency numerical damping:
541: \alpha_m = (2-\rho)/(1+\rho)
542: \alpha_f = 1/(1+\rho)
544: Input Parameters:
545: + ts - timestepping context
546: - radius - the desired spectral radius
548: Options Database Key:
549: . -ts_alpha_radius <radius> - set the desired spectral radius
551: Level: intermediate
553: .seealso: [](chapter_ts), `TS`, `TSALPHA2`, `TSAlpha2SetParams()`, `TSAlpha2GetParams()`
554: @*/
555: PetscErrorCode TSAlpha2SetRadius(TS ts, PetscReal radius)
556: {
560: PetscTryMethod(ts, "TSAlpha2SetRadius_C", (TS, PetscReal), (ts, radius));
561: return 0;
562: }
564: /*@
565: TSAlpha2SetParams - sets the algorithmic parameters for `TSALPHA2`
567: Logically Collective
569: Second-order accuracy can be obtained so long as:
570: \gamma = 1/2 + alpha_m - alpha_f
571: \beta = 1/4 (1 + alpha_m - alpha_f)^2
573: Unconditional stability requires:
574: \alpha_m >= \alpha_f >= 1/2
576: Input Parameters:
577: + ts - timestepping context
578: . \alpha_m - algorithmic parameter
579: . \alpha_f - algorithmic parameter
580: . \gamma - algorithmic parameter
581: - \beta - algorithmic parameter
583: Options Database Keys:
584: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
585: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
586: . -ts_alpha_gamma <gamma> - set gamma
587: - -ts_alpha_beta <beta> - set beta
589: Level: advanced
591: Note:
592: Use of this function is normally only required to hack `TSALPHA2` to
593: use a modified integration scheme. Users should call
594: `TSAlpha2SetRadius()` to set the desired spectral radius of the methods
595: (i.e. high-frequency damping) in order so select optimal values for
596: these parameters.
598: .seealso: [](chapter_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2GetParams()`
599: @*/
600: PetscErrorCode TSAlpha2SetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma, PetscReal beta)
601: {
607: PetscTryMethod(ts, "TSAlpha2SetParams_C", (TS, PetscReal, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma, beta));
608: return 0;
609: }
611: /*@
612: TSAlpha2GetParams - gets the algorithmic parameters for `TSALPHA2`
614: Not Collective
616: Input Parameter:
617: . ts - timestepping context
619: Output Parameters:
620: + \alpha_m - algorithmic parameter
621: . \alpha_f - algorithmic parameter
622: . \gamma - algorithmic parameter
623: - \beta - algorithmic parameter
625: Level: advanced
627: Note:
628: Use of this function is normally only required to hack `TSALPHA2` to
629: use a modified integration scheme. Users should call
630: `TSAlpha2SetRadius()` to set the high-frequency damping (i.e. spectral
631: radius of the method) in order so select optimal values for these
632: parameters.
634: .seealso: [](chapter_ts), `TS`, `TSALPHA2`, `TSAlpha2SetRadius()`, `TSAlpha2SetParams()`
635: @*/
636: PetscErrorCode TSAlpha2GetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma, PetscReal *beta)
637: {
643: PetscUseMethod(ts, "TSAlpha2GetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma, beta));
644: return 0;
645: }