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: /* We need to transfer X0 which will be copied into sol_prev */
38: static PetscErrorCode TSResizeRegister_Alpha(TS ts, PetscBool reg)
39: {
40: TS_Alpha *th = (TS_Alpha *)ts->data;
41: const char name[] = "ts:alpha:X0";
43: PetscFunctionBegin;
44: if (reg && th->vec_sol_prev) {
45: PetscCall(TSResizeRegisterVec(ts, name, th->X0));
46: } else if (!reg) {
47: PetscCall(TSResizeRetrieveVec(ts, name, &th->X0));
48: PetscCall(PetscObjectReference((PetscObject)th->X0));
49: }
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode TSAlpha_StageTime(TS ts)
54: {
55: TS_Alpha *th = (TS_Alpha *)ts->data;
56: PetscReal t = ts->ptime;
57: PetscReal dt = ts->time_step;
58: PetscReal Alpha_m = th->Alpha_m;
59: PetscReal Alpha_f = th->Alpha_f;
60: PetscReal Gamma = th->Gamma;
62: PetscFunctionBegin;
63: th->stage_time = t + Alpha_f * dt;
64: th->shift_V = Alpha_m / (Alpha_f * Gamma * dt);
65: th->scale_F = 1 / Alpha_f;
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode TSAlpha_StageVecs(TS ts, Vec X)
70: {
71: TS_Alpha *th = (TS_Alpha *)ts->data;
72: Vec X1 = X, V1 = th->V1;
73: Vec Xa = th->Xa, Va = th->Va;
74: Vec X0 = th->X0, V0 = th->V0;
75: PetscReal dt = ts->time_step;
76: PetscReal Alpha_m = th->Alpha_m;
77: PetscReal Alpha_f = th->Alpha_f;
78: PetscReal Gamma = th->Gamma;
80: PetscFunctionBegin;
81: /* V1 = 1/(Gamma*dT)*(X1-X0) + (1-1/Gamma)*V0 */
82: PetscCall(VecWAXPY(V1, -1.0, X0, X1));
83: PetscCall(VecAXPBY(V1, 1 - 1 / Gamma, 1 / (Gamma * dt), V0));
84: /* Xa = X0 + Alpha_f*(X1-X0) */
85: PetscCall(VecWAXPY(Xa, -1.0, X0, X1));
86: PetscCall(VecAYPX(Xa, Alpha_f, X0));
87: /* Va = V0 + Alpha_m*(V1-V0) */
88: PetscCall(VecWAXPY(Va, -1.0, V0, V1));
89: PetscCall(VecAYPX(Va, Alpha_m, V0));
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 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;
117: Vec X0 = ts->vec_sol, X1, X2 = th->X1;
118: PetscBool stageok;
120: PetscFunctionBegin;
121: PetscCall(VecDuplicate(X0, &X1));
123: /* Setup backward Euler with halved time step */
124: PetscCall(TSAlphaGetParams(ts, &alpha_m, &alpha_f, &gamma));
125: PetscCall(TSAlphaSetParams(ts, 1, 1, 1));
126: PetscCall(TSGetTimeStep(ts, &time_step));
127: ts->time_step = time_step / 2;
128: PetscCall(TSAlpha_StageTime(ts));
129: th->stage_time = ts->ptime;
130: PetscCall(VecZeroEntries(th->V0));
132: /* First BE step, (t0,X0) -> (t1,X1) */
133: th->stage_time += ts->time_step;
134: PetscCall(VecCopy(X0, th->X0));
135: PetscCall(TSPreStage(ts, th->stage_time));
136: PetscCall(VecCopy(th->X0, X1));
137: PetscCall(TSAlpha_SNESSolve(ts, NULL, X1));
138: PetscCall(TSPostStage(ts, th->stage_time, 0, &X1));
139: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
140: if (!stageok) goto finally;
142: /* Second BE step, (t1,X1) -> (t2,X2) */
143: th->stage_time += ts->time_step;
144: PetscCall(VecCopy(X1, th->X0));
145: PetscCall(TSPreStage(ts, th->stage_time));
146: PetscCall(VecCopy(th->X0, X2));
147: PetscCall(TSAlpha_SNESSolve(ts, NULL, X2));
148: PetscCall(TSPostStage(ts, th->stage_time, 0, &X2));
149: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X2, &stageok));
150: if (!stageok) goto finally;
152: /* Compute V0 ~ dX/dt at t0 with backward differences */
153: PetscCall(VecZeroEntries(th->V0));
154: PetscCall(VecAXPY(th->V0, -3 / ts->time_step, X0));
155: PetscCall(VecAXPY(th->V0, +4 / ts->time_step, X1));
156: PetscCall(VecAXPY(th->V0, -1 / ts->time_step, X2));
158: /* Rough, lower-order estimate LTE of the initial step */
159: if (th->vec_lte_work) {
160: PetscCall(VecZeroEntries(th->vec_lte_work));
161: PetscCall(VecAXPY(th->vec_lte_work, +2, X2));
162: PetscCall(VecAXPY(th->vec_lte_work, -4, X1));
163: PetscCall(VecAXPY(th->vec_lte_work, +2, X0));
164: }
166: finally:
167: /* Revert TSAlpha to the initial state (t0,X0) */
168: if (initok) *initok = stageok;
169: PetscCall(TSSetTimeStep(ts, time_step));
170: PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
171: PetscCall(VecCopy(ts->vec_sol, th->X0));
173: PetscCall(VecDestroy(&X1));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode TSStep_Alpha(TS ts)
178: {
179: TS_Alpha *th = (TS_Alpha *)ts->data;
180: PetscInt rejections = 0;
181: PetscBool stageok, accept = PETSC_TRUE;
182: PetscReal next_time_step = ts->time_step;
184: PetscFunctionBegin;
185: PetscCall(PetscCitationsRegister(citation, &cited));
187: if (!ts->steprollback) {
188: if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
189: PetscCall(VecCopy(ts->vec_sol, th->X0));
190: PetscCall(VecCopy(th->V1, th->V0));
191: }
193: th->status = TS_STEP_INCOMPLETE;
194: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
195: if (ts->steprestart) {
196: PetscCall(TSAlpha_Restart(ts, &stageok));
197: if (!stageok) goto reject_step;
198: }
200: PetscCall(TSAlpha_StageTime(ts));
201: PetscCall(VecCopy(th->X0, th->X1));
202: PetscCall(TSPreStage(ts, th->stage_time));
203: PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
204: PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
205: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
206: if (!stageok) goto reject_step;
208: th->status = TS_STEP_PENDING;
209: PetscCall(VecCopy(th->X1, ts->vec_sol));
210: PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
211: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
212: if (!accept) {
213: PetscCall(VecCopy(th->X0, ts->vec_sol));
214: ts->time_step = next_time_step;
215: goto reject_step;
216: }
218: ts->ptime += ts->time_step;
219: ts->time_step = next_time_step;
220: break;
222: reject_step:
223: ts->reject++;
224: accept = PETSC_FALSE;
225: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
226: ts->reason = TS_DIVERGED_STEP_REJECTED;
227: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
228: }
229: }
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
234: {
235: TS_Alpha *th = (TS_Alpha *)ts->data;
236: Vec X = th->X1; /* X = solution */
237: Vec Y = th->vec_lte_work; /* Y = X + LTE */
238: PetscReal wltea, wlter;
240: PetscFunctionBegin;
241: if (!th->vec_sol_prev) {
242: *wlte = -1;
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
245: if (!th->vec_lte_work) {
246: *wlte = -1;
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
249: if (ts->steprestart) {
250: /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */
251: PetscCall(VecAXPY(Y, 1, X));
252: } else {
253: /* Compute LTE using backward differences with non-constant time step */
254: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
255: PetscReal a = 1 + h_prev / h;
256: PetscScalar scal[3];
257: Vec vecs[3];
258: scal[0] = +1 / a;
259: scal[1] = -1 / (a - 1);
260: scal[2] = +1 / (a * (a - 1));
261: vecs[0] = th->X1;
262: vecs[1] = th->X0;
263: vecs[2] = th->vec_sol_prev;
264: PetscCall(VecCopy(X, Y));
265: PetscCall(VecMAXPY(Y, 3, scal, vecs));
266: }
267: PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
268: if (order) *order = 2;
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode TSRollBack_Alpha(TS ts)
273: {
274: TS_Alpha *th = (TS_Alpha *)ts->data;
276: PetscFunctionBegin;
277: PetscCall(VecCopy(th->X0, ts->vec_sol));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: static PetscErrorCode TSInterpolate_Alpha(TS ts, PetscReal t, Vec X)
282: {
283: TS_Alpha *th = (TS_Alpha *)ts->data;
284: PetscReal dt = t - ts->ptime;
286: PetscFunctionBegin;
287: PetscCall(VecCopy(ts->vec_sol, X));
288: PetscCall(VecAXPY(X, th->Gamma * dt, th->V1));
289: PetscCall(VecAXPY(X, (1 - th->Gamma) * dt, th->V0));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
294: {
295: TS_Alpha *th = (TS_Alpha *)ts->data;
296: PetscReal ta = th->stage_time;
297: Vec Xa = th->Xa, Va = th->Va;
299: PetscFunctionBegin;
300: PetscCall(TSAlpha_StageVecs(ts, X));
301: /* F = Function(ta,Xa,Va) */
302: PetscCall(TSComputeIFunction(ts, ta, Xa, Va, F, PETSC_FALSE));
303: PetscCall(VecScale(F, th->scale_F));
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, TS ts)
308: {
309: TS_Alpha *th = (TS_Alpha *)ts->data;
310: PetscReal ta = th->stage_time;
311: Vec Xa = th->Xa, Va = th->Va;
312: PetscReal dVdX = th->shift_V;
314: PetscFunctionBegin;
315: /* J,P = Jacobian(ta,Xa,Va) */
316: PetscCall(TSComputeIJacobian(ts, ta, Xa, Va, dVdX, J, P, PETSC_FALSE));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: static PetscErrorCode TSReset_Alpha(TS ts)
321: {
322: TS_Alpha *th = (TS_Alpha *)ts->data;
324: PetscFunctionBegin;
325: PetscCall(VecDestroy(&th->X0));
326: PetscCall(VecDestroy(&th->Xa));
327: PetscCall(VecDestroy(&th->X1));
328: PetscCall(VecDestroy(&th->V0));
329: PetscCall(VecDestroy(&th->Va));
330: PetscCall(VecDestroy(&th->V1));
331: PetscCall(VecDestroy(&th->vec_sol_prev));
332: PetscCall(VecDestroy(&th->vec_lte_work));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: static PetscErrorCode TSDestroy_Alpha(TS ts)
337: {
338: PetscFunctionBegin;
339: PetscCall(TSReset_Alpha(ts));
340: PetscCall(PetscFree(ts->data));
342: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", NULL));
343: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", NULL));
344: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", NULL));
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: static PetscErrorCode TSSetUp_Alpha(TS ts)
349: {
350: TS_Alpha *th = (TS_Alpha *)ts->data;
351: PetscBool match;
353: PetscFunctionBegin;
354: if (!th->X0) PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
355: PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
356: PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
357: PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
358: PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
359: PetscCall(VecDuplicate(ts->vec_sol, &th->V1));
361: PetscCall(TSGetAdapt(ts, &ts->adapt));
362: PetscCall(TSAdaptCandidatesClear(ts->adapt));
363: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
364: if (!match) {
365: PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
366: PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work));
367: }
369: PetscCall(TSGetSNES(ts, &ts->snes));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
374: {
375: TS_Alpha *th = (TS_Alpha *)ts->data;
377: PetscFunctionBegin;
378: PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
379: {
380: PetscBool flg;
381: PetscReal radius = 1;
382: PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlphaSetRadius", radius, &radius, &flg));
383: if (flg) PetscCall(TSAlphaSetRadius(ts, radius));
384: PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlphaSetParams", th->Alpha_m, &th->Alpha_m, NULL));
385: PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlphaSetParams", th->Alpha_f, &th->Alpha_f, NULL));
386: PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlphaSetParams", th->Gamma, &th->Gamma, NULL));
387: PetscCall(TSAlphaSetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma));
388: }
389: PetscOptionsHeadEnd();
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
394: {
395: TS_Alpha *th = (TS_Alpha *)ts->data;
396: PetscBool iascii;
398: PetscFunctionBegin;
399: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
400: 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));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts, PetscReal radius)
405: {
406: PetscReal alpha_m, alpha_f, gamma;
408: PetscFunctionBegin;
409: PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
410: alpha_m = (PetscReal)0.5 * (3 - radius) / (1 + radius);
411: alpha_f = 1 / (1 + radius);
412: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
413: PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static PetscErrorCode TSAlphaSetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
418: {
419: TS_Alpha *th = (TS_Alpha *)ts->data;
420: PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
421: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
423: PetscFunctionBegin;
424: th->Alpha_m = alpha_m;
425: th->Alpha_f = alpha_f;
426: th->Gamma = gamma;
427: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: static PetscErrorCode TSAlphaGetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
432: {
433: TS_Alpha *th = (TS_Alpha *)ts->data;
435: PetscFunctionBegin;
436: if (alpha_m) *alpha_m = th->Alpha_m;
437: if (alpha_f) *alpha_f = th->Alpha_f;
438: if (gamma) *gamma = th->Gamma;
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: /*MC
443: TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method {cite}`jansen_2000` {cite}`chung1993` for first-order systems
445: Level: beginner
447: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
448: M*/
449: PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts)
450: {
451: TS_Alpha *th;
453: PetscFunctionBegin;
454: ts->ops->reset = TSReset_Alpha;
455: ts->ops->destroy = TSDestroy_Alpha;
456: ts->ops->view = TSView_Alpha;
457: ts->ops->setup = TSSetUp_Alpha;
458: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
459: ts->ops->step = TSStep_Alpha;
460: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
461: ts->ops->rollback = TSRollBack_Alpha;
462: ts->ops->interpolate = TSInterpolate_Alpha;
463: ts->ops->resizeregister = TSResizeRegister_Alpha;
464: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
465: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
466: ts->default_adapt_type = TSADAPTNONE;
468: ts->usessnes = PETSC_TRUE;
470: PetscCall(PetscNew(&th));
471: ts->data = (void *)th;
473: th->Alpha_m = 0.5;
474: th->Alpha_f = 0.5;
475: th->Gamma = 0.5;
476: th->order = 2;
478: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", TSAlphaSetRadius_Alpha));
479: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", TSAlphaSetParams_Alpha));
480: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", TSAlphaGetParams_Alpha));
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: /*@
485: TSAlphaSetRadius - sets the desired spectral radius of the method for `TSALPHA`
486: (i.e. high-frequency numerical damping)
488: Logically Collective
490: Input Parameters:
491: + ts - timestepping context
492: - radius - the desired spectral radius
494: Options Database Key:
495: . -ts_alpha_radius <radius> - set alpha radius
497: Level: intermediate
499: Notes:
500: The algorithmic parameters $\alpha_m$ and $\alpha_f$ of the generalized-$\alpha$ method can
501: be computed in terms of a specified spectral radius $\rho$ in [0, 1] for infinite time step
502: in order to control high-frequency numerical damping\:
504: $$
505: \begin{align*}
506: \alpha_m = 0.5*(3-\rho)/(1+\rho) \\
507: \alpha_f = 1/(1+\rho)
508: \end{align*}
509: $$
511: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetParams()`, `TSAlphaGetParams()`
512: @*/
513: PetscErrorCode TSAlphaSetRadius(TS ts, PetscReal radius)
514: {
515: PetscFunctionBegin;
518: PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
519: PetscTryMethod(ts, "TSAlphaSetRadius_C", (TS, PetscReal), (ts, radius));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: /*@
524: TSAlphaSetParams - sets the algorithmic parameters for `TSALPHA`
526: Logically Collective
528: Input Parameters:
529: + ts - timestepping context
530: . alpha_m - algorithmic parameter
531: . alpha_f - algorithmic parameter
532: - gamma - algorithmic parameter
534: Options Database Keys:
535: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
536: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
537: - -ts_alpha_gamma <gamma> - set gamma
539: Level: advanced
541: Note:
542: Second-order accuracy can be obtained so long as\: $\gamma = 0.5 + \alpha_m - \alpha_f$
544: Unconditional stability requires\: $\alpha_m >= \alpha_f >= 0.5$
546: Backward Euler method is recovered with\: $\alpha_m = \alpha_f = \gamma = 1$
548: Use of this function is normally only required to hack `TSALPHA` to use a modified
549: integration scheme. Users should call `TSAlphaSetRadius()` to set the desired spectral radius
550: of the methods (i.e. high-frequency damping) in order so select optimal values for these
551: parameters.
553: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaGetParams()`
554: @*/
555: PetscErrorCode TSAlphaSetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
556: {
557: PetscFunctionBegin;
562: PetscTryMethod(ts, "TSAlphaSetParams_C", (TS, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma));
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: /*@
567: TSAlphaGetParams - gets the algorithmic parameters for `TSALPHA`
569: Not Collective
571: Input Parameter:
572: . ts - timestepping context
574: Output Parameters:
575: + alpha_m - algorithmic parameter
576: . alpha_f - algorithmic parameter
577: - gamma - algorithmic parameter
579: Level: advanced
581: Note:
582: Use of this function is normally only required to hack `TSALPHA` to use a modified
583: integration scheme. Users should call `TSAlphaSetRadius()` to set the high-frequency damping
584: (i.e. spectral radius of the method) in order so select optimal values for these parameters.
586: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
587: @*/
588: PetscErrorCode TSAlphaGetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
589: {
590: PetscFunctionBegin;
592: if (alpha_m) PetscAssertPointer(alpha_m, 2);
593: if (alpha_f) PetscAssertPointer(alpha_f, 3);
594: if (gamma) PetscAssertPointer(gamma, 4);
595: PetscUseMethod(ts, "TSAlphaGetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }