Actual source code: tssen.c

  1: #include <petsc/private/tsimpl.h>
  2: #include <petscdraw.h>

  4: PetscLogEvent TS_AdjointStep, TS_ForwardStep, TS_JacobianPEval;

  6: /* #define TSADJOINT_STAGE */

  8: /* ------------------------ Sensitivity Context ---------------------------*/

 10: /*@C
 11:   TSSetRHSJacobianP - Sets the function that computes the Jacobian of G w.r.t. the parameters P where U_t = G(U,P,t), as well as the location to store the matrix.

 13:   Logically Collective

 15:   Input Parameters:
 16: + ts - `TS` context obtained from `TSCreate()`
 17: . Amat - JacobianP matrix
 18: . func - function
 19: - ctx - [optional] user-defined function context

 21:   Calling sequence of func:
 22: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
 23: +   t - current timestep
 24: .   U - input vector (current ODE solution)
 25: .   A - output matrix
 26: -   ctx - [optional] user-defined function context

 28:   Level: intermediate

 30:   Note:
 31:     Amat has the same number of rows and the same row parallel layout as u, Amat has the same number of columns and parallel layout as p

 33: .seealso: [](chapter_ts), `TS`, `TSGetRHSJacobianP()`
 34: @*/
 35: PetscErrorCode TSSetRHSJacobianP(TS ts, Mat Amat, PetscErrorCode (*func)(TS, PetscReal, Vec, Mat, void *), void *ctx)
 36: {

 40:   ts->rhsjacobianp    = func;
 41:   ts->rhsjacobianpctx = ctx;
 42:   if (Amat) {
 43:     PetscObjectReference((PetscObject)Amat);
 44:     MatDestroy(&ts->Jacprhs);
 45:     ts->Jacprhs = Amat;
 46:   }
 47:   return 0;
 48: }

 50: /*@C
 51:   TSGetRHSJacobianP - Gets the function that computes the Jacobian of G w.r.t. the parameters P where U_t = G(U,P,t), as well as the location to store the matrix.

 53:   Logically Collective

 55:   Input Parameter:
 56: . ts - `TS` context obtained from `TSCreate()`

 58:   Output Parameters:
 59: + Amat - JacobianP matrix
 60: . func - function
 61: - ctx - [optional] user-defined function context

 63:   Calling sequence of func:
 64: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
 65: +   t - current timestep
 66: .   U - input vector (current ODE solution)
 67: .   A - output matrix
 68: -   ctx - [optional] user-defined function context

 70:   Level: intermediate

 72:   Note:
 73:     Amat has the same number of rows and the same row parallel layout as u, Amat has the same number of columns and parallel layout as p

 75: .seealso: [](chapter_ts), `TSSetRHSJacobianP()`, `TS`, `TSGetRHSJacobianP()`
 76: @*/
 77: PetscErrorCode TSGetRHSJacobianP(TS ts, Mat *Amat, PetscErrorCode (**func)(TS, PetscReal, Vec, Mat, void *), void **ctx)
 78: {
 79:   if (func) *func = ts->rhsjacobianp;
 80:   if (ctx) *ctx = ts->rhsjacobianpctx;
 81:   if (Amat) *Amat = ts->Jacprhs;
 82:   return 0;
 83: }

 85: /*@C
 86:   TSComputeRHSJacobianP - Runs the user-defined JacobianP function.

 88:   Collective

 90:   Input Parameters:
 91: . ts   - The `TS` context obtained from `TSCreate()`

 93:   Level: developer

 95: .seealso: [](chapter_ts), `TSSetRHSJacobianP()`, `TS`
 96: @*/
 97: PetscErrorCode TSComputeRHSJacobianP(TS ts, PetscReal t, Vec U, Mat Amat)
 98: {
 99:   if (!Amat) return 0;

103:   PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->rhsjacobianp)(ts, t, U, Amat, ts->rhsjacobianpctx));
104:   return 0;
105: }

107: /*@C
108:   TSSetIJacobianP - Sets the function that computes the Jacobian of F w.r.t. the parameters P where F(Udot,U,t) = G(U,P,t), as well as the location to store the matrix.

110:   Logically Collective

112:   Input Parameters:
113: + ts - `TS` context obtained from `TSCreate()`
114: . Amat - JacobianP matrix
115: . func - function
116: - ctx - [optional] user-defined function context

118:   Calling sequence of func:
119: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
120: +   t - current timestep
121: .   U - input vector (current ODE solution)
122: .   Udot - time derivative of state vector
123: .   shift - shift to apply, see note below
124: .   A - output matrix
125: -   ctx - [optional] user-defined function context

127:   Level: intermediate

129:   Note:
130:     Amat has the same number of rows and the same row parallel layout as u, Amat has the same number of columns and parallel layout as p

132: .seealso: [](chapter_ts), `TSSetRHSJacobianP()`, `TS`
133: @*/
134: PetscErrorCode TSSetIJacobianP(TS ts, Mat Amat, PetscErrorCode (*func)(TS, PetscReal, Vec, Vec, PetscReal, Mat, void *), void *ctx)
135: {

139:   ts->ijacobianp    = func;
140:   ts->ijacobianpctx = ctx;
141:   if (Amat) {
142:     PetscObjectReference((PetscObject)Amat);
143:     MatDestroy(&ts->Jacp);
144:     ts->Jacp = Amat;
145:   }
146:   return 0;
147: }

149: /*@C
150:   TSComputeIJacobianP - Runs the user-defined IJacobianP function.

152:   Collective

154:   Input Parameters:
155: + ts - the `TS` context
156: . t - current timestep
157: . U - state vector
158: . Udot - time derivative of state vector
159: . shift - shift to apply, see note below
160: - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate

162:   Output Parameters:
163: . A - Jacobian matrix

165:   Level: developer

167: .seealso: [](chapter_ts), `TS`, `TSSetIJacobianP()`
168: @*/
169: PetscErrorCode TSComputeIJacobianP(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat Amat, PetscBool imex)
170: {
171:   if (!Amat) return 0;

176:   PetscLogEventBegin(TS_JacobianPEval, ts, U, Amat, 0);
177:   if (ts->ijacobianp) PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->ijacobianp)(ts, t, U, Udot, shift, Amat, ts->ijacobianpctx));
178:   if (imex) {
179:     if (!ts->ijacobianp) { /* system was written as Udot = G(t,U) */
180:       PetscBool assembled;
181:       MatZeroEntries(Amat);
182:       MatAssembled(Amat, &assembled);
183:       if (!assembled) {
184:         MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
185:         MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
186:       }
187:     }
188:   } else {
189:     if (ts->rhsjacobianp) TSComputeRHSJacobianP(ts, t, U, ts->Jacprhs);
190:     if (ts->Jacprhs == Amat) { /* No IJacobian, so we only have the RHS matrix */
191:       MatScale(Amat, -1);
192:     } else if (ts->Jacprhs) { /* Both IJacobian and RHSJacobian */
193:       MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
194:       if (!ts->ijacobianp) { /* No IJacobianp provided, but we have a separate RHS matrix */
195:         MatZeroEntries(Amat);
196:       }
197:       MatAXPY(Amat, -1, ts->Jacprhs, axpy);
198:     }
199:   }
200:   PetscLogEventEnd(TS_JacobianPEval, ts, U, Amat, 0);
201:   return 0;
202: }

204: /*@C
205:     TSSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions

207:     Logically Collective

209:     Input Parameters:
210: +   ts - the `TS` context obtained from `TSCreate()`
211: .   numcost - number of gradients to be computed, this is the number of cost functions
212: .   costintegral - vector that stores the integral values
213: .   rf - routine for evaluating the integrand function
214: .   drduf - function that computes the gradients of the r's with respect to u
215: .   drdpf - function that computes the gradients of the r's with respect to p, can be NULL if parametric sensitivity is not desired (mu=NULL)
216: .   fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
217: -   ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

219:     Calling sequence of rf:
220: $   PetscErrorCode rf(TS ts,PetscReal t,Vec U,Vec F,void *ctx);

222:     Calling sequence of drduf:
223: $   PetscErroCode drduf(TS ts,PetscReal t,Vec U,Vec *dRdU,void *ctx);

225:     Calling sequence of drdpf:
226: $   PetscErroCode drdpf(TS ts,PetscReal t,Vec U,Vec *dRdP,void *ctx);

228:     Level: deprecated

230:     Note:
231:     For optimization there is usually a single cost function (numcost = 1). For sensitivities there may be multiple cost functions

233: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobianP()`, `TSGetCostGradients()`, `TSSetCostGradients()`
234: @*/
235: PetscErrorCode TSSetCostIntegrand(TS ts, PetscInt numcost, Vec costintegral, PetscErrorCode (*rf)(TS, PetscReal, Vec, Vec, void *), PetscErrorCode (*drduf)(TS, PetscReal, Vec, Vec *, void *), PetscErrorCode (*drdpf)(TS, PetscReal, Vec, Vec *, void *), PetscBool fwd, void *ctx)
236: {
240:   if (!ts->numcost) ts->numcost = numcost;

242:   if (costintegral) {
243:     PetscObjectReference((PetscObject)costintegral);
244:     VecDestroy(&ts->vec_costintegral);
245:     ts->vec_costintegral = costintegral;
246:   } else {
247:     if (!ts->vec_costintegral) { /* Create a seq vec if user does not provide one */
248:       VecCreateSeq(PETSC_COMM_SELF, numcost, &ts->vec_costintegral);
249:     } else {
250:       VecSet(ts->vec_costintegral, 0.0);
251:     }
252:   }
253:   if (!ts->vec_costintegrand) {
254:     VecDuplicate(ts->vec_costintegral, &ts->vec_costintegrand);
255:   } else {
256:     VecSet(ts->vec_costintegrand, 0.0);
257:   }
258:   ts->costintegralfwd  = fwd; /* Evaluate the cost integral in forward run if fwd is true */
259:   ts->costintegrand    = rf;
260:   ts->costintegrandctx = ctx;
261:   ts->drdufunction     = drduf;
262:   ts->drdpfunction     = drdpf;
263:   return 0;
264: }

266: /*@C
267:    TSGetCostIntegral - Returns the values of the integral term in the cost functions.
268:    It is valid to call the routine after a backward run.

270:    Not Collective

272:    Input Parameter:
273: .  ts - the `TS` context obtained from `TSCreate()`

275:    Output Parameter:
276: .  v - the vector containing the integrals for each cost function

278:    Level: intermediate

280: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, ``TSSetCostIntegrand()`
281: @*/
282: PetscErrorCode TSGetCostIntegral(TS ts, Vec *v)
283: {
284:   TS quadts;

288:   TSGetQuadratureTS(ts, NULL, &quadts);
289:   *v = quadts->vec_sol;
290:   return 0;
291: }

293: /*@C
294:    TSComputeCostIntegrand - Evaluates the integral function in the cost functions.

296:    Input Parameters:
297: +  ts - the `TS` context
298: .  t - current time
299: -  U - state vector, i.e. current solution

301:    Output Parameter:
302: .  Q - vector of size numcost to hold the outputs

304:    Level: deprecated

306:    Note:
307:    Most users should not need to explicitly call this routine, as it
308:    is used internally within the sensitivity analysis context.

310: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSSetCostIntegrand()`
311: @*/
312: PetscErrorCode TSComputeCostIntegrand(TS ts, PetscReal t, Vec U, Vec Q)
313: {

318:   PetscLogEventBegin(TS_FunctionEval, ts, U, Q, 0);
319:   if (ts->costintegrand) PetscCallBack("TS callback integrand in the cost function", (*ts->costintegrand)(ts, t, U, Q, ts->costintegrandctx));
320:   else VecZeroEntries(Q);
321:   PetscLogEventEnd(TS_FunctionEval, ts, U, Q, 0);
322:   return 0;
323: }

325: /*@C
326:   TSComputeDRDUFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobian()`

328:   Level: deprecated

330: @*/
331: PetscErrorCode TSComputeDRDUFunction(TS ts, PetscReal t, Vec U, Vec *DRDU)
332: {
333:   if (!DRDU) return 0;

337:   PetscCallBack("TS callback DRDU for sensitivity analysis", (*ts->drdufunction)(ts, t, U, DRDU, ts->costintegrandctx));
338:   return 0;
339: }

341: /*@C
342:   TSComputeDRDPFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobianP()`

344:   Level: deprecated

346: @*/
347: PetscErrorCode TSComputeDRDPFunction(TS ts, PetscReal t, Vec U, Vec *DRDP)
348: {
349:   if (!DRDP) return 0;

353:   PetscCallBack("TS callback DRDP for sensitivity analysis", (*ts->drdpfunction)(ts, t, U, DRDP, ts->costintegrandctx));
354:   return 0;
355: }

357: /*@C
358:   TSSetIHessianProduct - Sets the function that computes the vector-Hessian-vector product. The Hessian is the second-order derivative of F (IFunction) w.r.t. the state variable.

360:   Logically Collective

362:   Input Parameters:
363: + ts - `TS` context obtained from `TSCreate()`
364: . ihp1 - an array of vectors storing the result of vector-Hessian-vector product for F_UU
365: . hessianproductfunc1 - vector-Hessian-vector product function for F_UU
366: . ihp2 - an array of vectors storing the result of vector-Hessian-vector product for F_UP
367: . hessianproductfunc2 - vector-Hessian-vector product function for F_UP
368: . ihp3 - an array of vectors storing the result of vector-Hessian-vector product for F_PU
369: . hessianproductfunc3 - vector-Hessian-vector product function for F_PU
370: . ihp4 - an array of vectors storing the result of vector-Hessian-vector product for F_PP
371: - hessianproductfunc4 - vector-Hessian-vector product function for F_PP

373:   Calling sequence of ihessianproductfunc:
374: $ ihessianproductfunc (TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx);
375: +   t - current timestep
376: .   U - input vector (current ODE solution)
377: .   Vl - an array of input vectors to be left-multiplied with the Hessian
378: .   Vr - input vector to be right-multiplied with the Hessian
379: .   VHV - an array of output vectors for vector-Hessian-vector product
380: -   ctx - [optional] user-defined function context

382:   Level: intermediate

384:   Notes:
385:   The first Hessian function and the working array are required.
386:   As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
387:   $ Vl_n^T*F_UP*Vr
388:   where the vector Vl_n (n-th element in the array Vl) and Vr are of size N and M respectively, and the Hessian F_UP is of size N x N x M.
389:   Each entry of F_UP corresponds to the derivative
390:   $ F_UP[i][j][k] = \frac{\partial^2 F[i]}{\partial U[j] \partial P[k]}.
391:   The result of the vector-Hessian-vector product for Vl_n needs to be stored in vector VHV_n with the j-th entry being
392:   $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * F_UP[i][j][k] * Vr[k]}
393:   If the cost function is a scalar, there will be only one vector in Vl and VHV.

395: .seealso: [](chapter_ts), `TS`
396: @*/
397: PetscErrorCode TSSetIHessianProduct(TS ts, Vec *ihp1, PetscErrorCode (*ihessianproductfunc1)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *ihp2, PetscErrorCode (*ihessianproductfunc2)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *ihp3, PetscErrorCode (*ihessianproductfunc3)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *ihp4, PetscErrorCode (*ihessianproductfunc4)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), void *ctx)
398: {

402:   ts->ihessianproductctx = ctx;
403:   if (ihp1) ts->vecs_fuu = ihp1;
404:   if (ihp2) ts->vecs_fup = ihp2;
405:   if (ihp3) ts->vecs_fpu = ihp3;
406:   if (ihp4) ts->vecs_fpp = ihp4;
407:   ts->ihessianproduct_fuu = ihessianproductfunc1;
408:   ts->ihessianproduct_fup = ihessianproductfunc2;
409:   ts->ihessianproduct_fpu = ihessianproductfunc3;
410:   ts->ihessianproduct_fpp = ihessianproductfunc4;
411:   return 0;
412: }

414: /*@C
415:   TSComputeIHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Fuu.

417:   Collective

419:   Input Parameters:
420: . ts   - The `TS` context obtained from `TSCreate()`

422:   Level: developer

424:   Note:
425:   `TSComputeIHessianProductFunctionUU()` is typically used for sensitivity implementation,
426:   so most users would not generally call this routine themselves.

428: .seealso: [](chapter_ts), `TSSetIHessianProduct()`
429: @*/
430: PetscErrorCode TSComputeIHessianProductFunctionUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
431: {
432:   if (!VHV) return 0;

436:   if (ts->ihessianproduct_fuu) PetscCallBack("TS callback IHessianProduct 1 for sensitivity analysis", (*ts->ihessianproduct_fuu)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));

438:   /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
439:   if (ts->rhshessianproduct_guu) {
440:     PetscInt nadj;
441:     TSComputeRHSHessianProductFunctionUU(ts, t, U, Vl, Vr, VHV);
442:     for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
443:   }
444:   return 0;
445: }

447: /*@C
448:   TSComputeIHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Fup.

450:   Collective

452:   Input Parameters:
453: . ts   - The `TS` context obtained from `TSCreate()`

455:   Level: developer

457:   Note:
458:   `TSComputeIHessianProductFunctionUP()` is typically used for sensitivity implementation,
459:   so most users would not generally call this routine themselves.

461: .seealso: [](chapter_ts), `TSSetIHessianProduct()`
462: @*/
463: PetscErrorCode TSComputeIHessianProductFunctionUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
464: {
465:   if (!VHV) return 0;

469:   if (ts->ihessianproduct_fup) PetscCallBack("TS callback IHessianProduct 2 for sensitivity analysis", (*ts->ihessianproduct_fup)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));

471:   /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
472:   if (ts->rhshessianproduct_gup) {
473:     PetscInt nadj;
474:     TSComputeRHSHessianProductFunctionUP(ts, t, U, Vl, Vr, VHV);
475:     for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
476:   }
477:   return 0;
478: }

480: /*@C
481:   TSComputeIHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Fpu.

483:   Collective

485:   Input Parameters:
486: . ts   - The `TS` context obtained from `TSCreate()`

488:   Level: developer

490:   Note:
491:   `TSComputeIHessianProductFunctionPU()` is typically used for sensitivity implementation,
492:   so most users would not generally call this routine themselves.

494: .seealso: [](chapter_ts), `TSSetIHessianProduct()`
495: @*/
496: PetscErrorCode TSComputeIHessianProductFunctionPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
497: {
498:   if (!VHV) return 0;

502:   if (ts->ihessianproduct_fpu) PetscCallBack("TS callback IHessianProduct 3 for sensitivity analysis", (*ts->ihessianproduct_fpu)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));

504:   /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
505:   if (ts->rhshessianproduct_gpu) {
506:     PetscInt nadj;
507:     TSComputeRHSHessianProductFunctionPU(ts, t, U, Vl, Vr, VHV);
508:     for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
509:   }
510:   return 0;
511: }

513: /*@C
514:   TSComputeIHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Fpp.

516:   Collective

518:   Input Parameters:
519: . ts   - The `TS` context obtained from `TSCreate()`

521:   Level: developer

523:   Note:
524:   `TSComputeIHessianProductFunctionPP()` is typically used for sensitivity implementation,
525:   so most users would not generally call this routine themselves.

527: .seealso: [](chapter_ts), `TSSetIHessianProduct()`
528: @*/
529: PetscErrorCode TSComputeIHessianProductFunctionPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
530: {
531:   if (!VHV) return 0;

535:   if (ts->ihessianproduct_fpp) PetscCallBack("TS callback IHessianProduct 3 for sensitivity analysis", (*ts->ihessianproduct_fpp)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));

537:   /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
538:   if (ts->rhshessianproduct_gpp) {
539:     PetscInt nadj;
540:     TSComputeRHSHessianProductFunctionPP(ts, t, U, Vl, Vr, VHV);
541:     for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
542:   }
543:   return 0;
544: }

546: /*@C
547:   TSSetRHSHessianProduct - Sets the function that computes the vector-Hessian-vector product. The Hessian is the second-order derivative of G (RHSFunction) w.r.t. the state variable.

549:   Logically Collective

551:   Input Parameters:
552: + ts - `TS` context obtained from `TSCreate()`
553: . rhshp1 - an array of vectors storing the result of vector-Hessian-vector product for G_UU
554: . hessianproductfunc1 - vector-Hessian-vector product function for G_UU
555: . rhshp2 - an array of vectors storing the result of vector-Hessian-vector product for G_UP
556: . hessianproductfunc2 - vector-Hessian-vector product function for G_UP
557: . rhshp3 - an array of vectors storing the result of vector-Hessian-vector product for G_PU
558: . hessianproductfunc3 - vector-Hessian-vector product function for G_PU
559: . rhshp4 - an array of vectors storing the result of vector-Hessian-vector product for G_PP
560: - hessianproductfunc4 - vector-Hessian-vector product function for G_PP

562:   Calling sequence of ihessianproductfunc:
563: $ rhshessianproductfunc (TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx);
564: +   t - current timestep
565: .   U - input vector (current ODE solution)
566: .   Vl - an array of input vectors to be left-multiplied with the Hessian
567: .   Vr - input vector to be right-multiplied with the Hessian
568: .   VHV - an array of output vectors for vector-Hessian-vector product
569: -   ctx - [optional] user-defined function context

571:   Level: intermediate

573:   Notes:
574:   The first Hessian function and the working array are required.
575:   As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
576:   $ Vl_n^T*G_UP*Vr
577:   where the vector Vl_n (n-th element in the array Vl) and Vr are of size N and M respectively, and the Hessian G_UP is of size N x N x M.
578:   Each entry of G_UP corresponds to the derivative
579:   $ G_UP[i][j][k] = \frac{\partial^2 G[i]}{\partial U[j] \partial P[k]}.
580:   The result of the vector-Hessian-vector product for Vl_n needs to be stored in vector VHV_n with j-th entry being
581:   $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * G_UP[i][j][k] * Vr[k]}
582:   If the cost function is a scalar, there will be only one vector in Vl and VHV.

584: .seealso: `TS`
585: @*/
586: PetscErrorCode TSSetRHSHessianProduct(TS ts, Vec *rhshp1, PetscErrorCode (*rhshessianproductfunc1)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *rhshp2, PetscErrorCode (*rhshessianproductfunc2)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *rhshp3, PetscErrorCode (*rhshessianproductfunc3)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *rhshp4, PetscErrorCode (*rhshessianproductfunc4)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), void *ctx)
587: {

591:   ts->rhshessianproductctx = ctx;
592:   if (rhshp1) ts->vecs_guu = rhshp1;
593:   if (rhshp2) ts->vecs_gup = rhshp2;
594:   if (rhshp3) ts->vecs_gpu = rhshp3;
595:   if (rhshp4) ts->vecs_gpp = rhshp4;
596:   ts->rhshessianproduct_guu = rhshessianproductfunc1;
597:   ts->rhshessianproduct_gup = rhshessianproductfunc2;
598:   ts->rhshessianproduct_gpu = rhshessianproductfunc3;
599:   ts->rhshessianproduct_gpp = rhshessianproductfunc4;
600:   return 0;
601: }

603: /*@C
604:   TSComputeRHSHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Guu.

606:   Collective

608:   Input Parameters:
609: . ts   - The `TS` context obtained from `TSCreate()`

611:   Level: developer

613:   Note:
614:   `TSComputeRHSHessianProductFunctionUU()` is typically used for sensitivity implementation,
615:   so most users would not generally call this routine themselves.

617: .seealso: [](chapter_ts), `TS`, `TSSetRHSHessianProduct()`
618: @*/
619: PetscErrorCode TSComputeRHSHessianProductFunctionUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
620: {
621:   if (!VHV) return 0;

625:   PetscCallBack("TS callback RHSHessianProduct 1 for sensitivity analysis", (*ts->rhshessianproduct_guu)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
626:   return 0;
627: }

629: /*@C
630:   TSComputeRHSHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Gup.

632:   Collective

634:   Input Parameters:
635: . ts   - The `TS` context obtained from `TSCreate()`

637:   Level: developer

639:   Note:
640:   `TSComputeRHSHessianProductFunctionUP()` is typically used for sensitivity implementation,
641:   so most users would not generally call this routine themselves.

643: .seealso: [](chapter_ts), `TS`, `TSSetRHSHessianProduct()`
644: @*/
645: PetscErrorCode TSComputeRHSHessianProductFunctionUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
646: {
647:   if (!VHV) return 0;

651:   PetscCallBack("TS callback RHSHessianProduct 2 for sensitivity analysis", (*ts->rhshessianproduct_gup)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
652:   return 0;
653: }

655: /*@C
656:   TSComputeRHSHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Gpu.

658:   Collective

660:   Input Parameters:
661: . ts   - The `TS` context obtained from `TSCreate()`

663:   Level: developer

665:   Note:
666:   `TSComputeRHSHessianProductFunctionPU()` is typically used for sensitivity implementation,
667:   so most users would not generally call this routine themselves.

669: .seealso: [](chapter_ts), `TSSetRHSHessianProduct()`
670: @*/
671: PetscErrorCode TSComputeRHSHessianProductFunctionPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
672: {
673:   if (!VHV) return 0;

677:   PetscCallBack("TS callback RHSHessianProduct 3 for sensitivity analysis", (*ts->rhshessianproduct_gpu)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
678:   return 0;
679: }

681: /*@C
682:   TSComputeRHSHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Gpp.

684:   Collective

686:   Input Parameters:
687: . ts   - The `TS` context obtained from `TSCreate()`

689:   Level: developer

691:   Note:
692:   `TSComputeRHSHessianProductFunctionPP()` is typically used for sensitivity implementation,
693:   so most users would not generally call this routine themselves.

695: .seealso: [](chapter_ts), `TSSetRHSHessianProduct()`
696: @*/
697: PetscErrorCode TSComputeRHSHessianProductFunctionPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
698: {
699:   if (!VHV) return 0;

703:   PetscCallBack("TS callback RHSHessianProduct 3 for sensitivity analysis", (*ts->rhshessianproduct_gpp)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
704:   return 0;
705: }

707: /* --------------------------- Adjoint sensitivity ---------------------------*/

709: /*@
710:    TSSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial values and w.r.t. the problem parameters
711:       for use by the `TS` adjoint routines.

713:    Logically Collective

715:    Input Parameters:
716: +  ts - the `TS` context obtained from `TSCreate()`
717: .  numcost - number of gradients to be computed, this is the number of cost functions
718: .  lambda - gradients with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
719: -  mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters

721:    Level: beginner

723:    Notes:
724:     the entries in these vectors must be correctly initialized with the values lambda_i = df/dy|finaltime  mu_i = df/dp|finaltime

726:    After `TSAdjointSolve()` is called the lambda and the mu contain the computed sensitivities

728: .seealso: `TS`, `TSAdjointSolve()`, `TSGetCostGradients()`
729: @*/
730: PetscErrorCode TSSetCostGradients(TS ts, PetscInt numcost, Vec *lambda, Vec *mu)
731: {
734:   ts->vecs_sensi  = lambda;
735:   ts->vecs_sensip = mu;
737:   ts->numcost = numcost;
738:   return 0;
739: }

741: /*@
742:    TSGetCostGradients - Returns the gradients from the `TSAdjointSolve()`

744:    Not Collective, but the vectors returned are parallel if `TS` is parallel

746:    Input Parameter:
747: .  ts - the `TS` context obtained from `TSCreate()`

749:    Output Parameters:
750: +  numcost - size of returned arrays
751: .  lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
752: -  mu - vectors containing the gradients of the cost functions with respect to the problem parameters

754:    Level: intermediate

756: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSSetCostGradients()`
757: @*/
758: PetscErrorCode TSGetCostGradients(TS ts, PetscInt *numcost, Vec **lambda, Vec **mu)
759: {
761:   if (numcost) *numcost = ts->numcost;
762:   if (lambda) *lambda = ts->vecs_sensi;
763:   if (mu) *mu = ts->vecs_sensip;
764:   return 0;
765: }

767: /*@
768:    TSSetCostHessianProducts - Sets the initial value of the Hessian-vector products of the cost function w.r.t. initial values and w.r.t. the problem parameters
769:    for use by the `TS` adjoint routines.

771:    Logically Collective

773:    Input Parameters:
774: +  ts - the `TS` context obtained from `TSCreate()`
775: .  numcost - number of cost functions
776: .  lambda2 - Hessian-vector product with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
777: .  mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
778: -  dir - the direction vector that are multiplied with the Hessian of the cost functions

780:    Level: beginner

782:    Notes:
783:    Hessian of the cost function is completely different from Hessian of the ODE/DAE system

785:    For second-order adjoint, one needs to call this function and then `TSAdjointSetForward()` before `TSSolve()`.

787:    After `TSAdjointSolve()` is called, the lambda2 and the mu2 will contain the computed second-order adjoint sensitivities, and can be used to produce Hessian-vector product (not the full Hessian matrix). Users must provide a direction vector; it is usually generated by an optimization solver.

789:    Passing NULL for lambda2 disables the second-order calculation.

791: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSAdjointSetForward()`
792: @*/
793: PetscErrorCode TSSetCostHessianProducts(TS ts, PetscInt numcost, Vec *lambda2, Vec *mu2, Vec dir)
794: {
797:   ts->numcost      = numcost;
798:   ts->vecs_sensi2  = lambda2;
799:   ts->vecs_sensi2p = mu2;
800:   ts->vec_dir      = dir;
801:   return 0;
802: }

804: /*@
805:    TSGetCostHessianProducts - Returns the gradients from the `TSAdjointSolve()`

807:    Not Collective, but vectors returned are parallel if `TS` is parallel

809:    Input Parameter:
810: .  ts - the `TS` context obtained from `TSCreate()`

812:    Output Parameters:
813: +  numcost - number of cost functions
814: .  lambda2 - Hessian-vector product with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
815: .  mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
816: -  dir - the direction vector that are multiplied with the Hessian of the cost functions

818:    Level: intermediate

820: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSSetCostHessianProducts()`
821: @*/
822: PetscErrorCode TSGetCostHessianProducts(TS ts, PetscInt *numcost, Vec **lambda2, Vec **mu2, Vec *dir)
823: {
825:   if (numcost) *numcost = ts->numcost;
826:   if (lambda2) *lambda2 = ts->vecs_sensi2;
827:   if (mu2) *mu2 = ts->vecs_sensi2p;
828:   if (dir) *dir = ts->vec_dir;
829:   return 0;
830: }

832: /*@
833:   TSAdjointSetForward - Trigger the tangent linear solver and initialize the forward sensitivities

835:   Logically Collective

837:   Input Parameters:
838: +  ts - the `TS` context obtained from `TSCreate()`
839: -  didp - the derivative of initial values w.r.t. parameters

841:   Level: intermediate

843:   Notes:
844:   When computing sensitivies w.r.t. initial condition, set didp to NULL so that the solver will take it as an identity matrix mathematically.
845:   `TS` adjoint does not reset the tangent linear solver automatically, `TSAdjointResetForward()` should be called to reset the tangent linear solver.

847: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSSetCostHessianProducts()`, `TSAdjointResetForward()`
848: @*/
849: PetscErrorCode TSAdjointSetForward(TS ts, Mat didp)
850: {
851:   Mat          A;
852:   Vec          sp;
853:   PetscScalar *xarr;
854:   PetscInt     lsize;

856:   ts->forward_solve = PETSC_TRUE; /* turn on tangent linear mode */
859:   /* create a single-column dense matrix */
860:   VecGetLocalSize(ts->vec_sol, &lsize);
861:   MatCreateDense(PetscObjectComm((PetscObject)ts), lsize, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, &A);

863:   VecDuplicate(ts->vec_sol, &sp);
864:   MatDenseGetColumn(A, 0, &xarr);
865:   VecPlaceArray(sp, xarr);
866:   if (ts->vecs_sensi2p) { /* tangent linear variable initialized as 2*dIdP*dir */
867:     if (didp) {
868:       MatMult(didp, ts->vec_dir, sp);
869:       VecScale(sp, 2.);
870:     } else {
871:       VecZeroEntries(sp);
872:     }
873:   } else { /* tangent linear variable initialized as dir */
874:     VecCopy(ts->vec_dir, sp);
875:   }
876:   VecResetArray(sp);
877:   MatDenseRestoreColumn(A, &xarr);
878:   VecDestroy(&sp);

880:   TSForwardSetInitialSensitivities(ts, A); /* if didp is NULL, identity matrix is assumed */

882:   MatDestroy(&A);
883:   return 0;
884: }

886: /*@
887:   TSAdjointResetForward - Reset the tangent linear solver and destroy the tangent linear context

889:   Logically Collective

891:   Input Parameters:
892: .  ts - the `TS` context obtained from `TSCreate()`

894:   Level: intermediate

896: .seealso: [](chapter_ts), `TSAdjointSetForward()`
897: @*/
898: PetscErrorCode TSAdjointResetForward(TS ts)
899: {
900:   ts->forward_solve = PETSC_FALSE; /* turn off tangent linear mode */
901:   TSForwardReset(ts);
902:   return 0;
903: }

905: /*@
906:    TSAdjointSetUp - Sets up the internal data structures for the later use
907:    of an adjoint solver

909:    Collective

911:    Input Parameter:
912: .  ts - the `TS` context obtained from `TSCreate()`

914:    Level: advanced

916: .seealso: [](chapter_ts), `TSCreate()`, `TSAdjointStep()`, `TSSetCostGradients()`
917: @*/
918: PetscErrorCode TSAdjointSetUp(TS ts)
919: {
920:   TSTrajectory tj;
921:   PetscBool    match;

924:   if (ts->adjointsetupcalled) return 0;
927:   TSGetTrajectory(ts, &tj);
928:   PetscObjectTypeCompare((PetscObject)tj, TSTRAJECTORYBASIC, &match);
929:   if (match) {
930:     PetscBool solution_only;
931:     TSTrajectoryGetSolutionOnly(tj, &solution_only);
933:   }
934:   TSTrajectorySetUseHistory(tj, PETSC_FALSE); /* not use TSHistory */

936:   if (ts->quadraturets) { /* if there is integral in the cost function */
937:     VecDuplicate(ts->vecs_sensi[0], &ts->vec_drdu_col);
938:     if (ts->vecs_sensip) VecDuplicate(ts->vecs_sensip[0], &ts->vec_drdp_col);
939:   }

941:   PetscTryTypeMethod(ts, adjointsetup);
942:   ts->adjointsetupcalled = PETSC_TRUE;
943:   return 0;
944: }

946: /*@
947:   TSAdjointReset - Resets a `TS` adjoint context and removes any allocated `Vec`s and `Mat`s.

949:    Collective

951:    Input Parameter:
952: .  ts - the `TS` context obtained from `TSCreate()`

954:    Level: beginner

956: .seealso: [](chapter_ts), `TSCreate()`, `TSAdjointSetUp()`, `TSADestroy()`
957: @*/
958: PetscErrorCode TSAdjointReset(TS ts)
959: {
961:   PetscTryTypeMethod(ts, adjointreset);
962:   if (ts->quadraturets) { /* if there is integral in the cost function */
963:     VecDestroy(&ts->vec_drdu_col);
964:     if (ts->vecs_sensip) VecDestroy(&ts->vec_drdp_col);
965:   }
966:   ts->vecs_sensi         = NULL;
967:   ts->vecs_sensip        = NULL;
968:   ts->vecs_sensi2        = NULL;
969:   ts->vecs_sensi2p       = NULL;
970:   ts->vec_dir            = NULL;
971:   ts->adjointsetupcalled = PETSC_FALSE;
972:   return 0;
973: }

975: /*@
976:    TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time

978:    Logically Collective

980:    Input Parameters:
981: +  ts - the `TS` context obtained from `TSCreate()`
982: -  steps - number of steps to use

984:    Level: intermediate

986:    Notes:
987:     Normally one does not call this and `TSAdjointSolve()` integrates back to the original timestep. One can call this
988:           so as to integrate back to less than the original timestep

990: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TS`, `TSSetExactFinalTime()`
991: @*/
992: PetscErrorCode TSAdjointSetSteps(TS ts, PetscInt steps)
993: {
998:   ts->adjoint_max_steps = steps;
999:   return 0;
1000: }

1002: /*@C
1003:   TSAdjointSetRHSJacobian - Deprecated, use `TSSetRHSJacobianP()`

1005:   Level: deprecated

1007: @*/
1008: PetscErrorCode TSAdjointSetRHSJacobian(TS ts, Mat Amat, PetscErrorCode (*func)(TS, PetscReal, Vec, Mat, void *), void *ctx)
1009: {

1013:   ts->rhsjacobianp    = func;
1014:   ts->rhsjacobianpctx = ctx;
1015:   if (Amat) {
1016:     PetscObjectReference((PetscObject)Amat);
1017:     MatDestroy(&ts->Jacp);
1018:     ts->Jacp = Amat;
1019:   }
1020:   return 0;
1021: }

1023: /*@C
1024:   TSAdjointComputeRHSJacobian - Deprecated, use `TSComputeRHSJacobianP()`

1026:   Level: deprecated

1028: @*/
1029: PetscErrorCode TSAdjointComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat Amat)
1030: {

1035:   PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->rhsjacobianp)(ts, t, U, Amat, ts->rhsjacobianpctx));
1036:   return 0;
1037: }

1039: /*@
1040:   TSAdjointComputeDRDYFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobian()`

1042:   Level: deprecated

1044: @*/
1045: PetscErrorCode TSAdjointComputeDRDYFunction(TS ts, PetscReal t, Vec U, Vec *DRDU)
1046: {

1050:   PetscCallBack("TS callback DRDY for sensitivity analysis", (*ts->drdufunction)(ts, t, U, DRDU, ts->costintegrandctx));
1051:   return 0;
1052: }

1054: /*@
1055:   TSAdjointComputeDRDPFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobianP()`

1057:   Level: deprecated

1059: @*/
1060: PetscErrorCode TSAdjointComputeDRDPFunction(TS ts, PetscReal t, Vec U, Vec *DRDP)
1061: {

1065:   PetscCallBack("TS callback DRDP for sensitivity analysis", (*ts->drdpfunction)(ts, t, U, DRDP, ts->costintegrandctx));
1066:   return 0;
1067: }

1069: /*@C
1070:    TSAdjointMonitorSensi - monitors the first lambda sensitivity

1072:    Level: intermediate

1074: .seealso: [](chapter_ts), `TSAdjointMonitorSet()`
1075: @*/
1076: PetscErrorCode TSAdjointMonitorSensi(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscInt numcost, Vec *lambda, Vec *mu, PetscViewerAndFormat *vf)
1077: {
1078:   PetscViewer viewer = vf->viewer;

1081:   PetscViewerPushFormat(viewer, vf->format);
1082:   VecView(lambda[0], viewer);
1083:   PetscViewerPopFormat(viewer);
1084:   return 0;
1085: }

1087: /*@C
1088:    TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user

1090:    Collective

1092:    Input Parameters:
1093: +  ts - `TS` object you wish to monitor
1094: .  name - the monitor type one is seeking
1095: .  help - message indicating what monitoring is done
1096: .  manual - manual page for the monitor
1097: .  monitor - the monitor function
1098: -  monitorsetup - a function that is called once ONLY if the user selected this monitor that may set additional features of the `TS` or `PetscViewer` objects

1100:    Level: developer

1102: .seealso: [](chapter_ts), `PetscOptionsGetViewer()`, `PetscOptionsGetReal()`, `PetscOptionsHasName()`, `PetscOptionsGetString()`,
1103:           `PetscOptionsGetIntArray()`, `PetscOptionsGetRealArray()`, `PetscOptionsBool()`
1104:           `PetscOptionsInt()`, `PetscOptionsString()`, `PetscOptionsReal()`, `PetscOptionsBool()`,
1105:           `PetscOptionsName()`, `PetscOptionsBegin()`, `PetscOptionsEnd()`, `PetscOptionsHeadBegin()`,
1106:           `PetscOptionsStringArray()`, `PetscOptionsRealArray()`, `PetscOptionsScalar()`,
1107:           `PetscOptionsBoolGroupBegin()`, `PetscOptionsBoolGroup()`, `PetscOptionsBoolGroupEnd()`,
1108:           `PetscOptionsFList()`, `PetscOptionsEList()`
1109: @*/
1110: PetscErrorCode TSAdjointMonitorSetFromOptions(TS ts, const char name[], const char help[], const char manual[], PetscErrorCode (*monitor)(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, PetscViewerAndFormat *), PetscErrorCode (*monitorsetup)(TS, PetscViewerAndFormat *))
1111: {
1112:   PetscViewer       viewer;
1113:   PetscViewerFormat format;
1114:   PetscBool         flg;

1116:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, name, &viewer, &format, &flg);
1117:   if (flg) {
1118:     PetscViewerAndFormat *vf;
1119:     PetscViewerAndFormatCreate(viewer, format, &vf);
1120:     PetscObjectDereference((PetscObject)viewer);
1121:     if (monitorsetup) (*monitorsetup)(ts, vf);
1122:     TSAdjointMonitorSet(ts, (PetscErrorCode(*)(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, void *))monitor, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1123:   }
1124:   return 0;
1125: }

1127: /*@C
1128:    TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
1129:    timestep to display the iteration's  progress.

1131:    Logically Collective

1133:    Input Parameters:
1134: +  ts - the `TS` context obtained from `TSCreate()`
1135: .  adjointmonitor - monitoring routine
1136: .  adjointmctx - [optional] user-defined context for private data for the
1137:              monitor routine (use NULL if no context is desired)
1138: -  adjointmonitordestroy - [optional] routine that frees monitor context
1139:           (may be NULL)

1141:    Calling sequence of monitor:
1142: $    int adjointmonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *adjointmctx)

1144: +    ts - the `TS` context
1145: .    steps - iteration number (after the final time step the monitor routine is called with a step of -1, this is at the final time which may have
1146:                                been interpolated to)
1147: .    time - current time
1148: .    u - current iterate
1149: .    numcost - number of cost functionos
1150: .    lambda - sensitivities to initial conditions
1151: .    mu - sensitivities to parameters
1152: -    adjointmctx - [optional] adjoint monitoring context

1154:    Level: intermediate

1156:    Note:
1157:    This routine adds an additional monitor to the list of monitors that
1158:    already has been loaded.

1160:    Fortran Note:
1161:    Only a single monitor function can be set for each TS object

1163: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSAdjointMonitorCancel()`
1164: @*/
1165: PetscErrorCode TSAdjointMonitorSet(TS ts, PetscErrorCode (*adjointmonitor)(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, void *), void *adjointmctx, PetscErrorCode (*adjointmdestroy)(void **))
1166: {
1167:   PetscInt  i;
1168:   PetscBool identical;

1171:   for (i = 0; i < ts->numbermonitors; i++) {
1172:     PetscMonitorCompare((PetscErrorCode(*)(void))adjointmonitor, adjointmctx, adjointmdestroy, (PetscErrorCode(*)(void))ts->adjointmonitor[i], ts->adjointmonitorcontext[i], ts->adjointmonitordestroy[i], &identical);
1173:     if (identical) return 0;
1174:   }
1176:   ts->adjointmonitor[ts->numberadjointmonitors]          = adjointmonitor;
1177:   ts->adjointmonitordestroy[ts->numberadjointmonitors]   = adjointmdestroy;
1178:   ts->adjointmonitorcontext[ts->numberadjointmonitors++] = (void *)adjointmctx;
1179:   return 0;
1180: }

1182: /*@C
1183:    TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.

1185:    Logically Collective

1187:    Input Parameters:
1188: .  ts - the `TS` context obtained from `TSCreate()`

1190:    Notes:
1191:    There is no way to remove a single, specific monitor.

1193:    Level: intermediate

1195: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSAdjointMonitorSet()`
1196: @*/
1197: PetscErrorCode TSAdjointMonitorCancel(TS ts)
1198: {
1199:   PetscInt i;

1202:   for (i = 0; i < ts->numberadjointmonitors; i++) {
1203:     if (ts->adjointmonitordestroy[i]) (*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]);
1204:   }
1205:   ts->numberadjointmonitors = 0;
1206:   return 0;
1207: }

1209: /*@C
1210:    TSAdjointMonitorDefault - the default monitor of adjoint computations

1212:    Level: intermediate

1214: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSAdjointMonitorSet()`
1215: @*/
1216: PetscErrorCode TSAdjointMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscInt numcost, Vec *lambda, Vec *mu, PetscViewerAndFormat *vf)
1217: {
1218:   PetscViewer viewer = vf->viewer;

1221:   PetscViewerPushFormat(viewer, vf->format);
1222:   PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel);
1223:   PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " TS dt %g time %g%s", step, (double)ts->time_step, (double)ptime, ts->steprollback ? " (r)\n" : "\n");
1224:   PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel);
1225:   PetscViewerPopFormat(viewer);
1226:   return 0;
1227: }

1229: /*@C
1230:    TSAdjointMonitorDrawSensi - Monitors progress of the adjoint `TS` solvers by calling
1231:    `VecView()` for the sensitivities to initial states at each timestep

1233:    Collective

1235:    Input Parameters:
1236: +  ts - the `TS` context
1237: .  step - current time-step
1238: .  ptime - current time
1239: .  u - current state
1240: .  numcost - number of cost functions
1241: .  lambda - sensitivities to initial conditions
1242: .  mu - sensitivities to parameters
1243: -  dummy - either a viewer or NULL

1245:    Level: intermediate

1247: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSAdjointMonitorSet()`, `TSAdjointMonitorDefault()`, `VecView()`
1248: @*/
1249: PetscErrorCode TSAdjointMonitorDrawSensi(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscInt numcost, Vec *lambda, Vec *mu, void *dummy)
1250: {
1251:   TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
1252:   PetscDraw        draw;
1253:   PetscReal        xl, yl, xr, yr, h;
1254:   char             time[32];

1256:   if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return 0;

1258:   VecView(lambda[0], ictx->viewer);
1259:   PetscViewerDrawGetDraw(ictx->viewer, 0, &draw);
1260:   PetscSNPrintf(time, 32, "Timestep %d Time %g", (int)step, (double)ptime);
1261:   PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr);
1262:   h = yl + .95 * (yr - yl);
1263:   PetscDrawStringCentered(draw, .5 * (xl + xr), h, PETSC_DRAW_BLACK, time);
1264:   PetscDrawFlush(draw);
1265:   return 0;
1266: }

1268: /*
1269:    TSAdjointSetFromOptions - Sets various `TS` adjoint parameters from user options.

1271:    Collective

1273:    Input Parameter:
1274: .  ts - the `TS` context

1276:    Options Database Keys:
1277: +  -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
1278: .  -ts_adjoint_monitor - print information at each adjoint time step
1279: -  -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically

1281:    Level: developer

1283:    Note:
1284:     This is not normally called directly by users

1286: .seealso: [](chapter_ts), `TSSetSaveTrajectory()`, `TSTrajectorySetUp()`
1287: */
1288: PetscErrorCode TSAdjointSetFromOptions(TS ts, PetscOptionItems *PetscOptionsObject)
1289: {
1290:   PetscBool tflg, opt;

1293:   PetscOptionsHeadBegin(PetscOptionsObject, "TS Adjoint options");
1294:   tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
1295:   PetscOptionsBool("-ts_adjoint_solve", "Solve the adjoint problem immediately after solving the forward problem", "", tflg, &tflg, &opt);
1296:   if (opt) {
1297:     TSSetSaveTrajectory(ts);
1298:     ts->adjoint_solve = tflg;
1299:   }
1300:   TSAdjointMonitorSetFromOptions(ts, "-ts_adjoint_monitor", "Monitor adjoint timestep size", "TSAdjointMonitorDefault", TSAdjointMonitorDefault, NULL);
1301:   TSAdjointMonitorSetFromOptions(ts, "-ts_adjoint_monitor_sensi", "Monitor sensitivity in the adjoint computation", "TSAdjointMonitorSensi", TSAdjointMonitorSensi, NULL);
1302:   opt = PETSC_FALSE;
1303:   PetscOptionsName("-ts_adjoint_monitor_draw_sensi", "Monitor adjoint sensitivities (lambda only) graphically", "TSAdjointMonitorDrawSensi", &opt);
1304:   if (opt) {
1305:     TSMonitorDrawCtx ctx;
1306:     PetscInt         howoften = 1;

1308:     PetscOptionsInt("-ts_adjoint_monitor_draw_sensi", "Monitor adjoint sensitivities (lambda only) graphically", "TSAdjointMonitorDrawSensi", howoften, &howoften, NULL);
1309:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
1310:     TSAdjointMonitorSet(ts, TSAdjointMonitorDrawSensi, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
1311:   }
1312:   return 0;
1313: }

1315: /*@
1316:    TSAdjointStep - Steps one time step backward in the adjoint run

1318:    Collective

1320:    Input Parameter:
1321: .  ts - the `TS` context obtained from `TSCreate()`

1323:    Level: intermediate

1325: .seealso: [](chapter_ts), `TSAdjointSetUp()`, `TSAdjointSolve()`
1326: @*/
1327: PetscErrorCode TSAdjointStep(TS ts)
1328: {
1329:   DM dm;

1332:   TSGetDM(ts, &dm);
1333:   TSAdjointSetUp(ts);
1334:   ts->steps--; /* must decrease the step index before the adjoint step is taken. */

1336:   ts->reason     = TS_CONVERGED_ITERATING;
1337:   ts->ptime_prev = ts->ptime;
1338:   PetscLogEventBegin(TS_AdjointStep, ts, 0, 0, 0);
1339:   PetscUseTypeMethod(ts, adjointstep);
1340:   PetscLogEventEnd(TS_AdjointStep, ts, 0, 0, 0);
1341:   ts->adjoint_steps++;

1343:   if (ts->reason < 0) {
1345:   } else if (!ts->reason) {
1346:     if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1347:   }
1348:   return 0;
1349: }

1351: /*@
1352:    TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE

1354:    Collective
1355: `
1356:    Input Parameter:
1357: .  ts - the `TS` context obtained from `TSCreate()`

1359:    Options Database Key:
1360: . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial values

1362:    Level: intermediate

1364:    Notes:
1365:    This must be called after a call to `TSSolve()` that solves the forward problem

1367:    By default this will integrate back to the initial time, one can use `TSAdjointSetSteps()` to step back to a later time

1369: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSCreate()`, `TSSetCostGradients()`, `TSSetSolution()`, `TSAdjointStep()`
1370: @*/
1371: PetscErrorCode TSAdjointSolve(TS ts)
1372: {
1373:   static PetscBool cite = PETSC_FALSE;
1374: #if defined(TSADJOINT_STAGE)
1375:   PetscLogStage adjoint_stage;
1376: #endif

1379:   PetscCall(PetscCitationsRegister("@article{Zhang2022tsadjoint,\n"
1380:                                    "  title         = {{PETSc TSAdjoint: A Discrete Adjoint ODE Solver for First-Order and Second-Order Sensitivity Analysis}},\n"
1381:                                    "  author        = {Zhang, Hong and Constantinescu, Emil M.  and Smith, Barry F.},\n"
1382:                                    "  journal       = {SIAM Journal on Scientific Computing},\n"
1383:                                    "  volume        = {44},\n"
1384:                                    "  number        = {1},\n"
1385:                                    "  pages         = {C1-C24},\n"
1386:                                    "  doi           = {10.1137/21M140078X},\n"
1387:                                    "  year          = {2022}\n}\n",
1388:                                    &cite));
1389: #if defined(TSADJOINT_STAGE)
1390:   PetscLogStageRegister("TSAdjoint", &adjoint_stage);
1391:   PetscLogStagePush(adjoint_stage);
1392: #endif
1393:   TSAdjointSetUp(ts);

1395:   /* reset time step and iteration counters */
1396:   ts->adjoint_steps     = 0;
1397:   ts->ksp_its           = 0;
1398:   ts->snes_its          = 0;
1399:   ts->num_snes_failures = 0;
1400:   ts->reject            = 0;
1401:   ts->reason            = TS_CONVERGED_ITERATING;

1403:   if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->steps;
1404:   if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;

1406:   while (!ts->reason) {
1407:     TSTrajectoryGet(ts->trajectory, ts, ts->steps, &ts->ptime);
1408:     TSAdjointMonitor(ts, ts->steps, ts->ptime, ts->vec_sol, ts->numcost, ts->vecs_sensi, ts->vecs_sensip);
1409:     TSAdjointEventHandler(ts);
1410:     TSAdjointStep(ts);
1411:     if ((ts->vec_costintegral || ts->quadraturets) && !ts->costintegralfwd) TSAdjointCostIntegral(ts);
1412:   }
1413:   if (!ts->steps) {
1414:     TSTrajectoryGet(ts->trajectory, ts, ts->steps, &ts->ptime);
1415:     TSAdjointMonitor(ts, ts->steps, ts->ptime, ts->vec_sol, ts->numcost, ts->vecs_sensi, ts->vecs_sensip);
1416:   }
1417:   ts->solvetime = ts->ptime;
1418:   TSTrajectoryViewFromOptions(ts->trajectory, NULL, "-ts_trajectory_view");
1419:   VecViewFromOptions(ts->vecs_sensi[0], (PetscObject)ts, "-ts_adjoint_view_solution");
1420:   ts->adjoint_max_steps = 0;
1421: #if defined(TSADJOINT_STAGE)
1422:   PetscLogStagePop();
1423: #endif
1424:   return 0;
1425: }

1427: /*@C
1428:    TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using `TSAdjointMonitorSet()`

1430:    Collective

1432:    Input Parameters:
1433: +  ts - time stepping context obtained from `TSCreate()`
1434: .  step - step number that has just completed
1435: .  ptime - model time of the state
1436: .  u - state at the current model time
1437: .  numcost - number of cost functions (dimension of lambda  or mu)
1438: .  lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
1439: -  mu - vectors containing the gradients of the cost functions with respect to the problem parameters

1441:    Level: developer

1443:    Note:
1444:    `TSAdjointMonitor()` is typically used automatically within the time stepping implementations.
1445:    Users would almost never call this routine directly.

1447: .seealso: `TSAdjointMonitorSet()`, `TSAdjointSolve()`
1448: @*/
1449: PetscErrorCode TSAdjointMonitor(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscInt numcost, Vec *lambda, Vec *mu)
1450: {
1451:   PetscInt i, n = ts->numberadjointmonitors;

1455:   VecLockReadPush(u);
1456:   for (i = 0; i < n; i++) (*ts->adjointmonitor[i])(ts, step, ptime, u, numcost, lambda, mu, ts->adjointmonitorcontext[i]);
1457:   VecLockReadPop(u);
1458:   return 0;
1459: }

1461: /*@
1462:  TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.

1464:  Collective

1466:  Input Parameter:
1467:  .  ts - time stepping context

1469:  Level: advanced

1471:  Notes:
1472:  This function cannot be called until `TSAdjointStep()` has been completed.

1474:  .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSAdjointStep()`
1475:  @*/
1476: PetscErrorCode TSAdjointCostIntegral(TS ts)
1477: {
1479:   PetscUseTypeMethod(ts, adjointintegral);
1480:   return 0;
1481: }

1483: /* ------------------ Forward (tangent linear) sensitivity  ------------------*/

1485: /*@
1486:   TSForwardSetUp - Sets up the internal data structures for the later use
1487:   of forward sensitivity analysis

1489:   Collective

1491:   Input Parameter:
1492: . ts - the `TS` context obtained from `TSCreate()`

1494:   Level: advanced

1496: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSDestroy()`, `TSSetUp()`
1497: @*/
1498: PetscErrorCode TSForwardSetUp(TS ts)
1499: {
1501:   if (ts->forwardsetupcalled) return 0;
1502:   PetscTryTypeMethod(ts, forwardsetup);
1503:   VecDuplicate(ts->vec_sol, &ts->vec_sensip_col);
1504:   ts->forwardsetupcalled = PETSC_TRUE;
1505:   return 0;
1506: }

1508: /*@
1509:   TSForwardReset - Reset the internal data structures used by forward sensitivity analysis

1511:   Collective

1513:   Input Parameter:
1514: . ts - the `TS` context obtained from `TSCreate()`

1516:   Level: advanced

1518: .seealso: [](chapter_ts), `TSCreate()`, `TSDestroy()`, `TSForwardSetUp()`
1519: @*/
1520: PetscErrorCode TSForwardReset(TS ts)
1521: {
1522:   TS quadts = ts->quadraturets;

1525:   PetscTryTypeMethod(ts, forwardreset);
1526:   MatDestroy(&ts->mat_sensip);
1527:   if (quadts) MatDestroy(&quadts->mat_sensip);
1528:   VecDestroy(&ts->vec_sensip_col);
1529:   ts->forward_solve      = PETSC_FALSE;
1530:   ts->forwardsetupcalled = PETSC_FALSE;
1531:   return 0;
1532: }

1534: /*@
1535:   TSForwardSetIntegralGradients - Set the vectors holding forward sensitivities of the integral term.

1537:   Input Parameters:
1538: + ts - the `TS` context obtained from `TSCreate()`
1539: . numfwdint - number of integrals
1540: - vp - the vectors containing the gradients for each integral w.r.t. parameters

1542:   Level: deprecated

1544: .seealso: [](chapter_ts), `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1545: @*/
1546: PetscErrorCode TSForwardSetIntegralGradients(TS ts, PetscInt numfwdint, Vec *vp)
1547: {
1550:   if (!ts->numcost) ts->numcost = numfwdint;

1552:   ts->vecs_integral_sensip = vp;
1553:   return 0;
1554: }

1556: /*@
1557:   TSForwardGetIntegralGradients - Returns the forward sensitivities ofthe integral term.

1559:   Input Parameter:
1560: . ts - the `TS` context obtained from `TSCreate()`

1562:   Output Parameter:
1563: . vp - the vectors containing the gradients for each integral w.r.t. parameters

1565:   Level: deprecated

1567: .seealso: [](chapter_ts), `TSForwardSetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1568: @*/
1569: PetscErrorCode TSForwardGetIntegralGradients(TS ts, PetscInt *numfwdint, Vec **vp)
1570: {
1573:   if (numfwdint) *numfwdint = ts->numcost;
1574:   if (vp) *vp = ts->vecs_integral_sensip;
1575:   return 0;
1576: }

1578: /*@
1579:   TSForwardStep - Compute the forward sensitivity for one time step.

1581:   Collective

1583:   Input Parameter:
1584: . ts - time stepping context

1586:   Level: advanced

1588:   Notes:
1589:   This function cannot be called until `TSStep()` has been completed.

1591: .seealso: [](chapter_ts), `TSForwardSetSensitivities()`, `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardSetUp()`
1592: @*/
1593: PetscErrorCode TSForwardStep(TS ts)
1594: {
1596:   PetscLogEventBegin(TS_ForwardStep, ts, 0, 0, 0);
1597:   PetscUseTypeMethod(ts, forwardstep);
1598:   PetscLogEventEnd(TS_ForwardStep, ts, 0, 0, 0);
1600:   return 0;
1601: }

1603: /*@
1604:   TSForwardSetSensitivities - Sets the initial value of the trajectory sensitivities of solution  w.r.t. the problem parameters and initial values.

1606:   Logically Collective

1608:   Input Parameters:
1609: + ts - the `TS` context obtained from `TSCreate()`
1610: . nump - number of parameters
1611: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters

1613:   Level: beginner

1615:   Notes:
1616:   Forward sensitivity is also called 'trajectory sensitivity' in some fields such as power systems.
1617:   This function turns on a flag to trigger `TSSolve()` to compute forward sensitivities automatically.
1618:   You must call this function before `TSSolve()`.
1619:   The entries in the sensitivity matrix must be correctly initialized with the values S = dy/dp|startingtime.

1621: .seealso: [](chapter_ts), `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1622: @*/
1623: PetscErrorCode TSForwardSetSensitivities(TS ts, PetscInt nump, Mat Smat)
1624: {
1627:   ts->forward_solve = PETSC_TRUE;
1628:   if (nump == PETSC_DEFAULT) {
1629:     MatGetSize(Smat, NULL, &ts->num_parameters);
1630:   } else ts->num_parameters = nump;
1631:   PetscObjectReference((PetscObject)Smat);
1632:   MatDestroy(&ts->mat_sensip);
1633:   ts->mat_sensip = Smat;
1634:   return 0;
1635: }

1637: /*@
1638:   TSForwardGetSensitivities - Returns the trajectory sensitivities

1640:   Not Collective, but Smat returned is parallel if ts is parallel

1642:   Output Parameters:
1643: + ts - the `TS` context obtained from `TSCreate()`
1644: . nump - number of parameters
1645: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters

1647:   Level: intermediate

1649: .seealso: [](chapter_ts), `TSForwardSetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1650: @*/
1651: PetscErrorCode TSForwardGetSensitivities(TS ts, PetscInt *nump, Mat *Smat)
1652: {
1654:   if (nump) *nump = ts->num_parameters;
1655:   if (Smat) *Smat = ts->mat_sensip;
1656:   return 0;
1657: }

1659: /*@
1660:    TSForwardCostIntegral - Evaluate the cost integral in the forward run.

1662:    Collective

1664:    Input Parameter:
1665: .  ts - time stepping context

1667:    Level: advanced

1669:    Note:
1670:    This function cannot be called until `TSStep()` has been completed.

1672: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSAdjointCostIntegral()`
1673: @*/
1674: PetscErrorCode TSForwardCostIntegral(TS ts)
1675: {
1677:   PetscUseTypeMethod(ts, forwardintegral);
1678:   return 0;
1679: }

1681: /*@
1682:   TSForwardSetInitialSensitivities - Set initial values for tangent linear sensitivities

1684:   Collective

1686:   Input Parameters:
1687: + ts - the `TS` context obtained from `TSCreate()`
1688: - didp - parametric sensitivities of the initial condition

1690:   Level: intermediate

1692:   Notes:
1693:   `TSSolve()` allows users to pass the initial solution directly to `TS`. But the tangent linear variables cannot be initialized in this way.
1694:    This function is used to set initial values for tangent linear variables.

1696: .seealso: [](chapter_ts), `TS`, `TSForwardSetSensitivities()`
1697: @*/
1698: PetscErrorCode TSForwardSetInitialSensitivities(TS ts, Mat didp)
1699: {
1702:   if (!ts->mat_sensip) TSForwardSetSensitivities(ts, PETSC_DEFAULT, didp);
1703:   return 0;
1704: }

1706: /*@
1707:    TSForwardGetStages - Get the number of stages and the tangent linear sensitivities at the intermediate stages

1709:    Input Parameter:
1710: .  ts - the `TS` context obtained from `TSCreate()`

1712:    Output Parameters:
1713: +  ns - number of stages
1714: -  S - tangent linear sensitivities at the intermediate stages

1716:    Level: advanced

1718: .seealso: `TS`
1719: @*/
1720: PetscErrorCode TSForwardGetStages(TS ts, PetscInt *ns, Mat **S)
1721: {

1724:   if (!ts->ops->getstages) *S = NULL;
1725:   else PetscUseTypeMethod(ts, forwardgetstages, ns, S);
1726:   return 0;
1727: }

1729: /*@
1730:    TSCreateQuadratureTS - Create a sub-`TS` that evaluates integrals over time

1732:    Input Parameters:
1733: +  ts - the `TS` context obtained from `TSCreate()`
1734: -  fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run

1736:    Output Parameters:
1737: .  quadts - the child `TS` context

1739:    Level: intermediate

1741: .seealso: [](chapter_ts), `TSGetQuadratureTS()`
1742: @*/
1743: PetscErrorCode TSCreateQuadratureTS(TS ts, PetscBool fwd, TS *quadts)
1744: {
1745:   char prefix[128];

1749:   TSDestroy(&ts->quadraturets);
1750:   TSCreate(PetscObjectComm((PetscObject)ts), &ts->quadraturets);
1751:   PetscObjectIncrementTabLevel((PetscObject)ts->quadraturets, (PetscObject)ts, 1);
1752:   PetscSNPrintf(prefix, sizeof(prefix), "%squad_", ((PetscObject)ts)->prefix ? ((PetscObject)ts)->prefix : "");
1753:   TSSetOptionsPrefix(ts->quadraturets, prefix);
1754:   *quadts = ts->quadraturets;

1756:   if (ts->numcost) {
1757:     VecCreateSeq(PETSC_COMM_SELF, ts->numcost, &(*quadts)->vec_sol);
1758:   } else {
1759:     VecCreateSeq(PETSC_COMM_SELF, 1, &(*quadts)->vec_sol);
1760:   }
1761:   ts->costintegralfwd = fwd;
1762:   return 0;
1763: }

1765: /*@
1766:    TSGetQuadratureTS - Return the sub-`TS` that evaluates integrals over time

1768:    Input Parameter:
1769: .  ts - the `TS` context obtained from `TSCreate()`

1771:    Output Parameters:
1772: +  fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1773: -  quadts - the child `TS` context

1775:    Level: intermediate

1777: .seealso: [](chapter_ts), `TSCreateQuadratureTS()`
1778: @*/
1779: PetscErrorCode TSGetQuadratureTS(TS ts, PetscBool *fwd, TS *quadts)
1780: {
1782:   if (fwd) *fwd = ts->costintegralfwd;
1783:   if (quadts) *quadts = ts->quadraturets;
1784:   return 0;
1785: }

1787: /*@
1788:    TSComputeSNESJacobian - Compute the Jacobian needed for the `SNESSolve()` in `TS`

1790:    Collective

1792:    Input Parameters:
1793: +  ts - the `TS` context obtained from `TSCreate()`
1794: -  x - state vector

1796:    Output Parameters:
1797: +  J - Jacobian matrix
1798: -  Jpre - preconditioning matrix for J (may be same as J)

1800:    Level: developer

1802:    Note:
1803:    Uses finite differencing when `TS` Jacobian is not available.

1805: .seealso: `SNES`, `TS`, `SNESSetJacobian()`, TSSetRHSJacobian()`, `TSSetIJacobian()`
1806: @*/
1807: PetscErrorCode TSComputeSNESJacobian(TS ts, Vec x, Mat J, Mat Jpre)
1808: {
1809:   SNES snes                                          = ts->snes;
1810:   PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *) = NULL;

1812:   /*
1813:     Unlike implicit methods, explicit methods do not have SNESMatFDColoring in the snes object
1814:     because SNESSolve() has not been called yet; so querying SNESMatFDColoring does not work for
1815:     explicit methods. Instead, we check the Jacobian compute function directly to determine if FD
1816:     coloring is used.
1817:   */
1818:   SNESGetJacobian(snes, NULL, NULL, &jac, NULL);
1819:   if (jac == SNESComputeJacobianDefaultColor) {
1820:     Vec f;
1821:     SNESSetSolution(snes, x);
1822:     SNESGetFunction(snes, &f, NULL, NULL);
1823:     /* Force MatFDColoringApply to evaluate the SNES residual function for the base vector */
1824:     SNESComputeFunction(snes, x, f);
1825:   }
1826:   SNESComputeJacobian(snes, x, J, Jpre);
1827:   return 0;
1828: }