Actual source code: theta.c

petsc-master 2017-11-20
Report Typos and Errors
  1: /*
  2:   Code for timestepping with implicit Theta method
  3: */
  4:  #include <petsc/private/tsimpl.h>
  5:  #include <petscsnes.h>
  6:  #include <petscdm.h>
  7:  #include <petscmat.h>

  9: typedef struct {
 10:   /* context for time stepping */
 11:   PetscReal    stage_time;
 12:   Vec          X0,X,Xdot;                /* Storage for stages and time derivative */
 13:   Vec          affine;                   /* Affine vector needed for residual at beginning of step in endpoint formulation */
 14:   PetscReal    Theta;
 15:   PetscReal    ptime;
 16:   PetscReal    time_step;
 17:   PetscInt     order;
 18:   PetscBool    endpoint;
 19:   PetscBool    extrapolate;
 20:   TSStepStatus status;
 21:   Vec          VecCostIntegral0;         /* Backup for roll-backs due to events */

 23:   /* context for sensitivity analysis */
 24:   PetscInt     num_tlm;                  /* Total number of tangent linear equations */
 25:   Vec          *VecsDeltaLam;            /* Increment of the adjoint sensitivity w.r.t IC at stage */
 26:   Vec          *VecsDeltaMu;             /* Increment of the adjoint sensitivity w.r.t P at stage */
 27:   Vec          *VecsSensiTemp;           /* Vector to be timed with Jacobian transpose */
 28:   Vec          *VecsDeltaFwdSensi;       /* Increment of the forward sensitivity at stage */
 29:   Vec          VecIntegralSensiTemp;     /* Working vector for forward integral sensitivity */
 30:   Vec          VecIntegralSensipTemp;    /* Working vector for forward integral sensitivity */
 31:   Vec          *VecsFwdSensi0;           /* backup for roll-backs due to events */
 32:   Vec          *VecsIntegralSensi0;      /* backup for roll-backs due to events */
 33:   Vec          *VecsIntegralSensip0;     /* backup for roll-backs due to events */

 35:   /* context for error estimation */
 36:   Vec          vec_sol_prev;
 37:   Vec          vec_lte_work;
 38: } TS_Theta;

 40: static PetscErrorCode TSThetaGetX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
 41: {
 42:   TS_Theta       *th = (TS_Theta*)ts->data;

 46:   if (X0) {
 47:     if (dm && dm != ts->dm) {
 48:       DMGetNamedGlobalVector(dm,"TSTheta_X0",X0);
 49:     } else *X0 = ts->vec_sol;
 50:   }
 51:   if (Xdot) {
 52:     if (dm && dm != ts->dm) {
 53:       DMGetNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
 54:     } else *Xdot = th->Xdot;
 55:   }
 56:   return(0);
 57: }

 59: static PetscErrorCode TSThetaRestoreX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
 60: {

 64:   if (X0) {
 65:     if (dm && dm != ts->dm) {
 66:       DMRestoreNamedGlobalVector(dm,"TSTheta_X0",X0);
 67:     }
 68:   }
 69:   if (Xdot) {
 70:     if (dm && dm != ts->dm) {
 71:       DMRestoreNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
 72:     }
 73:   }
 74:   return(0);
 75: }

 77: static PetscErrorCode DMCoarsenHook_TSTheta(DM fine,DM coarse,void *ctx)
 78: {
 80:   return(0);
 81: }

 83: static PetscErrorCode DMRestrictHook_TSTheta(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
 84: {
 85:   TS             ts = (TS)ctx;
 87:   Vec            X0,Xdot,X0_c,Xdot_c;

 90:   TSThetaGetX0AndXdot(ts,fine,&X0,&Xdot);
 91:   TSThetaGetX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
 92:   MatRestrict(restrct,X0,X0_c);
 93:   MatRestrict(restrct,Xdot,Xdot_c);
 94:   VecPointwiseMult(X0_c,rscale,X0_c);
 95:   VecPointwiseMult(Xdot_c,rscale,Xdot_c);
 96:   TSThetaRestoreX0AndXdot(ts,fine,&X0,&Xdot);
 97:   TSThetaRestoreX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
 98:   return(0);
 99: }

101: static PetscErrorCode DMSubDomainHook_TSTheta(DM dm,DM subdm,void *ctx)
102: {
104:   return(0);
105: }

107: static PetscErrorCode DMSubDomainRestrictHook_TSTheta(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
108: {
109:   TS             ts = (TS)ctx;
111:   Vec            X0,Xdot,X0_sub,Xdot_sub;

114:   TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
115:   TSThetaGetX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);

117:   VecScatterBegin(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);
118:   VecScatterEnd(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);

120:   VecScatterBegin(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);
121:   VecScatterEnd(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);

123:   TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
124:   TSThetaRestoreX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);
125:   return(0);
126: }

128: static PetscErrorCode TSThetaEvaluateCostIntegral(TS ts)
129: {
130:   TS_Theta       *th = (TS_Theta*)ts->data;

134:   if (th->endpoint) {
135:     /* Evolve ts->vec_costintegral to compute integrals */
136:     if (th->Theta!=1.0) {
137:       TSComputeCostIntegrand(ts,th->ptime,th->X0,ts->vec_costintegrand);
138:       VecAXPY(ts->vec_costintegral,th->time_step*(1.0-th->Theta),ts->vec_costintegrand);
139:     }
140:     TSComputeCostIntegrand(ts,ts->ptime,ts->vec_sol,ts->vec_costintegrand);
141:     VecAXPY(ts->vec_costintegral,th->time_step*th->Theta,ts->vec_costintegrand);
142:   } else {
143:     TSComputeCostIntegrand(ts,th->stage_time,th->X,ts->vec_costintegrand);
144:     VecAXPY(ts->vec_costintegral,th->time_step,ts->vec_costintegrand);
145:   }
146:   return(0);
147: }

149: static PetscErrorCode TSForwardCostIntegral_Theta(TS ts)
150: {
151:   TS_Theta       *th = (TS_Theta*)ts->data;

155:   /* backup cost integral */
156:   VecCopy(ts->vec_costintegral,th->VecCostIntegral0);
157:   TSThetaEvaluateCostIntegral(ts);
158:   return(0);
159: }

161: static PetscErrorCode TSAdjointCostIntegral_Theta(TS ts)
162: {

166:   TSThetaEvaluateCostIntegral(ts);
167:   return(0);
168: }

170: static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x)
171: {
172:   PetscInt       nits,lits;

176:   SNESSolve(ts->snes,b,x);
177:   SNESGetIterationNumber(ts->snes,&nits);
178:   SNESGetLinearSolveIterations(ts->snes,&lits);
179:   ts->snes_its += nits; ts->ksp_its += lits;
180:   return(0);
181: }

183: static PetscErrorCode TSStep_Theta(TS ts)
184: {
185:   TS_Theta       *th = (TS_Theta*)ts->data;
186:   PetscInt       rejections = 0;
187:   PetscBool      stageok,accept = PETSC_TRUE;
188:   PetscReal      next_time_step = ts->time_step;

192:   if (!ts->steprollback) {
193:     if (th->vec_sol_prev) { VecCopy(th->X0,th->vec_sol_prev); }
194:     VecCopy(ts->vec_sol,th->X0);
195:   }

197:   th->status = TS_STEP_INCOMPLETE;
198:   while (!ts->reason && th->status != TS_STEP_COMPLETE) {

200:     PetscReal shift = 1/(th->Theta*ts->time_step);
201:     th->stage_time = ts->ptime + (th->endpoint ? (PetscReal)1 : th->Theta)*ts->time_step;

203:     VecCopy(th->X0,th->X);
204:     if (th->extrapolate && !ts->steprestart) {
205:       VecAXPY(th->X,1/shift,th->Xdot);
206:     }
207:     if (th->endpoint) { /* This formulation assumes linear time-independent mass matrix */
208:       if (!th->affine) {VecDuplicate(ts->vec_sol,&th->affine);}
209:       VecZeroEntries(th->Xdot);
210:       TSComputeIFunction(ts,ts->ptime,th->X0,th->Xdot,th->affine,PETSC_FALSE);
211:       VecScale(th->affine,(th->Theta-1)/th->Theta);
212:     } else if (th->affine) { /* Just in case th->endpoint is changed between calls to TSStep_Theta() */
213:       VecZeroEntries(th->affine);
214:     }
215:     TSPreStage(ts,th->stage_time);
216:     TS_SNESSolve(ts,th->affine,th->X);
217:     TSPostStage(ts,th->stage_time,0,&th->X);
218:     TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->X,&stageok);
219:     if (!stageok) goto reject_step;

221:     th->status = TS_STEP_PENDING;
222:     if (th->endpoint) {
223:       VecCopy(th->X,ts->vec_sol);
224:     } else {
225:       VecAXPBYPCZ(th->Xdot,-shift,shift,0,th->X0,th->X);
226:       VecAXPY(ts->vec_sol,ts->time_step,th->Xdot);
227:     }
228:     TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
229:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
230:     if (!accept) {
231:       VecCopy(th->X0,ts->vec_sol);
232:       ts->time_step = next_time_step;
233:       goto reject_step;
234:     }

236:     if (ts->forward_solve || ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */
237:       th->ptime     = ts->ptime;
238:       th->time_step = ts->time_step;
239:     }

241:     ts->ptime += ts->time_step;
242:     ts->time_step = next_time_step;
243:     break;

245:   reject_step:
246:     ts->reject++; accept = PETSC_FALSE;
247:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
248:       ts->reason = TS_DIVERGED_STEP_REJECTED;
249:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
250:     }
251:   }
252:   return(0);
253: }

255: static PetscErrorCode TSAdjointStep_Theta(TS ts)
256: {
257:   TS_Theta            *th = (TS_Theta*)ts->data;
258:   Vec                 *VecsDeltaLam = th->VecsDeltaLam,*VecsDeltaMu = th->VecsDeltaMu,*VecsSensiTemp = th->VecsSensiTemp;
259:   PetscInt            nadj;
260:   PetscErrorCode      ierr;
261:   Mat                 J,Jp;
262:   KSP                 ksp;
263:   PetscReal           shift;

266:   th->status = TS_STEP_INCOMPLETE;
267:   SNESGetKSP(ts->snes,&ksp);
268:   TSGetIJacobian(ts,&J,&Jp,NULL,NULL);

270:   /* If endpoint=1, th->ptime and th->X0 will be used; if endpoint=0, th->stage_time and th->X will be used. */
271:   th->stage_time = th->endpoint ? ts->ptime : (ts->ptime+(1.-th->Theta)*ts->time_step); /* time_step is negative*/
272:   th->ptime      = ts->ptime + ts->time_step;
273:   th->time_step  = -ts->time_step;

275:   /* Build RHS */
276:   if (ts->vec_costintegral) { /* Cost function has an integral term */
277:     if (th->endpoint) {
278:       TSAdjointComputeDRDYFunction(ts,th->stage_time,ts->vec_sol,ts->vecs_drdy);
279:     } else {
280:       TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
281:     }
282:   }
283:   for (nadj=0; nadj<ts->numcost; nadj++) {
284:     VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);
285:     VecScale(VecsSensiTemp[nadj],1./(th->Theta*th->time_step));
286:     if (ts->vec_costintegral) {
287:       VecAXPY(VecsSensiTemp[nadj],1.,ts->vecs_drdy[nadj]);
288:     }
289:   }

291:   /* Build LHS */
292:   shift = 1./(th->Theta*th->time_step);
293:   if (th->endpoint) {
294:     TSComputeIJacobian(ts,th->stage_time,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE);
295:   } else {
296:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
297:   }
298:   KSPSetOperators(ksp,J,Jp);

300:   /* Solve LHS X = RHS */
301:   for (nadj=0; nadj<ts->numcost; nadj++) {
302:     KSPSolveTranspose(ksp,VecsSensiTemp[nadj],VecsDeltaLam[nadj]);
303:   }

305:   /* Update sensitivities, and evaluate integrals if there is any */
306:   if(th->endpoint) { /* two-stage case */
307:     if (th->Theta!=1.) {
308:       shift = 1./((th->Theta-1.)*th->time_step);
309:       TSComputeIJacobian(ts,th->ptime,th->X0,th->Xdot,shift,J,Jp,PETSC_FALSE);
310:       if (ts->vec_costintegral) {
311:         TSAdjointComputeDRDYFunction(ts,th->ptime,th->X0,ts->vecs_drdy);
312:       }
313:       for (nadj=0; nadj<ts->numcost; nadj++) {
314:         MatMultTranspose(J,VecsDeltaLam[nadj],ts->vecs_sensi[nadj]);
315:         if (ts->vec_costintegral) {
316:           VecAXPY(ts->vecs_sensi[nadj],-1.,ts->vecs_drdy[nadj]);
317:         }
318:         VecScale(ts->vecs_sensi[nadj],1./shift);
319:       }
320:     } else { /* backward Euler */
321:       shift = 0.0;
322:       TSComputeIJacobian(ts,th->stage_time,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE); /* get -f_y */
323:       for (nadj=0; nadj<ts->numcost; nadj++) {
324:         MatMultTranspose(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj]);
325:         VecAXPY(ts->vecs_sensi[nadj],-th->time_step,VecsSensiTemp[nadj]);
326:         if (ts->vec_costintegral) {
327:           VecAXPY(ts->vecs_sensi[nadj],th->time_step,ts->vecs_drdy[nadj]);
328:         }
329:       }
330:     }

332:     if (ts->vecs_sensip) { /* sensitivities wrt parameters */
333:       TSAdjointComputeRHSJacobian(ts,th->stage_time,ts->vec_sol,ts->Jacp);
334:       for (nadj=0; nadj<ts->numcost; nadj++) {
335:         MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
336:         VecAXPY(ts->vecs_sensip[nadj],th->time_step*th->Theta,VecsDeltaMu[nadj]);
337:       }
338:       if (th->Theta!=1.) {
339:         TSAdjointComputeRHSJacobian(ts,th->ptime,th->X0,ts->Jacp);
340:         for (nadj=0; nadj<ts->numcost; nadj++) {
341:           MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
342:           VecAXPY(ts->vecs_sensip[nadj],th->time_step*(1.-th->Theta),VecsDeltaMu[nadj]);
343:         }
344:       }
345:       if (ts->vec_costintegral) {
346:         TSAdjointComputeDRDPFunction(ts,th->stage_time,ts->vec_sol,ts->vecs_drdp);
347:         for (nadj=0; nadj<ts->numcost; nadj++) {
348:           VecAXPY(ts->vecs_sensip[nadj],th->time_step*th->Theta,ts->vecs_drdp[nadj]);
349:         }
350:         if (th->Theta!=1.) {
351:           TSAdjointComputeDRDPFunction(ts,th->ptime,th->X0,ts->vecs_drdp);
352:           for (nadj=0; nadj<ts->numcost; nadj++) {
353:             VecAXPY(ts->vecs_sensip[nadj],th->time_step*(1.-th->Theta),ts->vecs_drdp[nadj]);
354:           }
355:         }
356:       }
357:     }
358:   } else { /* one-stage case */
359:     shift = 0.0;
360:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE); /* get -f_y */
361:     if (ts->vec_costintegral) {
362:       TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
363:     }
364:     for (nadj=0; nadj<ts->numcost; nadj++) {
365:       MatMultTranspose(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj]);
366:       VecAXPY(ts->vecs_sensi[nadj],-th->time_step,VecsSensiTemp[nadj]);
367:       if (ts->vec_costintegral) {
368:         VecAXPY(ts->vecs_sensi[nadj],th->time_step,ts->vecs_drdy[nadj]);
369:       }
370:     }
371:     if (ts->vecs_sensip) {
372:       TSAdjointComputeRHSJacobian(ts,th->stage_time,th->X,ts->Jacp);
373:       for (nadj=0; nadj<ts->numcost; nadj++) {
374:         MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
375:         VecAXPY(ts->vecs_sensip[nadj],th->time_step,VecsDeltaMu[nadj]);
376:       }
377:       if (ts->vec_costintegral) {
378:         TSAdjointComputeDRDPFunction(ts,th->stage_time,th->X,ts->vecs_drdp);
379:         for (nadj=0; nadj<ts->numcost; nadj++) {
380:           VecAXPY(ts->vecs_sensip[nadj],th->time_step,ts->vecs_drdp[nadj]);
381:         }
382:       }
383:     }
384:   }

386:   th->status = TS_STEP_COMPLETE;
387:   return(0);
388: }

390: static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X)
391: {
392:   TS_Theta       *th = (TS_Theta*)ts->data;
393:   PetscReal      dt  = t - ts->ptime;

397:   VecCopy(ts->vec_sol,th->X);
398:   if (th->endpoint) dt *= th->Theta;
399:   VecWAXPY(X,dt,th->Xdot,th->X);
400:   return(0);
401: }

403: static PetscErrorCode TSEvaluateWLTE_Theta(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
404: {
405:   TS_Theta       *th = (TS_Theta*)ts->data;
406:   Vec            X = ts->vec_sol;      /* X = solution */
407:   Vec            Y = th->vec_lte_work; /* Y = X + LTE  */
408:   PetscReal      wltea,wlter;

412:   if (!th->vec_sol_prev) {*wlte = -1; return(0);}
413:   /* Cannot compute LTE in first step or in restart after event */
414:   if (ts->steprestart) {*wlte = -1; return(0);}
415:   /* Compute LTE using backward differences with non-constant time step */
416:   {
417:     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
418:     PetscReal   a = 1 + h_prev/h;
419:     PetscScalar scal[3]; Vec vecs[3];
420:     scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
421:     vecs[0] = X;    vecs[1] = th->X0;   vecs[2] = th->vec_sol_prev;
422:     VecCopy(X,Y);
423:     VecMAXPY(Y,3,scal,vecs);
424:     TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
425:   }
426:   if (order) *order = 2;
427:   return(0);
428: }

430: static PetscErrorCode TSRollBack_Theta(TS ts)
431: {
432:   TS_Theta       *th = (TS_Theta*)ts->data;
433:   PetscInt       ntlm,ncost;

437:   VecCopy(th->X0,ts->vec_sol);
438:   if (ts->vec_costintegral && ts->costintegralfwd) {
439:     VecCopy(th->VecCostIntegral0,ts->vec_costintegral);
440:   }
441:   th->status = TS_STEP_INCOMPLETE;

443:   for (ntlm=0;ntlm<th->num_tlm;ntlm++) {
444:     VecCopy(th->VecsFwdSensi0[ntlm],ts->vecs_fwdsensipacked[ntlm]);
445:   }
446:   if (ts->vecs_integral_sensi) {
447:     for (ncost=0;ncost<ts->numcost;ncost++) {
448:       VecCopy(th->VecsIntegralSensi0[ncost],ts->vecs_integral_sensi[ncost]);
449:     }
450:   }
451:   if (ts->vecs_integral_sensip) {
452:     for (ncost=0;ncost<ts->numcost;ncost++) {
453:       VecCopy(th->VecsIntegralSensip0[ncost],ts->vecs_integral_sensip[ncost]);
454:     }
455:   }
456:   return(0);
457: }

459: static PetscErrorCode TSThetaIntegrandDerivative(TS ts,PetscInt numtlm,Vec DRncostDY,Vec* s,Vec DRncostDP,Vec VecIntegrandDerivative)
460: {
461:   PetscInt    ntlm,low,high;
462:   PetscScalar *v;


467:   VecGetOwnershipRange(VecIntegrandDerivative,&low,&high);
468:   VecGetArray(VecIntegrandDerivative,&v);
469:   for (ntlm=low; ntlm<high; ntlm++) {
470:     VecDot(DRncostDY,s[ntlm],&v[ntlm-low]);
471:   }
472:   VecRestoreArray(VecIntegrandDerivative,&v);
473:   if (DRncostDP) {
474:     VecAXPY(VecIntegrandDerivative,1.,DRncostDP);
475:   }
476:   return(0);
477: }

479: static PetscErrorCode TSForwardStep_Theta(TS ts)
480: {
481:   TS_Theta       *th = (TS_Theta*)ts->data;
482:   Vec            *VecsDeltaFwdSensi  = th->VecsDeltaFwdSensi;
483:   PetscInt       ntlm,ncost,np;
484:   KSP            ksp;
485:   Mat            J,Jp;
486:   PetscReal      shift;

490:   for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
491:     VecCopy(ts->vecs_fwdsensipacked[ntlm],th->VecsFwdSensi0[ntlm]);
492:   }
493:   for (ncost=0; ncost<ts->numcost; ncost++) {
494:     if (ts->vecs_integral_sensi) {
495:       VecCopy(ts->vecs_integral_sensi[ncost],th->VecsIntegralSensi0[ncost]);
496:     }
497:     if (ts->vecs_integral_sensip) {
498:       VecCopy(ts->vecs_integral_sensip[ncost],th->VecsIntegralSensip0[ncost]);
499:     }
500:   }

502:   SNESGetKSP(ts->snes,&ksp);
503:   TSGetIJacobian(ts,&J,&Jp,NULL,NULL);

505:   /* Build RHS */
506:   if (th->endpoint) { /* 2-stage method*/
507:     shift = 1./((th->Theta-1.)*th->time_step);
508:     TSComputeIJacobian(ts,th->ptime,th->X0,th->Xdot,shift,J,Jp,PETSC_FALSE);

510:     for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
511:       MatMult(J,ts->vecs_fwdsensipacked[ntlm],VecsDeltaFwdSensi[ntlm]);
512:       VecScale(VecsDeltaFwdSensi[ntlm],(th->Theta-1.)/th->Theta);
513:     }
514:     /* Add the f_p forcing terms */
515:     TSForwardComputeRHSJacobianP(ts,th->ptime,th->X0,ts->vecs_jacp);
516:     for (np=0; np<ts->num_parameters; np++) {
517:       VecAXPY(VecsDeltaFwdSensi[np+ts->num_initialvalues],(1.-th->Theta)/th->Theta,ts->vecs_jacp[np]);
518:     }
519:     TSForwardComputeRHSJacobianP(ts,th->stage_time,ts->vec_sol,ts->vecs_jacp);
520:     for (np=0; np<ts->num_parameters; np++) {
521:       VecAXPY(VecsDeltaFwdSensi[np+ts->num_initialvalues],1,ts->vecs_jacp[np]);
522:     }
523:   } else { /* 1-stage method */
524:     shift = 0.0;
525:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
526:     for (ntlm=0; ntlm<ts->num_parameters; ntlm++) {
527:       MatMult(J,ts->vecs_fwdsensipacked[ntlm],VecsDeltaFwdSensi[ntlm]);
528:       VecScale(VecsDeltaFwdSensi[ntlm],-1);
529:     }
530:     /* Add the f_p forcing terms */
531:     TSForwardComputeRHSJacobianP(ts,th->stage_time,th->X,ts->vecs_jacp);
532:     for (np=0; np<ts->num_parameters; np++) {
533:       VecAXPY(VecsDeltaFwdSensi[np+ts->num_initialvalues],1,ts->vecs_jacp[np]);
534:     }
535:   }

537:   /* Build LHS */
538:   shift = 1/(th->Theta*th->time_step);
539:   if (th->endpoint) {
540:     TSComputeIJacobian(ts,th->stage_time,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE);
541:   } else {
542:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
543:   }
544:   KSPSetOperators(ksp,J,Jp);

546:   /*
547:     Evaluate the first stage of integral gradients with the 2-stage method:
548:     drdy|t_n*S(t_n) + drdp|t_n
549:     This is done before the linear solve because the sensitivity variable S(t_n) will be propagated to S(t_{n+1})
550:     It is important to note that sensitivitesi to parameters (sensip) and sensitivities to initial values (sensi) are independent of each other, but integral sensip relies on sensip and integral sensi requires sensi
551:    */
552:   if (th->endpoint) { /* 2-stage method only */
553:     if (ts->vecs_integral_sensi || ts->vecs_integral_sensip) {
554:       TSAdjointComputeDRDYFunction(ts,th->ptime,th->X0,ts->vecs_drdy);
555:     }
556:     if (ts->vecs_integral_sensip) {
557:       TSAdjointComputeDRDPFunction(ts,th->ptime,th->X0,ts->vecs_drdp);
558:       for (ncost=0; ncost<ts->numcost; ncost++) {
559:         TSThetaIntegrandDerivative(ts,ts->num_parameters,ts->vecs_drdy[ncost],&ts->vecs_fwdsensipacked[ts->num_initialvalues],ts->vecs_drdp[ncost],th->VecIntegralSensipTemp);
560:         VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step*(1.-th->Theta),th->VecIntegralSensipTemp);
561:       }
562:     }
563:     if (ts->vecs_integral_sensi) {
564:       for (ncost=0; ncost<ts->numcost; ncost++) {
565:         TSThetaIntegrandDerivative(ts,ts->num_initialvalues,ts->vecs_drdy[ncost],ts->vecs_fwdsensipacked,NULL,th->VecIntegralSensiTemp);
566:         VecAXPY(ts->vecs_integral_sensi[ncost],th->time_step*(1.-th->Theta),th->VecIntegralSensiTemp);
567:       }
568:     }
569:   }

571:   /* Solve the tangent linear equation for forward sensitivities to parameters */
572:   for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
573:     if (th->endpoint) {
574:       KSPSolve(ksp,VecsDeltaFwdSensi[ntlm],ts->vecs_fwdsensipacked[ntlm]);
575:     } else {
576:       KSPSolve(ksp,VecsDeltaFwdSensi[ntlm],VecsDeltaFwdSensi[ntlm]);
577:       VecAXPY(ts->vecs_fwdsensipacked[ntlm],1.,VecsDeltaFwdSensi[ntlm]);
578:     }
579:   }
580:   /* Evaluate the second stage of integral gradients with the 2-stage method:
581:     drdy|t_{n+1}*S(t_{n+1}) + drdp|t_{n+1}
582:   */
583:   if (ts->vecs_integral_sensi || ts->vecs_integral_sensip) {
584:     TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
585:   }
586:   if (ts->vecs_integral_sensip) {
587:     TSAdjointComputeDRDPFunction(ts,th->stage_time,th->X,ts->vecs_drdp);
588:     for (ncost=0; ncost<ts->numcost; ncost++) {
589:       TSThetaIntegrandDerivative(ts,ts->num_parameters,ts->vecs_drdy[ncost],&ts->vecs_fwdsensipacked[ts->num_initialvalues],ts->vecs_drdp[ncost],th->VecIntegralSensipTemp);
590:       if (th->endpoint) {
591:         VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step*th->Theta,th->VecIntegralSensipTemp);
592:       } else {
593:         VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step,th->VecIntegralSensipTemp);
594:       }
595:     }
596:   }
597:   if (ts->vecs_integral_sensi) {
598:     for (ncost=0; ncost<ts->numcost; ncost++) {
599:       TSThetaIntegrandDerivative(ts,ts->num_initialvalues,ts->vecs_drdy[ncost],ts->vecs_fwdsensipacked,NULL,th->VecIntegralSensiTemp);
600:       if (th->endpoint) {
601:         VecAXPY(ts->vecs_integral_sensi[ncost],th->time_step*th->Theta,th->VecIntegralSensiTemp);
602:       } else {
603:         VecAXPY(ts->vecs_integral_sensi[ncost],th->time_step,th->VecIntegralSensiTemp);
604:       }
605:     }
606:   }

608:   if (!th->endpoint) { /* VecsDeltaFwdSensip are the increment for the stage values for the 1-stage method */
609:     for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
610:       VecAXPY(ts->vecs_fwdsensipacked[ntlm],(1.-th->Theta)/th->Theta,VecsDeltaFwdSensi[ntlm]);
611:     }
612:   }
613:   return(0);
614: }

616: /*------------------------------------------------------------*/
617: static PetscErrorCode TSReset_Theta(TS ts)
618: {
619:   TS_Theta       *th = (TS_Theta*)ts->data;

623:   VecDestroy(&th->X);
624:   VecDestroy(&th->Xdot);
625:   VecDestroy(&th->X0);
626:   VecDestroy(&th->affine);

628:   VecDestroy(&th->vec_sol_prev);
629:   VecDestroy(&th->vec_lte_work);

631:   VecDestroy(&th->VecCostIntegral0);
632:   if (ts->forward_solve) {
633:     VecDestroyVecs(th->num_tlm,&th->VecsDeltaFwdSensi);
634:     VecDestroyVecs(th->num_tlm,&th->VecsFwdSensi0);
635:     if (ts->vecs_integral_sensi) {
636:       VecDestroy(&th->VecIntegralSensiTemp);
637:       VecDestroyVecs(ts->numcost,&th->VecsIntegralSensi0);
638:     }
639:     if (ts->vecs_integral_sensip) {
640:       VecDestroy(&th->VecIntegralSensipTemp);
641:       VecDestroyVecs(ts->numcost,&th->VecsIntegralSensip0);
642:     }
643:   }
644:   VecDestroyVecs(ts->numcost,&th->VecsDeltaLam);
645:   VecDestroyVecs(ts->numcost,&th->VecsDeltaMu);
646:   VecDestroyVecs(ts->numcost,&th->VecsSensiTemp);

648:   return(0);
649: }

651: static PetscErrorCode TSDestroy_Theta(TS ts)
652: {

656:   TSReset_Theta(ts);
657:   if (ts->dm) {
658:     DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
659:     DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
660:   }
661:   PetscFree(ts->data);
662:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
663:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
664:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
665:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
666:   return(0);
667: }

669: /*
670:   This defines the nonlinear equation that is to be solved with SNES
671:   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
672: */
673: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
674: {
675:   TS_Theta       *th = (TS_Theta*)ts->data;
677:   Vec            X0,Xdot;
678:   DM             dm,dmsave;
679:   PetscReal      shift = 1/(th->Theta*ts->time_step);

682:   SNESGetDM(snes,&dm);
683:   /* When using the endpoint variant, this is actually 1/Theta * Xdot */
684:   TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
685:   VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);

687:   /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
688:   dmsave = ts->dm;
689:   ts->dm = dm;
690:   TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
691:   ts->dm = dmsave;
692:   TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
693:   return(0);
694: }

696: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat A,Mat B,TS ts)
697: {
698:   TS_Theta       *th = (TS_Theta*)ts->data;
700:   Vec            Xdot;
701:   DM             dm,dmsave;
702:   PetscReal      shift = 1/(th->Theta*ts->time_step);

705:   SNESGetDM(snes,&dm);
706:   /* Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
707:   TSThetaGetX0AndXdot(ts,dm,NULL,&Xdot);

709:   dmsave = ts->dm;
710:   ts->dm = dm;
711:   TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
712:   ts->dm = dmsave;
713:   TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
714:   return(0);
715: }

717: static PetscErrorCode TSForwardSetUp_Theta(TS ts)
718: {
719:   TS_Theta       *th = (TS_Theta*)ts->data;

723:   /* combine sensitivities to parameters and sensitivities to initial values into one array */
724:   th->num_tlm = ts->num_parameters + ts->num_initialvalues;
725:   VecDuplicateVecs(ts->vecs_fwdsensipacked[0],th->num_tlm,&th->VecsDeltaFwdSensi);
726:   if (ts->vecs_integral_sensi) {
727:     VecDuplicate(ts->vecs_integral_sensi[0],&th->VecIntegralSensiTemp);
728:   }
729:   if (ts->vecs_integral_sensip) {
730:     VecDuplicate(ts->vecs_integral_sensip[0],&th->VecIntegralSensipTemp);
731:   }
732:   /* backup sensitivity results for roll-backs */
733:   VecDuplicateVecs(ts->vecs_fwdsensipacked[0],th->num_tlm,&th->VecsFwdSensi0);
734:   if (ts->vecs_integral_sensi) {
735:     VecDuplicateVecs(ts->vecs_integral_sensi[0],ts->numcost,&th->VecsIntegralSensi0);
736:   }
737:   if (ts->vecs_integral_sensip) {
738:     VecDuplicateVecs(ts->vecs_integral_sensip[0],ts->numcost,&th->VecsIntegralSensip0);
739:   }
740:   return(0);
741: }

743: static PetscErrorCode TSSetUp_Theta(TS ts)
744: {
745:   TS_Theta       *th = (TS_Theta*)ts->data;
746:   PetscBool      match;

750:   if (!th->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */
751:     VecDuplicate(ts->vec_costintegral,&th->VecCostIntegral0);
752:   }
753:   if (!th->X) {
754:     VecDuplicate(ts->vec_sol,&th->X);
755:   }
756:   if (!th->Xdot) {
757:     VecDuplicate(ts->vec_sol,&th->Xdot);
758:   }
759:   if (!th->X0) {
760:     VecDuplicate(ts->vec_sol,&th->X0);
761:   }
762:   if (th->endpoint) {
763:     VecDuplicate(ts->vec_sol,&th->affine);
764:   }

766:   th->order = (th->Theta == 0.5) ? 2 : 1;

768:   TSGetDM(ts,&ts->dm);
769:   DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
770:   DMSubDomainHookAdd(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);

772:   TSGetAdapt(ts,&ts->adapt);
773:   TSAdaptCandidatesClear(ts->adapt);
774:   PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
775:   if (!match) {
776:     VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
777:     VecDuplicate(ts->vec_sol,&th->vec_lte_work);
778:   }
779:   TSGetSNES(ts,&ts->snes);
780:   return(0);
781: }

783: /*------------------------------------------------------------*/

