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 TSAlpha_StageTime(TS ts)
38: {
39: TS_Alpha *th = (TS_Alpha *)ts->data;
40: PetscReal t = ts->ptime;
41: PetscReal dt = ts->time_step;
42: PetscReal Alpha_m = th->Alpha_m;
43: PetscReal Alpha_f = th->Alpha_f;
44: PetscReal Gamma = th->Gamma;
46: PetscFunctionBegin;
47: th->stage_time = t + Alpha_f * dt;
48: th->shift_V = Alpha_m / (Alpha_f * Gamma * dt);
49: th->scale_F = 1 / Alpha_f;
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode TSAlpha_StageVecs(TS ts, Vec X)
54: {
55: TS_Alpha *th = (TS_Alpha *)ts->data;
56: Vec X1 = X, V1 = th->V1;
57: Vec Xa = th->Xa, Va = th->Va;
58: Vec X0 = th->X0, V0 = th->V0;
59: PetscReal dt = ts->time_step;
60: PetscReal Alpha_m = th->Alpha_m;
61: PetscReal Alpha_f = th->Alpha_f;
62: PetscReal Gamma = th->Gamma;
64: PetscFunctionBegin;
65: /* V1 = 1/(Gamma*dT)*(X1-X0) + (1-1/Gamma)*V0 */
66: PetscCall(VecWAXPY(V1, -1.0, X0, X1));
67: PetscCall(VecAXPBY(V1, 1 - 1 / Gamma, 1 / (Gamma * dt), V0));
68: /* Xa = X0 + Alpha_f*(X1-X0) */
69: PetscCall(VecWAXPY(Xa, -1.0, X0, X1));
70: PetscCall(VecAYPX(Xa, Alpha_f, X0));
71: /* Va = V0 + Alpha_m*(V1-V0) */
72: PetscCall(VecWAXPY(Va, -1.0, V0, V1));
73: PetscCall(VecAYPX(Va, Alpha_m, V0));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode TSAlpha_SNESSolve(TS ts, Vec b, Vec x)
78: {
79: PetscInt nits, lits;
81: PetscFunctionBegin;
82: PetscCall(SNESSolve(ts->snes, b, x));
83: PetscCall(SNESGetIterationNumber(ts->snes, &nits));
84: PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
85: ts->snes_its += nits;
86: ts->ksp_its += lits;
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: /*
91: Compute a consistent initial state for the generalized-alpha method.
92: - Solve two successive backward Euler steps with halved time step.
93: - Compute the initial time derivative using backward differences.
94: - If using adaptivity, estimate the LTE of the initial step.
95: */
96: static PetscErrorCode TSAlpha_Restart(TS ts, PetscBool *initok)
97: {
98: TS_Alpha *th = (TS_Alpha *)ts->data;
99: PetscReal time_step;
100: PetscReal alpha_m, alpha_f, gamma;
101: Vec X0 = ts->vec_sol, X1, X2 = th->X1;
102: PetscBool stageok;
104: PetscFunctionBegin;
105: PetscCall(VecDuplicate(X0, &X1));
107: /* Setup backward Euler with halved time step */
108: PetscCall(TSAlphaGetParams(ts, &alpha_m, &alpha_f, &gamma));
109: PetscCall(TSAlphaSetParams(ts, 1, 1, 1));
110: PetscCall(TSGetTimeStep(ts, &time_step));
111: ts->time_step = time_step / 2;
112: PetscCall(TSAlpha_StageTime(ts));
113: th->stage_time = ts->ptime;
114: PetscCall(VecZeroEntries(th->V0));
116: /* First BE step, (t0,X0) -> (t1,X1) */
117: th->stage_time += ts->time_step;
118: PetscCall(VecCopy(X0, th->X0));
119: PetscCall(TSPreStage(ts, th->stage_time));
120: PetscCall(VecCopy(th->X0, X1));
121: PetscCall(TSAlpha_SNESSolve(ts, NULL, X1));
122: PetscCall(TSPostStage(ts, th->stage_time, 0, &X1));
123: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X1, &stageok));
124: if (!stageok) goto finally;
126: /* Second BE step, (t1,X1) -> (t2,X2) */
127: th->stage_time += ts->time_step;
128: PetscCall(VecCopy(X1, th->X0));
129: PetscCall(TSPreStage(ts, th->stage_time));
130: PetscCall(VecCopy(th->X0, X2));
131: PetscCall(TSAlpha_SNESSolve(ts, NULL, X2));
132: PetscCall(TSPostStage(ts, th->stage_time, 0, &X2));
133: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, X2, &stageok));
134: if (!stageok) goto finally;
136: /* Compute V0 ~ dX/dt at t0 with backward differences */
137: PetscCall(VecZeroEntries(th->V0));
138: PetscCall(VecAXPY(th->V0, -3 / ts->time_step, X0));
139: PetscCall(VecAXPY(th->V0, +4 / ts->time_step, X1));
140: PetscCall(VecAXPY(th->V0, -1 / ts->time_step, X2));
142: /* Rough, lower-order estimate LTE of the initial step */
143: if (th->vec_lte_work) {
144: PetscCall(VecZeroEntries(th->vec_lte_work));
145: PetscCall(VecAXPY(th->vec_lte_work, +2, X2));
146: PetscCall(VecAXPY(th->vec_lte_work, -4, X1));
147: PetscCall(VecAXPY(th->vec_lte_work, +2, X0));
148: }
150: finally:
151: /* Revert TSAlpha to the initial state (t0,X0) */
152: if (initok) *initok = stageok;
153: PetscCall(TSSetTimeStep(ts, time_step));
154: PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
155: PetscCall(VecCopy(ts->vec_sol, th->X0));
157: PetscCall(VecDestroy(&X1));
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: static PetscErrorCode TSStep_Alpha(TS ts)
162: {
163: TS_Alpha *th = (TS_Alpha *)ts->data;
164: PetscInt rejections = 0;
165: PetscBool stageok, accept = PETSC_TRUE;
166: PetscReal next_time_step = ts->time_step;
168: PetscFunctionBegin;
169: PetscCall(PetscCitationsRegister(citation, &cited));
171: if (!ts->steprollback) {
172: if (th->vec_sol_prev) PetscCall(VecCopy(th->X0, th->vec_sol_prev));
173: PetscCall(VecCopy(ts->vec_sol, th->X0));
174: PetscCall(VecCopy(th->V1, th->V0));
175: }
177: th->status = TS_STEP_INCOMPLETE;
178: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
179: if (ts->steprestart) {
180: PetscCall(TSAlpha_Restart(ts, &stageok));
181: if (!stageok) goto reject_step;
182: }
184: PetscCall(TSAlpha_StageTime(ts));
185: PetscCall(VecCopy(th->X0, th->X1));
186: PetscCall(TSPreStage(ts, th->stage_time));
187: PetscCall(TSAlpha_SNESSolve(ts, NULL, th->X1));
188: PetscCall(TSPostStage(ts, th->stage_time, 0, &th->Xa));
189: PetscCall(TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->Xa, &stageok));
190: if (!stageok) goto reject_step;
192: th->status = TS_STEP_PENDING;
193: PetscCall(VecCopy(th->X1, ts->vec_sol));
194: PetscCall(TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
195: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
196: if (!accept) {
197: PetscCall(VecCopy(th->X0, ts->vec_sol));
198: ts->time_step = next_time_step;
199: goto reject_step;
200: }
202: ts->ptime += ts->time_step;
203: ts->time_step = next_time_step;
204: break;
206: reject_step:
207: ts->reject++;
208: accept = PETSC_FALSE;
209: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
210: ts->reason = TS_DIVERGED_STEP_REJECTED;
211: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
212: }
213: }
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
218: {
219: TS_Alpha *th = (TS_Alpha *)ts->data;
220: Vec X = th->X1; /* X = solution */
221: Vec Y = th->vec_lte_work; /* Y = X + LTE */
222: PetscReal wltea, wlter;
224: PetscFunctionBegin;
225: if (!th->vec_sol_prev) {
226: *wlte = -1;
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
229: if (!th->vec_lte_work) {
230: *wlte = -1;
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
233: if (ts->steprestart) {
234: /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */
235: PetscCall(VecAXPY(Y, 1, X));
236: } else {
237: /* Compute LTE using backward differences with non-constant time step */
238: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
239: PetscReal a = 1 + h_prev / h;
240: PetscScalar scal[3];
241: Vec vecs[3];
242: scal[0] = +1 / a;
243: scal[1] = -1 / (a - 1);
244: scal[2] = +1 / (a * (a - 1));
245: vecs[0] = th->X1;
246: vecs[1] = th->X0;
247: vecs[2] = th->vec_sol_prev;
248: PetscCall(VecCopy(X, Y));
249: PetscCall(VecMAXPY(Y, 3, scal, vecs));
250: }
251: PetscCall(TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter));
252: if (order) *order = 2;
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: static PetscErrorCode TSRollBack_Alpha(TS ts)
257: {
258: TS_Alpha *th = (TS_Alpha *)ts->data;
260: PetscFunctionBegin;
261: PetscCall(VecCopy(th->X0, ts->vec_sol));
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: static PetscErrorCode TSInterpolate_Alpha(TS ts, PetscReal t, Vec X)
266: {
267: TS_Alpha *th = (TS_Alpha *)ts->data;
268: PetscReal dt = t - ts->ptime;
270: PetscFunctionBegin;
271: PetscCall(VecCopy(ts->vec_sol, X));
272: PetscCall(VecAXPY(X, th->Gamma * dt, th->V1));
273: PetscCall(VecAXPY(X, (1 - th->Gamma) * dt, th->V0));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes, Vec X, Vec F, TS ts)
278: {
279: TS_Alpha *th = (TS_Alpha *)ts->data;
280: PetscReal ta = th->stage_time;
281: Vec Xa = th->Xa, Va = th->Va;
283: PetscFunctionBegin;
284: PetscCall(TSAlpha_StageVecs(ts, X));
285: /* F = Function(ta,Xa,Va) */
286: PetscCall(TSComputeIFunction(ts, ta, Xa, Va, F, PETSC_FALSE));
287: PetscCall(VecScale(F, th->scale_F));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes, PETSC_UNUSED Vec X, Mat J, Mat P, 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;
296: PetscReal dVdX = th->shift_V;
298: PetscFunctionBegin;
299: /* J,P = Jacobian(ta,Xa,Va) */
300: PetscCall(TSComputeIJacobian(ts, ta, Xa, Va, dVdX, J, P, PETSC_FALSE));
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: static PetscErrorCode TSReset_Alpha(TS ts)
305: {
306: TS_Alpha *th = (TS_Alpha *)ts->data;
308: PetscFunctionBegin;
309: PetscCall(VecDestroy(&th->X0));
310: PetscCall(VecDestroy(&th->Xa));
311: PetscCall(VecDestroy(&th->X1));
312: PetscCall(VecDestroy(&th->V0));
313: PetscCall(VecDestroy(&th->Va));
314: PetscCall(VecDestroy(&th->V1));
315: PetscCall(VecDestroy(&th->vec_sol_prev));
316: PetscCall(VecDestroy(&th->vec_lte_work));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: static PetscErrorCode TSDestroy_Alpha(TS ts)
321: {
322: PetscFunctionBegin;
323: PetscCall(TSReset_Alpha(ts));
324: PetscCall(PetscFree(ts->data));
326: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", NULL));
327: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", NULL));
328: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", NULL));
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: static PetscErrorCode TSSetUp_Alpha(TS ts)
333: {
334: TS_Alpha *th = (TS_Alpha *)ts->data;
335: PetscBool match;
337: PetscFunctionBegin;
338: PetscCall(VecDuplicate(ts->vec_sol, &th->X0));
339: PetscCall(VecDuplicate(ts->vec_sol, &th->Xa));
340: PetscCall(VecDuplicate(ts->vec_sol, &th->X1));
341: PetscCall(VecDuplicate(ts->vec_sol, &th->V0));
342: PetscCall(VecDuplicate(ts->vec_sol, &th->Va));
343: PetscCall(VecDuplicate(ts->vec_sol, &th->V1));
345: PetscCall(TSGetAdapt(ts, &ts->adapt));
346: PetscCall(TSAdaptCandidatesClear(ts->adapt));
347: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match));
348: if (!match) {
349: PetscCall(VecDuplicate(ts->vec_sol, &th->vec_sol_prev));
350: PetscCall(VecDuplicate(ts->vec_sol, &th->vec_lte_work));
351: }
353: PetscCall(TSGetSNES(ts, &ts->snes));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode TSSetFromOptions_Alpha(TS ts, PetscOptionItems *PetscOptionsObject)
358: {
359: TS_Alpha *th = (TS_Alpha *)ts->data;
361: PetscFunctionBegin;
362: PetscOptionsHeadBegin(PetscOptionsObject, "Generalized-Alpha ODE solver options");
363: {
364: PetscBool flg;
365: PetscReal radius = 1;
366: PetscCall(PetscOptionsReal("-ts_alpha_radius", "Spectral radius (high-frequency dissipation)", "TSAlphaSetRadius", radius, &radius, &flg));
367: if (flg) PetscCall(TSAlphaSetRadius(ts, radius));
368: PetscCall(PetscOptionsReal("-ts_alpha_alpha_m", "Algorithmic parameter alpha_m", "TSAlphaSetParams", th->Alpha_m, &th->Alpha_m, NULL));
369: PetscCall(PetscOptionsReal("-ts_alpha_alpha_f", "Algorithmic parameter alpha_f", "TSAlphaSetParams", th->Alpha_f, &th->Alpha_f, NULL));
370: PetscCall(PetscOptionsReal("-ts_alpha_gamma", "Algorithmic parameter gamma", "TSAlphaSetParams", th->Gamma, &th->Gamma, NULL));
371: PetscCall(TSAlphaSetParams(ts, th->Alpha_m, th->Alpha_f, th->Gamma));
372: }
373: PetscOptionsHeadEnd();
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: static PetscErrorCode TSView_Alpha(TS ts, PetscViewer viewer)
378: {
379: TS_Alpha *th = (TS_Alpha *)ts->data;
380: PetscBool iascii;
382: PetscFunctionBegin;
383: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
384: 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));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts, PetscReal radius)
389: {
390: PetscReal alpha_m, alpha_f, gamma;
392: PetscFunctionBegin;
393: PetscCheck(radius >= 0 && radius <= 1, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
394: alpha_m = (PetscReal)0.5 * (3 - radius) / (1 + radius);
395: alpha_f = 1 / (1 + radius);
396: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
397: PetscCall(TSAlphaSetParams(ts, alpha_m, alpha_f, gamma));
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: static PetscErrorCode TSAlphaSetParams_Alpha(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
402: {
403: TS_Alpha *th = (TS_Alpha *)ts->data;
404: PetscReal tol = 100 * PETSC_MACHINE_EPSILON;
405: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
407: PetscFunctionBegin;
408: th->Alpha_m = alpha_m;
409: th->Alpha_f = alpha_f;
410: th->Gamma = gamma;
411: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode TSAlphaGetParams_Alpha(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
416: {
417: TS_Alpha *th = (TS_Alpha *)ts->data;
419: PetscFunctionBegin;
420: if (alpha_m) *alpha_m = th->Alpha_m;
421: if (alpha_f) *alpha_f = th->Alpha_f;
422: if (gamma) *gamma = th->Gamma;
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: /*MC
427: TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method
428: for first-order systems
430: Level: beginner
432: References:
433: + * - K.E. Jansen, C.H. Whiting, G.M. Hulber, "A generalized-alpha
434: method for integrating the filtered Navier-Stokes equations with a
435: stabilized finite element method", Computer Methods in Applied
436: Mechanics and Engineering, 190, 305-319, 2000.
437: DOI: 10.1016/S0045-7825(00)00203-6.
438: - * - J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
439: Dynamics with Improved Numerical Dissipation: The Generalized-alpha
440: Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.
442: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetType()`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
443: M*/
444: PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts)
445: {
446: TS_Alpha *th;
448: PetscFunctionBegin;
449: ts->ops->reset = TSReset_Alpha;
450: ts->ops->destroy = TSDestroy_Alpha;
451: ts->ops->view = TSView_Alpha;
452: ts->ops->setup = TSSetUp_Alpha;
453: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
454: ts->ops->step = TSStep_Alpha;
455: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
456: ts->ops->rollback = TSRollBack_Alpha;
457: ts->ops->interpolate = TSInterpolate_Alpha;
458: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
459: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
460: ts->default_adapt_type = TSADAPTNONE;
462: ts->usessnes = PETSC_TRUE;
464: PetscCall(PetscNew(&th));
465: ts->data = (void *)th;
467: th->Alpha_m = 0.5;
468: th->Alpha_f = 0.5;
469: th->Gamma = 0.5;
470: th->order = 2;
472: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetRadius_C", TSAlphaSetRadius_Alpha));
473: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaSetParams_C", TSAlphaSetParams_Alpha));
474: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSAlphaGetParams_C", TSAlphaGetParams_Alpha));
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: /*@
479: TSAlphaSetRadius - sets the desired spectral radius of the method for `TSALPHA`
480: (i.e. high-frequency numerical damping)
482: Logically Collective
484: The algorithmic parameters \alpha_m and \alpha_f of the
485: generalized-\alpha method can be computed in terms of a specified
486: spectral radius \rho in [0,1] for infinite time step in order to
487: control high-frequency numerical damping:
488: \alpha_m = 0.5*(3-\rho)/(1+\rho)
489: \alpha_f = 1/(1+\rho)
491: Input Parameters:
492: + ts - timestepping context
493: - radius - the desired spectral radius
495: Options Database Key:
496: . -ts_alpha_radius <radius> - set alpha radius
498: Level: intermediate
500: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetParams()`, `TSAlphaGetParams()`
501: @*/
502: PetscErrorCode TSAlphaSetRadius(TS ts, PetscReal radius)
503: {
504: PetscFunctionBegin;
507: PetscCheck(radius >= 0 && radius <= 1, ((PetscObject)ts)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Radius %g not in range [0,1]", (double)radius);
508: PetscTryMethod(ts, "TSAlphaSetRadius_C", (TS, PetscReal), (ts, radius));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: /*@
513: TSAlphaSetParams - sets the algorithmic parameters for `TSALPHA`
515: Logically Collective
517: Second-order accuracy can be obtained so long as:
518: \gamma = 0.5 + alpha_m - alpha_f
520: Unconditional stability requires:
521: \alpha_m >= \alpha_f >= 0.5
523: Backward Euler method is recovered with:
524: \alpha_m = \alpha_f = gamma = 1
526: Input Parameters:
527: + ts - timestepping context
528: . alpha_m - algorithmic parameter
529: . alpha_f - algorithmic parameter
530: - gamma - algorithmic parameter
532: Options Database Keys:
533: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
534: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
535: - -ts_alpha_gamma <gamma> - set gamma
537: Level: advanced
539: Note:
540: Use of this function is normally only required to hack `TSALPHA` to
541: use a modified integration scheme. Users should call
542: `TSAlphaSetRadius()` to set the desired spectral radius of the methods
543: (i.e. high-frequency damping) in order so select optimal values for
544: these parameters.
546: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaGetParams()`
547: @*/
548: PetscErrorCode TSAlphaSetParams(TS ts, PetscReal alpha_m, PetscReal alpha_f, PetscReal gamma)
549: {
550: PetscFunctionBegin;
555: PetscTryMethod(ts, "TSAlphaSetParams_C", (TS, PetscReal, PetscReal, PetscReal), (ts, alpha_m, alpha_f, gamma));
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: /*@
560: TSAlphaGetParams - gets the algorithmic parameters for `TSALPHA`
562: Not Collective
564: Input Parameter:
565: . ts - timestepping context
567: Output Parameters:
568: + alpha_m - algorithmic parameter
569: . alpha_f - algorithmic parameter
570: - gamma - algorithmic parameter
572: Level: advanced
574: Note:
575: Use of this function is normally only required to hack `TSALPHA` to
576: use a modified integration scheme. Users should call
577: `TSAlphaSetRadius()` to set the high-frequency damping (i.e. spectral
578: radius of the method) in order so select optimal values for these
579: parameters.
581: .seealso: [](ch_ts), `TS`, `TSALPHA`, `TSAlphaSetRadius()`, `TSAlphaSetParams()`
582: @*/
583: PetscErrorCode TSAlphaGetParams(TS ts, PetscReal *alpha_m, PetscReal *alpha_f, PetscReal *gamma)
584: {
585: PetscFunctionBegin;
590: PetscUseMethod(ts, "TSAlphaGetParams_C", (TS, PetscReal *, PetscReal *, PetscReal *), (ts, alpha_m, alpha_f, gamma));
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }