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;
278: PetscFunctionBegin;
279: PetscCall(VecCopy(ts->vec_sol, X));
280: PetscCall(VecAXPY(X, th->Gamma * dt, th->V1));
281: PetscCall(VecAXPY(X, (1 - th->Gamma) * dt, th->V0));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
286: {
287: TS_Alpha *th = (TS_Alpha *)ts->data;
288: PetscReal ta = th->stage_time;
289: Vec Xa = th->Xa, Va = th->Va;
291: PetscFunctionBegin;
292: PetscCall(TSAlpha_StageVecs(ts, X));
293: /* F = Function(ta,Xa,Va) */
294: PetscCall(TSComputeIFunction(ts, ta, Xa, Va, F, PETSC_FALSE));
295: PetscCall(VecScale(F, th->scale_F));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
300: {
301: TS_Alpha *th = (TS_Alpha *)ts->data;
302: PetscReal ta = th->stage_time;
303: Vec Xa = th->Xa, Va = th->Va;
304: PetscReal dVdX = th->shift_V;
306: PetscFunctionBegin;
307: /* J,P = Jacobian(ta,Xa,Va) */
308: PetscCall(TSComputeIJacobian(ts, ta, Xa, Va, dVdX, J, P, PETSC_FALSE));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: static PetscErrorCode TSReset_Alpha(TS ts)
313: {
314: TS_Alpha *th = (TS_Alpha *)ts->data;
316: PetscFunctionBegin;
317: PetscCall(VecDestroy(&th->X0));
318: PetscCall(VecDestroy(&th->Xa));
319: PetscCall(VecDestroy(&th->X1));
320: PetscCall(VecDestroy(&th->V0));
321: PetscCall(VecDestroy(&th->Va));
322: PetscCall(VecDestroy(&th->V1));
323: PetscCall(VecDestroy(&th->vec_sol_prev));
324: PetscCall(VecDestroy(&th->vec_lte_work));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: static PetscErrorCode TSDestroy_Alpha(TS ts)
329: {
330: PetscFunctionBegin;
331: PetscCall(TSReset_Alpha(ts));
332: PetscCall(PetscFree(ts->data));
334: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", NULL));
335: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", NULL));
336: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", NULL));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: static PetscErrorCode TSSetUp_Alpha(TS ts)
341: {
342: TS_Alpha *th = (TS_Alpha *)ts->data;
343: PetscBool match;
345: PetscFunctionBegin;
346: if (!th->X0) PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
347: PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
348: PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
349: PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
350: PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
351: PetscCall(VecDuplicate(ts->vec_sol, &th->V1));
353: PetscCall(TSGetAdapt(ts, &ts->adapt));
354: PetscCall(TSAdaptCandidatesClear(ts->adapt));
355: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
356: if (!match) {
357: if (!th->vec_sol_prev) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
358: if (!th->vec_lte_work) PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work));
359: }
361: PetscCall(TSGetSNES(ts, &ts->snes));
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
366: {
367: TS_Alpha *th = (TS_Alpha *)ts->data;
369: PetscFunctionBegin;
370: PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
371: {
372: PetscBool flg;
373: PetscReal radius = 1;
374: PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlphaSetRadius", radius, &radius, &flg));
375: if (flg) PetscCall(TSAlphaSetRadius(ts, radius));
376: PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlphaSetParams", th->Alpha_m, &th->Alpha_m, NULL));
377: PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlphaSetParams", th->Alpha_f, &th->Alpha_f, NULL));
378: PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlphaSetParams", th->Gamma, &th->Gamma, NULL));
379: PetscCall(TSAlphaSetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma));
380: }
381: PetscOptionsHeadEnd();
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
386: {
387: TS_Alpha *th = (TS_Alpha *)ts->data;
388: PetscBool iascii;
390: PetscFunctionBegin;
391: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
392: 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));
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts, PetscReal radius)
397: {
398: PetscReal alpha_m, alpha_f, gamma;
400: PetscFunctionBegin;
401: PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
402: alpha_m = (PetscReal)0.5 * (3 - radius) / (1 + radius);
403: alpha_f = 1 / (1 + radius);
404: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
405: PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: static PetscErrorCode TSAlphaSetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
410: {
411: TS_Alpha *th = (TS_Alpha *)ts->data;
412: PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
413: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
415: PetscFunctionBegin;
416: th->Alpha_m = alpha_m;
417: th->Alpha_f = alpha_f;
418: th->Gamma = gamma;
419: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: static PetscErrorCode TSAlphaGetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
424: {
425: TS_Alpha *th = (TS_Alpha *)ts->data;
427: PetscFunctionBegin;
428: if (alpha_m) *alpha_m = th->Alpha_m;
429: if (alpha_f) *alpha_f = th->Alpha_f;
430: if (gamma) *gamma = th->Gamma;
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*MC
435: TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method {cite}`jansen_2000` {cite}`chung1993` for first-order systems
437: Level: beginner
439: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
440: M*/
441: PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts)
442: {
443: TS_Alpha *th;
445: PetscFunctionBegin;
446: ts->ops->reset = TSReset_Alpha;
447: ts->ops->destroy = TSDestroy_Alpha;
448: ts->ops->view = TSView_Alpha;
449: ts->ops->setup = TSSetUp_Alpha;
450: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
451: ts->ops->step = TSStep_Alpha;
452: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
453: ts->ops->interpolate = TSInterpolate_Alpha;
454: ts->ops->resizeregister = TSResizeRegister_Alpha;
455: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
456: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
457: ts->default_adapt_type = TSADAPTNONE;
459: ts->usessnes = PETSC_TRUE;
461: PetscCall(PetscNew(&th));
462: ts->data = (void *)th;
464: th->Alpha_m = 0.5;
465: th->Alpha_f = 0.5;
466: th->Gamma = 0.5;
467: th->order = 2;
469: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", TSAlphaSetRadius_Alpha));
470: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", TSAlphaSetParams_Alpha));
471: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", TSAlphaGetParams_Alpha));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: /*@
476: TSAlphaSetRadius - sets the desired spectral radius of the method for `TSALPHA`
477: (i.e. high-frequency numerical damping)
479: Logically Collective
481: Input Parameters:
482: + ts - timestepping context
483: - radius - the desired spectral radius
485: Options Database Key:
486: . -ts_alpha_radius <radius> - set alpha radius
488: Level: intermediate
490: Notes:
491: The algorithmic parameters $\alpha_m$ and $\alpha_f$ of the generalized-$\alpha$ method can
492: be computed in terms of a specified spectral radius $\rho$ in [0, 1] for infinite time step
493: in order to control high-frequency numerical damping\:
495: $$
496: \begin{align*}
497: \alpha_m = 0.5*(3-\rho)/(1+\rho) \\
498: \alpha_f = 1/(1+\rho)
499: \end{align*}
500: $$
502: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetParams()`, `TSAlphaGetParams()`
503: @*/
504: PetscErrorCode TSAlphaSetRadius(TS ts, PetscReal radius)
505: {
506: PetscFunctionBegin;
509: PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
510: PetscTryMethod(ts, "TSAlphaSetRadius_C", (TS, PetscReal), (ts, radius));
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: /*@
515: TSAlphaSetParams - sets the algorithmic parameters for `TSALPHA`
517: Logically Collective
519: Input Parameters:
520: + ts - timestepping context
521: . alpha_m - algorithmic parameter
522: . alpha_f - algorithmic parameter
523: - gamma - algorithmic parameter
525: Options Database Keys:
526: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
527: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
528: - -ts_alpha_gamma <gamma> - set gamma
530: Level: advanced
532: Note:
533: Second-order accuracy can be obtained so long as\: $\gamma = 0.5 + \alpha_m - \alpha_f$
535: Unconditional stability requires\: $\alpha_m >= \alpha_f >= 0.5$
537: Backward Euler method is recovered with\: $\alpha_m = \alpha_f = \gamma = 1$
539: Use of this function is normally only required to hack `TSALPHA` to use a modified
540: integration scheme. Users should call `TSAlphaSetRadius()` to set the desired spectral radius
541: of the methods (i.e. high-frequency damping) in order so select optimal values for these
542: parameters.
544: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaGetParams()`
545: @*/
546: PetscErrorCode TSAlphaSetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
547: {
548: PetscFunctionBegin;
553: PetscTryMethod(ts, "TSAlphaSetParams_C", (TS, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: /*@
558: TSAlphaGetParams - gets the algorithmic parameters for `TSALPHA`
560: Not Collective
562: Input Parameter:
563: . ts - timestepping context
565: Output Parameters:
566: + alpha_m - algorithmic parameter
567: . alpha_f - algorithmic parameter
568: - gamma - algorithmic parameter
570: Level: advanced
572: Note:
573: Use of this function is normally only required to hack `TSALPHA` to use a modified
574: integration scheme. Users should call `TSAlphaSetRadius()` to set the high-frequency damping
575: (i.e. spectral radius of the method) in order so select optimal values for these parameters.
577: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
578: @*/
579: PetscErrorCode TSAlphaGetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
580: {
581: PetscFunctionBegin;
583: if (alpha_m) PetscAssertPointer(alpha_m, 2);
584: if (alpha_f) PetscAssertPointer(alpha_f, 3);
585: if (gamma) PetscAssertPointer(gamma, 4);
586: PetscUseMethod(ts, "TSAlphaGetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }