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[] =
9: "@article{Jansen2000,\n"
10: " title = {A generalized-$\\alpha$ method for integrating the filtered {N}avier--{S}tokes equations with a stabilized finite element method},\n"
11: " author = {Kenneth E. Jansen and Christian H. Whiting and Gregory M. Hulbert},\n"
12: " journal = {Computer Methods in Applied Mechanics and Engineering},\n"
13: " volume = {190},\n"
14: " number = {3--4},\n"
15: " pages = {305--319},\n"
16: " year = {2000},\n"
17: " issn = {0045-7825},\n"
18: " doi = {http://dx.doi.org/10.1016/S0045-7825(00)00203-6}\n}\n";
20: typedef struct {
21: PetscReal stage_time;
22: PetscReal shift_V;
23: PetscReal scale_F;
24: Vec X0,Xa,X1;
25: Vec V0,Va,V1;
27: PetscReal Alpha_m;
28: PetscReal Alpha_f;
29: PetscReal Gamma;
30: PetscInt order;
32: Vec vec_sol_prev;
33: Vec vec_lte_work;
35: TSStepStatus status;
36: } TS_Alpha;
38: static PetscErrorCode TSAlpha_StageTime(TS ts)
39: {
40: TS_Alpha *th = (TS_Alpha*)ts->data;
41: PetscReal t = ts->ptime;
42: PetscReal dt = ts->time_step;
43: PetscReal Alpha_m = th->Alpha_m;
44: PetscReal Alpha_f = th->Alpha_f;
45: PetscReal Gamma = th->Gamma;
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: return 0;
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: /* V1 = 1/(Gamma*dT)*(X1-X0) + (1-1/Gamma)*V0 */
65: VecWAXPY(V1,-1.0,X0,X1);
66: VecAXPBY(V1,1-1/Gamma,1/(Gamma*dt),V0);
67: /* Xa = X0 + Alpha_f*(X1-X0) */
68: VecWAXPY(Xa,-1.0,X0,X1);
69: VecAYPX(Xa,Alpha_f,X0);
70: /* Va = V0 + Alpha_m*(V1-V0) */
71: VecWAXPY(Va,-1.0,V0,V1);
72: VecAYPX(Va,Alpha_m,V0);
73: return 0;
74: }
76: static PetscErrorCode TSAlpha_SNESSolve(TS ts,Vec b,Vec x)
77: {
78: PetscInt nits,lits;
80: SNESSolve(ts->snes,b,x);
81: SNESGetIterationNumber(ts->snes,&nits);
82: SNESGetLinearSolveIterations(ts->snes,&lits);
83: ts->snes_its += nits; 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) {
175: if (ts->steprestart) {
176: TSAlpha_Restart(ts,&stageok);
177: if (!stageok) goto reject_step;
178: }
180: TSAlpha_StageTime(ts);
181: VecCopy(th->X0,th->X1);
182: TSPreStage(ts,th->stage_time);
183: TSAlpha_SNESSolve(ts,NULL,th->X1);
184: TSPostStage(ts,th->stage_time,0,&th->Xa);
185: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->Xa,&stageok);
186: if (!stageok) goto reject_step;
188: th->status = TS_STEP_PENDING;
189: VecCopy(th->X1,ts->vec_sol);
190: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
191: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
192: if (!accept) {
193: VecCopy(th->X0,ts->vec_sol);
194: ts->time_step = next_time_step;
195: goto reject_step;
196: }
198: ts->ptime += ts->time_step;
199: ts->time_step = next_time_step;
200: break;
202: reject_step:
203: ts->reject++; 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=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
207: }
209: }
210: return 0;
211: }
213: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
214: {
215: TS_Alpha *th = (TS_Alpha*)ts->data;
216: Vec X = th->X1; /* X = solution */
217: Vec Y = th->vec_lte_work; /* Y = X + LTE */
218: PetscReal wltea,wlter;
220: if (!th->vec_sol_prev) {*wlte = -1; return 0;}
221: if (!th->vec_lte_work) {*wlte = -1; return 0;}
222: if (ts->steprestart) {
223: /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */
224: VecAXPY(Y,1,X);
225: } else {
226: /* Compute LTE using backward differences with non-constant time step */
227: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
228: PetscReal a = 1 + h_prev/h;
229: PetscScalar scal[3]; Vec vecs[3];
230: scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
231: vecs[0] = th->X1; vecs[1] = th->X0; vecs[2] = th->vec_sol_prev;
232: VecCopy(X,Y);
233: VecMAXPY(Y,3,scal,vecs);
234: }
235: TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
236: if (order) *order = 2;
237: return 0;
238: }
240: static PetscErrorCode TSRollBack_Alpha(TS ts)
241: {
242: TS_Alpha *th = (TS_Alpha*)ts->data;
244: VecCopy(th->X0,ts->vec_sol);
245: return 0;
246: }
248: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X)
249: {
250: TS_Alpha *th = (TS_Alpha*)ts->data;
251: PetscReal dt = t - ts->ptime;
253: VecCopy(ts->vec_sol,X);
254: VecAXPY(X,th->Gamma*dt,th->V1);
255: VecAXPY(X,(1-th->Gamma)*dt,th->V0);
256: return 0;
257: }
259: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
260: {
261: TS_Alpha *th = (TS_Alpha*)ts->data;
262: PetscReal ta = th->stage_time;
263: Vec Xa = th->Xa, Va = th->Va;
265: TSAlpha_StageVecs(ts,X);
266: /* F = Function(ta,Xa,Va) */
267: TSComputeIFunction(ts,ta,Xa,Va,F,PETSC_FALSE);
268: VecScale(F,th->scale_F);
269: return 0;
270: }
272: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,TS ts)
273: {
274: TS_Alpha *th = (TS_Alpha*)ts->data;
275: PetscReal ta = th->stage_time;
276: Vec Xa = th->Xa, Va = th->Va;
277: PetscReal dVdX = th->shift_V;
279: /* J,P = Jacobian(ta,Xa,Va) */
280: TSComputeIJacobian(ts,ta,Xa,Va,dVdX,J,P,PETSC_FALSE);
281: return 0;
282: }
284: static PetscErrorCode TSReset_Alpha(TS ts)
285: {
286: TS_Alpha *th = (TS_Alpha*)ts->data;
288: VecDestroy(&th->X0);
289: VecDestroy(&th->Xa);
290: VecDestroy(&th->X1);
291: VecDestroy(&th->V0);
292: VecDestroy(&th->Va);
293: VecDestroy(&th->V1);
294: VecDestroy(&th->vec_sol_prev);
295: VecDestroy(&th->vec_lte_work);
296: return 0;
297: }
299: static PetscErrorCode TSDestroy_Alpha(TS ts)
300: {
301: TSReset_Alpha(ts);
302: PetscFree(ts->data);
304: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetRadius_C",NULL);
305: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetParams_C",NULL);
306: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaGetParams_C",NULL);
307: return 0;
308: }
310: static PetscErrorCode TSSetUp_Alpha(TS ts)
311: {
312: TS_Alpha *th = (TS_Alpha*)ts->data;
313: PetscBool match;
315: VecDuplicate(ts->vec_sol,&th->X0);
316: VecDuplicate(ts->vec_sol,&th->Xa);
317: VecDuplicate(ts->vec_sol,&th->X1);
318: VecDuplicate(ts->vec_sol,&th->V0);
319: VecDuplicate(ts->vec_sol,&th->Va);
320: VecDuplicate(ts->vec_sol,&th->V1);
322: TSGetAdapt(ts,&ts->adapt);
323: TSAdaptCandidatesClear(ts->adapt);
324: PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
325: if (!match) {
326: VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
327: VecDuplicate(ts->vec_sol,&th->vec_lte_work);
328: }
330: TSGetSNES(ts,&ts->snes);
331: return 0;
332: }
334: static PetscErrorCode TSSetFromOptions_Alpha(PetscOptionItems *PetscOptionsObject,TS ts)
335: {
336: TS_Alpha *th = (TS_Alpha*)ts->data;
338: PetscOptionsHead(PetscOptionsObject,"Generalized-Alpha ODE solver options");
339: {
340: PetscBool flg;
341: PetscReal radius = 1;
342: PetscOptionsReal("-ts_alpha_radius","Spectral radius (high-frequency dissipation)","TSAlphaSetRadius",radius,&radius,&flg);
343: if (flg) TSAlphaSetRadius(ts,radius);
344: PetscOptionsReal("-ts_alpha_alpha_m","Algorithmic parameter alpha_m","TSAlphaSetParams",th->Alpha_m,&th->Alpha_m,NULL);
345: PetscOptionsReal("-ts_alpha_alpha_f","Algorithmic parameter alpha_f","TSAlphaSetParams",th->Alpha_f,&th->Alpha_f,NULL);
346: PetscOptionsReal("-ts_alpha_gamma","Algorithmic parameter gamma","TSAlphaSetParams",th->Gamma,&th->Gamma,NULL);
347: TSAlphaSetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma);
348: }
349: PetscOptionsTail();
350: return 0;
351: }
353: static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer)
354: {
355: TS_Alpha *th = (TS_Alpha*)ts->data;
356: PetscBool iascii;
358: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
359: if (iascii) {
360: PetscViewerASCIIPrintf(viewer," Alpha_m=%g, Alpha_f=%g, Gamma=%g\n",(double)th->Alpha_m,(double)th->Alpha_f,(double)th->Gamma);
361: }
362: return 0;
363: }
365: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts,PetscReal radius)
366: {
367: PetscReal alpha_m,alpha_f,gamma;
370: alpha_m = (PetscReal)0.5*(3-radius)/(1+radius);
371: alpha_f = 1/(1+radius);
372: gamma = (PetscReal)0.5 + alpha_m - alpha_f;
373: TSAlphaSetParams(ts,alpha_m,alpha_f,gamma);
374: return 0;
375: }
377: static PetscErrorCode TSAlphaSetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma)
378: {
379: TS_Alpha *th = (TS_Alpha*)ts->data;
380: PetscReal tol = 100*PETSC_MACHINE_EPSILON;
381: PetscReal res = ((PetscReal)0.5 + alpha_m - alpha_f) - gamma;
383: th->Alpha_m = alpha_m;
384: th->Alpha_f = alpha_f;
385: th->Gamma = gamma;
386: th->order = (PetscAbsReal(res) < tol) ? 2 : 1;
387: return 0;
388: }
390: static PetscErrorCode TSAlphaGetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma)
391: {
392: TS_Alpha *th = (TS_Alpha*)ts->data;
394: if (alpha_m) *alpha_m = th->Alpha_m;
395: if (alpha_f) *alpha_f = th->Alpha_f;
396: if (gamma) *gamma = th->Gamma;
397: return 0;
398: }
400: /*MC
401: TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method
402: for first-order systems
404: Level: beginner
406: References:
407: + * - K.E. Jansen, C.H. Whiting, G.M. Hulber, "A generalized-alpha
408: method for integrating the filtered Navier-Stokes equations with a
409: stabilized finite element method", Computer Methods in Applied
410: Mechanics and Engineering, 190, 305-319, 2000.
411: DOI: 10.1016/S0045-7825(00)00203-6.
412: - * - J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
413: Dynamics with Improved Numerical Dissipation: The Generalized-alpha
414: Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.
416: .seealso: TS, TSCreate(), TSSetType(), TSAlphaSetRadius(), TSAlphaSetParams()
417: M*/
418: PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts)
419: {
420: TS_Alpha *th;
422: ts->ops->reset = TSReset_Alpha;
423: ts->ops->destroy = TSDestroy_Alpha;
424: ts->ops->view = TSView_Alpha;
425: ts->ops->setup = TSSetUp_Alpha;
426: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
427: ts->ops->step = TSStep_Alpha;
428: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
429: ts->ops->rollback = TSRollBack_Alpha;
430: ts->ops->interpolate = TSInterpolate_Alpha;
431: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
432: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
433: ts->default_adapt_type = TSADAPTNONE;
435: ts->usessnes = PETSC_TRUE;
437: PetscNewLog(ts,&th);
438: ts->data = (void*)th;
440: th->Alpha_m = 0.5;
441: th->Alpha_f = 0.5;
442: th->Gamma = 0.5;
443: th->order = 2;
445: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetRadius_C",TSAlphaSetRadius_Alpha);
446: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetParams_C",TSAlphaSetParams_Alpha);
447: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaGetParams_C",TSAlphaGetParams_Alpha);
448: return 0;
449: }
451: /*@
452: TSAlphaSetRadius - sets the desired spectral radius of the method
453: (i.e. high-frequency numerical damping)
455: Logically Collective on TS
457: The algorithmic parameters \alpha_m and \alpha_f of the
458: generalized-\alpha method can be computed in terms of a specified
459: spectral radius \rho in [0,1] for infinite time step in order to
460: control high-frequency numerical damping:
461: \alpha_m = 0.5*(3-\rho)/(1+\rho)
462: \alpha_f = 1/(1+\rho)
464: Input Parameters:
465: + ts - timestepping context
466: - radius - the desired spectral radius
468: Options Database:
469: . -ts_alpha_radius <radius> - set alpha radius
471: Level: intermediate
473: .seealso: TSAlphaSetParams(), TSAlphaGetParams()
474: @*/
475: PetscErrorCode TSAlphaSetRadius(TS ts,PetscReal radius)
476: {
480: PetscTryMethod(ts,"TSAlphaSetRadius_C",(TS,PetscReal),(ts,radius));
481: return 0;
482: }
484: /*@
485: TSAlphaSetParams - sets the algorithmic parameters for TSALPHA
487: Logically Collective on TS
489: Second-order accuracy can be obtained so long as:
490: \gamma = 0.5 + alpha_m - alpha_f
492: Unconditional stability requires:
493: \alpha_m >= \alpha_f >= 0.5
495: Backward Euler method is recovered with:
496: \alpha_m = \alpha_f = gamma = 1
498: Input Parameters:
499: + ts - timestepping context
500: . \alpha_m - algorithmic parameter
501: . \alpha_f - algorithmic parameter
502: - \gamma - algorithmic parameter
504: Options Database:
505: + -ts_alpha_alpha_m <alpha_m> - set alpha_m
506: . -ts_alpha_alpha_f <alpha_f> - set alpha_f
507: - -ts_alpha_gamma <gamma> - set gamma
509: Note:
510: Use of this function is normally only required to hack TSALPHA to
511: use a modified integration scheme. Users should call
512: TSAlphaSetRadius() to set the desired spectral radius of the methods
513: (i.e. high-frequency damping) in order so select optimal values for
514: these parameters.
516: Level: advanced
518: .seealso: TSAlphaSetRadius(), TSAlphaGetParams()
519: @*/
520: PetscErrorCode TSAlphaSetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma)
521: {
526: PetscTryMethod(ts,"TSAlphaSetParams_C",(TS,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma));
527: return 0;
528: }
530: /*@
531: TSAlphaGetParams - gets the algorithmic parameters for TSALPHA
533: Not Collective
535: Input Parameter:
536: . ts - timestepping context
538: Output Parameters:
539: + \alpha_m - algorithmic parameter
540: . \alpha_f - algorithmic parameter
541: - \gamma - algorithmic parameter
543: Note:
544: Use of this function is normally only required to hack TSALPHA to
545: use a modified integration scheme. Users should call
546: TSAlphaSetRadius() to set the high-frequency damping (i.e. spectral
547: radius of the method) in order so select optimal values for these
548: parameters.
550: Level: advanced
552: .seealso: TSAlphaSetRadius(), TSAlphaSetParams()
553: @*/
554: PetscErrorCode TSAlphaGetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma)
555: {
560: PetscUseMethod(ts,"TSAlphaGetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma));
561: return 0;
562: }