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 / ts->time_step, X0));
156: PetscCall(VecAXPY(th->V0, +4 / ts->time_step, X1));
157: PetscCall(VecAXPY(th->V0, -1 / ts->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) */
169: if (initok) *initok = stageok;
170: PetscCall(TSSetTimeStep(ts, time_step));
171: PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
172: PetscCall(VecCopy(ts->vec_sol, th->X0));
174: PetscCall(VecDestroy(&X1));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode TSStep_Alpha(TS ts)
179: {
180: TS_Alpha *th = (TS_Alpha *)ts->data;
181: PetscInt rejections = 0;
182: PetscBool stageok, accept = PETSC_TRUE;
183: PetscReal next_time_step = ts->time_step;
185: PetscFunctionBegin;
186: PetscCall(PetscCitationsRegister(citation, &cited));
188: if (!ts->steprollback) {
189: if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
190: PetscCall(VecCopy(ts->vec_sol, th->X0));
191: PetscCall(VecCopy(th->V1, th->V0));
192: }
194: th->status = TS_STEP_INCOMPLETE;
195: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
196: if (ts->steprestart) {
197: PetscCall(TSAlpha_Restart(ts, &stageok));
198: if (!stageok) goto reject_step;
199: }
201: PetscCall(TSAlpha_StageTime(ts));
202: PetscCall(VecCopy(th->X0, th->X1));
203: PetscCall(TSPreStage(ts, th->stage_time));
204: PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
205: PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
206: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
207: if (!stageok) goto reject_step;
209: th->status = TS_STEP_PENDING;
210: PetscCall(VecCopy(th->X1, ts->vec_sol));
211: PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
212: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
213: if (!accept) {
214: PetscCall(VecCopy(th->X0, ts->vec_sol));
215: ts->time_step = next_time_step;
216: goto reject_step;
217: }
219: ts->ptime += ts->time_step;
220: ts->time_step = next_time_step;
221: break;
223: reject_step:
224: ts->reject++;
225: accept = PETSC_FALSE;
226: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
227: ts->reason = TS_DIVERGED_STEP_REJECTED;
228: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
229: }
230: }
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
235: {
236: TS_Alpha *th = (TS_Alpha *)ts->data;
237: Vec X = th->X1; /* X = solution */
238: Vec Y = th->vec_lte_work; /* Y = X + LTE */
239: PetscReal wltea, wlter;
241: PetscFunctionBegin;
242: if (!th->vec_sol_prev) {
243: *wlte = -1;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
246: if (!th->vec_lte_work) {
247: *wlte = -1;
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
250: if (ts->steprestart) {
251: /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */
252: PetscCall(VecAXPY(Y, 1, X));
253: } else {
254: /* Compute LTE using backward differences with non-constant time step */
255: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
256: PetscReal a = 1 + h_prev / h;
257: PetscScalar scal[3];
258: Vec vecs[3];
259: scal[0] = +1 / a;
260: scal[1] = -1 / (a - 1);
261: scal[2] = +1 / (a * (a - 1));
262: vecs[0] = th->X1;
263: vecs[1] = th->X0;
264: vecs[2] = th->vec_sol_prev;
265: PetscCall(VecCopy(X, Y));
266: PetscCall(VecMAXPY(Y, 3, scal, vecs));
267: }
268: PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
269: if (order) *order = 2;
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode TSInterpolate_Alpha(TS ts, PetscReal t, Vec X)
274: {
275: TS_Alpha *th = (TS_Alpha *)ts->data;
276: PetscReal dt = t - ts->ptime;
277: PetscReal Gamma = th->Gamma;
279: PetscFunctionBegin;
280: PetscCall(VecWAXPY(th->V1, -1.0, th->X0, ts->vec_sol));
281: PetscCall(VecAXPBY(th->V1, 1 - 1 / Gamma, 1 / (Gamma * ts->time_step), th->V0));
282: PetscCall(VecCopy(ts->vec_sol, X));
283: /* X = X + Gamma*dT*V1 */
284: PetscCall(VecAXPY(X, th->Gamma * dt, th->V1));
285: /* X = X + (1-Gamma)*dT*V0 */
286: PetscCall(VecAXPY(X, (1 - th->Gamma) * dt, th->V0));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
291: {
292: TS_Alpha *th = (TS_Alpha *)ts->data;
293: PetscReal ta = th->stage_time;
294: Vec Xa = th->Xa, Va = th->Va;
296: PetscFunctionBegin;
297: PetscCall(TSAlpha_StageVecs(ts, X));
298: /* F = Function(ta,Xa,Va) */
299: PetscCall(TSComputeIFunction(ts, ta, Xa, Va, F, PETSC_FALSE));
300: PetscCall(VecScale(F, th->scale_F));
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
305: {
306: TS_Alpha *th = (TS_Alpha *)ts->data;
307: PetscReal ta = th->stage_time;
308: Vec Xa = th->Xa, Va = th->Va;
309: PetscReal dVdX = th->shift_V;
311: PetscFunctionBegin;
312: /* J,P = Jacobian(ta,Xa,Va) */
313: PetscCall(TSComputeIJacobian(ts, ta, Xa, Va, dVdX, J, P, PETSC_FALSE));
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: static PetscErrorCode TSReset_Alpha(TS ts)
318: {
319: TS_Alpha *th = (TS_Alpha *)ts->data;
321: PetscFunctionBegin;
322: PetscCall(VecDestroy(&th->X0));
323: PetscCall(VecDestroy(&th->Xa));
324: PetscCall(VecDestroy(&th->X1));
325: PetscCall(VecDestroy(&th->V0));
326: PetscCall(VecDestroy(&th->Va));
327: PetscCall(VecDestroy(&th->V1));
328: PetscCall(VecDestroy(&th->vec_sol_prev));
329: PetscCall(VecDestroy(&th->vec_lte_work));
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode TSDestroy_Alpha(TS ts)
334: {
335: PetscFunctionBegin;
336: PetscCall(TSReset_Alpha(ts));
337: PetscCall(PetscFree(ts->data));
339: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", NULL));
340: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", NULL));
341: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", NULL));
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: static PetscErrorCode TSSetUp_Alpha(TS ts)
346: {
347: TS_Alpha *th = (TS_Alpha *)ts->data;
348: PetscBool match;
350: PetscFunctionBegin;
351: if (!th->X0) PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
352: PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
353: PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
354: PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
355: PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
356: PetscCall(VecDuplicate(ts->vec_sol, &th->V1));
358: PetscCall(TSGetAdapt(ts, &ts->adapt));
359: PetscCall(TSAdaptCandidatesClear(ts->adapt));
360: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
361: if (!match) {
362: if (!th->vec_sol_prev) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
363: if (!th->vec_lte_work) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work));
364: }
366: PetscCall(TSGetSNES(ts, &ts->snes));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
371: {
372: TS_Alpha *th = (TS_Alpha *)ts->data;
374: PetscFunctionBegin;
375: PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
376: {
377: PetscBool flg;
378: PetscReal radius = 1;
379: PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlphaSetRadius", radius, &radius, &flg));
380: if (flg) PetscCall(TSAlphaSetRadius(ts, radius));
381: PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlphaSetParams", th->Alpha_m, &th->Alpha_m, NULL));
382: PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlphaSetParams", th->Alpha_f, &th->Alpha_f, NULL));
383: PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlphaSetParams", th->Gamma, &th->Gamma, NULL));
384: PetscCall(TSAlphaSetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma));
385: }
386: PetscOptionsHeadEnd();
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
391: {
392: TS_Alpha *th = (TS_Alpha *)ts->data;
393: PetscBool iascii;
395: PetscFunctionBegin;
396: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
397: 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));
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts, PetscReal radius)
402: {
403: PetscReal alpha_m, alpha_f, gamma;
405: PetscFunctionBegin;
406: PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
407: alpha_m = (PetscReal)0.5 * (3 - radius) / (1 + radius);
408: alpha_f = 1 / (1 + radius);
409: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
410: PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode TSAlphaSetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
415: {
416: TS_Alpha *th = (TS_Alpha *)ts->data;
417: PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
418: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
420: PetscFunctionBegin;
421: th->Alpha_m = alpha_m;
422: th->Alpha_f = alpha_f;
423: th->Gamma = gamma;
424: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode TSAlphaGetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
429: {
430: TS_Alpha *th = (TS_Alpha *)ts->data;
432: PetscFunctionBegin;
433: if (alpha_m) *alpha_m = th->Alpha_m;
434: if (alpha_f) *alpha_f = th->Alpha_f;
435: if (gamma) *gamma = th->Gamma;
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: /*MC
440: TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method {cite}`jansen_2000` {cite}`chung1993` for first-order systems
442: Level: beginner
444: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
445: M*/
446: PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts)
447: {
448: TS_Alpha *th;
450: PetscFunctionBegin;
451: ts->ops->reset = TSReset_Alpha;
452: ts->ops->destroy = TSDestroy_Alpha;
453: ts->ops->view = TSView_Alpha;
454: ts->ops->setup = TSSetUp_Alpha;
455: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
456: ts->ops->step = TSStep_Alpha;
457: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
458: ts->ops->interpolate = TSInterpolate_Alpha;
459: ts->ops->resizeregister = TSResizeRegister_Alpha;
460: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
461: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
462: ts->default_adapt_type = TSADAPTNONE;
464: ts->usessnes = PETSC_TRUE;
466: PetscCall(PetscNew(&th));
467: ts->data = (void *)th;
469: th->Alpha_m = 0.5;
470: th->Alpha_f = 0.5;
471: th->Gamma = 0.5;
472: th->order = 2;
474: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", TSAlphaSetRadius_Alpha));
475: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", TSAlphaSetParams_Alpha));
476: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", TSAlphaGetParams_Alpha));
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: /*@
481: TSAlphaSetRadius - sets the desired spectral radius of the method for `TSALPHA`
482: (i.e. high-frequency numerical damping)
484: Logically Collective
486: Input Parameters:
487: + ts - timestepping context
488: - radius - the desired spectral radius
490: Options Database Key:
491: . -ts_alpha_radius <radius> - set alpha radius
493: Level: intermediate
495: Notes:
496: The algorithmic parameters $\alpha_m$ and $\alpha_f$ of the generalized-$\alpha$ method can
497: be computed in terms of a specified spectral radius $\rho$ in [0, 1] for infinite time step
498: in order to control high-frequency numerical damping\:
500: $$
501: \begin{align*}
502: \alpha_m = 0.5*(3-\rho)/(1+\rho) \\
503: \alpha_f = 1/(1+\rho)
504: \end{align*}
505: $$
507: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetParams()`, `TSAlphaGetParams()`
508: @*/
509: PetscErrorCode TSAlphaSetRadius(TS ts, PetscReal radius)
510: {
511: PetscFunctionBegin;
514: PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
515: PetscTryMethod(ts, "TSAlphaSetRadius_C", (TS, PetscReal), (ts, radius));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: /*@
520: TSAlphaSetParams - sets the algorithmic parameters for `TSALPHA`
522: Logically Collective
524: Input Parameters:
525: + ts - timestepping context
526: . alpha_m - algorithmic parameter
527: . alpha_f - algorithmic parameter
528: - gamma - algorithmic parameter
530: Options Database Keys:
531: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
532: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
533: - -ts_alpha_gamma <gamma> - set gamma
535: Level: advanced
537: Note:
538: Second-order accuracy can be obtained so long as\: $\gamma = 0.5 + \alpha_m - \alpha_f$
540: Unconditional stability requires\: $\alpha_m >= \alpha_f >= 0.5$
542: Backward Euler method is recovered with\: $\alpha_m = \alpha_f = \gamma = 1$
544: Use of this function is normally only required to hack `TSALPHA` to use a modified
545: integration scheme. Users should call `TSAlphaSetRadius()` to set the desired spectral radius
546: of the methods (i.e. high-frequency damping) in order so select optimal values for these
547: parameters.
549: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaGetParams()`
550: @*/
551: PetscErrorCode TSAlphaSetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
552: {
553: PetscFunctionBegin;
558: PetscTryMethod(ts, "TSAlphaSetParams_C", (TS, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma));
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: /*@
563: TSAlphaGetParams - gets the algorithmic parameters for `TSALPHA`
565: Not Collective
567: Input Parameter:
568: . ts - timestepping context
570: Output Parameters:
571: + alpha_m - algorithmic parameter
572: . alpha_f - algorithmic parameter
573: - gamma - algorithmic parameter
575: Level: advanced
577: Note:
578: Use of this function is normally only required to hack `TSALPHA` to use a modified
579: integration scheme. Users should call `TSAlphaSetRadius()` to set the high-frequency damping
580: (i.e. spectral radius of the method) in order so select optimal values for these parameters.
582: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
583: @*/
584: PetscErrorCode TSAlphaGetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
585: {
586: PetscFunctionBegin;
588: if (alpha_m) PetscAssertPointer(alpha_m, 2);
589: if (alpha_f) PetscAssertPointer(alpha_f, 3);
590: if (gamma) PetscAssertPointer(gamma, 4);
591: PetscUseMethod(ts, "TSAlphaGetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma));
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }