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