785: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
786: {
787:   TS_Theta       *th = (TS_Theta*)ts->data;

791:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam);
792:   if(ts->vecs_sensip) {
793:     VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&th->VecsDeltaMu);
794:   }
795:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsSensiTemp);
796:   return(0);
797: }

799: static PetscErrorCode TSSetFromOptions_Theta(PetscOptionItems *PetscOptionsObject,TS ts)
800: {
801:   TS_Theta       *th = (TS_Theta*)ts->data;

805:   PetscOptionsHead(PetscOptionsObject,"Theta ODE solver options");
806:   {
807:     PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,NULL);
808:     PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,NULL);
809:     PetscOptionsBool("-ts_theta_initial_guess_extrapolate","Extrapolate stage initial guess from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,NULL);
810:   }
811:   PetscOptionsTail();
812:   return(0);
813: }

815: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
816: {
817:   TS_Theta       *th = (TS_Theta*)ts->data;
818:   PetscBool      iascii;

822:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
823:   if (iascii) {
824:     PetscViewerASCIIPrintf(viewer,"  Theta=%g\n",(double)th->Theta);
825:     PetscViewerASCIIPrintf(viewer,"  Extrapolation=%s\n",th->extrapolate ? "yes" : "no");
826:   }
827:   return(0);
828: }

830: static PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
831: {
832:   TS_Theta *th = (TS_Theta*)ts->data;

835:   *theta = th->Theta;
836:   return(0);
837: }

839: static PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta)
840: {
841:   TS_Theta *th = (TS_Theta*)ts->data;

844:   if (theta <= 0 || 1 < theta) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Theta %g not in range (0,1]",(double)theta);
845:   th->Theta = theta;
846:   th->order = (th->Theta == 0.5) ? 2 : 1;
847:   return(0);
848: }

850: static PetscErrorCode TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
851: {
852:   TS_Theta *th = (TS_Theta*)ts->data;

855:   *endpoint = th->endpoint;
856:   return(0);
857: }

859: static PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
860: {
861:   TS_Theta *th = (TS_Theta*)ts->data;

864:   th->endpoint = flg;
865:   return(0);
866: }

868: #if defined(PETSC_HAVE_COMPLEX)
869: static PetscErrorCode TSComputeLinearStability_Theta(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
870: {
871:   PetscComplex   z   = xr + xi*PETSC_i,f;
872:   TS_Theta       *th = (TS_Theta*)ts->data;
873:   const PetscReal one = 1.0;

876:   f   = (one + (one - th->Theta)*z)/(one - th->Theta*z);
877:   *yr = PetscRealPartComplex(f);
878:   *yi = PetscImaginaryPartComplex(f);
879:   return(0);
880: }
881: #endif

883: static PetscErrorCode TSGetStages_Theta(TS ts,PetscInt *ns,Vec **Y)
884: {
885:   TS_Theta     *th = (TS_Theta*)ts->data;

888:   if (ns) *ns = 1;
889:   if (Y)  *Y  = th->endpoint ? &(th->X0) : &(th->X);
890:   return(0);
891: }

893: /* ------------------------------------------------------------ */
894: /*MC
895:       TSTHETA - DAE solver using the implicit Theta method

897:    Level: beginner

899:    Options Database:
900: +  -ts_theta_theta <Theta> - Location of stage (0<Theta<=1)
901: .  -ts_theta_endpoint <flag> - Use the endpoint (like Crank-Nicholson) instead of midpoint form of the Theta method
902: -  -ts_theta_initial_guess_extrapolate <flg> - Extrapolate stage initial guess from previous solution (sometimes unstable)

904:    Notes:
905: $  -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
906: $  -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
907: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)

909:    This method can be applied to DAE.

911:    This method is cast as a 1-stage implicit Runge-Kutta method.

913: .vb
914:   Theta | Theta
915:   -------------
916:         |  1
917: .ve

919:    For the default Theta=0.5, this is also known as the implicit midpoint rule.

921:    When the endpoint variant is chosen, the method becomes a 2-stage method with first stage explicit:

923: .vb
924:   0 | 0         0
925:   1 | 1-Theta   Theta
926:   -------------------
927:     | 1-Theta   Theta
928: .ve

930:    For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN).

932:    To apply a diagonally implicit RK method to DAE, the stage formula

934: $  Y_i = X + h sum_j a_ij Y'_j

936:    is interpreted as a formula for Y'_i in terms of Y_i and known values (Y'_j, j<i)

938: .seealso:  TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint()

940: M*/
941: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
942: {
943:   TS_Theta       *th;

947:   ts->ops->reset           = TSReset_Theta;
948:   ts->ops->destroy         = TSDestroy_Theta;
949:   ts->ops->view            = TSView_Theta;
950:   ts->ops->setup           = TSSetUp_Theta;
951:   ts->ops->adjointsetup    = TSAdjointSetUp_Theta;
952:   ts->ops->step            = TSStep_Theta;
953:   ts->ops->interpolate     = TSInterpolate_Theta;
954:   ts->ops->evaluatewlte    = TSEvaluateWLTE_Theta;
955:   ts->ops->rollback        = TSRollBack_Theta;
956:   ts->ops->setfromoptions  = TSSetFromOptions_Theta;
957:   ts->ops->snesfunction    = SNESTSFormFunction_Theta;
958:   ts->ops->snesjacobian    = SNESTSFormJacobian_Theta;
959: #if defined(PETSC_HAVE_COMPLEX)
960:   ts->ops->linearstability = TSComputeLinearStability_Theta;
961: #endif
962:   ts->ops->getstages       = TSGetStages_Theta;
963:   ts->ops->adjointstep     = TSAdjointStep_Theta;
964:   ts->ops->adjointintegral = TSAdjointCostIntegral_Theta;
965:   ts->ops->forwardintegral = TSForwardCostIntegral_Theta;
966:   ts->default_adapt_type   = TSADAPTNONE;
967:   ts->ops->forwardsetup    = TSForwardSetUp_Theta;
968:   ts->ops->forwardstep     = TSForwardStep_Theta;

970:   ts->usessnes = PETSC_TRUE;

972:   PetscNewLog(ts,&th);
973:   ts->data = (void*)th;

975:   th->VecsDeltaLam   = NULL;
976:   th->VecsDeltaMu    = NULL;
977:   th->VecsSensiTemp  = NULL;

979:   th->extrapolate = PETSC_FALSE;
980:   th->Theta       = 0.5;
981:   th->order       = 2;
982:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",TSThetaGetTheta_Theta);
983:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",TSThetaSetTheta_Theta);
984:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",TSThetaGetEndpoint_Theta);
985:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",TSThetaSetEndpoint_Theta);
986:   return(0);
987: }

989: /*@
990:   TSThetaGetTheta - Get the abscissa of the stage in (0,1].

992:   Not Collective

994:   Input Parameter:
995: .  ts - timestepping context

997:   Output Parameter:
998: .  theta - stage abscissa

1000:   Note:
1001:   Use of this function is normally only required to hack TSTHETA to use a modified integration scheme.

1003:   Level: Advanced

1005: .seealso: TSThetaSetTheta()
1006: @*/
1007: PetscErrorCode  TSThetaGetTheta(TS ts,PetscReal *theta)
1008: {

1014:   PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
1015:   return(0);
1016: }

1018: /*@
1019:   TSThetaSetTheta - Set the abscissa of the stage in (0,1].

1021:   Not Collective

1023:   Input Parameter:
1024: +  ts - timestepping context
1025: -  theta - stage abscissa

1027:   Options Database:
1028: .  -ts_theta_theta <theta>

1030:   Level: Intermediate

1032: .seealso: TSThetaGetTheta()
1033: @*/
1034: PetscErrorCode  TSThetaSetTheta(TS ts,PetscReal theta)
1035: {

1040:   PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
1041:   return(0);
1042: }

1044: /*@
1045:   TSThetaGetEndpoint - Gets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).

1047:   Not Collective

1049:   Input Parameter:
1050: .  ts - timestepping context

1052:   Output Parameter:
1053: .  endpoint - PETSC_TRUE when using the endpoint variant

1055:   Level: Advanced

1057: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
1058: @*/
1059: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
1060: {

1066:   PetscUseMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
1067:   return(0);
1068: }

1070: /*@
1071:   TSThetaSetEndpoint - Sets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).

1073:   Not Collective

1075:   Input Parameter:
1076: +  ts - timestepping context
1077: -  flg - PETSC_TRUE to use the endpoint variant

1079:   Options Database:
1080: .  -ts_theta_endpoint <flg>

1082:   Level: Intermediate

1084: .seealso: TSTHETA, TSCN
1085: @*/
1086: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
1087: {

1092:   PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
1093:   return(0);
1094: }

1096: /*
1097:  * TSBEULER and TSCN are straightforward specializations of TSTHETA.
1098:  * The creation functions for these specializations are below.
1099:  */

1101: static PetscErrorCode TSSetUp_BEuler(TS ts)
1102: {
1103:   TS_Theta       *th = (TS_Theta*)ts->data;

1107:   if (th->Theta != 1.0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_OPT_OVERWRITE,"Can not change the default value (1) of theta when using backward Euler\n");
1108:   if (th->endpoint) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_OPT_OVERWRITE,"Can not change to the endpoint form of the Theta methods when using backward Euler\n");
1109:   TSSetUp_Theta(ts);
1110:   return(0);
1111: }

1113: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
1114: {
1116:   return(0);
1117: }

1119: /*MC
1120:       TSBEULER - ODE solver using the implicit backward Euler method

1122:   Level: beginner

1124:   Notes:
1125:   TSBEULER is equivalent to TSTHETA with Theta=1.0

1127: $  -ts_type theta -ts_theta_theta 1.0

1129: .seealso:  TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA

1131: M*/
1132: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
1133: {

1137:   TSCreate_Theta(ts);
1138:   TSThetaSetTheta(ts,1.0);
1139:   TSThetaSetEndpoint(ts,PETSC_FALSE);
1140:   ts->ops->setup = TSSetUp_BEuler;
1141:   ts->ops->view  = TSView_BEuler;
1142:   return(0);
1143: }

1145: static PetscErrorCode TSSetUp_CN(TS ts)
1146: {
1147:   TS_Theta       *th = (TS_Theta*)ts->data;

1151:   if (th->Theta != 0.5) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_OPT_OVERWRITE,"Can not change the default value (0.5) of theta when using Crank-Nicolson\n");
1152:   if (!th->endpoint) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_OPT_OVERWRITE,"Can not change to the midpoint form of the Theta methods when using Crank-Nicolson\n");
1153:   TSSetUp_Theta(ts);
1154:   return(0);
1155: }

1157: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
1158: {
1160:   return(0);
1161: }

1163: /*MC
1164:       TSCN - ODE solver using the implicit Crank-Nicolson method.

1166:   Level: beginner

1168:   Notes:
1169:   TSCN is equivalent to TSTHETA with Theta=0.5 and the "endpoint" option set. I.e.

1171: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint

1173: .seealso:  TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA

1175: M*/
1176: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
1177: {

1181:   TSCreate_Theta(ts);
1182:   TSThetaSetTheta(ts,0.5);
1183:   TSThetaSetEndpoint(ts,PETSC_TRUE);
1184:   ts->ops->setup = TSSetUp_CN;
1185:   ts->ops->view  = TSView_CN;
1186:   return(0);
1187: }