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