Actual source code: theta.c

petsc-master 2018-06-22
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 multiplied with Jacobian transpose */
 28:   Mat          MatDeltaFwdSensip;        /* Increment of the forward sensitivity at stage */
 29:   Vec          VecDeltaFwdSensipTemp;    /* Working vector for holding one column of the sensitivity matrix */
 30:   Vec          VecDeltaFwdSensipTemp2;   /* Additional working vector for endpoint case */
 31:   Mat          MatFwdSensip0;            /* backup for roll-backs due to events */
 32:   Vec          VecIntegralSensipTemp;    /* Working vector for forward integral sensitivity */
 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 TSTheta_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:     TSTheta_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       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;
442:   if (ts->mat_sensip) {
443:     MatCopy(th->MatFwdSensip0,ts->mat_sensip,SAME_NONZERO_PATTERN);
444:   }
445:   if (ts->vecs_integral_sensip) {
446:     for (ncost=0;ncost<ts->numcost;ncost++) {
447:       VecCopy(th->VecsIntegralSensip0[ncost],ts->vecs_integral_sensip[ncost]);
448:     }
449:   }
450:   return(0);
451: }

453: static PetscErrorCode TSForwardStep_Theta(TS ts)
454: {
455:   TS_Theta       *th = (TS_Theta*)ts->data;
456:   Mat            MatDeltaFwdSensip = th->MatDeltaFwdSensip;
457:   Vec            VecDeltaFwdSensipTemp = th->VecDeltaFwdSensipTemp,VecDeltaFwdSensipTemp2 = th->VecDeltaFwdSensipTemp2;
458:   PetscInt       ncost,ntlm;
459:   KSP            ksp;
460:   Mat            J,Jp;
461:   PetscReal      shift;
462:   PetscScalar    *barr,*xarr;

466:   MatCopy(ts->mat_sensip,th->MatFwdSensip0,SAME_NONZERO_PATTERN);

468:   for (ncost=0; ncost<ts->numcost; ncost++) {
469:     if (ts->vecs_integral_sensip) {
470:       VecCopy(ts->vecs_integral_sensip[ncost],th->VecsIntegralSensip0[ncost]);
471:     }
472:   }

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

477:   /* Build RHS */
478:   if (th->endpoint) { /* 2-stage method*/
479:     shift = 1./((th->Theta-1.)*th->time_step);
480:     TSComputeIJacobian(ts,th->ptime,th->X0,th->Xdot,shift,J,Jp,PETSC_FALSE);
481:     MatMatMult(J,ts->mat_sensip,MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatDeltaFwdSensip);
482:     MatScale(MatDeltaFwdSensip,(th->Theta-1.)/th->Theta);

484:     /* Add the f_p forcing terms */
485:     TSAdjointComputeRHSJacobian(ts,th->ptime,th->X0,ts->Jacp);
486:     MatAXPY(MatDeltaFwdSensip,(1.-th->Theta)/th->Theta,ts->Jacp,SUBSET_NONZERO_PATTERN);
487:     TSAdjointComputeRHSJacobian(ts,th->stage_time,ts->vec_sol,ts->Jacp);
488:     MatAXPY(MatDeltaFwdSensip,1.,ts->Jacp,SUBSET_NONZERO_PATTERN);
489:   } else { /* 1-stage method */
490:     shift = 0.0;
491:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
492:     MatMatMult(J,ts->mat_sensip,MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatDeltaFwdSensip);
493:     MatScale(MatDeltaFwdSensip,-1.);

495:     /* Add the f_p forcing terms */
496:     TSAdjointComputeRHSJacobian(ts,th->stage_time,th->X,ts->Jacp);
497:     MatAXPY(MatDeltaFwdSensip,1.,ts->Jacp,SUBSET_NONZERO_PATTERN);
498:   }

500:   /* Build LHS */
501:   shift = 1/(th->Theta*th->time_step);
502:   if (th->endpoint) {
503:     TSComputeIJacobian(ts,th->stage_time,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE);
504:   } else {
505:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
506:   }
507:   KSPSetOperators(ksp,J,Jp);

509:   /*
510:     Evaluate the first stage of integral gradients with the 2-stage method:
511:     drdy|t_n*S(t_n) + drdp|t_n
512:     This is done before the linear solve because the sensitivity variable S(t_n) will be propagated to S(t_{n+1})
513:   */
514:   if (th->endpoint) { /* 2-stage method only */
515:     if (ts->vecs_integral_sensip) {
516:       TSAdjointComputeDRDYFunction(ts,th->ptime,th->X0,ts->vecs_drdy);
517:       TSAdjointComputeDRDPFunction(ts,th->ptime,th->X0,ts->vecs_drdp);
518:       for (ncost=0; ncost<ts->numcost; ncost++) {
519:         MatMultTranspose(ts->mat_sensip,ts->vecs_drdy[ncost],th->VecIntegralSensipTemp);
520:         VecAXPY(th->VecIntegralSensipTemp,1,ts->vecs_drdp[ncost]);
521:         VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step*(1.-th->Theta),th->VecIntegralSensipTemp);
522:       }
523:     }
524:   }

526:   /* Solve the tangent linear equation for forward sensitivities to parameters */
527:   for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
528:     MatDenseGetColumn(MatDeltaFwdSensip,ntlm,&barr);
529:     VecPlaceArray(VecDeltaFwdSensipTemp,barr);
530:     if (th->endpoint) {
531:       MatDenseGetColumn(ts->mat_sensip,ntlm,&xarr);
532:       VecPlaceArray(VecDeltaFwdSensipTemp2,xarr);
533:       KSPSolve(ksp,VecDeltaFwdSensipTemp,VecDeltaFwdSensipTemp2);
534:       VecResetArray(VecDeltaFwdSensipTemp2);
535:       MatDenseRestoreColumn(ts->mat_sensip,&xarr);
536:     } else {
537:       KSPSolve(ksp,VecDeltaFwdSensipTemp,VecDeltaFwdSensipTemp);
538:     }
539:     VecResetArray(VecDeltaFwdSensipTemp);
540:     MatDenseRestoreColumn(MatDeltaFwdSensip,&barr);
541:   }

543:   /*
544:     Evaluate the second stage of integral gradients with the 2-stage method:
545:     drdy|t_{n+1}*S(t_{n+1}) + drdp|t_{n+1}
546:   */
547:   if (ts->vecs_integral_sensip) {
548:     if (!th->endpoint) {
549:       MatAXPY(ts->mat_sensip,1,MatDeltaFwdSensip,SAME_NONZERO_PATTERN);
550:       TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
551:       TSAdjointComputeDRDPFunction(ts,th->stage_time,th->X,ts->vecs_drdp);
552:       for (ncost=0; ncost<ts->numcost; ncost++) {
553:         MatMultTranspose(ts->mat_sensip,ts->vecs_drdy[ncost],th->VecIntegralSensipTemp);
554:         VecAXPY(th->VecIntegralSensipTemp,1,ts->vecs_drdp[ncost]);
555:         VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step,th->VecIntegralSensipTemp);
556:       }
557:       MatAXPY(ts->mat_sensip,(1.-th->Theta)/th->Theta,MatDeltaFwdSensip,SAME_NONZERO_PATTERN);
558:     } else {
559:       TSAdjointComputeDRDYFunction(ts,th->stage_time,ts->vec_sol,ts->vecs_drdy);
560:       TSAdjointComputeDRDPFunction(ts,th->stage_time,ts->vec_sol,ts->vecs_drdp);
561:       for (ncost=0; ncost<ts->numcost; ncost++) {
562:         MatMultTranspose(ts->mat_sensip,ts->vecs_drdy[ncost],th->VecIntegralSensipTemp);
563:         VecAXPY(th->VecIntegralSensipTemp,1,ts->vecs_drdp[ncost]);
564:         VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step*th->Theta,th->VecIntegralSensipTemp);
565:       }
566:     }
567:   } else {
568:     if (!th->endpoint) {
569:       MatAXPY(ts->mat_sensip,1./th->Theta,MatDeltaFwdSensip,SAME_NONZERO_PATTERN);
570:     }
571:   }
572:   return(0);
573: }

575: /*------------------------------------------------------------*/
576: static PetscErrorCode TSReset_Theta(TS ts)
577: {
578:   TS_Theta       *th = (TS_Theta*)ts->data;

582:   VecDestroy(&th->X);
583:   VecDestroy(&th->Xdot);
584:   VecDestroy(&th->X0);
585:   VecDestroy(&th->affine);

587:   VecDestroy(&th->vec_sol_prev);
588:   VecDestroy(&th->vec_lte_work);

590:   VecDestroy(&th->VecCostIntegral0);
591:   if (ts->forward_solve) {
592:     if (ts->vecs_integral_sensip) {
593:       VecDestroy(&th->VecIntegralSensipTemp);
594:       VecDestroyVecs(ts->numcost,&th->VecsIntegralSensip0);
595:     }
596:     VecDestroy(&th->VecDeltaFwdSensipTemp);
597:     if (th->endpoint) {
598:       VecDestroy(&th->VecDeltaFwdSensipTemp2);
599:     }
600:     MatDestroy(&th->MatDeltaFwdSensip);
601:     MatDestroy(&th->MatFwdSensip0);
602:   }
603:   VecDestroyVecs(ts->numcost,&th->VecsDeltaLam);
604:   VecDestroyVecs(ts->numcost,&th->VecsDeltaMu);
605:   VecDestroyVecs(ts->numcost,&th->VecsSensiTemp);

607:   return(0);
608: }

