Actual source code: theta.c
petsc-3.3-p7 2013-05-11
1: /*
2: Code for timestepping with implicit Theta method
3: */
4: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
5: #include <petscsnesfas.h>
7: typedef struct {
8: Vec X,Xdot; /* Storage for one stage */
9: Vec affine; /* Affine vector needed for residual at beginning of step */
10: PetscBool extrapolate;
11: PetscBool endpoint;
12: PetscReal Theta;
13: PetscReal shift;
14: PetscReal stage_time;
15: } TS_Theta;
19: static PetscErrorCode TSThetaGetX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
20: {
21: TS_Theta *th = (TS_Theta*)ts->data;
25: if (X0) {
26: if (dm && dm != ts->dm) {
27: PetscObjectQuery((PetscObject)dm,"TSTheta_X0",(PetscObject*)X0);
28: if (!*X0) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_ARG_INCOMP,"TSTheta_X0 has not been composed with DM from SNES");
29: } else *X0 = ts->vec_sol;
30: }
31: if (Xdot) {
32: if (dm && dm != ts->dm) {
33: PetscObjectQuery((PetscObject)dm,"TSTheta_Xdot",(PetscObject*)Xdot);
34: if (!*Xdot) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_ARG_INCOMP,"TSTheta_Xdot has not been composed with DM from SNES");
35: } else *Xdot = th->Xdot;
36: }
37: return(0);
38: }
42: static PetscErrorCode DMCoarsenHook_TSTheta(DM fine,DM coarse,void *ctx)
43: {
44: Vec X0,Xdot;
48: DMCreateGlobalVector(coarse,&X0);
49: DMCreateGlobalVector(coarse,&Xdot);
50: /* Oh noes, this would create a loop because the Vec holds a reference to the DM.
51: Making a PetscContainer to hold these Vecs would make the following call succeed, but would create a reference loop.
52: Need to decide on a way to break the reference counting loop.
53: */
54: PetscObjectCompose((PetscObject)coarse,"TSTheta_X0",(PetscObject)X0);
55: PetscObjectCompose((PetscObject)coarse,"TSTheta_Xdot",(PetscObject)Xdot);
56: VecDestroy(&X0);
57: VecDestroy(&Xdot);
58: return(0);
59: }
63: static PetscErrorCode DMRestrictHook_TSTheta(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
64: {
65: TS ts = (TS)ctx;
67: Vec X0,Xdot,X0_c,Xdot_c;
70: TSThetaGetX0AndXdot(ts,fine,&X0,&Xdot);
71: TSThetaGetX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
72: MatRestrict(restrct,X0,X0_c);
73: MatRestrict(restrct,Xdot,Xdot_c);
74: VecPointwiseMult(X0_c,rscale,X0_c);
75: VecPointwiseMult(Xdot_c,rscale,Xdot_c);
76: return(0);
77: }
81: static PetscErrorCode TSStep_Theta(TS ts)
82: {
83: TS_Theta *th = (TS_Theta*)ts->data;
84: PetscInt its,lits;
85: PetscReal next_time_step;
86: SNESConvergedReason snesreason;
87: PetscErrorCode ierr;
90: next_time_step = ts->time_step;
91: th->stage_time = ts->ptime + (th->endpoint ? 1. : th->Theta)*ts->time_step;
92: th->shift = 1./(th->Theta*ts->time_step);
93: TSPreStep(ts);
94: TSPreStage(ts,th->stage_time);
96: if (th->endpoint) { /* This formulation assumes linear time-independent mass matrix */
97: VecZeroEntries(th->Xdot);
98: if (!th->affine) {VecDuplicate(ts->vec_sol,&th->affine);}
99: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,th->Xdot,th->affine,PETSC_FALSE);
100: VecScale(th->affine,(th->Theta-1.)/th->Theta);
101: }
102: if (th->extrapolate) {
103: VecWAXPY(th->X,1./th->shift,th->Xdot,ts->vec_sol);
104: } else {
105: VecCopy(ts->vec_sol,th->X);
106: }
107: SNESSolve(ts->snes,th->affine,th->X);
108: SNESGetIterationNumber(ts->snes,&its);
109: SNESGetLinearSolveIterations(ts->snes,&lits);
110: SNESGetConvergedReason(ts->snes,&snesreason);
111: ts->snes_its += its; ts->ksp_its += lits;
112: if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
113: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
114: PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
115: return(0);
116: }
117: if (th->endpoint) {
118: VecCopy(th->X,ts->vec_sol);
119: } else {
120: VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,ts->vec_sol,th->X);
121: VecAXPY(ts->vec_sol,ts->time_step,th->Xdot);
122: }
123: ts->ptime += ts->time_step;
124: ts->time_step = next_time_step;
125: ts->steps++;
126: return(0);
127: }
131: static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X)
132: {
133: TS_Theta *th = (TS_Theta*)ts->data;
134: PetscReal alpha = t - ts->ptime;
138: VecCopy(ts->vec_sol,th->X);
139: if (th->endpoint) alpha *= th->Theta;
140: VecWAXPY(X,alpha,th->Xdot,th->X);
141: return(0);
142: }
144: /*------------------------------------------------------------*/
147: static PetscErrorCode TSReset_Theta(TS ts)
148: {
149: TS_Theta *th = (TS_Theta*)ts->data;
150: PetscErrorCode ierr;
153: VecDestroy(&th->X);
154: VecDestroy(&th->Xdot);
155: VecDestroy(&th->affine);
156: return(0);
157: }
161: static PetscErrorCode TSDestroy_Theta(TS ts)
162: {
163: PetscErrorCode ierr;
166: TSReset_Theta(ts);
167: PetscFree(ts->data);
168: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","",PETSC_NULL);
169: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","",PETSC_NULL);
170: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetEndpoint_C","",PETSC_NULL);
171: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetEndpoint_C","",PETSC_NULL);
172: return(0);
173: }
175: /*
176: This defines the nonlinear equation that is to be solved with SNES
177: G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
178: */
181: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
182: {
183: TS_Theta *th = (TS_Theta*)ts->data;
185: Vec X0,Xdot;
186: DM dm,dmsave;
189: SNESGetDM(snes,&dm);
190: /* When using the endpoint variant, this is actually 1/Theta * Xdot */
191: TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
192: VecAXPBYPCZ(Xdot,-th->shift,th->shift,0,X0,x);
194: /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
195: dmsave = ts->dm;
196: ts->dm = dm;
197: TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
198: ts->dm = dmsave;
199: return(0);
200: }
204: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,TS ts)
205: {
206: TS_Theta *th = (TS_Theta*)ts->data;
208: Vec Xdot;
209: DM dm,dmsave;
212: SNESGetDM(snes,&dm);
214: /* th->Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
215: TSThetaGetX0AndXdot(ts,dm,PETSC_NULL,&Xdot);
217: dmsave = ts->dm;
218: ts->dm = dm;
219: TSComputeIJacobian(ts,th->stage_time,x,Xdot,th->shift,A,B,str,PETSC_FALSE);
220: ts->dm = dmsave;
221: return(0);
222: }
226: static PetscErrorCode TSSetUp_Theta(TS ts)
227: {
228: TS_Theta *th = (TS_Theta*)ts->data;
230: SNES snes;
231: DM dm;
234: VecDuplicate(ts->vec_sol,&th->X);
235: VecDuplicate(ts->vec_sol,&th->Xdot);
236: TSGetSNES(ts,&snes);
237: TSGetDM(ts,&dm);
238: if (dm) {
239: DMCoarsenHookAdd(dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
240: }
241: return(0);
242: }
243: /*------------------------------------------------------------*/
247: static PetscErrorCode TSSetFromOptions_Theta(TS ts)
248: {
249: TS_Theta *th = (TS_Theta*)ts->data;
253: PetscOptionsHead("Theta ODE solver options");
254: {
255: PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,PETSC_NULL);
256: PetscOptionsBool("-ts_theta_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,PETSC_NULL);
257: PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,PETSC_NULL);
258: SNESSetFromOptions(ts->snes);
259: }
260: PetscOptionsTail();
261: return(0);
262: }
266: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
267: {
268: TS_Theta *th = (TS_Theta*)ts->data;
269: PetscBool iascii;
270: PetscErrorCode ierr;
273: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
274: if (iascii) {
275: PetscViewerASCIIPrintf(viewer," Theta=%G\n",th->Theta);
276: PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate?"yes":"no");
277: }
278: SNESView(ts->snes,viewer);
279: return(0);
280: }
282: EXTERN_C_BEGIN
285: PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
286: {
287: TS_Theta *th = (TS_Theta*)ts->data;
290: *theta = th->Theta;
291: return(0);
292: }
296: PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta)
297: {
298: TS_Theta *th = (TS_Theta*)ts->data;
301: if (theta <= 0 || 1 < theta) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Theta %G not in range (0,1]",theta);
302: th->Theta = theta;
303: return(0);
304: }
308: PetscErrorCode TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
309: {
310: TS_Theta *th = (TS_Theta*)ts->data;
313: *endpoint = th->endpoint;
314: return(0);
315: }
319: PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
320: {
321: TS_Theta *th = (TS_Theta*)ts->data;
324: th->endpoint = flg;
325: return(0);
326: }
327: EXTERN_C_END
329: /* ------------------------------------------------------------ */
330: /*MC
331: TSTHETA - DAE solver using the implicit Theta method
333: Level: beginner
335: Notes:
336: This method can be applied to DAE.
338: This method is cast as a 1-stage implicit Runge-Kutta method.
340: .vb
341: Theta | Theta
342: -------------
343: | 1
344: .ve
346: For the default Theta=0.5, this is also known as the implicit midpoint rule.
348: When the endpoint variant is chosen, the method becomes a 2-stage method with first stage explicit:
350: .vb
351: 0 | 0 0
352: 1 | 1-Theta Theta
353: -------------------
354: | 1-Theta Theta
355: .ve
357: For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN).
359: To apply a diagonally implicit RK method to DAE, the stage formula
361: $ Y_i = X + h sum_j a_ij Y'_j
363: is interpreted as a formula for Y'_i in terms of Y_i and known stuff (Y'_j, j<i)
365: .seealso: TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint()
367: M*/
368: EXTERN_C_BEGIN
371: PetscErrorCode TSCreate_Theta(TS ts)
372: {
373: TS_Theta *th;
377: ts->ops->reset = TSReset_Theta;
378: ts->ops->destroy = TSDestroy_Theta;
379: ts->ops->view = TSView_Theta;
380: ts->ops->setup = TSSetUp_Theta;
381: ts->ops->step = TSStep_Theta;
382: ts->ops->interpolate = TSInterpolate_Theta;
383: ts->ops->setfromoptions = TSSetFromOptions_Theta;
384: ts->ops->snesfunction = SNESTSFormFunction_Theta;
385: ts->ops->snesjacobian = SNESTSFormJacobian_Theta;
387: PetscNewLog(ts,TS_Theta,&th);
388: ts->data = (void*)th;
390: th->extrapolate = PETSC_FALSE;
391: th->Theta = 0.5;
393: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetTheta_C","TSThetaGetTheta_Theta",TSThetaGetTheta_Theta);
394: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetTheta_C","TSThetaSetTheta_Theta",TSThetaSetTheta_Theta);
395: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaGetEndpoint_C","TSThetaGetEndpoint_Theta",TSThetaGetEndpoint_Theta);
396: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSThetaSetEndpoint_C","TSThetaSetEndpoint_Theta",TSThetaSetEndpoint_Theta);
397: return(0);
398: }
399: EXTERN_C_END
403: /*@
404: TSThetaGetTheta - Get the abscissa of the stage in (0,1].
406: Not Collective
408: Input Parameter:
409: . ts - timestepping context
411: Output Parameter:
412: . theta - stage abscissa
414: Note:
415: Use of this function is normally only required to hack TSTHETA to use a modified integration scheme.
417: Level: Advanced
419: .seealso: TSThetaSetTheta()
420: @*/
421: PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta)
422: {
428: PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
429: return(0);
430: }
434: /*@
435: TSThetaSetTheta - Set the abscissa of the stage in (0,1].
437: Not Collective
439: Input Parameter:
440: + ts - timestepping context
441: - theta - stage abscissa
443: Options Database:
444: . -ts_theta_theta <theta>
446: Level: Intermediate
448: .seealso: TSThetaGetTheta()
449: @*/
450: PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta)
451: {
456: PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
457: return(0);
458: }
462: /*@
463: TSThetaGetEndpoint - Gets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
465: Not Collective
467: Input Parameter:
468: . ts - timestepping context
470: Output Parameter:
471: . endpoint - PETSC_TRUE when using the endpoint variant
473: Level: Advanced
475: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
476: @*/
477: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
478: {
484: PetscTryMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
485: return(0);
486: }
490: /*@
491: TSThetaSetEndpoint - Sets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
493: Not Collective
495: Input Parameter:
496: + ts - timestepping context
497: - flg - PETSC_TRUE to use the endpoint variant
499: Options Database:
500: . -ts_theta_endpoint <flg>
502: Level: Intermediate
504: .seealso: TSTHETA, TSCN
505: @*/
506: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
507: {
512: PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
513: return(0);
514: }
516: /*
517: * TSBEULER and TSCN are straightforward specializations of TSTHETA.
518: * The creation functions for these specializations are below.
519: */
523: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
524: {
528: SNESView(ts->snes,viewer);
529: return(0);
530: }
532: /*MC
533: TSBEULER - ODE solver using the implicit backward Euler method
535: Level: beginner
537: .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA
539: M*/
540: EXTERN_C_BEGIN
543: PetscErrorCode TSCreate_BEuler(TS ts)
544: {
548: TSCreate_Theta(ts);
549: TSThetaSetTheta(ts,1.0);
550: ts->ops->view = TSView_BEuler;
551: return(0);
552: }
553: EXTERN_C_END
557: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
558: {
562: SNESView(ts->snes,viewer);
563: return(0);
564: }
566: /*MC
567: TSCN - ODE solver using the implicit Crank-Nicolson method.
569: Level: beginner
571: Notes:
572: TSCN is equivalent to TSTHETA with Theta=0.5 and the "endpoint" option set. I.e.
574: $ -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint
576: .seealso: TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA
578: M*/
579: EXTERN_C_BEGIN
582: PetscErrorCode TSCreate_CN(TS ts)
583: {
587: TSCreate_Theta(ts);
588: TSThetaSetTheta(ts,0.5);
589: TSThetaSetEndpoint(ts,PETSC_TRUE);
590: ts->ops->view = TSView_CN;
591: return(0);
592: }
593: EXTERN_C_END