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;
48: th->stage_time = t + Alpha_f*dt;
49: th->shift_V = Alpha_m/(Alpha_f*Gamma*dt);
50: th->scale_F = 1/Alpha_f;
51: return(0);
52: }
54: static PetscErrorCode TSAlpha_StageVecs(TS ts,Vec X)
55: {
56: TS_Alpha *th = (TS_Alpha*)ts->data;
57: Vec X1 = X, V1 = th->V1;
58: Vec Xa = th->Xa, Va = th->Va;
59: Vec X0 = th->X0, V0 = th->V0;
60: PetscReal dt = ts->time_step;
61: PetscReal Alpha_m = th->Alpha_m;
62: PetscReal Alpha_f = th->Alpha_f;
63: PetscReal Gamma = th->Gamma;
67: /* V1 = 1/(Gamma*dT)*(X1-X0) + (1-1/Gamma)*V0 */
68: VecWAXPY(V1,-1.0,X0,X1);
69: VecAXPBY(V1,1-1/Gamma,1/(Gamma*dt),V0);
70: /* Xa = X0 + Alpha_f*(X1-X0) */
71: VecWAXPY(Xa,-1.0,X0,X1);
72: VecAYPX(Xa,Alpha_f,X0);
73: /* Va = V0 + Alpha_m*(V1-V0) */
74: VecWAXPY(Va,-1.0,V0,V1);
75: VecAYPX(Va,Alpha_m,V0);
76: return(0);
77: }
79: static PetscErrorCode TSAlpha_SNESSolve(TS ts,Vec b,Vec x)
80: {
81: PetscInt nits,lits;
85: SNESSolve(ts->snes,b,x);
86: SNESGetIterationNumber(ts->snes,&nits);
87: SNESGetLinearSolveIterations(ts->snes,&lits);
88: ts->snes_its += nits; ts->ksp_its += lits;
89: return(0);
90: }
92: /*
93: Compute a consistent initial state for the generalized-alpha method.
94: - Solve two successive backward Euler steps with halved time step.
95: - Compute the initial time derivative using backward differences.
96: - If using adaptivity, estimate the LTE of the initial step.
97: */
98: static PetscErrorCode TSAlpha_Restart(TS ts,PetscBool *initok)
99: {
100: TS_Alpha *th = (TS_Alpha*)ts->data;
101: PetscReal time_step;
102: PetscReal alpha_m,alpha_f,gamma;
103: Vec X0 = ts->vec_sol, X1, X2 = th->X1;
104: PetscBool stageok;
108: VecDuplicate(X0,&X1);
110: /* Setup backward Euler with halved time step */
111: TSAlphaGetParams(ts,&alpha_m,&alpha_f,&gamma);
112: TSAlphaSetParams(ts,1,1,1);
113: TSGetTimeStep(ts,&time_step);
114: ts->time_step = time_step/2;
115: TSAlpha_StageTime(ts);
116: th->stage_time = ts->ptime;
117: VecZeroEntries(th->V0);
119: /* First BE step, (t0,X0) -> (t1,X1) */
120: th->stage_time += ts->time_step;
121: VecCopy(X0,th->X0);
122: TSPreStage(ts,th->stage_time);
123: VecCopy(th->X0,X1);
124: TSAlpha_SNESSolve(ts,NULL,X1);
125: TSPostStage(ts,th->stage_time,0,&X1);
126: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X1,&stageok);
127: if (!stageok) goto finally;
129: /* Second BE step, (t1,X1) -> (t2,X2) */
130: th->stage_time += ts->time_step;
131: VecCopy(X1,th->X0);
132: TSPreStage(ts,th->stage_time);
133: VecCopy(th->X0,X2);
134: TSAlpha_SNESSolve(ts,NULL,X2);
135: TSPostStage(ts,th->stage_time,0,&X2);
136: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,X2,&stageok);
137: if (!stageok) goto finally;
139: /* Compute V0 ~ dX/dt at t0 with backward differences */
140: VecZeroEntries(th->V0);
141: VecAXPY(th->V0,-3/ts->time_step,X0);
142: VecAXPY(th->V0,+4/ts->time_step,X1);
143: VecAXPY(th->V0,-1/ts->time_step,X2);
145: /* Rough, lower-order estimate LTE of the initial step */
146: if (th->vec_lte_work) {
147: VecZeroEntries(th->vec_lte_work);
148: VecAXPY(th->vec_lte_work,+2,X2);
149: VecAXPY(th->vec_lte_work,-4,X1);
150: VecAXPY(th->vec_lte_work,+2,X0);
151: }
153: finally:
154: /* Revert TSAlpha to the initial state (t0,X0) */
155: if (initok) *initok = stageok;
156: TSSetTimeStep(ts,time_step);
157: TSAlphaSetParams(ts,alpha_m,alpha_f,gamma);
158: VecCopy(ts->vec_sol,th->X0);
160: VecDestroy(&X1);
161: return(0);
162: }
164: static PetscErrorCode TSStep_Alpha(TS ts)
165: {
166: TS_Alpha *th = (TS_Alpha*)ts->data;
167: PetscInt rejections = 0;
168: PetscBool stageok,accept = PETSC_TRUE;
169: PetscReal next_time_step = ts->time_step;
173: PetscCitationsRegister(citation,&cited);
175: if (!ts->steprollback) {
176: if (th->vec_sol_prev) { VecCopy(th->X0,th->vec_sol_prev); }
177: VecCopy(ts->vec_sol,th->X0);
178: VecCopy(th->V1,th->V0);
179: }
181: th->status = TS_STEP_INCOMPLETE;
182: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
184: if (ts->steprestart) {
185: TSAlpha_Restart(ts,&stageok);
186: if (!stageok) goto reject_step;
187: }
189: TSAlpha_StageTime(ts);
190: VecCopy(th->X0,th->X1);
191: TSPreStage(ts,th->stage_time);
192: TSAlpha_SNESSolve(ts,NULL,th->X1);
193: TSPostStage(ts,th->stage_time,0,&th->Xa);
194: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->Xa,&stageok);
195: if (!stageok) goto reject_step;
197: th->status = TS_STEP_PENDING;
198: VecCopy(th->X1,ts->vec_sol);
199: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
200: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
201: if (!accept) {
202: VecCopy(th->X0,ts->vec_sol);
203: ts->time_step = next_time_step;
204: goto reject_step;
205: }
207: ts->ptime += ts->time_step;
208: ts->time_step = next_time_step;
209: break;
211: reject_step:
212: ts->reject++; accept = PETSC_FALSE;
213: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
214: ts->reason = TS_DIVERGED_STEP_REJECTED;
215: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
216: }
218: }
219: return(0);
220: }
222: static PetscErrorCode TSEvaluateWLTE_Alpha(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
223: {
224: TS_Alpha *th = (TS_Alpha*)ts->data;
225: Vec X = th->X1; /* X = solution */
226: Vec Y = th->vec_lte_work; /* Y = X + LTE */
227: PetscReal wltea,wlter;
231: if (!th->vec_sol_prev) {*wlte = -1; return(0);}
232: if (!th->vec_lte_work) {*wlte = -1; return(0);}
233: if (ts->steprestart) {
234: /* th->vec_lte_work is set to the LTE in TSAlpha_Restart() */
235: 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]; Vec vecs[3];
241: scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
242: vecs[0] = th->X1; vecs[1] = th->X0; vecs[2] = th->vec_sol_prev;
243: VecCopy(X,Y);
244: VecMAXPY(Y,3,scal,vecs);
245: }
246: TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
247: if (order) *order = 2;
248: return(0);
249: }
251: static PetscErrorCode TSRollBack_Alpha(TS ts)
252: {
253: TS_Alpha *th = (TS_Alpha*)ts->data;
257: VecCopy(th->X0,ts->vec_sol);
258: return(0);
259: }
261: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X)
262: {
263: TS_Alpha *th = (TS_Alpha*)ts->data;
264: PetscReal dt = t - ts->ptime;
268: VecCopy(ts->vec_sol,X);
269: VecAXPY(X,th->Gamma*dt,th->V1);
270: VecAXPY(X,(1-th->Gamma)*dt,th->V0);
271: return(0);
272: }
274: static PetscErrorCode SNESTSFormFunction_Alpha(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts)
275: {
276: TS_Alpha *th = (TS_Alpha*)ts->data;
277: PetscReal ta = th->stage_time;
278: Vec Xa = th->Xa, Va = th->Va;
282: TSAlpha_StageVecs(ts,X);
283: /* F = Function(ta,Xa,Va) */
284: TSComputeIFunction(ts,ta,Xa,Va,F,PETSC_FALSE);
285: VecScale(F,th->scale_F);
286: return(0);
287: }
289: static PetscErrorCode SNESTSFormJacobian_Alpha(PETSC_UNUSED SNES snes,PETSC_UNUSED Vec X,Mat J,Mat P,TS ts)
290: {
291: TS_Alpha *th = (TS_Alpha*)ts->data;
292: PetscReal ta = th->stage_time;
293: Vec Xa = th->Xa, Va = th->Va;
294: PetscReal dVdX = th->shift_V;
298: /* J,P = Jacobian(ta,Xa,Va) */
299: TSComputeIJacobian(ts,ta,Xa,Va,dVdX,J,P,PETSC_FALSE);
300: return(0);
301: }
303: static PetscErrorCode TSReset_Alpha(TS ts)
304: {
305: TS_Alpha *th = (TS_Alpha*)ts->data;
309: VecDestroy(&th->X0);
310: VecDestroy(&th->Xa);
311: VecDestroy(&th->X1);
312: VecDestroy(&th->V0);
313: VecDestroy(&th->Va);
314: VecDestroy(&th->V1);
315: VecDestroy(&th->vec_sol_prev);
316: VecDestroy(&th->vec_lte_work);
317: return(0);
318: }
320: static PetscErrorCode TSDestroy_Alpha(TS ts)
321: {
325: TSReset_Alpha(ts);
326: PetscFree(ts->data);
328: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetRadius_C",NULL);
329: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetParams_C",NULL);
330: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaGetParams_C",NULL);
331: return(0);
332: }
334: static PetscErrorCode TSSetUp_Alpha(TS ts)
335: {
336: TS_Alpha *th = (TS_Alpha*)ts->data;
337: PetscBool match;
341: VecDuplicate(ts->vec_sol,&th->X0);
342: VecDuplicate(ts->vec_sol,&th->Xa);
343: VecDuplicate(ts->vec_sol,&th->X1);
344: VecDuplicate(ts->vec_sol,&th->V0);
345: VecDuplicate(ts->vec_sol,&th->Va);
346: VecDuplicate(ts->vec_sol,&th->V1);
348: TSGetAdapt(ts,&ts->adapt);
349: TSAdaptCandidatesClear(ts->adapt);
350: PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
351: if (!match) {
352: VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
353: VecDuplicate(ts->vec_sol,&th->vec_lte_work);
354: }
356: TSGetSNES(ts,&ts->snes);
357: return(0);
358: }
360: static PetscErrorCode TSSetFromOptions_Alpha(PetscOptionItems *PetscOptionsObject,TS ts)
361: {
362: TS_Alpha *th = (TS_Alpha*)ts->data;
366: PetscOptionsHead(PetscOptionsObject,"Generalized-Alpha ODE solver options");
367: {
368: PetscBool flg;
369: PetscReal radius = 1;
370: PetscOptionsReal("-ts_alpha_radius","Spectral radius (high-frequency dissipation)","TSAlphaSetRadius",radius,&radius,&flg);
371: if (flg) {TSAlphaSetRadius(ts,radius);}
372: PetscOptionsReal("-ts_alpha_alpha_m","Algorithmic parameter alpha_m","TSAlphaSetParams",th->Alpha_m,&th->Alpha_m,NULL);
373: PetscOptionsReal("-ts_alpha_alpha_f","Algorithmic parameter alpha_f","TSAlphaSetParams",th->Alpha_f,&th->Alpha_f,NULL);
374: PetscOptionsReal("-ts_alpha_gamma","Algorithmic parameter gamma","TSAlphaSetParams",th->Gamma,&th->Gamma,NULL);
375: TSAlphaSetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma);
376: }
377: PetscOptionsTail();
378: return(0);
379: }
381: static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer)
382: {
383: TS_Alpha *th = (TS_Alpha*)ts->data;
384: PetscBool iascii;
388: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
389: if (iascii) {
390: PetscViewerASCIIPrintf(viewer," Alpha_m=%g, Alpha_f=%g, Gamma=%g\n",(double)th->Alpha_m,(double)th->Alpha_f,(double)th->Gamma);
391: }
392: return(0);
393: }
395: static PetscErrorCode TSAlphaSetRadius_Alpha(TS ts,PetscReal radius)
396: {
397: PetscReal alpha_m,alpha_f,gamma;
401: if (radius < 0 || radius > 1) SETERRQ1(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: TSAlphaSetParams(ts,alpha_m,alpha_f,gamma);
406: return(0);
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;
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: return(0);
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;
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: return(0);
432: }
434: /*MC
435: TSALPHA - ODE/DAE solver using the implicit Generalized-Alpha method
436: for first-order systems
438: Level: beginner
440: References:
441: K.E. Jansen, C.H. Whiting, G.M. Hulber, "A generalized-alpha
442: method for integrating the filtered Navier-Stokes equations with a
443: stabilized finite element method", Computer Methods in Applied
444: Mechanics and Engineering, 190, 305-319, 2000.
445: DOI: 10.1016/S0045-7825(00)00203-6.
447: J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
448: Dynamics with Improved Numerical Dissipation: The Generalized-alpha
449: Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.
451: .seealso: TS, TSCreate(), TSSetType(), TSAlphaSetRadius(), TSAlphaSetParams()
452: M*/
453: PETSC_EXTERN PetscErrorCode TSCreate_Alpha(TS ts)
454: {
455: TS_Alpha *th;
459: ts->ops->reset = TSReset_Alpha;
460: ts->ops->destroy = TSDestroy_Alpha;
461: ts->ops->view = TSView_Alpha;
462: ts->ops->setup = TSSetUp_Alpha;
463: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
464: ts->ops->step = TSStep_Alpha;
465: ts->ops->evaluatewlte = TSEvaluateWLTE_Alpha;
466: ts->ops->rollback = TSRollBack_Alpha;
467: ts->ops->interpolate = TSInterpolate_Alpha;
468: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
469: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
470: ts->default_adapt_type = TSADAPTNONE;
472: ts->usessnes = PETSC_TRUE;
474: PetscNewLog(ts,&th);
475: ts->data = (void*)th;
477: th->Alpha_m = 0.5;
478: th->Alpha_f = 0.5;
479: th->Gamma = 0.5;
480: th->order = 2;
482: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetRadius_C",TSAlphaSetRadius_Alpha);
483: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaSetParams_C",TSAlphaSetParams_Alpha);
484: PetscObjectComposeFunction((PetscObject)ts,"TSAlphaGetParams_C",TSAlphaGetParams_Alpha);
485: return(0);
486: }
488: /*@
489: TSAlphaSetRadius - sets the desired spectral radius of the method
490: (i.e. high-frequency numerical damping)
492: Logically Collective on TS
494: The algorithmic parameters \alpha_m and \alpha_f of the
495: generalized-\alpha method can be computed in terms of a specified
496: spectral radius \rho in [0,1] for infinite time step in order to
497: control high-frequency numerical damping:
498: \alpha_m = 0.5*(3-\rho)/(1+\rho)
499: \alpha_f = 1/(1+\rho)
501: Input Parameter:
502: + ts - timestepping context
503: - radius - the desired spectral radius
505: Options Database:
506: . -ts_alpha_radius <radius>
508: Level: intermediate
510: .seealso: TSAlphaSetParams(), TSAlphaGetParams()
511: @*/
512: PetscErrorCode TSAlphaSetRadius(TS ts,PetscReal radius)
513: {
519: if (radius < 0 || radius > 1) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Radius %g not in range [0,1]",(double)radius);
520: PetscTryMethod(ts,"TSAlphaSetRadius_C",(TS,PetscReal),(ts,radius));
521: return(0);
522: }
524: /*@
525: TSAlphaSetParams - sets the algorithmic parameters for TSALPHA
527: Logically Collective on TS
529: Second-order accuracy can be obtained so long as:
530: \gamma = 0.5 + alpha_m - alpha_f
532: Unconditional stability requires:
533: \alpha_m >= \alpha_f >= 0.5
535: Backward Euler method is recovered with:
536: \alpha_m = \alpha_f = gamma = 1
538: Input Parameter:
539: + ts - timestepping context
540: . \alpha_m - algorithmic paramenter
541: . \alpha_f - algorithmic paramenter
542: - \gamma - algorithmic paramenter
544: Options Database:
545: + -ts_alpha_alpha_m <alpha_m>
546: . -ts_alpha_alpha_f <alpha_f>
547: - -ts_alpha_gamma <gamma>
549: Note:
550: Use of this function is normally only required to hack TSALPHA to
551: use a modified integration scheme. Users should call
552: TSAlphaSetRadius() to set the desired spectral radius of the methods
553: (i.e. high-frequency damping) in order so select optimal values for
554: these parameters.
556: Level: advanced
558: .seealso: TSAlphaSetRadius(), TSAlphaGetParams()
559: @*/
560: PetscErrorCode TSAlphaSetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma)
561: {
569: PetscTryMethod(ts,"TSAlphaSetParams_C",(TS,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma));
570: return(0);
571: }
573: /*@
574: TSAlphaGetParams - gets the algorithmic parameters for TSALPHA
576: Not Collective
578: Input Parameter:
579: . ts - timestepping context
581: Output Parameters:
582: + \alpha_m - algorithmic parameter
583: . \alpha_f - algorithmic parameter
584: - \gamma - algorithmic parameter
586: Note:
587: Use of this function is normally only required to hack TSALPHA to
588: use a modified integration scheme. Users should call
589: TSAlphaSetRadius() to set the high-frequency damping (i.e. spectral
590: radius of the method) in order so select optimal values for these
591: parameters.
593: Level: advanced
595: .seealso: TSAlphaSetRadius(), TSAlphaSetParams()
596: @*/
597: PetscErrorCode TSAlphaGetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma)
598: {
606: PetscUseMethod(ts,"TSAlphaGetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma));
607: return(0);
608: }