Actual source code: theta.c

petsc-master 2017-09-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:   PetscFree(ts->data);
658:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
659:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
660:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
661:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
662:   return(0);
663: }

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

678:   SNESGetDM(snes,&dm);
679:   /* When using the endpoint variant, this is actually 1/Theta * Xdot */
680:   TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
681:   VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);

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

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

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

705:   dmsave = ts->dm;
706:   ts->dm = dm;
707:   TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
708:   ts->dm = dmsave;
709:   TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
710:   return(0);
711: }

713: static PetscErrorCode TSForwardSetUp_Theta(TS ts)
714: {
715:   TS_Theta       *th = (TS_Theta*)ts->data;

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

739: static PetscErrorCode TSSetUp_Theta(TS ts)
740: {
741:   TS_Theta       *th = (TS_Theta*)ts->data;
742:   PetscBool      match;

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

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

764:   TSGetDM(ts,&ts->dm);
765:   DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
766:   DMSubDomainHookAdd(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);

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

779: /*------------------------------------------------------------*/

781: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
782: {
783:   TS_Theta       *th = (TS_Theta*)ts->data;

787:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam);
788:   if(ts->vecs_sensip) {
789:     VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&th->VecsDeltaMu);
790:   }
791:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsSensiTemp);
792:   return(0);
793: }

795: static PetscErrorCode TSSetFromOptions_Theta(PetscOptionItems *PetscOptionsObject,TS ts)
796: {
797:   TS_Theta       *th = (TS_Theta*)ts->data;

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

811: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
812: {
813:   TS_Theta       *th = (TS_Theta*)ts->data;
814:   PetscBool      iascii;

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

826: static PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
827: {
828:   TS_Theta *th = (TS_Theta*)ts->data;

831:   *theta = th->Theta;
832:   return(0);
833: }

835: static PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta)
836: {
837:   TS_Theta *th = (TS_Theta*)ts->data;

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

846: static PetscErrorCode TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
847: {
848:   TS_Theta *th = (TS_Theta*)ts->data;

851:   *endpoint = th->endpoint;
852:   return(0);
853: }

855: static PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
856: {
857:   TS_Theta *th = (TS_Theta*)ts->data;

860:   th->endpoint = flg;
861:   return(0);
862: }

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

872:   f   = (one + (one - th->Theta)*z)/(one - th->Theta*z);
873:   *yr = PetscRealPartComplex(f);
874:   *yi = PetscImaginaryPartComplex(f);
875:   return(0);
876: }
877: #endif

879: static PetscErrorCode TSGetStages_Theta(TS ts,PetscInt *ns,Vec **Y)
880: {
881:   TS_Theta     *th = (TS_Theta*)ts->data;

884:   if (ns) *ns = 1;
885:   if (Y)  *Y  = th->endpoint ? &(th->X0) : &(th->X);
886:   return(0);
887: }

889: /* ------------------------------------------------------------ */
890: /*MC
891:       TSTHETA - DAE solver using the implicit Theta method

893:    Level: beginner

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

900:    Notes:
901: $  -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
902: $  -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
903: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)

905:    This method can be applied to DAE.

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

909: .vb
910:   Theta | Theta
911:   -------------
912:         |  1
913: .ve

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

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

919: .vb
920:   0 | 0         0
921:   1 | 1-Theta   Theta
922:   -------------------
923:     | 1-Theta   Theta
924: .ve

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

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

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

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

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

936: M*/
937: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
938: {
939:   TS_Theta       *th;

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

966:   ts->usessnes = PETSC_TRUE;

968:   PetscNewLog(ts,&th);
969:   ts->data = (void*)th;

971:   th->VecsDeltaLam   = NULL;
972:   th->VecsDeltaMu    = NULL;
973:   th->VecsSensiTemp  = NULL;

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

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

988:   Not Collective

990:   Input Parameter:
991: .  ts - timestepping context

993:   Output Parameter:
994: .  theta - stage abscissa

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

999:   Level: Advanced

1001: .seealso: TSThetaSetTheta()
1002: @*/
1003: PetscErrorCode  TSThetaGetTheta(TS ts,PetscReal *theta)
1004: {

1010:   PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
1011:   return(0);
1012: }

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

1017:   Not Collective

1019:   Input Parameter:
1020: +  ts - timestepping context
1021: -  theta - stage abscissa

1023:   Options Database:
1024: .  -ts_theta_theta <theta>

1026:   Level: Intermediate

1028: .seealso: TSThetaGetTheta()
1029: @*/
1030: PetscErrorCode  TSThetaSetTheta(TS ts,PetscReal theta)
1031: {

1036:   PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
1037:   return(0);
1038: }

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

1043:   Not Collective

1045:   Input Parameter:
1046: .  ts - timestepping context

1048:   Output Parameter:
1049: .  endpoint - PETSC_TRUE when using the endpoint variant

1051:   Level: Advanced

1053: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
1054: @*/
1055: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
1056: {

1062:   PetscUseMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
1063:   return(0);
1064: }

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

1069:   Not Collective

1071:   Input Parameter:
1072: +  ts - timestepping context
1073: -  flg - PETSC_TRUE to use the endpoint variant

1075:   Options Database:
1076: .  -ts_theta_endpoint <flg>

1078:   Level: Intermediate

1080: .seealso: TSTHETA, TSCN
1081: @*/
1082: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
1083: {

1088:   PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
1089:   return(0);
1090: }

1092: /*
1093:  * TSBEULER and TSCN are straightforward specializations of TSTHETA.
1094:  * The creation functions for these specializations are below.
1095:  */

1097: static PetscErrorCode TSSetUp_BEuler(TS ts)
1098: {
1099:   TS_Theta       *th = (TS_Theta*)ts->data;

1103:   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");
1104:   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");
1105:   TSSetUp_Theta(ts);
1106:   return(0);
1107: }

1109: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
1110: {
1112:   return(0);
1113: }

1115: /*MC
1116:       TSBEULER - ODE solver using the implicit backward Euler method

1118:   Level: beginner

1120:   Notes:
1121:   TSBEULER is equivalent to TSTHETA with Theta=1.0

1123: $  -ts_type theta -ts_theta_theta 1.0

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

1127: M*/
1128: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
1129: {

1133:   TSCreate_Theta(ts);
1134:   TSThetaSetTheta(ts,1.0);
1135:   TSThetaSetEndpoint(ts,PETSC_FALSE);
1136:   ts->ops->setup = TSSetUp_BEuler;
1137:   ts->ops->view  = TSView_BEuler;
1138:   return(0);
1139: }

1141: static PetscErrorCode TSSetUp_CN(TS ts)
1142: {
1143:   TS_Theta       *th = (TS_Theta*)ts->data;

1147:   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");
1148:   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");
1149:   TSSetUp_Theta(ts);
1150:   return(0);
1151: }

1153: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
1154: {
1156:   return(0);
1157: }

1159: /*MC
1160:       TSCN - ODE solver using the implicit Crank-Nicolson method.

1162:   Level: beginner

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

1167: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint

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

1171: M*/
1172: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
1173: {

1177:   TSCreate_Theta(ts);
1178:   TSThetaSetTheta(ts,0.5);
1179:   TSThetaSetEndpoint(ts,PETSC_TRUE);
1180:   ts->ops->setup = TSSetUp_CN;
1181:   ts->ops->view  = TSView_CN;
1182:   return(0);
1183: }