Actual source code: alpha1.c
1: /*
2: Code for timestepping with implicit generalized-\alpha method
3: for first order systems.
4: */
5: #include <petsc/private/tsimpl.h>
7: static PetscBool cited = PETSC_FALSE;
8: static const char citation[] = "@article{Jansen2000,\n"
9: " title = {A generalized-$\\alpha$ method for integrating the filtered {N}avier--{S}tokes equations with a stabilized finite element method},\n"
10: " author = {Kenneth E. Jansen and Christian H. Whiting and Gregory M. Hulbert},\n"
11: " journal = {Computer Methods in Applied Mechanics and Engineering},\n"
12: " volume = {190},\n"
13: " number = {3--4},\n"
14: " pages = {305--319},\n"
15: " year = {2000},\n"
16: " issn = {0045-7825},\n"
17: " doi = {http://dx.doi.org/10.1016/S0045-7825(00)00203-6}\n}\n";
19: typedef struct {
20: PetscReal stage_time;
21: PetscReal shift_V;
22: PetscReal scale_F;
23: Vec X0, Xa, X1;
24: Vec V0, Va, V1;
26: PetscReal Alpha_m;
27: PetscReal Alpha_f;
28: PetscReal Gamma;
29: PetscInt order;
31: Vec vec_sol_prev;
32: Vec vec_lte_work;
34: TSStepStatus status;
35: } TS_Alpha;
37: static PetscErrorCode TSResizeRegister_Alpha(TS ts, PetscBool reg)
38: {
39: TS_Alpha *th = (TS_Alpha *)ts->data;
41: PetscFunctionBegin;
42: if (reg) {
43: PetscCall(TSResizeRegisterVec(ts, "ts:theta:sol_prev", th->vec_sol_prev));
44: PetscCall(TSResizeRegisterVec(ts, "ts:theta:X0", th->X0));
45: } else {
46: PetscCall(TSResizeRetrieveVec(ts, "ts:theta:sol_prev", &th->vec_sol_prev));
47: PetscCall(PetscObjectReference((PetscObject)th->vec_sol_prev));
48: PetscCall(TSResizeRetrieveVec(ts, "ts:theta:X0", &th->X0));
49: PetscCall(PetscObjectReference((PetscObject)th->X0));
50: }
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode TSAlpha_StageTime(TS ts)
55: {
56: TS_Alpha *th = (TS_Alpha *)ts->data;
57: PetscReal t = ts->ptime;
58: PetscReal dt = ts->time_step;
59: PetscReal Alpha_m = th->Alpha_m;
60: PetscReal Alpha_f = th->Alpha_f;
61: PetscReal Gamma = th->Gamma;
63: PetscFunctionBegin;
64: th->stage_time = t + Alpha_f * dt;
65: th->shift_V = Alpha_m / (Alpha_f * Gamma * dt);
66: th->scale_F = 1 / Alpha_f;
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: static PetscErrorCode TSAlpha_StageVecs(TS ts, Vec X)
71: {
72: TS_Alpha *th = (TS_Alpha *)ts->data;
73: Vec X1 = X, V1 = th->V1;
74: Vec Xa = th->Xa, Va = th->Va;
75: Vec X0 = th->X0, V0 = th->V0;
76: PetscReal dt = ts->time_step;
77: PetscReal Alpha_m = th->Alpha_m;
78: PetscReal Alpha_f = th->Alpha_f;
79: PetscReal Gamma = th->Gamma;
81: PetscFunctionBegin;
82: /* V1 = 1/(Gamma*dT)*(X1-X0) + (1-1/Gamma)*V0 */
83: PetscCall(VecWAXPY(V1, -1.0, X0, X1));
84: PetscCall(VecAXPBY(V1, 1 - 1 / Gamma, 1 / (Gamma * dt), V0));
85: /* Xa = X0 + Alpha_f*(X1-X0) */
86: PetscCall(VecWAXPY(Xa, -1.0, X0, X1));
87: PetscCall(VecAYPX(Xa, Alpha_f, X0));
88: /* Va = V0 + Alpha_m*(V1-V0) */
89: PetscCall(VecWAXPY(Va, -1.0, V0, V1));
90: PetscCall(VecAYPX(Va, Alpha_m, V0));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
95: {
96: PetscInt nits, lits;
98: PetscFunctionBegin;
99: PetscCall(SNESSolve(ts->snes, b, x));
100: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
101: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
102: ts->snes_its += nits;
103: ts->ksp_its += lits;
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: /*
108: Compute a consistent initial state for the generalized-alpha method.
109: - Solve two successive backward Euler steps with halved time step.
110: - Compute the initial time derivative using backward differences.
111: - If using adaptivity, estimate the LTE of the initial step.
112: */
113: static PetscErrorCode TSAlpha_Restart(TS ts, PetscBool *initok)
114: {
115: TS_Alpha *th = (TS_Alpha *)ts->data;
116: PetscReal time_step;
117: PetscReal alpha_m, alpha_f, gamma;
118: Vec X0 = ts->vec_sol, X1, X2 = th->X1;
119: PetscBool stageok;
121: PetscFunctionBegin;
122: PetscCall(VecDuplicate(X0, &X1));
124: /* Setup backward Euler with halved time step */
125: PetscCall(TSAlphaGetParams(ts, &alpha_m, &alpha_f, &gamma));
126: PetscCall(TSAlphaSetParams(ts, 1, 1, 1));
127: PetscCall(TSGetTimeStep(ts, &time_step));
128: ts->time_step = time_step / 2;
129: PetscCall(TSAlpha_StageTime(ts));
130: th->stage_time = ts->ptime;
131: PetscCall(VecZeroEntries(th->V0));
133: /* First BE step, (t0,X0) -> (t1,X1) */
134: th->stage_time += ts->time_step;
135: PetscCall(VecCopy(X0, th->X0));
136: PetscCall(TSPreStage(ts, th->stage_time));
137: PetscCall(VecCopy(th->X0, X1));
138: PetscCall(TSAlpha_SNESSolve(ts, NULL, X1));
139: PetscCall(TSPostStage(ts, th->stage_time, 0, &X1));
140: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
141: if (!stageok) goto finally;
143: /* Second BE step, (t1,X1) -> (t2,X2) */
144: th->stage_time += ts->time_step;
145: PetscCall(VecCopy(X1, th->X0));
146: PetscCall(TSPreStage(ts, th->stage_time));
147: PetscCall(VecCopy(th->X0, X2));
148: PetscCall(TSAlpha_SNESSolve(ts, NULL, X2));
149: PetscCall(TSPostStage(ts, th->stage_time, 0, &X2));
150: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X2, &stageok));
151: if (!stageok) goto finally;
153: /* Compute V0 ~ dX/dt at t0 with backward differences */
154: PetscCall(VecZeroEntries(th->V0));
155: PetscCall(VecAXPY(th->V0, -3 / time_step, X0));
156: PetscCall(VecAXPY(th->V0, +4 / time_step, X1));
157: PetscCall(VecAXPY(th->V0, -1 / time_step, X2));
159: /* Rough, lower-order estimate LTE of the initial step */
160: if (th->vec_lte_work) {
161: PetscCall(VecZeroEntries(th->vec_lte_work));
162: PetscCall(VecAXPY(th->vec_lte_work, +2, X2));
163: PetscCall(VecAXPY(th->vec_lte_work, -4, X1));
164: PetscCall(VecAXPY(th->vec_lte_work, +2, X0));
165: }
167: finally:
168: /* Revert TSAlpha to the initial state (t0,X0), but retain
169: potential time step reduction by factor after failed solve. */
170: if (initok) *initok = stageok;
171: PetscCall(TSSetTimeStep(ts, 2 * ts->time_step));
172: PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
173: PetscCall(VecCopy(ts->vec_sol, th->X0));
175: PetscCall(VecDestroy(&X1));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: static PetscErrorCode TSStep_Alpha(TS ts)
180: {
181: TS_Alpha *th = (TS_Alpha *)ts->data;
182: PetscInt rejections = 0;
183: PetscBool stageok, accept = PETSC_TRUE;
184: PetscReal next_time_step = ts->time_step;
186: PetscFunctionBegin;
187: PetscCall(PetscCitationsRegister(citation, &cited));
189: if (!ts->steprollback) {
190: if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
191: PetscCall(VecCopy(ts->vec_sol, th->X0));
192: PetscCall(VecCopy(th->V1, th->V0));
193: }
195: th->status = TS_STEP_INCOMPLETE;
196: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
197: if (ts->steprestart) {
198: PetscCall(TSAlpha_Restart(ts, &stageok));
199: if (!stageok) goto reject_step;
200: }
202: PetscCall(TSAlpha_StageTime(ts));
203: PetscCall(VecCopy(th->X0, th->X1));
204: PetscCall(TSPreStage(ts, th->stage_time));
205: PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
206: PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
207: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
208: if (!stageok) goto reject_step;
210: th->status = TS_STEP_PENDING;
211: PetscCall(VecCopy(th->X1, ts->vec_sol));
212: PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
213: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
214: if (!accept) {
215: PetscCall(VecCopy(th->X0, ts->vec_sol));
216: ts->time_step = next_time_step;
217: goto reject_step;
218: }
220: ts->ptime += ts->time_step;
221: ts->time_step = next_time_step;
222: break;
224: reject_step:
225: ts->reject++;
226: accept = PETSC_FALSE;
227: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
228: ts->reason = TS_DIVERGED_STEP_REJECTED;
229: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
230: }
231: }
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
236: {
237: TS_Alpha *th = (TS_Alpha *)ts->data;
238: Vec X = th->X1; /* X = solution */
239: Vec Y = th->vec_lte_work; /* Y = X + LTE */
240: PetscReal wltea, wlter;
242: PetscFunctionBegin;
243: if (!th->vec_sol_prev) {
244: *wlte = -1;
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
247: if (!th->vec_lte_work) {
248: *wlte = -1;
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
251: if (ts->steprestart) {
252: /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */
253: PetscCall(VecAXPY(Y, 1, X));
254: } else {
255: /* Compute LTE using backward differences with non-constant time step */
256: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
257: PetscReal a = 1 + h_prev / h;
258: PetscScalar scal[3];
259: Vec vecs[3];
260: scal[0] = +1 / a;
261: scal[1] = -1 / (a - 1);
262: scal[2] = +1 / (a * (a - 1));
263: vecs[0] = th->X1;
264: vecs[1] = th->X0;
265: vecs[2] = th->vec_sol_prev;
266: PetscCall(VecCopy(X, Y));
267: PetscCall(VecMAXPY(Y, 3, scal, vecs));
268: }
269: PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
270: if (order) *order = 2;
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: static PetscErrorCode TSInterpolate_Alpha(TS ts, PetscReal t, Vec X)
275: {
276: TS_Alpha *th = (TS_Alpha *)ts->data;
277: PetscReal dt = t - ts->ptime;
278: PetscReal Gamma = th->Gamma;
280: PetscFunctionBegin;
281: PetscCall(VecWAXPY(th->V1, -1.0, th->X0, ts->vec_sol));
282: PetscCall(VecAXPBY(th->V1, 1 - 1 / Gamma, 1 / (Gamma * ts->time_step), th->V0));
283: PetscCall(VecCopy(ts->vec_sol, X));
284: /* X = X + Gamma*dT*V1 */
285: PetscCall(VecAXPY(X, th->Gamma * dt, th->V1));
286: /* X = X + (1-Gamma)*dT*V0 */
287: PetscCall(VecAXPY(X, (1 - th->Gamma) * dt, th->V0));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
292: {
293: TS_Alpha *th = (TS_Alpha *)ts->data;
294: PetscReal ta = th->stage_time;
295: Vec Xa = th->Xa, Va = th->Va;
297: PetscFunctionBegin;
298: PetscCall(TSAlpha_StageVecs(ts, X));
299: /* F = Function(ta,Xa,Va) */
300: PetscCall(TSComputeIFunction(ts, ta, Xa, Va, F, PETSC_FALSE));
301: PetscCall(VecScale(F, th->scale_F));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
306: {
307: TS_Alpha *th = (TS_Alpha *)ts->data;
308: PetscReal ta = th->stage_time;
309: Vec Xa = th->Xa, Va = th->Va;
310: PetscReal dVdX = th->shift_V;
312: PetscFunctionBegin;
313: /* J,P = Jacobian(ta,Xa,Va) */
314: PetscCall(TSComputeIJacobian(ts, ta, Xa, Va, dVdX, J, P, PETSC_FALSE));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: static PetscErrorCode TSReset_Alpha(TS ts)
319: {
320: TS_Alpha *th = (TS_Alpha *)ts->data;
322: PetscFunctionBegin;
323: PetscCall(VecDestroy(&th->X0));
324: PetscCall(VecDestroy(&th->Xa));
325: PetscCall(VecDestroy(&th->X1));
326: PetscCall(VecDestroy(&th->V0));
327: PetscCall(VecDestroy(&th->Va));
328: PetscCall(VecDestroy(&th->V1));
329: PetscCall(VecDestroy(&th->vec_sol_prev));
330: PetscCall(VecDestroy(&th->vec_lte_work));
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: static PetscErrorCode TSDestroy_Alpha(TS ts)
335: {
336: PetscFunctionBegin;
337: PetscCall(TSReset_Alpha(ts));
338: PetscCall(PetscFree(ts->data));
340: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", NULL));
341: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", NULL));
342: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", NULL));
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: static PetscErrorCode TSSetUp_Alpha(TS ts)
347: {
348: TS_Alpha *th = (TS_Alpha *)ts->data;
349: PetscBool match;
351: PetscFunctionBegin;
352: if (!th->X0) PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
353: PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
354: PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
355: PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
356: PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
357: PetscCall(VecDuplicate(ts->vec_sol, &th->V1));
359: PetscCall(TSGetAdapt(ts, &ts->adapt));
360: PetscCall(TSAdaptCandidatesClear(ts->adapt));
361: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
362: if (!match) {
363: if (!th->vec_sol_prev) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
364: if (!th->vec_lte_work) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work));
365: }
367: PetscCall(TSGetSNES(ts, &ts->snes));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
372: {
373: TS_Alpha *th = (TS_Alpha *)ts->data;
375: PetscFunctionBegin;
376: PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
377: {
378: PetscBool flg;
379: PetscReal radius = 1;
380: PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlphaSetRadius", radius, &radius, &flg));
381: if (flg) PetscCall(TSAlphaSetRadius(ts, radius));
382: PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlphaSetParams", th->Alpha_m, &th->Alpha_m, NULL));
383: PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlphaSetParams", th->Alpha_f, &th->Alpha_f, NULL));
384: PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlphaSetParams", th->Gamma, &th->Gamma, NULL));
385: PetscCall(TSAlphaSetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma));
386: }
387: PetscOptionsHeadEnd();
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
392: {
393: TS_Alpha *th = (TS_Alpha *)ts->data;
394: PetscBool iascii;
396: PetscFunctionBegin;
397: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
398: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " Alpha_m=%g, Alpha_f=%g, Gamma=%g\n", (double)th->Alpha_m, (double)th->Alpha_f, (double)th->Gamma));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts, PetscReal radius)
403: {
404: PetscReal alpha_m, alpha_f, gamma;
406: PetscFunctionBegin;
407: PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
408: alpha_m = (PetscReal)0.5 * (3 - radius) / (1 + radius);
409: alpha_f = 1 / (1 + radius);
410: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
411: PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode TSAlphaSetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
416: {
417: TS_Alpha *th = (TS_Alpha *)ts->data;
418: PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
419: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
421: PetscFunctionBegin;
422: th->Alpha_m = alpha_m;
423: th->Alpha_f = alpha_f;
424: th->Gamma = gamma;
425: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: static PetscErrorCode TSAlphaGetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
430: {
431: TS_Alpha *th = (TS_Alpha *)ts->data;
433: PetscFunctionBegin;
434: if (alpha_m) *alpha_m = th->Alpha_m;
435: if (alpha_f) *alpha_f = th->Alpha_f;
436: if (gamma) *gamma = th->Gamma;
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
440: /*MC
441: TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method {cite}`jansen_2000` {cite}`chung1993` for first-order systems
443: Level: beginner
445: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
446: M*/
447: PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts)
448: {
449: TS_Alpha *th;
451: PetscFunctionBegin;
452: ts->ops->reset = TSReset_Alpha;
453: ts->ops->destroy = TSDestroy_Alpha;
454: ts->ops->view = TSView_Alpha;
455: ts->ops->setup = TSSetUp_Alpha;
456: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
457: ts->ops->step = TSStep_Alpha;
458: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
459: ts->ops->interpolate = TSInterpolate_Alpha;
460: ts->ops->resizeregister = TSResizeRegister_Alpha;
461: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
462: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
463: ts->default_adapt_type = TSADAPTNONE;
465: ts->usessnes = PETSC_TRUE;
467: PetscCall(PetscNew(&th));
468: ts->data = (void *)th;
470: th->Alpha_m = 0.5;
471: th->Alpha_f = 0.5;
472: th->Gamma = 0.5;
473: th->order = 2;
475: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", TSAlphaSetRadius_Alpha));
476: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", TSAlphaSetParams_Alpha));
477: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", TSAlphaGetParams_Alpha));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@
482: TSAlphaSetRadius - sets the desired spectral radius of the method for `TSALPHA`
483: (i.e. high-frequency numerical damping)
485: Logically Collective
487: Input Parameters:
488: + ts - timestepping context
489: - radius - the desired spectral radius
491: Options Database Key:
492: . -ts_alpha_radius <radius> - set alpha radius
494: Level: intermediate
496: Notes:
497: The algorithmic parameters $\alpha_m$ and $\alpha_f$ of the generalized-$\alpha$ method can
498: be computed in terms of a specified spectral radius $\rho$ in [0, 1] for infinite time step
499: in order to control high-frequency numerical damping\:
501: $$
502: \begin{align*}
503: \alpha_m = 0.5*(3-\rho)/(1+\rho) \\
504: \alpha_f = 1/(1+\rho)
505: \end{align*}
506: $$
508: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetParams()`, `TSAlphaGetParams()`
509: @*/
510: PetscErrorCode TSAlphaSetRadius(TS ts, PetscReal radius)
511: {
512: PetscFunctionBegin;
515: PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
516: PetscTryMethod(ts, "TSAlphaSetRadius_C", (TS, PetscReal), (ts, radius));
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: /*@
521: TSAlphaSetParams - sets the algorithmic parameters for `TSALPHA`
523: Logically Collective
525: Input Parameters:
526: + ts - timestepping context
527: . alpha_m - algorithmic parameter
528: . alpha_f - algorithmic parameter
529: - gamma - algorithmic parameter
531: Options Database Keys:
532: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
533: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
534: - -ts_alpha_gamma <gamma> - set gamma
536: Level: advanced
538: Note:
539: Second-order accuracy can be obtained so long as\: $\gamma = 0.5 + \alpha_m - \alpha_f$
541: Unconditional stability requires\: $\alpha_m >= \alpha_f >= 0.5$
543: Backward Euler method is recovered with\: $\alpha_m = \alpha_f = \gamma = 1$
545: Use of this function is normally only required to hack `TSALPHA` to use a modified
546: integration scheme. Users should call `TSAlphaSetRadius()` to set the desired spectral radius
547: of the methods (i.e. high-frequency damping) in order so select optimal values for these
548: parameters.
550: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaGetParams()`
551: @*/
552: PetscErrorCode TSAlphaSetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
553: {
554: PetscFunctionBegin;
559: PetscTryMethod(ts, "TSAlphaSetParams_C", (TS, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma));
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: /*@
564: TSAlphaGetParams - gets the algorithmic parameters for `TSALPHA`
566: Not Collective
568: Input Parameter:
569: . ts - timestepping context
571: Output Parameters:
572: + alpha_m - algorithmic parameter
573: . alpha_f - algorithmic parameter
574: - gamma - algorithmic parameter
576: Level: advanced
578: Note:
579: Use of this function is normally only required to hack `TSALPHA` to use a modified
580: integration scheme. Users should call `TSAlphaSetRadius()` to set the high-frequency damping
581: (i.e. spectral radius of the method) in order so select optimal values for these parameters.
583: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
584: @*/
585: PetscErrorCode TSAlphaGetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
586: {
587: PetscFunctionBegin;
589: if (alpha_m) PetscAssertPointer(alpha_m, 2);
590: if (alpha_f) PetscAssertPointer(alpha_f, 3);
591: if (gamma) PetscAssertPointer(gamma, 4);
592: PetscUseMethod(ts, "TSAlphaGetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }