Actual source code: theta.c

  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          Stages[2];   /* Storage for stage solutions */
 13:   Vec          X0, X, Xdot; /* Storage for u^n, u^n + dt a_{11} k_1, and time derivative u^{n+1}_t */
 14:   Vec          affine;      /* Affine vector needed for residual at beginning of step in endpoint formulation */
 15:   PetscReal    Theta;
 16:   PetscReal    shift; /* Shift parameter for SNES Jacobian, used by forward, TLM and adjoint */
 17:   PetscInt     order;
 18:   PetscBool    endpoint;
 19:   PetscBool    extrapolate;
 20:   TSStepStatus status;
 21:   Vec          VecCostIntegral0; /* Backup for roll-backs due to events, used by cost integral */
 22:   PetscReal    ptime0;           /* Backup for ts->ptime, the start time of current time step, used by TLM and cost integral */
 23:   PetscReal    time_step0;       /* Backup for ts->timestep, the step size of current time step, used by TLM and cost integral*/

 25:   /* context for sensitivity analysis */
 26:   PetscInt num_tlm;               /* Total number of tangent linear equations */
 27:   Vec     *VecsDeltaLam;          /* Increment of the adjoint sensitivity w.r.t IC at stage */
 28:   Vec     *VecsDeltaMu;           /* Increment of the adjoint sensitivity w.r.t P at stage */
 29:   Vec     *VecsSensiTemp;         /* Vector to be multiplied with Jacobian transpose */
 30:   Mat      MatFwdStages[2];       /* TLM Stages */
 31:   Mat      MatDeltaFwdSensip;     /* Increment of the forward sensitivity at stage */
 32:   Vec      VecDeltaFwdSensipCol;  /* Working vector for holding one column of the sensitivity matrix */
 33:   Mat      MatFwdSensip0;         /* backup for roll-backs due to events */
 34:   Mat      MatIntegralSensipTemp; /* Working vector for forward integral sensitivity */
 35:   Mat      MatIntegralSensip0;    /* backup for roll-backs due to events */
 36:   Vec     *VecsDeltaLam2;         /* Increment of the 2nd-order adjoint sensitivity w.r.t IC at stage */
 37:   Vec     *VecsDeltaMu2;          /* Increment of the 2nd-order adjoint sensitivity w.r.t P at stage */
 38:   Vec     *VecsSensi2Temp;        /* Working vectors that holds the residual for the second-order adjoint */
 39:   Vec     *VecsAffine;            /* Working vectors to store residuals */
 40:   /* context for error estimation */
 41:   Vec vec_sol_prev;
 42:   Vec vec_lte_work;
 43: } TS_Theta;

 45: static PetscErrorCode TSThetaGetX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
 46: {
 47:   TS_Theta *th = (TS_Theta *)ts->data;

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

 62: static PetscErrorCode TSThetaRestoreX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
 63: {
 64:   if (X0) {
 65:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSTheta_X0", X0);
 66:   }
 67:   if (Xdot) {
 68:     if (dm && dm != ts->dm) DMRestoreNamedGlobalVector(dm, "TSTheta_Xdot", Xdot);
 69:   }
 70:   return 0;
 71: }

 73: static PetscErrorCode DMCoarsenHook_TSTheta(DM fine, DM coarse, void *ctx)
 74: {
 75:   return 0;
 76: }

 78: static PetscErrorCode DMRestrictHook_TSTheta(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
 79: {
 80:   TS  ts = (TS)ctx;
 81:   Vec X0, Xdot, X0_c, Xdot_c;

 83:   TSThetaGetX0AndXdot(ts, fine, &X0, &Xdot);
 84:   TSThetaGetX0AndXdot(ts, coarse, &X0_c, &Xdot_c);
 85:   MatRestrict(restrct, X0, X0_c);
 86:   MatRestrict(restrct, Xdot, Xdot_c);
 87:   VecPointwiseMult(X0_c, rscale, X0_c);
 88:   VecPointwiseMult(Xdot_c, rscale, Xdot_c);
 89:   TSThetaRestoreX0AndXdot(ts, fine, &X0, &Xdot);
 90:   TSThetaRestoreX0AndXdot(ts, coarse, &X0_c, &Xdot_c);
 91:   return 0;
 92: }

 94: static PetscErrorCode DMSubDomainHook_TSTheta(DM dm, DM subdm, void *ctx)
 95: {
 96:   return 0;
 97: }

 99: static PetscErrorCode DMSubDomainRestrictHook_TSTheta(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
100: {
101:   TS  ts = (TS)ctx;
102:   Vec X0, Xdot, X0_sub, Xdot_sub;

104:   TSThetaGetX0AndXdot(ts, dm, &X0, &Xdot);
105:   TSThetaGetX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);

107:   VecScatterBegin(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);
108:   VecScatterEnd(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);

110:   VecScatterBegin(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);
111:   VecScatterEnd(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);

113:   TSThetaRestoreX0AndXdot(ts, dm, &X0, &Xdot);
114:   TSThetaRestoreX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);
115:   return 0;
116: }

118: static PetscErrorCode TSThetaEvaluateCostIntegral(TS ts)
119: {
120:   TS_Theta *th     = (TS_Theta *)ts->data;
121:   TS        quadts = ts->quadraturets;

123:   if (th->endpoint) {
124:     /* Evolve ts->vec_costintegral to compute integrals */
125:     if (th->Theta != 1.0) {
126:       TSComputeRHSFunction(quadts, th->ptime0, th->X0, ts->vec_costintegrand);
127:       VecAXPY(quadts->vec_sol, th->time_step0 * (1.0 - th->Theta), ts->vec_costintegrand);
128:     }
129:     TSComputeRHSFunction(quadts, ts->ptime, ts->vec_sol, ts->vec_costintegrand);
130:     VecAXPY(quadts->vec_sol, th->time_step0 * th->Theta, ts->vec_costintegrand);
131:   } else {
132:     TSComputeRHSFunction(quadts, th->stage_time, th->X, ts->vec_costintegrand);
133:     VecAXPY(quadts->vec_sol, th->time_step0, ts->vec_costintegrand);
134:   }
135:   return 0;
136: }

138: static PetscErrorCode TSForwardCostIntegral_Theta(TS ts)
139: {
140:   TS_Theta *th     = (TS_Theta *)ts->data;
141:   TS        quadts = ts->quadraturets;

143:   /* backup cost integral */
144:   VecCopy(quadts->vec_sol, th->VecCostIntegral0);
145:   TSThetaEvaluateCostIntegral(ts);
146:   return 0;
147: }

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

153:   /* Like TSForwardCostIntegral(), the adjoint cost integral evaluation relies on ptime0 and time_step0. */
154:   th->ptime0     = ts->ptime + ts->time_step;
155:   th->time_step0 = -ts->time_step;
156:   TSThetaEvaluateCostIntegral(ts);
157:   return 0;
158: }

160: static PetscErrorCode TSTheta_SNESSolve(TS ts, Vec b, Vec x)
161: {
162:   PetscInt nits, lits;

164:   SNESSolve(ts->snes, b, x);
165:   SNESGetIterationNumber(ts->snes, &nits);
166:   SNESGetLinearSolveIterations(ts->snes, &lits);
167:   ts->snes_its += nits;
168:   ts->ksp_its += lits;
169:   return 0;
170: }

172: static PetscErrorCode TSStep_Theta(TS ts)
173: {
174:   TS_Theta *th         = (TS_Theta *)ts->data;
175:   PetscInt  rejections = 0;
176:   PetscBool stageok, accept = PETSC_TRUE;
177:   PetscReal next_time_step = ts->time_step;

179:   if (!ts->steprollback) {
180:     if (th->vec_sol_prev) VecCopy(th->X0, th->vec_sol_prev);
181:     VecCopy(ts->vec_sol, th->X0);
182:   }

184:   th->status = TS_STEP_INCOMPLETE;
185:   while (!ts->reason && th->status != TS_STEP_COMPLETE) {
186:     th->shift      = 1 / (th->Theta * ts->time_step);
187:     th->stage_time = ts->ptime + (th->endpoint ? (PetscReal)1 : th->Theta) * ts->time_step;
188:     VecCopy(th->X0, th->X);
189:     if (th->extrapolate && !ts->steprestart) VecAXPY(th->X, 1 / th->shift, th->Xdot);
190:     if (th->endpoint) { /* This formulation assumes linear time-independent mass matrix */
191:       if (!th->affine) VecDuplicate(ts->vec_sol, &th->affine);
192:       VecZeroEntries(th->Xdot);
193:       TSComputeIFunction(ts, ts->ptime, th->X0, th->Xdot, th->affine, PETSC_FALSE);
194:       VecScale(th->affine, (th->Theta - 1) / th->Theta);
195:     } else if (th->affine) { /* Just in case th->endpoint is changed between calls to TSStep_Theta() */
196:       VecZeroEntries(th->affine);
197:     }
198:     TSPreStage(ts, th->stage_time);
199:     TSTheta_SNESSolve(ts, th->affine, th->X);
200:     TSPostStage(ts, th->stage_time, 0, &th->X);
201:     TSAdaptCheckStage(ts->adapt, ts, th->stage_time, th->X, &stageok);
202:     if (!stageok) goto reject_step;

204:     th->status = TS_STEP_PENDING;
205:     if (th->endpoint) {
206:       VecCopy(th->X, ts->vec_sol);
207:     } else {
208:       VecAXPBYPCZ(th->Xdot, -th->shift, th->shift, 0, th->X0, th->X);
209:       VecAXPY(ts->vec_sol, ts->time_step, th->Xdot);
210:     }
211:     TSAdaptChoose(ts->adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
212:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
213:     if (!accept) {
214:       VecCopy(th->X0, ts->vec_sol);
215:       ts->time_step = next_time_step;
216:       goto reject_step;
217:     }

219:     if (ts->forward_solve || ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */
220:       th->ptime0     = ts->ptime;
221:       th->time_step0 = ts->time_step;
222:     }
223:     ts->ptime += ts->time_step;
224:     ts->time_step = next_time_step;
225:     break;

227:   reject_step:
228:     ts->reject++;
229:     accept = PETSC_FALSE;
230:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
231:       ts->reason = TS_DIVERGED_STEP_REJECTED;
232:       PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections);
233:     }
234:   }
235:   return 0;
236: }

238: static PetscErrorCode TSAdjointStepBEuler_Private(TS ts)
239: {
240:   TS_Theta      *th           = (TS_Theta *)ts->data;
241:   TS             quadts       = ts->quadraturets;
242:   Vec           *VecsDeltaLam = th->VecsDeltaLam, *VecsDeltaMu = th->VecsDeltaMu, *VecsSensiTemp = th->VecsSensiTemp;
243:   Vec           *VecsDeltaLam2 = th->VecsDeltaLam2, *VecsDeltaMu2 = th->VecsDeltaMu2, *VecsSensi2Temp = th->VecsSensi2Temp;
244:   PetscInt       nadj;
245:   Mat            J, Jpre, quadJ = NULL, quadJp = NULL;
246:   KSP            ksp;
247:   PetscScalar   *xarr;
248:   TSEquationType eqtype;
249:   PetscBool      isexplicitode = PETSC_FALSE;
250:   PetscReal      adjoint_time_step;

252:   TSGetEquationType(ts, &eqtype);
253:   if (eqtype == TS_EQ_ODE_EXPLICIT) {
254:     isexplicitode = PETSC_TRUE;
255:     VecsDeltaLam  = ts->vecs_sensi;
256:     VecsDeltaLam2 = ts->vecs_sensi2;
257:   }
258:   th->status = TS_STEP_INCOMPLETE;
259:   SNESGetKSP(ts->snes, &ksp);
260:   TSGetIJacobian(ts, &J, &Jpre, NULL, NULL);
261:   if (quadts) {
262:     TSGetRHSJacobian(quadts, &quadJ, NULL, NULL, NULL);
263:     TSGetRHSJacobianP(quadts, &quadJp, NULL, NULL);
264:   }

266:   th->stage_time    = ts->ptime;
267:   adjoint_time_step = -ts->time_step; /* always positive since time_step is negative */

269:   /* Build RHS for first-order adjoint lambda_{n+1}/h + r_u^T(n+1) */
270:   if (quadts) TSComputeRHSJacobian(quadts, th->stage_time, ts->vec_sol, quadJ, NULL);

272:   for (nadj = 0; nadj < ts->numcost; nadj++) {
273:     VecCopy(ts->vecs_sensi[nadj], VecsSensiTemp[nadj]);
274:     VecScale(VecsSensiTemp[nadj], 1. / adjoint_time_step); /* lambda_{n+1}/h */
275:     if (quadJ) {
276:       MatDenseGetColumn(quadJ, nadj, &xarr);
277:       VecPlaceArray(ts->vec_drdu_col, xarr);
278:       VecAXPY(VecsSensiTemp[nadj], 1., ts->vec_drdu_col);
279:       VecResetArray(ts->vec_drdu_col);
280:       MatDenseRestoreColumn(quadJ, &xarr);
281:     }
282:   }

284:   /* Build LHS for first-order adjoint */
285:   th->shift = 1. / adjoint_time_step;
286:   TSComputeSNESJacobian(ts, ts->vec_sol, J, Jpre);
287:   KSPSetOperators(ksp, J, Jpre);

289:   /* Solve stage equation LHS*lambda_s = RHS for first-order adjoint */
290:   for (nadj = 0; nadj < ts->numcost; nadj++) {
291:     KSPConvergedReason kspreason;
292:     KSPSolveTranspose(ksp, VecsSensiTemp[nadj], VecsDeltaLam[nadj]);
293:     KSPGetConvergedReason(ksp, &kspreason);
294:     if (kspreason < 0) {
295:       ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
296:       PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj);
297:     }
298:   }

300:   if (ts->vecs_sensi2) { /* U_{n+1} */
301:     /* Get w1 at t_{n+1} from TLM matrix */
302:     MatDenseGetColumn(ts->mat_sensip, 0, &xarr);
303:     VecPlaceArray(ts->vec_sensip_col, xarr);
304:     /* lambda_s^T F_UU w_1 */
305:     TSComputeIHessianProductFunctionUU(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_sensip_col, ts->vecs_fuu);
306:     /* lambda_s^T F_UP w_2 */
307:     TSComputeIHessianProductFunctionUP(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_dir, ts->vecs_fup);
308:     for (nadj = 0; nadj < ts->numcost; nadj++) { /* compute the residual */
309:       VecCopy(ts->vecs_sensi2[nadj], VecsSensi2Temp[nadj]);
310:       VecScale(VecsSensi2Temp[nadj], 1. / adjoint_time_step);
311:       VecAXPY(VecsSensi2Temp[nadj], -1., ts->vecs_fuu[nadj]);
312:       if (ts->vecs_fup) VecAXPY(VecsSensi2Temp[nadj], -1., ts->vecs_fup[nadj]);
313:     }
314:     /* Solve stage equation LHS X = RHS for second-order adjoint */
315:     for (nadj = 0; nadj < ts->numcost; nadj++) {
316:       KSPConvergedReason kspreason;
317:       KSPSolveTranspose(ksp, VecsSensi2Temp[nadj], VecsDeltaLam2[nadj]);
318:       KSPGetConvergedReason(ksp, &kspreason);
319:       if (kspreason < 0) {
320:         ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
321:         PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 2nd-order adjoint solve\n", ts->steps, nadj);
322:       }
323:     }
324:   }

326:   /* Update sensitivities, and evaluate integrals if there is any */
327:   if (!isexplicitode) {
328:     th->shift = 0.0;
329:     TSComputeSNESJacobian(ts, ts->vec_sol, J, Jpre);
330:     KSPSetOperators(ksp, J, Jpre);
331:     for (nadj = 0; nadj < ts->numcost; nadj++) {
332:       /* Add f_U \lambda_s to the original RHS */
333:       VecScale(VecsSensiTemp[nadj], -1.);
334:       MatMultTransposeAdd(J, VecsDeltaLam[nadj], VecsSensiTemp[nadj], VecsSensiTemp[nadj]);
335:       VecScale(VecsSensiTemp[nadj], -adjoint_time_step);
336:       VecCopy(VecsSensiTemp[nadj], ts->vecs_sensi[nadj]);
337:       if (ts->vecs_sensi2) {
338:         MatMultTransposeAdd(J, VecsDeltaLam2[nadj], VecsSensi2Temp[nadj], VecsSensi2Temp[nadj]);
339:         VecScale(VecsSensi2Temp[nadj], -adjoint_time_step);
340:         VecCopy(VecsSensi2Temp[nadj], ts->vecs_sensi2[nadj]);
341:       }
342:     }
343:   }
344:   if (ts->vecs_sensip) {
345:     VecAXPBYPCZ(th->Xdot, -1. / adjoint_time_step, 1.0 / adjoint_time_step, 0, th->X0, ts->vec_sol);
346:     TSComputeIJacobianP(ts, th->stage_time, ts->vec_sol, th->Xdot, 1. / adjoint_time_step, ts->Jacp, PETSC_FALSE); /* get -f_p */
347:     if (quadts) TSComputeRHSJacobianP(quadts, th->stage_time, ts->vec_sol, quadJp);
348:     if (ts->vecs_sensi2p) {
349:       /* lambda_s^T F_PU w_1 */
350:       TSComputeIHessianProductFunctionPU(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_sensip_col, ts->vecs_fpu);
351:       /* lambda_s^T F_PP w_2 */
352:       TSComputeIHessianProductFunctionPP(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_dir, ts->vecs_fpp);
353:     }

355:     for (nadj = 0; nadj < ts->numcost; nadj++) {
356:       MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj], VecsDeltaMu[nadj]);
357:       VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step, VecsDeltaMu[nadj]);
358:       if (quadJp) {
359:         MatDenseGetColumn(quadJp, nadj, &xarr);
360:         VecPlaceArray(ts->vec_drdp_col, xarr);
361:         VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, ts->vec_drdp_col);
362:         VecResetArray(ts->vec_drdp_col);
363:         MatDenseRestoreColumn(quadJp, &xarr);
364:       }
365:       if (ts->vecs_sensi2p) {
366:         MatMultTranspose(ts->Jacp, VecsDeltaLam2[nadj], VecsDeltaMu2[nadj]);
367:         VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step, VecsDeltaMu2[nadj]);
368:         if (ts->vecs_fpu) VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step, ts->vecs_fpu[nadj]);
369:         if (ts->vecs_fpp) VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step, ts->vecs_fpp[nadj]);
370:       }
371:     }
372:   }

374:   if (ts->vecs_sensi2) {
375:     VecResetArray(ts->vec_sensip_col);
376:     MatDenseRestoreColumn(ts->mat_sensip, &xarr);
377:   }
378:   th->status = TS_STEP_COMPLETE;
379:   return 0;
380: }

382: static PetscErrorCode TSAdjointStep_Theta(TS ts)
383: {
384:   TS_Theta    *th           = (TS_Theta *)ts->data;
385:   TS           quadts       = ts->quadraturets;
386:   Vec         *VecsDeltaLam = th->VecsDeltaLam, *VecsDeltaMu = th->VecsDeltaMu, *VecsSensiTemp = th->VecsSensiTemp;
387:   Vec         *VecsDeltaLam2 = th->VecsDeltaLam2, *VecsDeltaMu2 = th->VecsDeltaMu2, *VecsSensi2Temp = th->VecsSensi2Temp;
388:   PetscInt     nadj;
389:   Mat          J, Jpre, quadJ = NULL, quadJp = NULL;
390:   KSP          ksp;
391:   PetscScalar *xarr;
392:   PetscReal    adjoint_time_step;
393:   PetscReal    adjoint_ptime; /* end time of the adjoint time step (ts->ptime is the start time, ususally ts->ptime is larger than adjoint_ptime) */

395:   if (th->Theta == 1.) {
396:     TSAdjointStepBEuler_Private(ts);
397:     return 0;
398:   }
399:   th->status = TS_STEP_INCOMPLETE;
400:   SNESGetKSP(ts->snes, &ksp);
401:   TSGetIJacobian(ts, &J, &Jpre, NULL, NULL);
402:   if (quadts) {
403:     TSGetRHSJacobian(quadts, &quadJ, NULL, NULL, NULL);
404:     TSGetRHSJacobianP(quadts, &quadJp, NULL, NULL);
405:   }
406:   /* If endpoint=1, th->ptime and th->X0 will be used; if endpoint=0, th->stage_time and th->X will be used. */
407:   th->stage_time    = th->endpoint ? ts->ptime : (ts->ptime + (1. - th->Theta) * ts->time_step);
408:   adjoint_ptime     = ts->ptime + ts->time_step;
409:   adjoint_time_step = -ts->time_step; /* always positive since time_step is negative */

411:   if (!th->endpoint) {
412:     /* recover th->X0 using vec_sol and the stage value th->X */
413:     VecAXPBYPCZ(th->X0, 1.0 / (1.0 - th->Theta), th->Theta / (th->Theta - 1.0), 0, th->X, ts->vec_sol);
414:   }

416:   /* Build RHS for first-order adjoint */
417:   /* Cost function has an integral term */
418:   if (quadts) {
419:     if (th->endpoint) {
420:       TSComputeRHSJacobian(quadts, th->stage_time, ts->vec_sol, quadJ, NULL);
421:     } else {
422:       TSComputeRHSJacobian(quadts, th->stage_time, th->X, quadJ, NULL);
423:     }
424:   }

426:   for (nadj = 0; nadj < ts->numcost; nadj++) {
427:     VecCopy(ts->vecs_sensi[nadj], VecsSensiTemp[nadj]);
428:     VecScale(VecsSensiTemp[nadj], 1. / (th->Theta * adjoint_time_step));
429:     if (quadJ) {
430:       MatDenseGetColumn(quadJ, nadj, &xarr);
431:       VecPlaceArray(ts->vec_drdu_col, xarr);
432:       VecAXPY(VecsSensiTemp[nadj], 1., ts->vec_drdu_col);
433:       VecResetArray(ts->vec_drdu_col);
434:       MatDenseRestoreColumn(quadJ, &xarr);
435:     }
436:   }

438:   /* Build LHS for first-order adjoint */
439:   th->shift = 1. / (th->Theta * adjoint_time_step);
440:   if (th->endpoint) {
441:     TSComputeSNESJacobian(ts, ts->vec_sol, J, Jpre);
442:   } else {
443:     TSComputeSNESJacobian(ts, th->X, J, Jpre);
444:   }
445:   KSPSetOperators(ksp, J, Jpre);

447:   /* Solve stage equation LHS*lambda_s = RHS for first-order adjoint */
448:   for (nadj = 0; nadj < ts->numcost; nadj++) {
449:     KSPConvergedReason kspreason;
450:     KSPSolveTranspose(ksp, VecsSensiTemp[nadj], VecsDeltaLam[nadj]);
451:     KSPGetConvergedReason(ksp, &kspreason);
452:     if (kspreason < 0) {
453:       ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
454:       PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n", ts->steps, nadj);
455:     }
456:   }

458:   /* Second-order adjoint */
459:   if (ts->vecs_sensi2) { /* U_{n+1} */
461:     /* Get w1 at t_{n+1} from TLM matrix */
462:     MatDenseGetColumn(ts->mat_sensip, 0, &xarr);
463:     VecPlaceArray(ts->vec_sensip_col, xarr);
464:     /* lambda_s^T F_UU w_1 */
465:     TSComputeIHessianProductFunctionUU(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_sensip_col, ts->vecs_fuu);
466:     VecResetArray(ts->vec_sensip_col);
467:     MatDenseRestoreColumn(ts->mat_sensip, &xarr);
468:     /* lambda_s^T F_UP w_2 */
469:     TSComputeIHessianProductFunctionUP(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_dir, ts->vecs_fup);
470:     for (nadj = 0; nadj < ts->numcost; nadj++) { /* compute the residual */
471:       VecCopy(ts->vecs_sensi2[nadj], VecsSensi2Temp[nadj]);
472:       VecScale(VecsSensi2Temp[nadj], th->shift);
473:       VecAXPY(VecsSensi2Temp[nadj], -1., ts->vecs_fuu[nadj]);
474:       if (ts->vecs_fup) VecAXPY(VecsSensi2Temp[nadj], -1., ts->vecs_fup[nadj]);
475:     }
476:     /* Solve stage equation LHS X = RHS for second-order adjoint */
477:     for (nadj = 0; nadj < ts->numcost; nadj++) {
478:       KSPConvergedReason kspreason;
479:       KSPSolveTranspose(ksp, VecsSensi2Temp[nadj], VecsDeltaLam2[nadj]);
480:       KSPGetConvergedReason(ksp, &kspreason);
481:       if (kspreason < 0) {
482:         ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
483:         PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th cost function, transposed linear solve fails, stopping 2nd-order adjoint solve\n", ts->steps, nadj);
484:       }
485:     }
486:   }

488:   /* Update sensitivities, and evaluate integrals if there is any */
489:   if (th->endpoint) { /* two-stage Theta methods with th->Theta!=1, th->Theta==1 leads to BEuler */
490:     th->shift      = 1. / ((th->Theta - 1.) * adjoint_time_step);
491:     th->stage_time = adjoint_ptime;
492:     TSComputeSNESJacobian(ts, th->X0, J, Jpre);
493:     KSPSetOperators(ksp, J, Jpre);
494:     /* R_U at t_n */
495:     if (quadts) TSComputeRHSJacobian(quadts, adjoint_ptime, th->X0, quadJ, NULL);
496:     for (nadj = 0; nadj < ts->numcost; nadj++) {
497:       MatMultTranspose(J, VecsDeltaLam[nadj], ts->vecs_sensi[nadj]);
498:       if (quadJ) {
499:         MatDenseGetColumn(quadJ, nadj, &xarr);
500:         VecPlaceArray(ts->vec_drdu_col, xarr);
501:         VecAXPY(ts->vecs_sensi[nadj], -1., ts->vec_drdu_col);
502:         VecResetArray(ts->vec_drdu_col);
503:         MatDenseRestoreColumn(quadJ, &xarr);
504:       }
505:       VecScale(ts->vecs_sensi[nadj], 1. / th->shift);
506:     }

508:     /* Second-order adjoint */
509:     if (ts->vecs_sensi2) { /* U_n */
510:       /* Get w1 at t_n from TLM matrix */
511:       MatDenseGetColumn(th->MatFwdSensip0, 0, &xarr);
512:       VecPlaceArray(ts->vec_sensip_col, xarr);
513:       /* lambda_s^T F_UU w_1 */
514:       TSComputeIHessianProductFunctionUU(ts, adjoint_ptime, th->X0, VecsDeltaLam, ts->vec_sensip_col, ts->vecs_fuu);
515:       VecResetArray(ts->vec_sensip_col);
516:       MatDenseRestoreColumn(th->MatFwdSensip0, &xarr);
517:       /* lambda_s^T F_UU w_2 */
518:       TSComputeIHessianProductFunctionUP(ts, adjoint_ptime, th->X0, VecsDeltaLam, ts->vec_dir, ts->vecs_fup);
519:       for (nadj = 0; nadj < ts->numcost; nadj++) {
520:         /* M^T Lambda_s + h(1-theta) F_U^T Lambda_s + h(1-theta) lambda_s^T F_UU w_1 + lambda_s^T F_UP w_2  */
521:         MatMultTranspose(J, VecsDeltaLam2[nadj], ts->vecs_sensi2[nadj]);
522:         VecAXPY(ts->vecs_sensi2[nadj], 1., ts->vecs_fuu[nadj]);
523:         if (ts->vecs_fup) VecAXPY(ts->vecs_sensi2[nadj], 1., ts->vecs_fup[nadj]);
524:         VecScale(ts->vecs_sensi2[nadj], 1. / th->shift);
525:       }
526:     }

528:     th->stage_time = ts->ptime; /* recover the old value */

530:     if (ts->vecs_sensip) { /* sensitivities wrt parameters */
531:       /* U_{n+1} */
532:       th->shift = 1.0 / (adjoint_time_step * th->Theta);
533:       VecAXPBYPCZ(th->Xdot, -th->shift, th->shift, 0, th->X0, ts->vec_sol);
534:       TSComputeIJacobianP(ts, th->stage_time, ts->vec_sol, th->Xdot, -1. / (th->Theta * adjoint_time_step), ts->Jacp, PETSC_FALSE);
535:       if (quadts) TSComputeRHSJacobianP(quadts, th->stage_time, ts->vec_sol, quadJp);
536:       for (nadj = 0; nadj < ts->numcost; nadj++) {
537:         MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj], VecsDeltaMu[nadj]);
538:         VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step * th->Theta, VecsDeltaMu[nadj]);
539:         if (quadJp) {
540:           MatDenseGetColumn(quadJp, nadj, &xarr);
541:           VecPlaceArray(ts->vec_drdp_col, xarr);
542:           VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step * th->Theta, ts->vec_drdp_col);
543:           VecResetArray(ts->vec_drdp_col);
544:           MatDenseRestoreColumn(quadJp, &xarr);
545:         }
546:       }
547:       if (ts->vecs_sensi2p) { /* second-order */
548:         /* Get w1 at t_{n+1} from TLM matrix */
549:         MatDenseGetColumn(ts->mat_sensip, 0, &xarr);
550:         VecPlaceArray(ts->vec_sensip_col, xarr);
551:         /* lambda_s^T F_PU w_1 */
552:         TSComputeIHessianProductFunctionPU(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_sensip_col, ts->vecs_fpu);
553:         VecResetArray(ts->vec_sensip_col);
554:         MatDenseRestoreColumn(ts->mat_sensip, &xarr);

556:         /* lambda_s^T F_PP w_2 */
557:         TSComputeIHessianProductFunctionPP(ts, th->stage_time, ts->vec_sol, VecsDeltaLam, ts->vec_dir, ts->vecs_fpp);
558:         for (nadj = 0; nadj < ts->numcost; nadj++) {
559:           /* Mu2 <- Mu2 + h theta F_P^T Lambda_s + h theta (lambda_s^T F_UU w_1 + lambda_s^T F_UP w_2)  */
560:           MatMultTranspose(ts->Jacp, VecsDeltaLam2[nadj], VecsDeltaMu2[nadj]);
561:           VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * th->Theta, VecsDeltaMu2[nadj]);
562:           if (ts->vecs_fpu) VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * th->Theta, ts->vecs_fpu[nadj]);
563:           if (ts->vecs_fpp) VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * th->Theta, ts->vecs_fpp[nadj]);
564:         }
565:       }

567:       /* U_s */
568:       VecZeroEntries(th->Xdot);
569:       TSComputeIJacobianP(ts, adjoint_ptime, th->X0, th->Xdot, 1. / ((th->Theta - 1.0) * adjoint_time_step), ts->Jacp, PETSC_FALSE);
570:       if (quadts) TSComputeRHSJacobianP(quadts, adjoint_ptime, th->X0, quadJp);
571:       for (nadj = 0; nadj < ts->numcost; nadj++) {
572:         MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj], VecsDeltaMu[nadj]);
573:         VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step * (1.0 - th->Theta), VecsDeltaMu[nadj]);
574:         if (quadJp) {
575:           MatDenseGetColumn(quadJp, nadj, &xarr);
576:           VecPlaceArray(ts->vec_drdp_col, xarr);
577:           VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step * (1.0 - th->Theta), ts->vec_drdp_col);
578:           VecResetArray(ts->vec_drdp_col);
579:           MatDenseRestoreColumn(quadJp, &xarr);
580:         }
581:         if (ts->vecs_sensi2p) { /* second-order */
582:           /* Get w1 at t_n from TLM matrix */
583:           MatDenseGetColumn(th->MatFwdSensip0, 0, &xarr);
584:           VecPlaceArray(ts->vec_sensip_col, xarr);
585:           /* lambda_s^T F_PU w_1 */
586:           TSComputeIHessianProductFunctionPU(ts, adjoint_ptime, th->X0, VecsDeltaLam, ts->vec_sensip_col, ts->vecs_fpu);
587:           VecResetArray(ts->vec_sensip_col);
588:           MatDenseRestoreColumn(th->MatFwdSensip0, &xarr);
589:           /* lambda_s^T F_PP w_2 */
590:           TSComputeIHessianProductFunctionPP(ts, adjoint_ptime, th->X0, VecsDeltaLam, ts->vec_dir, ts->vecs_fpp);
591:           for (nadj = 0; nadj < ts->numcost; nadj++) {
592:             /* Mu2 <- Mu2 + h(1-theta) F_P^T Lambda_s + h(1-theta) (lambda_s^T F_UU w_1 + lambda_s^T F_UP w_2) */
593:             MatMultTranspose(ts->Jacp, VecsDeltaLam2[nadj], VecsDeltaMu2[nadj]);
594:             VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * (1.0 - th->Theta), VecsDeltaMu2[nadj]);
595:             if (ts->vecs_fpu) VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * (1.0 - th->Theta), ts->vecs_fpu[nadj]);
596:             if (ts->vecs_fpp) VecAXPY(ts->vecs_sensi2p[nadj], -adjoint_time_step * (1.0 - th->Theta), ts->vecs_fpp[nadj]);
597:           }
598:         }
599:       }
600:     }
601:   } else { /* one-stage case */
602:     th->shift = 0.0;
603:     TSComputeSNESJacobian(ts, th->X, J, Jpre); /* get -f_y */
604:     KSPSetOperators(ksp, J, Jpre);
605:     if (quadts) TSComputeRHSJacobian(quadts, th->stage_time, th->X, quadJ, NULL);
606:     for (nadj = 0; nadj < ts->numcost; nadj++) {
607:       MatMultTranspose(J, VecsDeltaLam[nadj], VecsSensiTemp[nadj]);
608:       VecAXPY(ts->vecs_sensi[nadj], -adjoint_time_step, VecsSensiTemp[nadj]);
609:       if (quadJ) {
610:         MatDenseGetColumn(quadJ, nadj, &xarr);
611:         VecPlaceArray(ts->vec_drdu_col, xarr);
612:         VecAXPY(ts->vecs_sensi[nadj], adjoint_time_step, ts->vec_drdu_col);
613:         VecResetArray(ts->vec_drdu_col);
614:         MatDenseRestoreColumn(quadJ, &xarr);
615:       }
616:     }
617:     if (ts->vecs_sensip) {
618:       th->shift = 1.0 / (adjoint_time_step * th->Theta);
619:       VecAXPBYPCZ(th->Xdot, -th->shift, th->shift, 0, th->X0, th->X);
620:       TSComputeIJacobianP(ts, th->stage_time, th->X, th->Xdot, th->shift, ts->Jacp, PETSC_FALSE);
621:       if (quadts) TSComputeRHSJacobianP(quadts, th->stage_time, th->X, quadJp);
622:       for (nadj = 0; nadj < ts->numcost; nadj++) {
623:         MatMultTranspose(ts->Jacp, VecsDeltaLam[nadj], VecsDeltaMu[nadj]);
624:         VecAXPY(ts->vecs_sensip[nadj], -adjoint_time_step, VecsDeltaMu[nadj]);
625:         if (quadJp) {
626:           MatDenseGetColumn(quadJp, nadj, &xarr);
627:           VecPlaceArray(ts->vec_drdp_col, xarr);
628:           VecAXPY(ts->vecs_sensip[nadj], adjoint_time_step, ts->vec_drdp_col);
629:           VecResetArray(ts->vec_drdp_col);
630:           MatDenseRestoreColumn(quadJp, &xarr);
631:         }
632:       }
633:     }
634:   }

636:   th->status = TS_STEP_COMPLETE;
637:   return 0;
638: }

640: static PetscErrorCode TSInterpolate_Theta(TS ts, PetscReal t, Vec X)
641: {
642:   TS_Theta *th = (TS_Theta *)ts->data;
643:   PetscReal dt = t - ts->ptime;

645:   VecCopy(ts->vec_sol, th->X);
646:   if (th->endpoint) dt *= th->Theta;
647:   VecWAXPY(X, dt, th->Xdot, th->X);
648:   return 0;
649: }

651: static PetscErrorCode TSEvaluateWLTE_Theta(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
652: {
653:   TS_Theta *th = (TS_Theta *)ts->data;
654:   Vec       X  = ts->vec_sol;      /* X = solution */
655:   Vec       Y  = th->vec_lte_work; /* Y = X + LTE  */
656:   PetscReal wltea, wlter;

658:   if (!th->vec_sol_prev) {
659:     *wlte = -1;
660:     return 0;
661:   }
662:   /* Cannot compute LTE in first step or in restart after event */
663:   if (ts->steprestart) {
664:     *wlte = -1;
665:     return 0;
666:   }
667:   /* Compute LTE using backward differences with non-constant time step */
668:   {
669:     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
670:     PetscReal   a = 1 + h_prev / h;
671:     PetscScalar scal[3];
672:     Vec         vecs[3];
673:     scal[0] = +1 / a;
674:     scal[1] = -1 / (a - 1);
675:     scal[2] = +1 / (a * (a - 1));
676:     vecs[0] = X;
677:     vecs[1] = th->X0;
678:     vecs[2] = th->vec_sol_prev;
679:     VecCopy(X, Y);
680:     VecMAXPY(Y, 3, scal, vecs);
681:     TSErrorWeightedNorm(ts, X, Y, wnormtype, wlte, &wltea, &wlter);
682:   }
683:   if (order) *order = 2;
684:   return 0;
685: }

687: static PetscErrorCode TSRollBack_Theta(TS ts)
688: {
689:   TS_Theta *th     = (TS_Theta *)ts->data;
690:   TS        quadts = ts->quadraturets;

692:   VecCopy(th->X0, ts->vec_sol);
693:   if (quadts && ts->costintegralfwd) VecCopy(th->VecCostIntegral0, quadts->vec_sol);
694:   th->status = TS_STEP_INCOMPLETE;
695:   if (ts->mat_sensip) MatCopy(th->MatFwdSensip0, ts->mat_sensip, SAME_NONZERO_PATTERN);
696:   if (quadts && quadts->mat_sensip) MatCopy(th->MatIntegralSensip0, quadts->mat_sensip, SAME_NONZERO_PATTERN);
697:   return 0;
698: }

700: static PetscErrorCode TSForwardStep_Theta(TS ts)
701: {
702:   TS_Theta    *th                   = (TS_Theta *)ts->data;
703:   TS           quadts               = ts->quadraturets;
704:   Mat          MatDeltaFwdSensip    = th->MatDeltaFwdSensip;
705:   Vec          VecDeltaFwdSensipCol = th->VecDeltaFwdSensipCol;
706:   PetscInt     ntlm;
707:   KSP          ksp;
708:   Mat          J, Jpre, quadJ = NULL, quadJp = NULL;
709:   PetscScalar *barr, *xarr;
710:   PetscReal    previous_shift;

712:   previous_shift = th->shift;
713:   MatCopy(ts->mat_sensip, th->MatFwdSensip0, SAME_NONZERO_PATTERN);

715:   if (quadts && quadts->mat_sensip) MatCopy(quadts->mat_sensip, th->MatIntegralSensip0, SAME_NONZERO_PATTERN);
716:   SNESGetKSP(ts->snes, &ksp);
717:   TSGetIJacobian(ts, &J, &Jpre, NULL, NULL);
718:   if (quadts) {
719:     TSGetRHSJacobian(quadts, &quadJ, NULL, NULL, NULL);
720:     TSGetRHSJacobianP(quadts, &quadJp, NULL, NULL);
721:   }

723:   /* Build RHS */
724:   if (th->endpoint) { /* 2-stage method*/
725:     th->shift = 1. / ((th->Theta - 1.) * th->time_step0);
726:     TSComputeIJacobian(ts, th->ptime0, th->X0, th->Xdot, th->shift, J, Jpre, PETSC_FALSE);
727:     MatMatMult(J, ts->mat_sensip, MAT_REUSE_MATRIX, PETSC_DEFAULT, &MatDeltaFwdSensip);
728:     MatScale(MatDeltaFwdSensip, (th->Theta - 1.) / th->Theta);

730:     /* Add the f_p forcing terms */
731:     if (ts->Jacp) {
732:       VecZeroEntries(th->Xdot);
733:       TSComputeIJacobianP(ts, th->ptime0, th->X0, th->Xdot, th->shift, ts->Jacp, PETSC_FALSE);
734:       MatAXPY(MatDeltaFwdSensip, (th->Theta - 1.) / th->Theta, ts->Jacp, SUBSET_NONZERO_PATTERN);
735:       th->shift = previous_shift;
736:       VecAXPBYPCZ(th->Xdot, -th->shift, th->shift, 0, th->X0, ts->vec_sol);
737:       TSComputeIJacobianP(ts, th->stage_time, ts->vec_sol, th->Xdot, th->shift, ts->Jacp, PETSC_FALSE);
738:       MatAXPY(MatDeltaFwdSensip, -1., ts->Jacp, SUBSET_NONZERO_PATTERN);
739:     }
740:   } else { /* 1-stage method */
741:     th->shift = 0.0;
742:     TSComputeIJacobian(ts, th->stage_time, th->X, th->Xdot, th->shift, J, Jpre, PETSC_FALSE);
743:     MatMatMult(J, ts->mat_sensip, MAT_REUSE_MATRIX, PETSC_DEFAULT, &MatDeltaFwdSensip);
744:     MatScale(MatDeltaFwdSensip, -1.);

746:     /* Add the f_p forcing terms */
747:     if (ts->Jacp) {
748:       th->shift = previous_shift;
749:       VecAXPBYPCZ(th->Xdot, -th->shift, th->shift, 0, th->X0, th->X);
750:       TSComputeIJacobianP(ts, th->stage_time, th->X, th->Xdot, th->shift, ts->Jacp, PETSC_FALSE);
751:       MatAXPY(MatDeltaFwdSensip, -1., ts->Jacp, SUBSET_NONZERO_PATTERN);
752:     }
753:   }

755:   /* Build LHS */
756:   th->shift = previous_shift; /* recover the previous shift used in TSStep_Theta() */
757:   if (th->endpoint) {
758:     TSComputeIJacobian(ts, th->stage_time, ts->vec_sol, th->Xdot, th->shift, J, Jpre, PETSC_FALSE);
759:   } else {
760:     TSComputeIJacobian(ts, th->stage_time, th->X, th->Xdot, th->shift, J, Jpre, PETSC_FALSE);
761:   }
762:   KSPSetOperators(ksp, J, Jpre);

764:   /*
765:     Evaluate the first stage of integral gradients with the 2-stage method:
766:     drdu|t_n*S(t_n) + drdp|t_n
767:     This is done before the linear solve because the sensitivity variable S(t_n) will be propagated to S(t_{n+1})
768:   */
769:   if (th->endpoint) { /* 2-stage method only */
770:     if (quadts && quadts->mat_sensip) {
771:       TSComputeRHSJacobian(quadts, th->ptime0, th->X0, quadJ, NULL);
772:       TSComputeRHSJacobianP(quadts, th->ptime0, th->X0, quadJp);
773:       MatTransposeMatMult(ts->mat_sensip, quadJ, MAT_REUSE_MATRIX, PETSC_DEFAULT, &th->MatIntegralSensipTemp);
774:       MatAXPY(th->MatIntegralSensipTemp, 1, quadJp, SAME_NONZERO_PATTERN);
775:       MatAXPY(quadts->mat_sensip, th->time_step0 * (1. - th->Theta), th->MatIntegralSensipTemp, SAME_NONZERO_PATTERN);
776:     }
777:   }

779:   /* Solve the tangent linear equation for forward sensitivities to parameters */
780:   for (ntlm = 0; ntlm < th->num_tlm; ntlm++) {
781:     KSPConvergedReason kspreason;
782:     MatDenseGetColumn(MatDeltaFwdSensip, ntlm, &barr);
783:     VecPlaceArray(VecDeltaFwdSensipCol, barr);
784:     if (th->endpoint) {
785:       MatDenseGetColumn(ts->mat_sensip, ntlm, &xarr);
786:       VecPlaceArray(ts->vec_sensip_col, xarr);
787:       KSPSolve(ksp, VecDeltaFwdSensipCol, ts->vec_sensip_col);
788:       VecResetArray(ts->vec_sensip_col);
789:       MatDenseRestoreColumn(ts->mat_sensip, &xarr);
790:     } else {
791:       KSPSolve(ksp, VecDeltaFwdSensipCol, VecDeltaFwdSensipCol);
792:     }
793:     KSPGetConvergedReason(ksp, &kspreason);
794:     if (kspreason < 0) {
795:       ts->reason = TSFORWARD_DIVERGED_LINEAR_SOLVE;
796:       PetscInfo(ts, "Step=%" PetscInt_FMT ", %" PetscInt_FMT "th tangent linear solve, linear solve fails, stopping tangent linear solve\n", ts->steps, ntlm);
797:     }
798:     VecResetArray(VecDeltaFwdSensipCol);
799:     MatDenseRestoreColumn(MatDeltaFwdSensip, &barr);
800:   }

802:   /*
803:     Evaluate the second stage of integral gradients with the 2-stage method:
804:     drdu|t_{n+1}*S(t_{n+1}) + drdp|t_{n+1}
805:   */
806:   if (quadts && quadts->mat_sensip) {
807:     if (!th->endpoint) {
808:       MatAXPY(ts->mat_sensip, 1, MatDeltaFwdSensip, SAME_NONZERO_PATTERN); /* stage sensitivity */
809:       TSComputeRHSJacobian(quadts, th->stage_time, th->X, quadJ, NULL);
810:       TSComputeRHSJacobianP(quadts, th->stage_time, th->X, quadJp);
811:       MatTransposeMatMult(ts->mat_sensip, quadJ, MAT_REUSE_MATRIX, PETSC_DEFAULT, &th->MatIntegralSensipTemp);
812:       MatAXPY(th->MatIntegralSensipTemp, 1, quadJp, SAME_NONZERO_PATTERN);
813:       MatAXPY(quadts->mat_sensip, th->time_step0, th->MatIntegralSensipTemp, SAME_NONZERO_PATTERN);
814:       MatAXPY(ts->mat_sensip, (1. - th->Theta) / th->Theta, MatDeltaFwdSensip, SAME_NONZERO_PATTERN);
815:     } else {
816:       TSComputeRHSJacobian(quadts, th->stage_time, ts->vec_sol, quadJ, NULL);
817:       TSComputeRHSJacobianP(quadts, th->stage_time, ts->vec_sol, quadJp);
818:       MatTransposeMatMult(ts->mat_sensip, quadJ, MAT_REUSE_MATRIX, PETSC_DEFAULT, &th->MatIntegralSensipTemp);
819:       MatAXPY(th->MatIntegralSensipTemp, 1, quadJp, SAME_NONZERO_PATTERN);
820:       MatAXPY(quadts->mat_sensip, th->time_step0 * th->Theta, th->MatIntegralSensipTemp, SAME_NONZERO_PATTERN);
821:     }
822:   } else {
823:     if (!th->endpoint) MatAXPY(ts->mat_sensip, 1. / th->Theta, MatDeltaFwdSensip, SAME_NONZERO_PATTERN);
824:   }
825:   return 0;
826: }

828: static PetscErrorCode TSForwardGetStages_Theta(TS ts, PetscInt *ns, Mat *stagesensip[])
829: {
830:   TS_Theta *th = (TS_Theta *)ts->data;

832:   if (ns) {
833:     if (!th->endpoint && th->Theta != 1.0) *ns = 1; /* midpoint form */
834:     else *ns = 2;                                   /* endpoint form */
835:   }
836:   if (stagesensip) {
837:     if (!th->endpoint && th->Theta != 1.0) {
838:       th->MatFwdStages[0] = th->MatDeltaFwdSensip;
839:     } else {
840:       th->MatFwdStages[0] = th->MatFwdSensip0;
841:       th->MatFwdStages[1] = ts->mat_sensip; /* stiffly accurate */
842:     }
843:     *stagesensip = th->MatFwdStages;
844:   }
845:   return 0;
846: }

848: /*------------------------------------------------------------*/
849: static PetscErrorCode TSReset_Theta(TS ts)
850: {
851:   TS_Theta *th = (TS_Theta *)ts->data;

853:   VecDestroy(&th->X);
854:   VecDestroy(&th->Xdot);
855:   VecDestroy(&th->X0);
856:   VecDestroy(&th->affine);

858:   VecDestroy(&th->vec_sol_prev);
859:   VecDestroy(&th->vec_lte_work);

861:   VecDestroy(&th->VecCostIntegral0);
862:   return 0;
863: }

865: static PetscErrorCode TSAdjointReset_Theta(TS ts)
866: {
867:   TS_Theta *th = (TS_Theta *)ts->data;

869:   VecDestroyVecs(ts->numcost, &th->VecsDeltaLam);
870:   VecDestroyVecs(ts->numcost, &th->VecsDeltaMu);
871:   VecDestroyVecs(ts->numcost, &th->VecsDeltaLam2);
872:   VecDestroyVecs(ts->numcost, &th->VecsDeltaMu2);
873:   VecDestroyVecs(ts->numcost, &th->VecsSensiTemp);
874:   VecDestroyVecs(ts->numcost, &th->VecsSensi2Temp);
875:   return 0;
876: }

878: static PetscErrorCode TSDestroy_Theta(TS ts)
879: {
880:   TSReset_Theta(ts);
881:   if (ts->dm) {
882:     DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSTheta, DMRestrictHook_TSTheta, ts);
883:     DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSTheta, DMSubDomainRestrictHook_TSTheta, ts);
884:   }
885:   PetscFree(ts->data);
886:   PetscObjectComposeFunction((PetscObject)ts, "TSThetaGetTheta_C", NULL);
887:   PetscObjectComposeFunction((PetscObject)ts, "TSThetaSetTheta_C", NULL);
888:   PetscObjectComposeFunction((PetscObject)ts, "TSThetaGetEndpoint_C", NULL);
889:   PetscObjectComposeFunction((PetscObject)ts, "TSThetaSetEndpoint_C", NULL);
890:   return 0;
891: }

893: /*
894:   This defines the nonlinear equation that is to be solved with SNES
895:   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0

897:   Note that U here is the stage argument. This means that U = U_{n+1} only if endpoint = true,
898:   otherwise U = theta U_{n+1} + (1 - theta) U0, which for the case of implicit midpoint is
899:   U = (U_{n+1} + U0)/2
900: */
901: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes, Vec x, Vec y, TS ts)
902: {
903:   TS_Theta *th = (TS_Theta *)ts->data;
904:   Vec       X0, Xdot;
905:   DM        dm, dmsave;
906:   PetscReal shift = th->shift;

908:   SNESGetDM(snes, &dm);
909:   /* When using the endpoint variant, this is actually 1/Theta * Xdot */
910:   TSThetaGetX0AndXdot(ts, dm, &X0, &Xdot);
911:   if (x != X0) {
912:     VecAXPBYPCZ(Xdot, -shift, shift, 0, X0, x);
913:   } else {
914:     VecZeroEntries(Xdot);
915:   }
916:   /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
917:   dmsave = ts->dm;
918:   ts->dm = dm;
919:   TSComputeIFunction(ts, th->stage_time, x, Xdot, y, PETSC_FALSE);
920:   ts->dm = dmsave;
921:   TSThetaRestoreX0AndXdot(ts, dm, &X0, &Xdot);
922:   return 0;
923: }

925: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes, Vec x, Mat A, Mat B, TS ts)
926: {
927:   TS_Theta *th = (TS_Theta *)ts->data;
928:   Vec       Xdot;
929:   DM        dm, dmsave;
930:   PetscReal shift = th->shift;

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

936:   dmsave = ts->dm;
937:   ts->dm = dm;
938:   TSComputeIJacobian(ts, th->stage_time, x, Xdot, shift, A, B, PETSC_FALSE);
939:   ts->dm = dmsave;
940:   TSThetaRestoreX0AndXdot(ts, dm, NULL, &Xdot);
941:   return 0;
942: }

944: static PetscErrorCode TSForwardSetUp_Theta(TS ts)
945: {
946:   TS_Theta *th     = (TS_Theta *)ts->data;
947:   TS        quadts = ts->quadraturets;

949:   /* combine sensitivities to parameters and sensitivities to initial values into one array */
950:   th->num_tlm = ts->num_parameters;
951:   MatDuplicate(ts->mat_sensip, MAT_DO_NOT_COPY_VALUES, &th->MatDeltaFwdSensip);
952:   if (quadts && quadts->mat_sensip) {
953:     MatDuplicate(quadts->mat_sensip, MAT_DO_NOT_COPY_VALUES, &th->MatIntegralSensipTemp);
954:     MatDuplicate(quadts->mat_sensip, MAT_DO_NOT_COPY_VALUES, &th->MatIntegralSensip0);
955:   }
956:   /* backup sensitivity results for roll-backs */
957:   MatDuplicate(ts->mat_sensip, MAT_DO_NOT_COPY_VALUES, &th->MatFwdSensip0);

959:   VecDuplicate(ts->vec_sol, &th->VecDeltaFwdSensipCol);
960:   return 0;
961: }

963: static PetscErrorCode TSForwardReset_Theta(TS ts)
964: {
965:   TS_Theta *th     = (TS_Theta *)ts->data;
966:   TS        quadts = ts->quadraturets;

968:   if (quadts && quadts->mat_sensip) {
969:     MatDestroy(&th->MatIntegralSensipTemp);
970:     MatDestroy(&th->MatIntegralSensip0);
971:   }
972:   VecDestroy(&th->VecDeltaFwdSensipCol);
973:   MatDestroy(&th->MatDeltaFwdSensip);
974:   MatDestroy(&th->MatFwdSensip0);
975:   return 0;
976: }

978: static PetscErrorCode TSSetUp_Theta(TS ts)
979: {
980:   TS_Theta *th     = (TS_Theta *)ts->data;
981:   TS        quadts = ts->quadraturets;
982:   PetscBool match;

984:   if (!th->VecCostIntegral0 && quadts && ts->costintegralfwd) { /* back up cost integral */
985:     VecDuplicate(quadts->vec_sol, &th->VecCostIntegral0);
986:   }
987:   if (!th->X) VecDuplicate(ts->vec_sol, &th->X);
988:   if (!th->Xdot) VecDuplicate(ts->vec_sol, &th->Xdot);
989:   if (!th->X0) VecDuplicate(ts->vec_sol, &th->X0);
990:   if (th->endpoint) VecDuplicate(ts->vec_sol, &th->affine);

992:   th->order = (th->Theta == 0.5) ? 2 : 1;
993:   th->shift = 1 / (th->Theta * ts->time_step);

995:   TSGetDM(ts, &ts->dm);
996:   DMCoarsenHookAdd(ts->dm, DMCoarsenHook_TSTheta, DMRestrictHook_TSTheta, ts);
997:   DMSubDomainHookAdd(ts->dm, DMSubDomainHook_TSTheta, DMSubDomainRestrictHook_TSTheta, ts);

999:   TSGetAdapt(ts, &ts->adapt);
1000:   TSAdaptCandidatesClear(ts->adapt);
1001:   PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &match);
1002:   if (!match) {
1003:     VecDuplicate(ts->vec_sol, &th->vec_sol_prev);
1004:     VecDuplicate(ts->vec_sol, &th->vec_lte_work);
1005:   }
1006:   TSGetSNES(ts, &ts->snes);

1008:   ts->stifflyaccurate = (!th->endpoint && th->Theta != 1.0) ? PETSC_FALSE : PETSC_TRUE;
1009:   return 0;
1010: }

1012: /*------------------------------------------------------------*/

1014: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
1015: {
1016:   TS_Theta *th = (TS_Theta *)ts->data;

1018:   VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &th->VecsDeltaLam);
1019:   VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &th->VecsSensiTemp);
1020:   if (ts->vecs_sensip) VecDuplicateVecs(ts->vecs_sensip[0], ts->numcost, &th->VecsDeltaMu);
1021:   if (ts->vecs_sensi2) {
1022:     VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &th->VecsDeltaLam2);
1023:     VecDuplicateVecs(ts->vecs_sensi2[0], ts->numcost, &th->VecsSensi2Temp);
1024:     /* hack ts to make implicit TS solver work when users provide only explicit versions of callbacks (RHSFunction,RHSJacobian,RHSHessian etc.) */
1025:     if (!ts->ihessianproduct_fuu) ts->vecs_fuu = ts->vecs_guu;
1026:     if (!ts->ihessianproduct_fup) ts->vecs_fup = ts->vecs_gup;
1027:   }
1028:   if (ts->vecs_sensi2p) {
1029:     VecDuplicateVecs(ts->vecs_sensi2p[0], ts->numcost, &th->VecsDeltaMu2);
1030:     /* hack ts to make implicit TS solver work when users provide only explicit versions of callbacks (RHSFunction,RHSJacobian,RHSHessian etc.) */
1031:     if (!ts->ihessianproduct_fpu) ts->vecs_fpu = ts->vecs_gpu;
1032:     if (!ts->ihessianproduct_fpp) ts->vecs_fpp = ts->vecs_gpp;
1033:   }
1034:   return 0;
1035: }

1037: static PetscErrorCode TSSetFromOptions_Theta(TS ts, PetscOptionItems *PetscOptionsObject)
1038: {
1039:   TS_Theta *th = (TS_Theta *)ts->data;

1041:   PetscOptionsHeadBegin(PetscOptionsObject, "Theta ODE solver options");
1042:   {
1043:     PetscOptionsReal("-ts_theta_theta", "Location of stage (0<Theta<=1)", "TSThetaSetTheta", th->Theta, &th->Theta, NULL);
1044:     PetscOptionsBool("-ts_theta_endpoint", "Use the endpoint instead of midpoint form of the Theta method", "TSThetaSetEndpoint", th->endpoint, &th->endpoint, NULL);
1045:     PetscOptionsBool("-ts_theta_initial_guess_extrapolate", "Extrapolate stage initial guess from previous solution (sometimes unstable)", "TSThetaSetExtrapolate", th->extrapolate, &th->extrapolate, NULL);
1046:   }
1047:   PetscOptionsHeadEnd();
1048:   return 0;
1049: }

1051: static PetscErrorCode TSView_Theta(TS ts, PetscViewer viewer)
1052: {
1053:   TS_Theta *th = (TS_Theta *)ts->data;
1054:   PetscBool iascii;

1056:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1057:   if (iascii) {
1058:     PetscViewerASCIIPrintf(viewer, "  Theta=%g\n", (double)th->Theta);
1059:     PetscViewerASCIIPrintf(viewer, "  Extrapolation=%s\n", th->extrapolate ? "yes" : "no");
1060:   }
1061:   return 0;
1062: }

1064: static PetscErrorCode TSThetaGetTheta_Theta(TS ts, PetscReal *theta)
1065: {
1066:   TS_Theta *th = (TS_Theta *)ts->data;

1068:   *theta = th->Theta;
1069:   return 0;
1070: }

1072: static PetscErrorCode TSThetaSetTheta_Theta(TS ts, PetscReal theta)
1073: {
1074:   TS_Theta *th = (TS_Theta *)ts->data;

1077:   th->Theta = theta;
1078:   th->order = (th->Theta == 0.5) ? 2 : 1;
1079:   return 0;
1080: }

1082: static PetscErrorCode TSThetaGetEndpoint_Theta(TS ts, PetscBool *endpoint)
1083: {
1084:   TS_Theta *th = (TS_Theta *)ts->data;

1086:   *endpoint = th->endpoint;
1087:   return 0;
1088: }

1090: static PetscErrorCode TSThetaSetEndpoint_Theta(TS ts, PetscBool flg)
1091: {
1092:   TS_Theta *th = (TS_Theta *)ts->data;

1094:   th->endpoint = flg;
1095:   return 0;
1096: }

1098: #if defined(PETSC_HAVE_COMPLEX)
1099: static PetscErrorCode TSComputeLinearStability_Theta(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
1100: {
1101:   PetscComplex    z   = xr + xi * PETSC_i, f;
1102:   TS_Theta       *th  = (TS_Theta *)ts->data;
1103:   const PetscReal one = 1.0;

1105:   f   = (one + (one - th->Theta) * z) / (one - th->Theta * z);
1106:   *yr = PetscRealPartComplex(f);
1107:   *yi = PetscImaginaryPartComplex(f);
1108:   return 0;
1109: }
1110: #endif

1112: static PetscErrorCode TSGetStages_Theta(TS ts, PetscInt *ns, Vec *Y[])
1113: {
1114:   TS_Theta *th = (TS_Theta *)ts->data;

1116:   if (ns) {
1117:     if (!th->endpoint && th->Theta != 1.0) *ns = 1; /* midpoint form */
1118:     else *ns = 2;                                   /* endpoint form */
1119:   }
1120:   if (Y) {
1121:     if (!th->endpoint && th->Theta != 1.0) {
1122:       th->Stages[0] = th->X;
1123:     } else {
1124:       th->Stages[0] = th->X0;
1125:       th->Stages[1] = ts->vec_sol; /* stiffly accurate */
1126:     }
1127:     *Y = th->Stages;
1128:   }
1129:   return 0;
1130: }

1132: /* ------------------------------------------------------------ */
1133: /*MC
1134:       TSTHETA - DAE solver using the implicit Theta method

1136:    Level: beginner

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

1143:    Notes:
1144: $  -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
1145: $  -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
1146: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)

1148:    The endpoint variant of the Theta method and backward Euler can be applied to DAE. The midpoint variant is not suitable for DAEs because it is not stiffly accurate.

1150:    The midpoint variant is cast as a 1-stage implicit Runge-Kutta method.

1152: .vb
1153:   Theta | Theta
1154:   -------------
1155:         |  1
1156: .ve

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

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

1162: .vb
1163:   0 | 0         0
1164:   1 | 1-Theta   Theta
1165:   -------------------
1166:     | 1-Theta   Theta
1167: .ve

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

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

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

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

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

1179: M*/
1180: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
1181: {
1182:   TS_Theta *th;

1184:   ts->ops->reset          = TSReset_Theta;
1185:   ts->ops->adjointreset   = TSAdjointReset_Theta;
1186:   ts->ops->destroy        = TSDestroy_Theta;
1187:   ts->ops->view           = TSView_Theta;
1188:   ts->ops->setup          = TSSetUp_Theta;
1189:   ts->ops->adjointsetup   = TSAdjointSetUp_Theta;
1190:   ts->ops->adjointreset   = TSAdjointReset_Theta;
1191:   ts->ops->step           = TSStep_Theta;
1192:   ts->ops->interpolate    = TSInterpolate_Theta;
1193:   ts->ops->evaluatewlte   = TSEvaluateWLTE_Theta;
1194:   ts->ops->rollback       = TSRollBack_Theta;
1195:   ts->ops->setfromoptions = TSSetFromOptions_Theta;
1196:   ts->ops->snesfunction   = SNESTSFormFunction_Theta;
1197:   ts->ops->snesjacobian   = SNESTSFormJacobian_Theta;
1198: #if defined(PETSC_HAVE_COMPLEX)
1199:   ts->ops->linearstability = TSComputeLinearStability_Theta;
1200: #endif
1201:   ts->ops->getstages       = TSGetStages_Theta;
1202:   ts->ops->adjointstep     = TSAdjointStep_Theta;
1203:   ts->ops->adjointintegral = TSAdjointCostIntegral_Theta;
1204:   ts->ops->forwardintegral = TSForwardCostIntegral_Theta;
1205:   ts->default_adapt_type   = TSADAPTNONE;

1207:   ts->ops->forwardsetup     = TSForwardSetUp_Theta;
1208:   ts->ops->forwardreset     = TSForwardReset_Theta;
1209:   ts->ops->forwardstep      = TSForwardStep_Theta;
1210:   ts->ops->forwardgetstages = TSForwardGetStages_Theta;

1212:   ts->usessnes = PETSC_TRUE;

1214:   PetscNew(&th);
1215:   ts->data = (void *)th;

1217:   th->VecsDeltaLam   = NULL;
1218:   th->VecsDeltaMu    = NULL;
1219:   th->VecsSensiTemp  = NULL;
1220:   th->VecsSensi2Temp = NULL;

1222:   th->extrapolate = PETSC_FALSE;
1223:   th->Theta       = 0.5;
1224:   th->order       = 2;
1225:   PetscObjectComposeFunction((PetscObject)ts, "TSThetaGetTheta_C", TSThetaGetTheta_Theta);
1226:   PetscObjectComposeFunction((PetscObject)ts, "TSThetaSetTheta_C", TSThetaSetTheta_Theta);
1227:   PetscObjectComposeFunction((PetscObject)ts, "TSThetaGetEndpoint_C", TSThetaGetEndpoint_Theta);
1228:   PetscObjectComposeFunction((PetscObject)ts, "TSThetaSetEndpoint_C", TSThetaSetEndpoint_Theta);
1229:   return 0;
1230: }

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

1235:   Not Collective

1237:   Input Parameter:
1238: .  ts - timestepping context

1240:   Output Parameter:
1241: .  theta - stage abscissa

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

1246:   Level: Advanced

1248: .seealso: `TSThetaSetTheta()`
1249: @*/
1250: PetscErrorCode TSThetaGetTheta(TS ts, PetscReal *theta)
1251: {
1254:   PetscUseMethod(ts, "TSThetaGetTheta_C", (TS, PetscReal *), (ts, theta));
1255:   return 0;
1256: }

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

1261:   Not Collective

1263:   Input Parameters:
1264: +  ts - timestepping context
1265: -  theta - stage abscissa

1267:   Options Database:
1268: .  -ts_theta_theta <theta> - set theta

1270:   Level: Intermediate

1272: .seealso: `TSThetaGetTheta()`
1273: @*/
1274: PetscErrorCode TSThetaSetTheta(TS ts, PetscReal theta)
1275: {
1277:   PetscTryMethod(ts, "TSThetaSetTheta_C", (TS, PetscReal), (ts, theta));
1278:   return 0;
1279: }

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

1284:   Not Collective

1286:   Input Parameter:
1287: .  ts - timestepping context

1289:   Output Parameter:
1290: .  endpoint - PETSC_TRUE when using the endpoint variant

1292:   Level: Advanced

1294: .seealso: `TSThetaSetEndpoint()`, `TSTHETA`, `TSCN`
1295: @*/
1296: PetscErrorCode TSThetaGetEndpoint(TS ts, PetscBool *endpoint)
1297: {
1300:   PetscUseMethod(ts, "TSThetaGetEndpoint_C", (TS, PetscBool *), (ts, endpoint));
1301:   return 0;
1302: }

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

1307:   Not Collective

1309:   Input Parameters:
1310: +  ts - timestepping context
1311: -  flg - PETSC_TRUE to use the endpoint variant

1313:   Options Database:
1314: .  -ts_theta_endpoint <flg> - use the endpoint variant

1316:   Level: Intermediate

1318: .seealso: `TSTHETA`, `TSCN`
1319: @*/
1320: PetscErrorCode TSThetaSetEndpoint(TS ts, PetscBool flg)
1321: {
1323:   PetscTryMethod(ts, "TSThetaSetEndpoint_C", (TS, PetscBool), (ts, flg));
1324:   return 0;
1325: }

1327: /*
1328:  * TSBEULER and TSCN are straightforward specializations of TSTHETA.
1329:  * The creation functions for these specializations are below.
1330:  */

1332: static PetscErrorCode TSSetUp_BEuler(TS ts)
1333: {
1334:   TS_Theta *th = (TS_Theta *)ts->data;

1338:   TSSetUp_Theta(ts);
1339:   return 0;
1340: }

1342: static PetscErrorCode TSView_BEuler(TS ts, PetscViewer viewer)
1343: {
1344:   return 0;
1345: }

1347: /*MC
1348:       TSBEULER - ODE solver using the implicit backward Euler method

1350:   Level: beginner

1352:   Notes:
1353:   TSBEULER is equivalent to TSTHETA with Theta=1.0

1355: $  -ts_type theta -ts_theta_theta 1.0

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

1359: M*/
1360: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
1361: {
1362:   TSCreate_Theta(ts);
1363:   TSThetaSetTheta(ts, 1.0);
1364:   TSThetaSetEndpoint(ts, PETSC_FALSE);
1365:   ts->ops->setup = TSSetUp_BEuler;
1366:   ts->ops->view  = TSView_BEuler;
1367:   return 0;
1368: }

1370: static PetscErrorCode TSSetUp_CN(TS ts)
1371: {
1372:   TS_Theta *th = (TS_Theta *)ts->data;

1376:   TSSetUp_Theta(ts);
1377:   return 0;
1378: }

1380: static PetscErrorCode TSView_CN(TS ts, PetscViewer viewer)
1381: {
1382:   return 0;
1383: }

1385: /*MC
1386:       TSCN - ODE solver using the implicit Crank-Nicolson method.

1388:   Level: beginner

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

1393: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint

1395: .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSBEULER`, `TSTHETA`

1397: M*/
1398: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
1399: {
1400:   TSCreate_Theta(ts);
1401:   TSThetaSetTheta(ts, 0.5);
1402:   TSThetaSetEndpoint(ts, PETSC_TRUE);
1403:   ts->ops->setup = TSSetUp_CN;
1404:   ts->ops->view  = TSView_CN;
1405:   return 0;
1406: }