610: static PetscErrorCode TSDestroy_Theta(TS ts)
611: {

615:   TSReset_Theta(ts);
616:   if (ts->dm) {
617:     DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
618:     DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
619:   }
620:   PetscFree(ts->data);
621:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
622:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
623:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
624:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
625:   return(0);
626: }

628: /*
629:   This defines the nonlinear equation that is to be solved with SNES
630:   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
631: */
632: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
633: {
634:   TS_Theta       *th = (TS_Theta*)ts->data;
636:   Vec            X0,Xdot;
637:   DM             dm,dmsave;
638:   PetscReal      shift = 1/(th->Theta*ts->time_step);

641:   SNESGetDM(snes,&dm);
642:   /* When using the endpoint variant, this is actually 1/Theta * Xdot */
643:   TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
644:   VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);

646:   /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
647:   dmsave = ts->dm;
648:   ts->dm = dm;
649:   TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
650:   ts->dm = dmsave;
651:   TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
652:   return(0);
653: }

655: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat A,Mat B,TS ts)
656: {
657:   TS_Theta       *th = (TS_Theta*)ts->data;
659:   Vec            Xdot;
660:   DM             dm,dmsave;
661:   PetscReal      shift = 1/(th->Theta*ts->time_step);

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

668:   dmsave = ts->dm;
669:   ts->dm = dm;
670:   TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
671:   ts->dm = dmsave;
672:   TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
673:   return(0);
674: }

676: static PetscErrorCode TSForwardSetUp_Theta(TS ts)
677: {
678:   TS_Theta       *th = (TS_Theta*)ts->data;

682:   /* combine sensitivities to parameters and sensitivities to initial values into one array */
683:   th->num_tlm = ts->num_parameters;
684:   MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatDeltaFwdSensip);
685:   if (ts->vecs_integral_sensip) {
686:     VecDuplicate(ts->vecs_integral_sensip[0],&th->VecIntegralSensipTemp);
687:   }
688:   /* backup sensitivity results for roll-backs */
689:   MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatFwdSensip0);

691:   if (ts->vecs_integral_sensip) {
692:     VecDuplicateVecs(ts->vecs_integral_sensip[0],ts->numcost,&th->VecsIntegralSensip0);
693:   }
694:   VecDuplicate(ts->vec_sol,&th->VecDeltaFwdSensipTemp);
695:   if (th->endpoint) {
696:     VecDuplicate(ts->vec_sol,&th->VecDeltaFwdSensipTemp2);
697:   }
698:   return(0);
699: }

701: static PetscErrorCode TSSetUp_Theta(TS ts)
702: {
703:   TS_Theta       *th = (TS_Theta*)ts->data;
704:   PetscBool      match;

708:   if (!th->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */
709:     VecDuplicate(ts->vec_costintegral,&th->VecCostIntegral0);
710:   }
711:   if (!th->X) {
712:     VecDuplicate(ts->vec_sol,&th->X);
713:   }
714:   if (!th->Xdot) {
715:     VecDuplicate(ts->vec_sol,&th->Xdot);
716:   }
717:   if (!th->X0) {
718:     VecDuplicate(ts->vec_sol,&th->X0);
719:   }
720:   if (th->endpoint) {
721:     VecDuplicate(ts->vec_sol,&th->affine);
722:   }

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

726:   TSGetDM(ts,&ts->dm);
727:   DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
728:   DMSubDomainHookAdd(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);

730:   TSGetAdapt(ts,&ts->adapt);
731:   TSAdaptCandidatesClear(ts->adapt);
732:   PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
733:   if (!match) {
734:     VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
735:     VecDuplicate(ts->vec_sol,&th->vec_lte_work);
736:   }
737:   TSGetSNES(ts,&ts->snes);
738:   return(0);
739: }

741: /*------------------------------------------------------------*/

743: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
744: {
745:   TS_Theta       *th = (TS_Theta*)ts->data;

749:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam);
750:   if(ts->vecs_sensip) {
751:     VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&th->VecsDeltaMu);
752:   }
753:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsSensiTemp);
754:   return(0);
755: }

757: static PetscErrorCode TSSetFromOptions_Theta(PetscOptionItems *PetscOptionsObject,TS ts)
758: {
759:   TS_Theta       *th = (TS_Theta*)ts->data;

763:   PetscOptionsHead(PetscOptionsObject,"Theta ODE solver options");
764:   {
765:     PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,NULL);
766:     PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,NULL);
767:     PetscOptionsBool("-ts_theta_initial_guess_extrapolate","Extrapolate stage initial guess from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,NULL);
768:   }
769:   PetscOptionsTail();
770:   return(0);
771: }

773: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
774: {
775:   TS_Theta       *th = (TS_Theta*)ts->data;
776:   PetscBool      iascii;

780:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
781:   if (iascii) {
782:     PetscViewerASCIIPrintf(viewer,"  Theta=%g\n",(double)th->Theta);
783:     PetscViewerASCIIPrintf(viewer,"  Extrapolation=%s\n",th->extrapolate ? "yes" : "no");
784:   }
785:   return(0);
786: }

788: static PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
789: {
790:   TS_Theta *th = (TS_Theta*)ts->data;

793:   *theta = th->Theta;
794:   return(0);
795: }

797: static PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta)
798: {
799:   TS_Theta *th = (TS_Theta*)ts->data;

802:   if (theta <= 0 || 1 < theta) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Theta %g not in range (0,1]",(double)theta);
803:   th->Theta = theta;
804:   th->order = (th->Theta == 0.5) ? 2 : 1;
805:   return(0);
806: }

808: static PetscErrorCode TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
809: {
810:   TS_Theta *th = (TS_Theta*)ts->data;

813:   *endpoint = th->endpoint;
814:   return(0);
815: }

817: static PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
818: {
819:   TS_Theta *th = (TS_Theta*)ts->data;

822:   th->endpoint = flg;
823:   return(0);
824: }

826: #if defined(PETSC_HAVE_COMPLEX)
827: static PetscErrorCode TSComputeLinearStability_Theta(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
828: {
829:   PetscComplex   z   = xr + xi*PETSC_i,f;
830:   TS_Theta       *th = (TS_Theta*)ts->data;
831:   const PetscReal one = 1.0;

834:   f   = (one + (one - th->Theta)*z)/(one - th->Theta*z);
835:   *yr = PetscRealPartComplex(f);
836:   *yi = PetscImaginaryPartComplex(f);
837:   return(0);
838: }
839: #endif

841: static PetscErrorCode TSGetStages_Theta(TS ts,PetscInt *ns,Vec **Y)
842: {
843:   TS_Theta     *th = (TS_Theta*)ts->data;

846:   if (ns) *ns = 1;
847:   if (Y)  *Y  = th->endpoint ? &(th->X0) : &(th->X);
848:   return(0);
849: }

851: /* ------------------------------------------------------------ */
852: /*MC
853:       TSTHETA - DAE solver using the implicit Theta method

855:    Level: beginner

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

862:    Notes:
863: $  -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
864: $  -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
865: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)

867:    This method can be applied to DAE.

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

871: .vb
872:   Theta | Theta
873:   -------------
874:         |  1
875: .ve

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

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

881: .vb
882:   0 | 0         0
883:   1 | 1-Theta   Theta
884:   -------------------
885:     | 1-Theta   Theta
886: .ve

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

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

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

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

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

898: M*/
899: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
900: {
901:   TS_Theta       *th;

905:   ts->ops->reset           = TSReset_Theta;
906:   ts->ops->destroy         = TSDestroy_Theta;
907:   ts->ops->view            = TSView_Theta;
908:   ts->ops->setup           = TSSetUp_Theta;
909:   ts->ops->adjointsetup    = TSAdjointSetUp_Theta;
910:   ts->ops->step            = TSStep_Theta;
911:   ts->ops->interpolate     = TSInterpolate_Theta;
912:   ts->ops->evaluatewlte    = TSEvaluateWLTE_Theta;
913:   ts->ops->rollback        = TSRollBack_Theta;
914:   ts->ops->setfromoptions  = TSSetFromOptions_Theta;
915:   ts->ops->snesfunction    = SNESTSFormFunction_Theta;
916:   ts->ops->snesjacobian    = SNESTSFormJacobian_Theta;
917: #if defined(PETSC_HAVE_COMPLEX)
918:   ts->ops->linearstability = TSComputeLinearStability_Theta;
919: #endif
920:   ts->ops->getstages       = TSGetStages_Theta;
921:   ts->ops->adjointstep     = TSAdjointStep_Theta;
922:   ts->ops->adjointintegral = TSAdjointCostIntegral_Theta;
923:   ts->ops->forwardintegral = TSForwardCostIntegral_Theta;
924:   ts->default_adapt_type   = TSADAPTNONE;
925:   ts->ops->forwardsetup    = TSForwardSetUp_Theta;
926:   ts->ops->forwardstep     = TSForwardStep_Theta;

928:   ts->usessnes = PETSC_TRUE;

930:   PetscNewLog(ts,&th);
931:   ts->data = (void*)th;

933:   th->VecsDeltaLam   = NULL;
934:   th->VecsDeltaMu    = NULL;
935:   th->VecsSensiTemp  = NULL;

937:   th->extrapolate = PETSC_FALSE;
938:   th->Theta       = 0.5;
939:   th->order       = 2;
940:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",TSThetaGetTheta_Theta);
941:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",TSThetaSetTheta_Theta);
942:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",TSThetaGetEndpoint_Theta);
943:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",TSThetaSetEndpoint_Theta);
944:   return(0);
945: }

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

950:   Not Collective

952:   Input Parameter:
953: .  ts - timestepping context

955:   Output Parameter:
956: .  theta - stage abscissa

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

961:   Level: Advanced

963: .seealso: TSThetaSetTheta()
964: @*/
965: PetscErrorCode  TSThetaGetTheta(TS ts,PetscReal *theta)
966: {

972:   PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
973:   return(0);
974: }

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

979:   Not Collective

981:   Input Parameter:
982: +  ts - timestepping context
983: -  theta - stage abscissa

985:   Options Database:
986: .  -ts_theta_theta <theta>

988:   Level: Intermediate

990: .seealso: TSThetaGetTheta()
991: @*/
992: PetscErrorCode  TSThetaSetTheta(TS ts,PetscReal theta)
993: {

998:   PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
999:   return(0);
1000: }

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

1005:   Not Collective

1007:   Input Parameter:
1008: .  ts - timestepping context

1010:   Output Parameter:
1011: .  endpoint - PETSC_TRUE when using the endpoint variant

1013:   Level: Advanced

1015: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
1016: @*/
1017: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
1018: {

1024:   PetscUseMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
1025:   return(0);
1026: }

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

1031:   Not Collective

1033:   Input Parameter:
1034: +  ts - timestepping context
1035: -  flg - PETSC_TRUE to use the endpoint variant

1037:   Options Database:
1038: .  -ts_theta_endpoint <flg>

1040:   Level: Intermediate

1042: .seealso: TSTHETA, TSCN
1043: @*/
1044: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
1045: {

1050:   PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
1051:   return(0);
1052: }

1054: /*
1055:  * TSBEULER and TSCN are straightforward specializations of TSTHETA.
1056:  * The creation functions for these specializations are below.
1057:  */

1059: static PetscErrorCode TSSetUp_BEuler(TS ts)
1060: {
1061:   TS_Theta       *th = (TS_Theta*)ts->data;

1065:   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");
1066:   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");
1067:   TSSetUp_Theta(ts);
1068:   return(0);
1069: }

1071: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
1072: {
1074:   return(0);
1075: }

1077: /*MC
1078:       TSBEULER - ODE solver using the implicit backward Euler method

1080:   Level: beginner

1082:   Notes:
1083:   TSBEULER is equivalent to TSTHETA with Theta=1.0

1085: $  -ts_type theta -ts_theta_theta 1.0

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

1089: M*/
1090: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
1091: {

1095:   TSCreate_Theta(ts);
1096:   TSThetaSetTheta(ts,1.0);
1097:   TSThetaSetEndpoint(ts,PETSC_FALSE);
1098:   ts->ops->setup = TSSetUp_BEuler;
1099:   ts->ops->view  = TSView_BEuler;
1100:   return(0);
1101: }

1103: static PetscErrorCode TSSetUp_CN(TS ts)
1104: {
1105:   TS_Theta       *th = (TS_Theta*)ts->data;

1109:   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");
1110:   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");
1111:   TSSetUp_Theta(ts);
1112:   return(0);
1113: }

1115: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
1116: {
1118:   return(0);
1119: }

1121: /*MC
1122:       TSCN - ODE solver using the implicit Crank-Nicolson method.

1124:   Level: beginner

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

1129: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint

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

1133: M*/
1134: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
1135: {

1139:   TSCreate_Theta(ts);
1140:   TSThetaSetTheta(ts,0.5);
1141:   TSThetaSetEndpoint(ts,PETSC_TRUE);
1142:   ts->ops->setup = TSSetUp_CN;
1143:   ts->ops->view  = TSView_CN;
1144:   return(0);
1145: }