Actual source code: alpha.c
petsc-3.3-p7 2013-05-11
1: /*
2: Code for timestepping with implicit generalized-\alpha method
3: for first order systems.
4: */
5: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
7: typedef PetscErrorCode (*TSAlphaAdaptFunction)(TS,PetscReal,Vec,Vec,PetscReal*,PetscBool*,void*);
9: typedef struct {
10: Vec X0,Xa,X1;
11: Vec V0,Va,V1;
12: Vec R,E;
13: PetscReal Alpha_m;
14: PetscReal Alpha_f;
15: PetscReal Gamma;
16: PetscReal stage_time;
17: PetscReal shift;
19: TSAlphaAdaptFunction adapt;
20: void *adaptctx;
21: PetscReal rtol;
22: PetscReal atol;
23: PetscReal rho;
24: PetscReal scale_min;
25: PetscReal scale_max;
26: PetscReal dt_min;
27: PetscReal dt_max;
28: } TS_Alpha;
32: static PetscErrorCode TSStep_Alpha(TS ts)
33: {
34: TS_Alpha *th = (TS_Alpha*)ts->data;
35: PetscInt its,lits,reject;
36: PetscReal next_time_step;
37: SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING;
38: PetscErrorCode ierr;
41: if (ts->steps == 0) {
42: VecSet(th->V0,0.0);
43: } else {
44: VecCopy(th->V1,th->V0);
45: }
46: VecCopy(ts->vec_sol,th->X0);
47: next_time_step = ts->time_step;
48: for (reject=0; reject<ts->max_reject; reject++,ts->reject++) {
49: ts->time_step = next_time_step;
50: th->stage_time = ts->ptime + th->Alpha_f*ts->time_step;
51: th->shift = th->Alpha_m/(th->Alpha_f*th->Gamma*ts->time_step);
52: TSPreStep(ts);
53: TSPreStage(ts,th->stage_time);
54: /* predictor */
55: VecCopy(th->X0,th->X1);
56: /* solve R(X,V) = 0 */
57: SNESSolve(ts->snes,PETSC_NULL,th->X1);
58: /* V1 = (1-1/Gamma)*V0 + 1/(Gamma*dT)*(X1-X0) */
59: VecWAXPY(th->V1,-1,th->X0,th->X1);
60: VecAXPBY(th->V1,1-1/th->Gamma,1/(th->Gamma*ts->time_step),th->V0);
61: /* nonlinear solve convergence */
62: SNESGetConvergedReason(ts->snes,&snesreason);
63: if (snesreason < 0 && !th->adapt) break;
64: SNESGetIterationNumber(ts->snes,&its);
65: SNESGetLinearSolveIterations(ts->snes,&lits);
66: ts->snes_its += its; ts->ksp_its += lits;
67: PetscInfo3(ts,"step=%D, nonlinear solve iterations=%D, linear solve iterations=%D\n",ts->steps,its,lits);
68: /* time step adaptativity */
69: if (!th->adapt) break;
70: else {
71: PetscReal t1 = ts->ptime + ts->time_step;
72: PetscBool stepok = (reject==0) ? PETSC_TRUE : PETSC_FALSE;
73: th->adapt(ts,t1,th->X1,th->V1,&next_time_step,&stepok,th->adaptctx);
74: PetscInfo5(ts,"Step %D (t=%G,dt=%G) %s, next dt=%G\n",ts->steps,ts->ptime,ts->time_step,stepok?"accepted":"rejected",next_time_step);
75: if (stepok) break;
76: }
77: }
78: if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
79: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
80: PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
81: return(0);
82: }
83: if (reject >= ts->max_reject) {
84: ts->reason = TS_DIVERGED_STEP_REJECTED;
85: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
86: return(0);
87: }
88: VecCopy(th->X1,ts->vec_sol);
89: ts->ptime += ts->time_step;
90: ts->time_step = next_time_step;
91: ts->steps++;
92: return(0);
93: }
97: static PetscErrorCode TSInterpolate_Alpha(TS ts,PetscReal t,Vec X)
98: {
99: TS_Alpha *th = (TS_Alpha*)ts->data;
100: PetscReal dt = t - ts->ptime;
104: VecCopy(ts->vec_sol,X);
105: VecAXPY(X,th->Gamma*dt,th->V1);
106: VecAXPY(X,(1-th->Gamma)*dt,th->V0);
107: return(0);
108: }
110: /*------------------------------------------------------------*/
113: static PetscErrorCode TSReset_Alpha(TS ts)
114: {
115: TS_Alpha *th = (TS_Alpha*)ts->data;
116: PetscErrorCode ierr;
119: VecDestroy(&th->X0);
120: VecDestroy(&th->Xa);
121: VecDestroy(&th->X1);
122: VecDestroy(&th->V0);
123: VecDestroy(&th->Va);
124: VecDestroy(&th->V1);
125: VecDestroy(&th->E);
126: return(0);
127: }
131: static PetscErrorCode TSDestroy_Alpha(TS ts)
132: {
133: PetscErrorCode ierr;
136: TSReset_Alpha(ts);
137: PetscFree(ts->data);
139: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetRadius_C","",PETSC_NULL);
140: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetParams_C","",PETSC_NULL);
141: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaGetParams_C","",PETSC_NULL);
142: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetAdapt_C","",PETSC_NULL);
143: return(0);
144: }
148: static PetscErrorCode SNESTSFormFunction_Alpha(SNES snes,Vec x,Vec y,TS ts)
149: {
150: TS_Alpha *th = (TS_Alpha*)ts->data;
151: Vec X0 = th->X0, V0 = th->V0;
152: Vec X1 = x, V1 = th->V1, R = y;
156: /* V1 = (1-1/Gamma)*V0 + 1/(Gamma*dT)*(X1-X0) */
157: VecWAXPY(V1,-1,X0,X1);
158: VecAXPBY(V1,1-1/th->Gamma,1/(th->Gamma*ts->time_step),V0);
159: /* Xa = X0 + Alpha_f*(X1-X0) */
160: VecWAXPY(th->Xa,-1,X0,X1);
161: VecAYPX(th->Xa,th->Alpha_f,X0);
162: /* Va = V0 + Alpha_m*(V1-V0) */
163: VecWAXPY(th->Va,-1,V0,V1);
164: VecAYPX(th->Va,th->Alpha_m,V0);
165: /* F = Function(ta,Xa,Va) */
166: TSComputeIFunction(ts,th->stage_time,th->Xa,th->Va,R,PETSC_FALSE);
167: VecScale(R,1/th->Alpha_f);
168: return(0);
169: }
173: static PetscErrorCode SNESTSFormJacobian_Alpha(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,TS ts)
174: {
175: TS_Alpha *th = (TS_Alpha*)ts->data;
179: /* A,B = Jacobian(ta,Xa,Va) */
180: TSComputeIJacobian(ts,th->stage_time,th->Xa,th->Va,th->shift,A,B,str,PETSC_FALSE);
181: return(0);
182: }
186: static PetscErrorCode TSSetUp_Alpha(TS ts)
187: {
188: TS_Alpha *th = (TS_Alpha*)ts->data;
192: VecDuplicate(ts->vec_sol,&th->X0);
193: VecDuplicate(ts->vec_sol,&th->Xa);
194: VecDuplicate(ts->vec_sol,&th->X1);
195: VecDuplicate(ts->vec_sol,&th->V0);
196: VecDuplicate(ts->vec_sol,&th->Va);
197: VecDuplicate(ts->vec_sol,&th->V1);
198: return(0);
199: }
203: static PetscErrorCode TSSetFromOptions_Alpha(TS ts)
204: {
205: TS_Alpha *th = (TS_Alpha*)ts->data;
209: PetscOptionsHead("Alpha ODE solver options");
210: {
211: PetscBool flag, adapt = PETSC_FALSE;
212: PetscReal radius = 1.0;
213: PetscOptionsReal("-ts_alpha_radius","spectral radius","TSAlphaSetRadius",radius,&radius,&flag);
214: if (flag) { TSAlphaSetRadius(ts,radius); }
215: PetscOptionsReal("-ts_alpha_alpha_m","algoritmic parameter alpha_m","TSAlphaSetParams",th->Alpha_m,&th->Alpha_m,PETSC_NULL);
216: PetscOptionsReal("-ts_alpha_alpha_f","algoritmic parameter alpha_f","TSAlphaSetParams",th->Alpha_f,&th->Alpha_f,PETSC_NULL);
217: PetscOptionsReal("-ts_alpha_gamma","algoritmic parameter gamma","TSAlphaSetParams",th->Gamma,&th->Gamma,PETSC_NULL);
218: TSAlphaSetParams(ts,th->Alpha_m,th->Alpha_f,th->Gamma);
220: PetscOptionsBool("-ts_alpha_adapt","default time step adaptativity","TSAlphaSetAdapt",adapt,&adapt,&flag);
221: if (flag) { TSAlphaSetAdapt(ts,adapt?TSAlphaAdaptDefault:PETSC_NULL,PETSC_NULL); }
222: PetscOptionsReal("-ts_alpha_adapt_rtol","relative tolerance for dt adaptativity","",th->rtol,&th->rtol,PETSC_NULL);
223: PetscOptionsReal("-ts_alpha_adapt_atol","absolute tolerance for dt adaptativity","",th->atol,&th->atol,PETSC_NULL);
224: PetscOptionsReal("-ts_alpha_adapt_min","minimum dt scale","",th->scale_min,&th->scale_min,PETSC_NULL);
225: PetscOptionsReal("-ts_alpha_adapt_max","maximum dt scale","",th->scale_max,&th->scale_max,PETSC_NULL);
226: PetscOptionsReal("-ts_alpha_adapt_dt_min","minimum dt","",th->dt_min,&th->dt_min,PETSC_NULL);
227: PetscOptionsReal("-ts_alpha_adapt_dt_max","maximum dt","",th->dt_max,&th->dt_max,PETSC_NULL);
228: SNESSetFromOptions(ts->snes);
229: }
230: PetscOptionsTail();
231: return(0);
232: }
236: static PetscErrorCode TSView_Alpha(TS ts,PetscViewer viewer)
237: {
238: TS_Alpha *th = (TS_Alpha*)ts->data;
239: PetscBool iascii;
240: PetscErrorCode ierr;
243: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
244: if (iascii) {
245: PetscViewerASCIIPrintf(viewer," Alpha_m=%G, Alpha_f=%G, Gamma=%G\n",th->Alpha_m,th->Alpha_f,th->Gamma);
246: }
247: SNESView(ts->snes,viewer);
248: return(0);
249: }
251: /*------------------------------------------------------------*/
253: EXTERN_C_BEGIN
257: PetscErrorCode TSAlphaSetRadius_Alpha(TS ts,PetscReal radius)
258: {
259: TS_Alpha *th = (TS_Alpha*)ts->data;
262: if (radius < 0 || radius > 1) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Radius %G not in range [0,1]",radius);
263: th->Alpha_m = 0.5*(3-radius)/(1+radius);
264: th->Alpha_f = 1/(1+radius);
265: th->Gamma = 0.5 + th->Alpha_m - th->Alpha_f;
266: return(0);
267: }
271: PetscErrorCode TSAlphaSetParams_Alpha(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma)
272: {
273: TS_Alpha *th = (TS_Alpha*)ts->data;
276: th->Alpha_m = alpha_m;
277: th->Alpha_f = alpha_f;
278: th->Gamma = gamma;
279: return(0);
280: }
284: PetscErrorCode TSAlphaGetParams_Alpha(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma)
285: {
286: TS_Alpha *th = (TS_Alpha*)ts->data;
289: if (alpha_m) *alpha_m = th->Alpha_m;
290: if (alpha_f) *alpha_f = th->Alpha_f;
291: if (gamma) *gamma = th->Gamma;
292: return(0);
293: }
297: PetscErrorCode TSAlphaSetAdapt_Alpha(TS ts,TSAlphaAdaptFunction adapt,void *ctx)
298: {
299: TS_Alpha *th = (TS_Alpha*)ts->data;
302: th->adapt = adapt;
303: th->adaptctx = ctx;
304: return(0);
305: }
307: EXTERN_C_END
309: /* ------------------------------------------------------------ */
310: /*MC
311: TSALPHA - DAE solver using the implicit Generalized-Alpha method
313: Level: beginner
315: References:
316: K.E. Jansen, C.H. Whiting, G.M. Hulber, "A generalized-alpha
317: method for integrating the filtered Navier-Stokes equations with a
318: stabilized finite element method", Computer Methods in Applied
319: Mechanics and Engineering, 190, 305-319, 2000.
320: DOI: 10.1016/S0045-7825(00)00203-6.
322: J. Chung, G.M.Hubert. "A Time Integration Algorithm for Structural
323: Dynamics with Improved Numerical Dissipation: The Generalized-alpha
324: Method" ASME Journal of Applied Mechanics, 60, 371:375, 1993.
326: .seealso: TSCreate(), TS, TSSetType()
328: M*/
329: EXTERN_C_BEGIN
332: PetscErrorCode TSCreate_Alpha(TS ts)
333: {
334: TS_Alpha *th;
338: ts->ops->reset = TSReset_Alpha;
339: ts->ops->destroy = TSDestroy_Alpha;
340: ts->ops->view = TSView_Alpha;
341: ts->ops->setup = TSSetUp_Alpha;
342: ts->ops->step = TSStep_Alpha;
343: ts->ops->interpolate = TSInterpolate_Alpha;
344: ts->ops->setfromoptions = TSSetFromOptions_Alpha;
345: ts->ops->snesfunction = SNESTSFormFunction_Alpha;
346: ts->ops->snesjacobian = SNESTSFormJacobian_Alpha;
348: PetscNewLog(ts,TS_Alpha,&th);
349: ts->data = (void*)th;
351: th->Alpha_m = 0.5;
352: th->Alpha_f = 0.5;
353: th->Gamma = 0.5;
355: th->rtol = 1e-3;
356: th->atol = 1e-3;
357: th->rho = 0.9;
358: th->scale_min = 0.1;
359: th->scale_max = 5.0;
360: th->dt_min = 0.0;
361: th->dt_max = PETSC_MAX_REAL;
363: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetAdapt_C","TSAlphaSetAdapt_Alpha",TSAlphaSetAdapt_Alpha);
364: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetRadius_C","TSAlphaSetRadius_Alpha",TSAlphaSetRadius_Alpha);
365: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaSetParams_C","TSAlphaSetParams_Alpha",TSAlphaSetParams_Alpha);
366: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSAlphaGetParams_C","TSAlphaGetParams_Alpha",TSAlphaGetParams_Alpha);
368: return(0);
369: }
370: EXTERN_C_END
374: /*@C
375: TSAlphaSetAdapt - sets the time step adaptativity and acceptance test routine
377: This function allows to accept/reject a step and select the
378: next time step to use.
380: Not Collective
382: Input Parameter:
383: + ts - timestepping context
384: . adapt - user-defined adapt routine
385: - ctx - [optional] user-defined context for private data for the
386: adapt routine (may be PETSC_NULL)
388: Calling sequence of adapt:
389: $ adapt (TS ts,PetscReal t,Vec X,Vec Xdot,
390: $ PetscReal *next_dt,PetscBool *accepted,void *ctx);
392: Level: intermediate
394: @*/
395: PetscErrorCode TSAlphaSetAdapt(TS ts,TSAlphaAdaptFunction adapt,void *ctx)
396: {
400: PetscTryMethod(ts,"TSAlphaSetAdapt_C",(TS,TSAlphaAdaptFunction,void*),(ts,adapt,ctx));
401: return(0);
402: }
406: PetscErrorCode TSAlphaAdaptDefault(TS ts,PetscReal t,Vec X,Vec Xdot, PetscReal *nextdt,PetscBool *ok,void *ctx)
407: {
408: TS_Alpha *th;
409: SNESConvergedReason snesreason;
410: PetscReal dt,normX,normE,Emax,scale;
411: PetscErrorCode ierr;
415: #if PETSC_USE_DEBUG
416: {
417: PetscBool match;
418: PetscObjectTypeCompare((PetscObject)ts,TSALPHA,&match);
419: if (!match) SETERRQ(((PetscObject)ts)->comm,1,"Only for TSALPHA");
420: }
421: #endif
422: th = (TS_Alpha*)ts->data;
424: SNESGetConvergedReason(ts->snes,&snesreason);
425: if (snesreason < 0) {
426: *ok = PETSC_FALSE;
427: *nextdt *= th->scale_min;
428: goto finally;
429: }
431: /* first-order aproximation to the local error */
432: /* E = (X0 + dt*Xdot) - X */
433: TSGetTimeStep(ts,&dt);
434: if (!th->E) {VecDuplicate(th->X0,&th->E);}
435: VecWAXPY(th->E,dt,Xdot,th->X0);
436: VecAXPY(th->E,-1,X);
437: VecNorm(th->E,NORM_2,&normE);
438: /* compute maximum allowable error */
439: VecNorm(X,NORM_2,&normX);
440: if (normX == 0) {VecNorm(th->X0,NORM_2,&normX);}
441: Emax = th->rtol * normX + th->atol;
442: /* compute next time step */
443: if (normE > 0) {
444: scale = th->rho * PetscRealPart(PetscSqrtScalar((PetscScalar)(Emax/normE)));
445: scale = PetscMax(scale,th->scale_min);
446: scale = PetscMin(scale,th->scale_max);
447: if (!(*ok))
448: scale = PetscMin(1.0,scale);
449: *nextdt *= scale;
450: }
451: /* accept or reject step */
452: if (normE <= Emax)
453: *ok = PETSC_TRUE;
454: else
455: *ok = PETSC_FALSE;
457: finally:
458: *nextdt = PetscMax(*nextdt,th->dt_min);
459: *nextdt = PetscMin(*nextdt,th->dt_max);
460: return(0);
461: }
465: /*@
466: TSAlphaSetRadius - sets the desired spectral radius of the method
467: (i.e. high-frequency numerical damping)
469: Logically Collective on TS
471: The algorithmic parameters \alpha_m and \alpha_f of the
472: generalized-\alpha method can be computed in terms of a specified
473: spectral radius \rho in [0,1] for infinite time step in order to
474: control high-frequency numerical damping:
475: alpha_m = 0.5*(3-\rho)/(1+\rho)
476: alpha_f = 1/(1+\rho)
478: Input Parameter:
479: + ts - timestepping context
480: - radius - the desired spectral radius
482: Options Database:
483: . -ts_alpha_radius <radius>
485: Level: intermediate
487: .seealso: TSAlphaSetParams(), TSAlphaGetParams()
488: @*/
489: PetscErrorCode TSAlphaSetRadius(TS ts,PetscReal radius)
490: {
495: PetscTryMethod(ts,"TSAlphaSetRadius_C",(TS,PetscReal),(ts,radius));
496: return(0);
497: }
501: /*@
502: TSAlphaSetParams - sets the algorithmic parameters for TSALPHA
504: Not Collective
506: Second-order accuracy can be obtained so long as:
507: \gamma = 0.5 + alpha_m - alpha_f
509: Unconditional stability requires:
510: \alpha_m >= \alpha_f >= 0.5
512: Backward Euler method is recovered when:
513: \alpha_m = \alpha_f = gamma = 1
516: Input Parameter:
517: + ts - timestepping context
518: . \alpha_m - algorithmic paramenter
519: . \alpha_f - algorithmic paramenter
520: - \gamma - algorithmic paramenter
522: Options Database:
523: + -ts_alpha_alpha_m <alpha_m>
524: . -ts_alpha_alpha_f <alpha_f>
525: - -ts_alpha_gamma <gamma>
527: Note:
528: Use of this function is normally only required to hack TSALPHA to
529: use a modified integration scheme. Users should call
530: TSAlphaSetRadius() to set the desired spectral radius of the methods
531: (i.e. high-frequency damping) in order so select optimal values for
532: these parameters.
534: Level: advanced
536: .seealso: TSAlphaSetRadius(), TSAlphaGetParams()
537: @*/
538: PetscErrorCode TSAlphaSetParams(TS ts,PetscReal alpha_m,PetscReal alpha_f,PetscReal gamma)
539: {
544: PetscTryMethod(ts,"TSAlphaSetParams_C",(TS,PetscReal,PetscReal,PetscReal),(ts,alpha_m,alpha_f,gamma));
545: return(0);
546: }
550: /*@
551: TSAlphaGetParams - gets the algorithmic parameters for TSALPHA
553: Not Collective
555: Input Parameter:
556: + ts - timestepping context
557: . \alpha_m - algorithmic parameter
558: . \alpha_f - algorithmic parameter
559: - \gamma - algorithmic parameter
561: Note:
562: Use of this function is normally only required to hack TSALPHA to
563: use a modified integration scheme. Users should call
564: TSAlphaSetRadius() to set the high-frequency damping (i.e. spectral
565: radius of the method) in order so select optimal values for these
566: parameters.
568: Level: advanced
570: .seealso: TSAlphaSetRadius(), TSAlphaSetParams()
571: @*/
572: PetscErrorCode TSAlphaGetParams(TS ts,PetscReal *alpha_m,PetscReal *alpha_f,PetscReal *gamma)
573: {
581: PetscUseMethod(ts,"TSAlphaGetParams_C",(TS,PetscReal*,PetscReal*,PetscReal*),(ts,alpha_m,alpha_f,gamma));
582: return(0);
583: }