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] function context
21: Level: intermediate
23: Note:
24: `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`
26: .seealso: [](ch_ts), `TS`, `TSRHSJacobianPFn`, `TSGetRHSJacobianP()`
27: @*/
28: PetscErrorCode TSSetRHSJacobianP(TS ts, Mat Amat, TSRHSJacobianPFn *func, PetscCtx ctx)
29: {
30: PetscFunctionBegin;
34: ts->rhsjacobianp = func;
35: ts->rhsjacobianpctx = ctx;
36: if (Amat) {
37: PetscCall(PetscObjectReference((PetscObject)Amat));
38: PetscCall(MatDestroy(&ts->Jacprhs));
39: ts->Jacprhs = Amat;
40: }
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: /*@C
45: 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.
47: Logically Collective
49: Input Parameter:
50: . ts - `TS` context obtained from `TSCreate()`
52: Output Parameters:
53: + Amat - JacobianP matrix
54: . func - function
55: - ctx - [optional] function context
57: Level: intermediate
59: Note:
60: `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`
62: .seealso: [](ch_ts), `TSSetRHSJacobianP()`, `TS`, `TSRHSJacobianPFn`
63: @*/
64: PetscErrorCode TSGetRHSJacobianP(TS ts, Mat *Amat, TSRHSJacobianPFn **func, PetscCtxRt ctx)
65: {
66: PetscFunctionBegin;
67: if (func) *func = ts->rhsjacobianp;
68: if (ctx) *(void **)ctx = ts->rhsjacobianpctx;
69: if (Amat) *Amat = ts->Jacprhs;
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*@C
74: TSComputeRHSJacobianP - Runs the user-defined JacobianP function.
76: Collective
78: Input Parameters:
79: + ts - The `TS` context obtained from `TSCreate()`
80: . t - the time
81: - U - the solution at which to compute the Jacobian
83: Output Parameter:
84: . Amat - the computed Jacobian
86: Level: developer
88: .seealso: [](ch_ts), `TSSetRHSJacobianP()`, `TS`
89: @*/
90: PetscErrorCode TSComputeRHSJacobianP(TS ts, PetscReal t, Vec U, Mat Amat)
91: {
92: PetscFunctionBegin;
93: if (!Amat) PetscFunctionReturn(PETSC_SUCCESS);
97: if (ts->rhsjacobianp) PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->rhsjacobianp)(ts, t, U, Amat, ts->rhsjacobianpctx));
98: else {
99: PetscBool assembled;
100: PetscCall(MatZeroEntries(Amat));
101: PetscCall(MatAssembled(Amat, &assembled));
102: if (!assembled) {
103: PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY));
104: PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY));
105: }
106: }
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: /*@C
111: TSSetIJacobianP - Sets the function that computes the Jacobian of $F$ w.r.t. the parameters $p$ where $F(Udot,U,p,t) = G(U,p,t)$, as well as the location to store the matrix.
113: Logically Collective
115: Input Parameters:
116: + ts - `TS` context obtained from `TSCreate()`
117: . Amat - JacobianP matrix
118: . func - function
119: - ctx - [optional] function context
121: Calling sequence of `func`:
122: + ts - the `TS` context
123: . t - current timestep
124: . U - input vector (current ODE solution)
125: . Udot - time derivative of state vector
126: . shift - shift to apply, see the note in `TSSetIJacobian()`
127: . A - output matrix
128: - ctx - [optional] function context
130: Level: intermediate
132: Note:
133: `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`
135: .seealso: [](ch_ts), `TSSetRHSJacobianP()`, `TS`
136: @*/
137: PetscErrorCode TSSetIJacobianP(TS ts, Mat Amat, PetscErrorCode (*func)(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, PetscCtx ctx), PetscCtx ctx)
138: {
139: PetscFunctionBegin;
143: ts->ijacobianp = func;
144: ts->ijacobianpctx = ctx;
145: if (Amat) {
146: PetscCall(PetscObjectReference((PetscObject)Amat));
147: PetscCall(MatDestroy(&ts->Jacp));
148: ts->Jacp = Amat;
149: }
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: /*@C
154: TSGetIJacobianP - Gets the function that computes the Jacobian of $ F$ w.r.t. the parameters $p$ where $F(Udot,U,p,t) = G(U,p,t) $, as well as the location to store the matrix.
156: Logically Collective
158: Input Parameter:
159: . ts - `TS` context obtained from `TSCreate()`
161: Output Parameters:
162: + Amat - JacobianP matrix
163: . func - the function that computes the JacobianP
164: - ctx - [optional] function context
166: Calling sequence of `func`:
167: + ts - the `TS` context
168: . t - current timestep
169: . U - input vector (current ODE solution)
170: . Udot - time derivative of state vector
171: . shift - shift to apply, see the note in `TSSetIJacobian()`
172: . A - output matrix
173: - ctx - [optional] function context
175: Level: intermediate
177: Note:
178: `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`
180: .seealso: [](ch_ts), `TSSetRHSJacobianP()`, `TS`, `TSSetIJacobianP()`, `TSGetRHSJacobianP()`
181: @*/
182: PetscErrorCode TSGetIJacobianP(TS ts, Mat *Amat, PetscErrorCode (**func)(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, PetscCtx ctx), PetscCtxRt ctx)
183: {
184: PetscFunctionBegin;
187: if (func) *func = ts->ijacobianp;
188: if (ctx) *(void **)ctx = ts->ijacobianpctx;
189: if (Amat) *Amat = ts->Jacp;
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: /*@
194: TSComputeIJacobianP - Runs the user-defined IJacobianP function.
196: Collective
198: Input Parameters:
199: + ts - the `TS` context
200: . t - current timestep
201: . U - state vector
202: . Udot - time derivative of state vector
203: . shift - shift to apply, see note below
204: - imex - flag indicates if the method is IMEX so that the `RHSJacobianP` should be kept separate
206: Output Parameter:
207: . Amat - Jacobian matrix
209: Level: developer
211: .seealso: [](ch_ts), `TS`, `TSSetIJacobianP()`
212: @*/
213: PetscErrorCode TSComputeIJacobianP(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat Amat, PetscBool imex)
214: {
215: PetscFunctionBegin;
216: if (!Amat) PetscFunctionReturn(PETSC_SUCCESS);
221: PetscCall(PetscLogEventBegin(TS_JacobianPEval, ts, U, Amat, 0));
222: if (ts->ijacobianp) PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->ijacobianp)(ts, t, U, Udot, shift, Amat, ts->ijacobianpctx));
223: if (imex) {
224: if (!ts->ijacobianp) { /* system was written as Udot = G(t,U) */
225: PetscBool assembled;
226: PetscCall(MatZeroEntries(Amat));
227: PetscCall(MatAssembled(Amat, &assembled));
228: if (!assembled) {
229: PetscCall(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY));
230: PetscCall(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY));
231: }
232: }
233: } else {
234: if (ts->rhsjacobianp) PetscCall(TSComputeRHSJacobianP(ts, t, U, ts->Jacprhs));
235: if (ts->Jacprhs == Amat) { /* No IJacobian, so we only have the RHS matrix */
236: PetscCall(MatScale(Amat, -1));
237: } else if (ts->Jacprhs) { /* Both IJacobian and RHSJacobian */
238: MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
239: if (!ts->ijacobianp) { /* No IJacobianp provided, but we have a separate RHS matrix */
240: PetscCall(MatZeroEntries(Amat));
241: }
242: PetscCall(MatAXPY(Amat, -1, ts->Jacprhs, axpy));
243: }
244: }
245: PetscCall(PetscLogEventEnd(TS_JacobianPEval, ts, U, Amat, 0));
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: /*@C
250: TSSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions
252: Logically Collective
254: Input Parameters:
255: + ts - the `TS` context obtained from `TSCreate()`
256: . numcost - number of gradients to be computed, this is the number of cost functions
257: . costintegral - vector that stores the integral values
258: . rf - routine for evaluating the integrand function
259: . drduf - function that computes the gradients of the `r` with respect to `u`
260: . drdpf - function that computes the gradients of the `r` with respect to p, can be `NULL` if parametric sensitivity is not desired (`mu` = `NULL`)
261: . fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
262: - ctx - [optional] application context for private data for the function evaluation routine (may be `NULL`)
264: Calling sequence of `rf`:
265: + ts - the integrator
266: . t - the time
267: . U - the solution
268: . F - the computed value of the function
269: - ctx - the application context
271: Calling sequence of `drduf`:
272: + ts - the integrator
273: . t - the time
274: . U - the solution
275: . dRdU - the computed gradients of the `r` with respect to `u`
276: - ctx - the application context
278: Calling sequence of `drdpf`:
279: + ts - the integrator
280: . t - the time
281: . U - the solution
282: . dRdP - the computed gradients of the `r` with respect to `p`
283: - ctx - the application context
285: Level: deprecated
287: Notes:
288: For optimization there is usually a single cost function (numcost = 1). For sensitivities there may be multiple cost functions
290: Use `TSCreateQuadratureTS()` and `TSForwardSetSensitivities()` instead
292: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobianP()`, `TSGetCostGradients()`, `TSSetCostGradients()`,
293: `TSCreateQuadratureTS()`, `TSForwardSetSensitivities()`
294: @*/
295: PetscErrorCode TSSetCostIntegrand(TS ts, PetscInt numcost, Vec costintegral, PetscErrorCode (*rf)(TS ts, PetscReal t, Vec U, Vec F, PetscCtx ctx), PetscErrorCode (*drduf)(TS ts, PetscReal t, Vec U, Vec *dRdU, PetscCtx ctx), PetscErrorCode (*drdpf)(TS ts, PetscReal t, Vec U, Vec *dRdP, PetscCtx ctx), PetscBool fwd, PetscCtx ctx)
296: {
297: PetscFunctionBegin;
300: PetscCheck(!ts->numcost || ts->numcost == numcost, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "The number of cost functions (2nd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostGradients() or TSForwardSetIntegralGradients()");
301: if (!ts->numcost) ts->numcost = numcost;
303: if (costintegral) {
304: PetscCall(PetscObjectReference((PetscObject)costintegral));
305: PetscCall(VecDestroy(&ts->vec_costintegral));
306: ts->vec_costintegral = costintegral;
307: } else {
308: if (!ts->vec_costintegral) { /* Create a seq vec if user does not provide one */
309: PetscCall(VecCreateSeq(PETSC_COMM_SELF, numcost, &ts->vec_costintegral));
310: } else {
311: PetscCall(VecSet(ts->vec_costintegral, 0.0));
312: }
313: }
314: if (!ts->vec_costintegrand) {
315: PetscCall(VecDuplicate(ts->vec_costintegral, &ts->vec_costintegrand));
316: } else {
317: PetscCall(VecSet(ts->vec_costintegrand, 0.0));
318: }
319: ts->costintegralfwd = fwd; /* Evaluate the cost integral in forward run if fwd is true */
320: ts->costintegrand = rf;
321: ts->costintegrandctx = ctx;
322: ts->drdufunction = drduf;
323: ts->drdpfunction = drdpf;
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /*@
328: TSGetCostIntegral - Returns the values of the integral term in the cost functions.
329: It is valid to call the routine after a backward run.
331: Not Collective
333: Input Parameter:
334: . ts - the `TS` context obtained from `TSCreate()`
336: Output Parameter:
337: . v - the vector containing the integrals for each cost function
339: Level: intermediate
341: .seealso: [](ch_ts), `TS`, `TSAdjointSolve()`, `TSSetCostIntegrand()`
342: @*/
343: PetscErrorCode TSGetCostIntegral(TS ts, Vec *v)
344: {
345: TS quadts;
347: PetscFunctionBegin;
349: PetscAssertPointer(v, 2);
350: PetscCall(TSGetQuadratureTS(ts, NULL, &quadts));
351: *v = quadts->vec_sol;
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@
356: TSComputeCostIntegrand - Evaluates the integral function in the cost functions.
358: Input Parameters:
359: + ts - the `TS` context
360: . t - current time
361: - U - state vector, i.e. current solution
363: Output Parameter:
364: . Q - vector of size numcost to hold the outputs
366: Level: deprecated
368: Note:
369: Most users should not need to explicitly call this routine, as it
370: is used internally within the sensitivity analysis context.
372: .seealso: [](ch_ts), `TS`, `TSAdjointSolve()`, `TSSetCostIntegrand()`
373: @*/
374: PetscErrorCode TSComputeCostIntegrand(TS ts, PetscReal t, Vec U, Vec Q)
375: {
376: PetscFunctionBegin;
381: PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, Q, 0));
382: if (ts->costintegrand) PetscCallBack("TS callback integrand in the cost function", (*ts->costintegrand)(ts, t, U, Q, ts->costintegrandctx));
383: else PetscCall(VecZeroEntries(Q));
384: PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, Q, 0));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: // PetscClangLinter pragma disable: -fdoc-*
389: /*@
390: TSComputeDRDUFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobian()`
392: Level: deprecated
394: @*/
395: PetscErrorCode TSComputeDRDUFunction(TS ts, PetscReal t, Vec U, Vec *DRDU)
396: {
397: PetscFunctionBegin;
398: if (!DRDU) PetscFunctionReturn(PETSC_SUCCESS);
402: PetscCallBack("TS callback DRDU for sensitivity analysis", (*ts->drdufunction)(ts, t, U, DRDU, ts->costintegrandctx));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: // PetscClangLinter pragma disable: -fdoc-*
407: /*@
408: TSComputeDRDPFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobianP()`
410: Level: deprecated
412: @*/
413: PetscErrorCode TSComputeDRDPFunction(TS ts, PetscReal t, Vec U, Vec *DRDP)
414: {
415: PetscFunctionBegin;
416: if (!DRDP) PetscFunctionReturn(PETSC_SUCCESS);
420: PetscCallBack("TS callback DRDP for sensitivity analysis", (*ts->drdpfunction)(ts, t, U, DRDP, ts->costintegrandctx));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
425: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
426: /*@C
427: 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.
429: Logically Collective
431: Input Parameters:
432: + ts - `TS` context obtained from `TSCreate()`
433: . ihp1 - an array of vectors storing the result of vector-Hessian-vector product for $F_{UU}$
434: . ihessianproductfunc1 - vector-Hessian-vector product function for $F_{UU}$
435: . ihp2 - an array of vectors storing the result of vector-Hessian-vector product for $F_{UP}$
436: . ihessianproductfunc2 - vector-Hessian-vector product function for $F_{UP}$
437: . ihp3 - an array of vectors storing the result of vector-Hessian-vector product for $F_{PU}$
438: . ihessianproductfunc3 - vector-Hessian-vector product function for $F_{PU}$
439: . ihp4 - an array of vectors storing the result of vector-Hessian-vector product for $F_{PP}$
440: . ihessianproductfunc4 - vector-Hessian-vector product function for $F_{PP}$
441: - ctx - [optional] function context
443: Calling sequence of `ihessianproductfunc1`:
444: + ts - the `TS` context
445: . t - current timestep
446: . U - input vector (current ODE solution)
447: . Vl - an array of input vectors to be left-multiplied with the Hessian
448: . Vr - input vector to be right-multiplied with the Hessian
449: . VHV - an array of output vectors for vector-Hessian-vector product
450: - ctx - [optional] function context
452: Level: intermediate
454: Notes:
455: All other functions have the same calling sequence as `ihessianproductfunc1`, so their
456: descriptions are omitted for brevity.
458: The first Hessian function and the working array are required.
459: As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
460: $Vl_n^T*F_UP*Vr$
461: 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.$
462: Each entry of $F_{UP}$ corresponds to the derivative
463: $ F_UP[i][j][k] = \frac{\partial^2 F[i]}{\partial U[j] \partial P[k]}.$
464: 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
465: $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * F_UP[i][j][k] * Vr[k]}$
466: If the cost function is a scalar, there will be only one vector in `Vl` and `VHV`.
468: .seealso: [](ch_ts), `TS`
469: @*/
470: PetscErrorCode TSSetIHessianProduct(TS ts, Vec ihp1[], PetscErrorCode (*ihessianproductfunc1)(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[], PetscCtx ctx), Vec ihp2[], PetscErrorCode (*ihessianproductfunc2)(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[], PetscCtx ctx), Vec ihp3[], PetscErrorCode (*ihessianproductfunc3)(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[], PetscCtx ctx), Vec ihp4[], PetscErrorCode (*ihessianproductfunc4)(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[], PetscCtx ctx), PetscCtx ctx)
471: {
472: PetscFunctionBegin;
474: PetscAssertPointer(ihp1, 2);
476: ts->ihessianproductctx = ctx;
477: if (ihp1) ts->vecs_fuu = ihp1;
478: if (ihp2) ts->vecs_fup = ihp2;
479: if (ihp3) ts->vecs_fpu = ihp3;
480: if (ihp4) ts->vecs_fpp = ihp4;
481: ts->ihessianproduct_fuu = ihessianproductfunc1;
482: ts->ihessianproduct_fup = ihessianproductfunc2;
483: ts->ihessianproduct_fpu = ihessianproductfunc3;
484: ts->ihessianproduct_fpp = ihessianproductfunc4;
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: /*@
489: TSComputeIHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Fuu.
491: Collective
493: Input Parameters:
494: + ts - The `TS` context obtained from `TSCreate()`
495: . t - the time
496: . U - the solution at which to compute the Hessian product
497: . Vl - the array of input vectors to be multiplied with the Hessian from the left
498: - Vr - the input vector to be multiplied with the Hessian from the right
500: Output Parameter:
501: . VHV - the array of output vectors that store the Hessian product
503: Level: developer
505: Note:
506: `TSComputeIHessianProductFunctionUU()` is typically used for sensitivity implementation,
507: so most users would not generally call this routine themselves.
509: .seealso: [](ch_ts), `TSSetIHessianProduct()`
510: @*/
511: PetscErrorCode TSComputeIHessianProductFunctionUU(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[])
512: {
513: PetscFunctionBegin;
514: if (!VHV) PetscFunctionReturn(PETSC_SUCCESS);
518: if (ts->ihessianproduct_fuu) PetscCallBack("TS callback IHessianProduct 1 for sensitivity analysis", (*ts->ihessianproduct_fuu)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
520: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
521: if (ts->rhshessianproduct_guu) {
522: PetscInt nadj;
523: PetscCall(TSComputeRHSHessianProductFunctionUU(ts, t, U, Vl, Vr, VHV));
524: for (nadj = 0; nadj < ts->numcost; nadj++) PetscCall(VecScale(VHV[nadj], -1));
525: }
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /*@
530: TSComputeIHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Fup.
532: Collective
534: Input Parameters:
535: + ts - The `TS` context obtained from `TSCreate()`
536: . t - the time
537: . U - the solution at which to compute the Hessian product
538: . Vl - the array of input vectors to be multiplied with the Hessian from the left
539: - Vr - the input vector to be multiplied with the Hessian from the right
541: Output Parameter:
542: . VHV - the array of output vectors that store the Hessian product
544: Level: developer
546: Note:
547: `TSComputeIHessianProductFunctionUP()` is typically used for sensitivity implementation,
548: so most users would not generally call this routine themselves.
550: .seealso: [](ch_ts), `TSSetIHessianProduct()`
551: @*/
552: PetscErrorCode TSComputeIHessianProductFunctionUP(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[])
553: {
554: PetscFunctionBegin;
555: if (!VHV) PetscFunctionReturn(PETSC_SUCCESS);
559: if (ts->ihessianproduct_fup) PetscCallBack("TS callback IHessianProduct 2 for sensitivity analysis", (*ts->ihessianproduct_fup)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
561: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
562: if (ts->rhshessianproduct_gup) {
563: PetscInt nadj;
564: PetscCall(TSComputeRHSHessianProductFunctionUP(ts, t, U, Vl, Vr, VHV));
565: for (nadj = 0; nadj < ts->numcost; nadj++) PetscCall(VecScale(VHV[nadj], -1));
566: }
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: /*@
571: TSComputeIHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Fpu.
573: Collective
575: Input Parameters:
576: + ts - The `TS` context obtained from `TSCreate()`
577: . t - the time
578: . U - the solution at which to compute the Hessian product
579: . Vl - the array of input vectors to be multiplied with the Hessian from the left
580: - Vr - the input vector to be multiplied with the Hessian from the right
582: Output Parameter:
583: . VHV - the array of output vectors that store the Hessian product
585: Level: developer
587: Note:
588: `TSComputeIHessianProductFunctionPU()` is typically used for sensitivity implementation,
589: so most users would not generally call this routine themselves.
591: .seealso: [](ch_ts), `TSSetIHessianProduct()`
592: @*/
593: PetscErrorCode TSComputeIHessianProductFunctionPU(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[])
594: {
595: PetscFunctionBegin;
596: if (!VHV) PetscFunctionReturn(PETSC_SUCCESS);
600: if (ts->ihessianproduct_fpu) PetscCallBack("TS callback IHessianProduct 3 for sensitivity analysis", (*ts->ihessianproduct_fpu)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
602: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
603: if (ts->rhshessianproduct_gpu) {
604: PetscInt nadj;
605: PetscCall(TSComputeRHSHessianProductFunctionPU(ts, t, U, Vl, Vr, VHV));
606: for (nadj = 0; nadj < ts->numcost; nadj++) PetscCall(VecScale(VHV[nadj], -1));
607: }
608: PetscFunctionReturn(PETSC_SUCCESS);
609: }
611: /*@C
612: TSComputeIHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Fpp.
614: Collective
616: Input Parameters:
617: + ts - The `TS` context obtained from `TSCreate()`
618: . t - the time
619: . U - the solution at which to compute the Hessian product
620: . Vl - the array of input vectors to be multiplied with the Hessian from the left
621: - Vr - the input vector to be multiplied with the Hessian from the right
623: Output Parameter:
624: . VHV - the array of output vectors that store the Hessian product
626: Level: developer
628: Note:
629: `TSComputeIHessianProductFunctionPP()` is typically used for sensitivity implementation,
630: so most users would not generally call this routine themselves.
632: .seealso: [](ch_ts), `TSSetIHessianProduct()`
633: @*/
634: PetscErrorCode TSComputeIHessianProductFunctionPP(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[])
635: {
636: PetscFunctionBegin;
637: if (!VHV) PetscFunctionReturn(PETSC_SUCCESS);
641: if (ts->ihessianproduct_fpp) PetscCallBack("TS callback IHessianProduct 3 for sensitivity analysis", (*ts->ihessianproduct_fpp)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
643: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
644: if (ts->rhshessianproduct_gpp) {
645: PetscInt nadj;
646: PetscCall(TSComputeRHSHessianProductFunctionPP(ts, t, U, Vl, Vr, VHV));
647: for (nadj = 0; nadj < ts->numcost; nadj++) PetscCall(VecScale(VHV[nadj], -1));
648: }
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
653: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
654: /*@C
655: TSSetRHSHessianProduct - Sets the function that computes the vector-Hessian-vector
656: product. The Hessian is the second-order derivative of `G` (RHSFunction) w.r.t. the state
657: variable.
659: Logically Collective
661: Input Parameters:
662: + ts - `TS` context obtained from `TSCreate()`
663: . rhshp1 - an array of vectors storing the result of vector-Hessian-vector product for $G_{UU}$
664: . rhshessianproductfunc1 - vector-Hessian-vector product function for $G_{UU}$
665: . rhshp2 - an array of vectors storing the result of vector-Hessian-vector product for $G_{UP}$
666: . rhshessianproductfunc2 - vector-Hessian-vector product function for $G_{UP}$
667: . rhshp3 - an array of vectors storing the result of vector-Hessian-vector product for $G_{PU}$
668: . rhshessianproductfunc3 - vector-Hessian-vector product function for $G_{PU}$
669: . rhshp4 - an array of vectors storing the result of vector-Hessian-vector product for $G_{PP}$
670: . rhshessianproductfunc4 - vector-Hessian-vector product function for $G_{PP}$
671: - ctx - [optional] function context
673: Calling sequence of `rhshessianproductfunc1`:
674: + ts - the `TS` context
675: . t - current timestep
676: . U - input vector (current ODE solution)
677: . Vl - an array of input vectors to be left-multiplied with the Hessian
678: . Vr - input vector to be right-multiplied with the Hessian
679: . VHV - an array of output vectors for vector-Hessian-vector product
680: - ctx - [optional] function context
682: Level: intermediate
684: Notes:
685: All other functions have the same calling sequence as `rhshessianproductfunc1`, so their
686: descriptions are omitted for brevity.
688: The first Hessian function and the working array are required.
690: As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
691: $ Vl_n^T*G_UP*Vr$
692: 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$.
693: Each entry of $G_{UP}$ corresponds to the derivative
694: $ G_UP[i][j][k] = \frac{\partial^2 G[i]}{\partial U[j] \partial P[k]}.$
695: The result of the vector-Hessian-vector product for $Vl_n$ needs to be stored in vector $VHV_n$ with j-th entry being
696: $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * G_UP[i][j][k] * Vr[k]}$
697: If the cost function is a scalar, there will be only one vector in $Vl$ and $VHV$.
699: .seealso: `TS`, `TSAdjoint`
700: @*/
701: PetscErrorCode TSSetRHSHessianProduct(TS ts, Vec rhshp1[], PetscErrorCode (*rhshessianproductfunc1)(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[], PetscCtx ctx), Vec rhshp2[], PetscErrorCode (*rhshessianproductfunc2)(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[], PetscCtx ctx), Vec rhshp3[], PetscErrorCode (*rhshessianproductfunc3)(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[], PetscCtx ctx), Vec rhshp4[], PetscErrorCode (*rhshessianproductfunc4)(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[], PetscCtx ctx), PetscCtx ctx)
702: {
703: PetscFunctionBegin;
705: PetscAssertPointer(rhshp1, 2);
707: ts->rhshessianproductctx = ctx;
708: if (rhshp1) ts->vecs_guu = rhshp1;
709: if (rhshp2) ts->vecs_gup = rhshp2;
710: if (rhshp3) ts->vecs_gpu = rhshp3;
711: if (rhshp4) ts->vecs_gpp = rhshp4;
712: ts->rhshessianproduct_guu = rhshessianproductfunc1;
713: ts->rhshessianproduct_gup = rhshessianproductfunc2;
714: ts->rhshessianproduct_gpu = rhshessianproductfunc3;
715: ts->rhshessianproduct_gpp = rhshessianproductfunc4;
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: /*@
720: TSComputeRHSHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for $G_{uu}$.
722: Collective
724: Input Parameters:
725: + ts - The `TS` context obtained from `TSCreate()`
726: . t - the time
727: . U - the solution at which to compute the Hessian product
728: . Vl - the array of input vectors to be multiplied with the Hessian from the left
729: - Vr - the input vector to be multiplied with the Hessian from the right
731: Output Parameter:
732: . VHV - the array of output vectors that store the Hessian product
734: Level: developer
736: Note:
737: `TSComputeRHSHessianProductFunctionUU()` is typically used for sensitivity implementation,
738: so most users would not generally call this routine themselves.
740: .seealso: [](ch_ts), `TS`, `TSSetRHSHessianProduct()`
741: @*/
742: PetscErrorCode TSComputeRHSHessianProductFunctionUU(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[])
743: {
744: PetscFunctionBegin;
745: if (!VHV) PetscFunctionReturn(PETSC_SUCCESS);
749: PetscCallBack("TS callback RHSHessianProduct 1 for sensitivity analysis", (*ts->rhshessianproduct_guu)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: /*@
754: TSComputeRHSHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for $G_{up}$.
756: Collective
758: Input Parameters:
759: + ts - The `TS` context obtained from `TSCreate()`
760: . t - the time
761: . U - the solution at which to compute the Hessian product
762: . Vl - the array of input vectors to be multiplied with the Hessian from the left
763: - Vr - the input vector to be multiplied with the Hessian from the right
765: Output Parameter:
766: . VHV - the array of output vectors that store the Hessian product
768: Level: developer
770: Note:
771: `TSComputeRHSHessianProductFunctionUP()` is typically used for sensitivity implementation,
772: so most users would not generally call this routine themselves.
774: .seealso: [](ch_ts), `TS`, `TSSetRHSHessianProduct()`
775: @*/
776: PetscErrorCode TSComputeRHSHessianProductFunctionUP(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[])
777: {
778: PetscFunctionBegin;
779: if (!VHV) PetscFunctionReturn(PETSC_SUCCESS);
783: PetscCallBack("TS callback RHSHessianProduct 2 for sensitivity analysis", (*ts->rhshessianproduct_gup)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
787: /*@
788: TSComputeRHSHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for $G_{pu}.$
790: Collective
792: Input Parameters:
793: + ts - The `TS` context obtained from `TSCreate()`
794: . t - the time
795: . U - the solution at which to compute the Hessian product
796: . Vl - the array of input vectors to be multiplied with the Hessian from the left
797: - Vr - the input vector to be multiplied with the Hessian from the right
799: Output Parameter:
800: . VHV - the array of output vectors that store the Hessian product
802: Level: developer
804: Note:
805: `TSComputeRHSHessianProductFunctionPU()` is typically used for sensitivity implementation,
806: so most users would not generally call this routine themselves.
808: .seealso: [](ch_ts), `TSSetRHSHessianProduct()`
809: @*/
810: PetscErrorCode TSComputeRHSHessianProductFunctionPU(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[])
811: {
812: PetscFunctionBegin;
813: if (!VHV) PetscFunctionReturn(PETSC_SUCCESS);
817: PetscCallBack("TS callback RHSHessianProduct 3 for sensitivity analysis", (*ts->rhshessianproduct_gpu)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: /*@
822: TSComputeRHSHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for $G_{pp}$.
824: Collective
826: Input Parameters:
827: + ts - The `TS` context obtained from `TSCreate()`
828: . t - the time
829: . U - the solution at which to compute the Hessian product
830: . Vl - the array of input vectors to be multiplied with the Hessian from the left
831: - Vr - the input vector to be multiplied with the Hessian from the right
833: Output Parameter:
834: . VHV - the array of output vectors that store the Hessian product
836: Level: developer
838: Note:
839: `TSComputeRHSHessianProductFunctionPP()` is typically used for sensitivity implementation,
840: so most users would not generally call this routine themselves.
842: .seealso: [](ch_ts), `TSSetRHSHessianProduct()`
843: @*/
844: PetscErrorCode TSComputeRHSHessianProductFunctionPP(TS ts, PetscReal t, Vec U, Vec Vl[], Vec Vr, Vec VHV[])
845: {
846: PetscFunctionBegin;
847: if (!VHV) PetscFunctionReturn(PETSC_SUCCESS);
851: PetscCallBack("TS callback RHSHessianProduct 3 for sensitivity analysis", (*ts->rhshessianproduct_gpp)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: /* --------------------------- Adjoint sensitivity ---------------------------*/
857: /*@
858: TSSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial values and w.r.t. the problem parameters
859: for use by the `TS` adjoint routines.
861: Logically Collective
863: Input Parameters:
864: + ts - the `TS` context obtained from `TSCreate()`
865: . numcost - number of gradients to be computed, this is the number of cost functions
866: . 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
867: - mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
869: Level: beginner
871: Notes:
872: the entries in these vectors must be correctly initialized with the values lambda_i = df/dy|finaltime mu_i = df/dp|finaltime
874: After `TSAdjointSolve()` is called the lambda and the mu contain the computed sensitivities
876: .seealso: `TS`, `TSAdjointSolve()`, `TSGetCostGradients()`
877: @*/
878: PetscErrorCode TSSetCostGradients(TS ts, PetscInt numcost, Vec lambda[], Vec mu[])
879: {
880: PetscFunctionBegin;
882: PetscAssertPointer(lambda, 3);
883: ts->vecs_sensi = lambda;
884: ts->vecs_sensip = mu;
885: PetscCheck(!ts->numcost || ts->numcost == numcost, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "The number of cost functions (2nd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand");
886: ts->numcost = numcost;
887: PetscFunctionReturn(PETSC_SUCCESS);
888: }
890: /*@
891: TSGetCostGradients - Returns the gradients from the `TSAdjointSolve()`
893: Not Collective, but the vectors returned are parallel if `TS` is parallel
895: Input Parameter:
896: . ts - the `TS` context obtained from `TSCreate()`
898: Output Parameters:
899: + numcost - size of returned arrays
900: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
901: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
903: Level: intermediate
905: .seealso: [](ch_ts), `TS`, `TSAdjointSolve()`, `TSSetCostGradients()`
906: @*/
907: PetscErrorCode TSGetCostGradients(TS ts, PetscInt *numcost, Vec *lambda[], Vec *mu[])
908: {
909: PetscFunctionBegin;
911: if (numcost) *numcost = ts->numcost;
912: if (lambda) *lambda = ts->vecs_sensi;
913: if (mu) *mu = ts->vecs_sensip;
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: /*@
918: 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
919: for use by the `TS` adjoint routines.
921: Logically Collective
923: Input Parameters:
924: + ts - the `TS` context obtained from `TSCreate()`
925: . numcost - number of cost functions
926: . 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
927: . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
928: - dir - the direction vector that are multiplied with the Hessian of the cost functions
930: Level: beginner
932: Notes:
933: Hessian of the cost function is completely different from Hessian of the ODE/DAE system
935: For second-order adjoint, one needs to call this function and then `TSAdjointSetForward()` before `TSSolve()`.
937: 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.
939: Passing `NULL` for `lambda2` disables the second-order calculation.
941: .seealso: [](ch_ts), `TS`, `TSAdjointSolve()`, `TSAdjointSetForward()`
942: @*/
943: PetscErrorCode TSSetCostHessianProducts(TS ts, PetscInt numcost, Vec lambda2[], Vec mu2[], Vec dir)
944: {
945: PetscFunctionBegin;
947: PetscCheck(!ts->numcost || ts->numcost == numcost, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "The number of cost functions (2nd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand");
948: ts->numcost = numcost;
949: ts->vecs_sensi2 = lambda2;
950: ts->vecs_sensi2p = mu2;
951: ts->vec_dir = dir;
952: PetscFunctionReturn(PETSC_SUCCESS);
953: }
955: /*@
956: TSGetCostHessianProducts - Returns the gradients from the `TSAdjointSolve()`
958: Not Collective, but vectors returned are parallel if `TS` is parallel
960: Input Parameter:
961: . ts - the `TS` context obtained from `TSCreate()`
963: Output Parameters:
964: + numcost - number of cost functions
965: . 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
966: . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
967: - dir - the direction vector that are multiplied with the Hessian of the cost functions
969: Level: intermediate
971: .seealso: [](ch_ts), `TSAdjointSolve()`, `TSSetCostHessianProducts()`
972: @*/
973: PetscErrorCode TSGetCostHessianProducts(TS ts, PetscInt *numcost, Vec *lambda2[], Vec *mu2[], Vec *dir)
974: {
975: PetscFunctionBegin;
977: if (numcost) *numcost = ts->numcost;
978: if (lambda2) *lambda2 = ts->vecs_sensi2;
979: if (mu2) *mu2 = ts->vecs_sensi2p;
980: if (dir) *dir = ts->vec_dir;
981: PetscFunctionReturn(PETSC_SUCCESS);
982: }
984: /*@
985: TSAdjointSetForward - Trigger the tangent linear solver and initialize the forward sensitivities
987: Logically Collective
989: Input Parameters:
990: + ts - the `TS` context obtained from `TSCreate()`
991: - didp - the derivative of initial values w.r.t. parameters
993: Level: intermediate
995: Notes:
996: When computing sensitivities w.r.t. initial condition, set didp to `NULL` so that the solver will take it as an identity matrix mathematically.
997: `TSAdjoint` does not reset the tangent linear solver automatically, `TSAdjointResetForward()` should be called to reset the tangent linear solver.
999: .seealso: [](ch_ts), `TSAdjointSolve()`, `TSSetCostHessianProducts()`, `TSAdjointResetForward()`
1000: @*/
1001: PetscErrorCode TSAdjointSetForward(TS ts, Mat didp)
1002: {
1003: Mat A;
1004: Vec sp;
1005: PetscScalar *xarr;
1006: PetscInt lsize;
1008: PetscFunctionBegin;
1009: ts->forward_solve = PETSC_TRUE; /* turn on tangent linear mode */
1010: PetscCheck(ts->vecs_sensi2, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetCostHessianProducts() first");
1011: PetscCheck(ts->vec_dir, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Directional vector is missing. Call TSSetCostHessianProducts() to set it.");
1012: /* create a single-column dense matrix */
1013: PetscCall(VecGetLocalSize(ts->vec_sol, &lsize));
1014: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ts), lsize, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, &A));
1016: PetscCall(VecDuplicate(ts->vec_sol, &sp));
1017: PetscCall(MatDenseGetColumn(A, 0, &xarr));
1018: PetscCall(VecPlaceArray(sp, xarr));
1019: if (ts->vecs_sensi2p) { /* tangent linear variable initialized as 2*dIdP*dir */
1020: if (didp) {
1021: PetscCall(MatMult(didp, ts->vec_dir, sp));
1022: PetscCall(VecScale(sp, 2.));
1023: } else {
1024: PetscCall(VecZeroEntries(sp));
1025: }
1026: } else { /* tangent linear variable initialized as dir */
1027: PetscCall(VecCopy(ts->vec_dir, sp));
1028: }
1029: PetscCall(VecResetArray(sp));
1030: PetscCall(MatDenseRestoreColumn(A, &xarr));
1031: PetscCall(VecDestroy(&sp));
1033: PetscCall(TSForwardSetInitialSensitivities(ts, A)); /* if didp is NULL, identity matrix is assumed */
1035: PetscCall(MatDestroy(&A));
1036: PetscFunctionReturn(PETSC_SUCCESS);
1037: }
1039: /*@
1040: TSAdjointResetForward - Reset the tangent linear solver and destroy the tangent linear context
1042: Logically Collective
1044: Input Parameter:
1045: . ts - the `TS` context obtained from `TSCreate()`
1047: Level: intermediate
1049: .seealso: [](ch_ts), `TSAdjointSetForward()`
1050: @*/
1051: PetscErrorCode TSAdjointResetForward(TS ts)
1052: {
1053: PetscFunctionBegin;
1054: ts->forward_solve = PETSC_FALSE; /* turn off tangent linear mode */
1055: PetscCall(TSForwardReset(ts));
1056: PetscFunctionReturn(PETSC_SUCCESS);
1057: }
1059: /*@
1060: TSAdjointSetUp - Sets up the internal data structures for the later use
1061: of an adjoint solver
1063: Collective
1065: Input Parameter:
1066: . ts - the `TS` context obtained from `TSCreate()`
1068: Level: advanced
1070: .seealso: [](ch_ts), `TSCreate()`, `TSAdjointStep()`, `TSSetCostGradients()`
1071: @*/
1072: PetscErrorCode TSAdjointSetUp(TS ts)
1073: {
1074: TSTrajectory tj;
1075: PetscBool match;
1077: PetscFunctionBegin;
1079: if (ts->adjointsetupcalled) PetscFunctionReturn(PETSC_SUCCESS);
1080: PetscCheck(ts->vecs_sensi, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetCostGradients() first");
1081: PetscCheck(!ts->vecs_sensip || ts->Jacp || ts->Jacprhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetRHSJacobianP() or TSSetIJacobianP() first");
1082: PetscCall(TSGetTrajectory(ts, &tj));
1083: PetscCall(PetscObjectTypeCompare((PetscObject)tj, TSTRAJECTORYBASIC, &match));
1084: if (match) {
1085: PetscBool solution_only;
1086: PetscCall(TSTrajectoryGetSolutionOnly(tj, &solution_only));
1087: PetscCheck(!solution_only, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "TSAdjoint cannot use the solution-only mode when choosing the Basic TSTrajectory type. Turn it off with -ts_trajectory_solution_only 0");
1088: }
1089: PetscCall(TSTrajectorySetUseHistory(tj, PETSC_FALSE)); /* not use TSHistory */
1091: if (ts->quadraturets) { /* if there is integral in the cost function */
1092: PetscCall(VecDuplicate(ts->vecs_sensi[0], &ts->vec_drdu_col));
1093: if (ts->vecs_sensip) PetscCall(VecDuplicate(ts->vecs_sensip[0], &ts->vec_drdp_col));
1094: }
1096: PetscTryTypeMethod(ts, adjointsetup);
1097: ts->adjointsetupcalled = PETSC_TRUE;
1098: PetscFunctionReturn(PETSC_SUCCESS);
1099: }
1101: /*@
1102: TSAdjointReset - Resets a `TS` adjoint context and removes any allocated `Vec`s and `Mat`s.
1104: Collective
1106: Input Parameter:
1107: . ts - the `TS` context obtained from `TSCreate()`
1109: Level: beginner
1111: .seealso: [](ch_ts), `TSCreate()`, `TSAdjointSetUp()`, `TSDestroy()`
1112: @*/
1113: PetscErrorCode TSAdjointReset(TS ts)
1114: {
1115: PetscFunctionBegin;
1117: PetscTryTypeMethod(ts, adjointreset);
1118: if (ts->quadraturets) { /* if there is integral in the cost function */
1119: PetscCall(VecDestroy(&ts->vec_drdu_col));
1120: if (ts->vecs_sensip) PetscCall(VecDestroy(&ts->vec_drdp_col));
1121: }
1122: ts->vecs_sensi = NULL;
1123: ts->vecs_sensip = NULL;
1124: ts->vecs_sensi2 = NULL;
1125: ts->vecs_sensi2p = NULL;
1126: ts->vec_dir = NULL;
1127: ts->adjointsetupcalled = PETSC_FALSE;
1128: PetscFunctionReturn(PETSC_SUCCESS);
1129: }
1131: /*@
1132: TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time
1134: Logically Collective
1136: Input Parameters:
1137: + ts - the `TS` context obtained from `TSCreate()`
1138: - steps - number of steps to use
1140: Level: intermediate
1142: Notes:
1143: Normally one does not call this and `TSAdjointSolve()` integrates back to the original timestep. One can call this
1144: so as to integrate back to less than the original timestep
1146: .seealso: [](ch_ts), `TSAdjointSolve()`, `TS`, `TSSetExactFinalTime()`
1147: @*/
1148: PetscErrorCode TSAdjointSetSteps(TS ts, PetscInt steps)
1149: {
1150: PetscFunctionBegin;
1153: PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot step back a negative number of steps");
1154: PetscCheck(steps <= ts->steps, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot step back more than the total number of forward steps");
1155: ts->adjoint_max_steps = steps;
1156: PetscFunctionReturn(PETSC_SUCCESS);
1157: }
1159: // PetscClangLinter pragma disable: -fdoc-*
1160: /*@C
1161: TSAdjointSetRHSJacobian - Deprecated, use `TSSetRHSJacobianP()`
1163: Level: deprecated
1164: @*/
1165: PetscErrorCode TSAdjointSetRHSJacobian(TS ts, Mat Amat, PetscErrorCode (*func)(TS, PetscReal, Vec, Mat, void *), PetscCtx ctx)
1166: {
1167: PetscFunctionBegin;
1171: ts->rhsjacobianp = func;
1172: ts->rhsjacobianpctx = ctx;
1173: if (Amat) {
1174: PetscCall(PetscObjectReference((PetscObject)Amat));
1175: PetscCall(MatDestroy(&ts->Jacp));
1176: ts->Jacp = Amat;
1177: }
1178: PetscFunctionReturn(PETSC_SUCCESS);
1179: }
1181: // PetscClangLinter pragma disable: -fdoc-*
1182: /*@C
1183: TSAdjointComputeRHSJacobian - Deprecated, use `TSComputeRHSJacobianP()`
1185: Level: deprecated
1186: @*/
1187: PetscErrorCode TSAdjointComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat Amat)
1188: {
1189: PetscFunctionBegin;
1194: PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->rhsjacobianp)(ts, t, U, Amat, ts->rhsjacobianpctx));
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: // PetscClangLinter pragma disable: -fdoc-*
1199: /*@
1200: TSAdjointComputeDRDYFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobian()`
1202: Level: deprecated
1203: @*/
1204: PetscErrorCode TSAdjointComputeDRDYFunction(TS ts, PetscReal t, Vec U, Vec *DRDU)
1205: {
1206: PetscFunctionBegin;
1210: PetscCallBack("TS callback DRDY for sensitivity analysis", (*ts->drdufunction)(ts, t, U, DRDU, ts->costintegrandctx));
1211: PetscFunctionReturn(PETSC_SUCCESS);
1212: }
1214: // PetscClangLinter pragma disable: -fdoc-*
1215: /*@
1216: TSAdjointComputeDRDPFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobianP()`
1218: Level: deprecated
1219: @*/
1220: PetscErrorCode TSAdjointComputeDRDPFunction(TS ts, PetscReal t, Vec U, Vec *DRDP)
1221: {
1222: PetscFunctionBegin;
1226: PetscCallBack("TS callback DRDP for sensitivity analysis", (*ts->drdpfunction)(ts, t, U, DRDP, ts->costintegrandctx));
1227: PetscFunctionReturn(PETSC_SUCCESS);
1228: }
1230: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
1231: /*@C
1232: TSAdjointMonitorSensi - monitors the first lambda sensitivity
1234: Level: intermediate
1236: .seealso: [](ch_ts), `TSAdjointMonitorSet()`
1237: @*/
1238: static PetscErrorCode TSAdjointMonitorSensi(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscInt numcost, Vec *lambda, Vec *mu, PetscViewerAndFormat *vf)
1239: {
1240: PetscViewer viewer = vf->viewer;
1242: PetscFunctionBegin;
1244: PetscCall(PetscViewerPushFormat(viewer, vf->format));
1245: PetscCall(VecView(lambda[0], viewer));
1246: PetscCall(PetscViewerPopFormat(viewer));
1247: PetscFunctionReturn(PETSC_SUCCESS);
1248: }
1250: /*@C
1251: TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
1253: Collective
1255: Input Parameters:
1256: + ts - `TS` object you wish to monitor
1257: . name - the monitor type one is seeking
1258: . help - message indicating what monitoring is done
1259: . manual - manual page for the monitor
1260: . monitor - the monitor function, its context must be a `PetscViewerAndFormat`
1261: - 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
1263: Calling sequence of `monitor`:
1264: + ts - the `TS` context
1265: . step - iteration number (after the final time step the monitor routine is called with
1266: a step of -1, this is at the final time which may have been interpolated to)
1267: . time - current time
1268: . u - current iterate
1269: . numcost - number of cost functions
1270: . lambda - sensitivities to initial conditions
1271: . mu - sensitivities to parameters
1272: - vf - the `PetscViewer` and format the monitor is using
1274: Calling sequence of `monitorsetup`:
1275: + ts - the `TS` object being monitored
1276: - vf - the `PetscViewer` and format the monitor is using
1278: Level: developer
1280: .seealso: [](ch_ts), `PetscOptionsCreateViewer()`, `PetscOptionsGetReal()`, `PetscOptionsHasName()`, `PetscOptionsGetString()`,
1281: `PetscOptionsGetIntArray()`, `PetscOptionsGetRealArray()`, `PetscOptionsBool()`
1282: `PetscOptionsInt()`, `PetscOptionsString()`, `PetscOptionsReal()`,
1283: `PetscOptionsName()`, `PetscOptionsBegin()`, `PetscOptionsEnd()`, `PetscOptionsHeadBegin()`,
1284: `PetscOptionsStringArray()`, `PetscOptionsRealArray()`, `PetscOptionsScalar()`,
1285: `PetscOptionsBoolGroupBegin()`, `PetscOptionsBoolGroup()`, `PetscOptionsBoolGroupEnd()`,
1286: `PetscOptionsFList()`, `PetscOptionsEList()`, `PetscViewerAndFormat`
1287: @*/
1288: PetscErrorCode TSAdjointMonitorSetFromOptions(TS ts, const char name[], const char help[], const char manual[], PetscErrorCode (*monitor)(TS ts, PetscInt step, PetscReal time, Vec u, PetscInt numcost, Vec *lambda, Vec *mu, PetscViewerAndFormat *vf), PetscErrorCode (*monitorsetup)(TS ts, PetscViewerAndFormat *vf))
1289: {
1290: PetscViewer viewer;
1291: PetscViewerFormat format;
1292: PetscBool flg;
1294: PetscFunctionBegin;
1295: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, name, &viewer, &format, &flg));
1296: if (flg) {
1297: PetscViewerAndFormat *vf;
1298: PetscCall(PetscViewerAndFormatCreate(viewer, format, &vf));
1299: PetscCall(PetscViewerDestroy(&viewer));
1300: if (monitorsetup) PetscCall((*monitorsetup)(ts, vf));
1301: PetscCall(TSAdjointMonitorSet(ts, (PetscErrorCode (*)(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, PetscCtx))monitor, vf, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy));
1302: }
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: /*@C
1307: TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
1308: timestep to display the iteration's progress.
1310: Logically Collective
1312: Input Parameters:
1313: + ts - the `TS` context obtained from `TSCreate()`
1314: . adjointmonitor - monitoring routine
1315: . adjointmctx - [optional] context for private data for the monitor routine (use `NULL` if no context is desired)
1316: - adjointmdestroy - [optional] routine that frees monitor context (may be `NULL`), see `PetscCtxDestroyFn` for its calling sequence
1318: Calling sequence of `adjointmonitor`:
1319: + ts - the `TS` context
1320: . steps - iteration number (after the final time step the monitor routine is called with
1321: a step of -1, this is at the final time which may have been interpolated to)
1322: . time - current time
1323: . u - current iterate
1324: . numcost - number of cost functions
1325: . lambda - sensitivities to initial conditions
1326: . mu - sensitivities to parameters
1327: - adjointmctx - [optional] adjoint monitoring context
1329: Level: intermediate
1331: Note:
1332: This routine adds an additional monitor to the list of monitors that
1333: already has been loaded.
1335: Fortran Notes:
1336: Only a single monitor function can be set for each `TS` object
1338: .seealso: [](ch_ts), `TS`, `TSAdjointSolve()`, `TSAdjointMonitorCancel()`, `PetscCtxDestroyFn`
1339: @*/
1340: PetscErrorCode TSAdjointMonitorSet(TS ts, PetscErrorCode (*adjointmonitor)(TS ts, PetscInt steps, PetscReal time, Vec u, PetscInt numcost, Vec *lambda, Vec *mu, PetscCtx adjointmctx), PetscCtx adjointmctx, PetscCtxDestroyFn *adjointmdestroy)
1341: {
1342: PetscFunctionBegin;
1344: for (PetscInt i = 0; i < ts->numbermonitors; i++) {
1345: PetscBool identical;
1347: PetscCall(PetscMonitorCompare((PetscErrorCode (*)(void))(PetscVoidFn *)adjointmonitor, adjointmctx, adjointmdestroy, (PetscErrorCode (*)(void))(PetscVoidFn *)ts->adjointmonitor[i], ts->adjointmonitorcontext[i], ts->adjointmonitordestroy[i], &identical));
1348: if (identical) PetscFunctionReturn(PETSC_SUCCESS);
1349: }
1350: PetscCheck(ts->numberadjointmonitors < MAXTSMONITORS, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many adjoint monitors set");
1351: ts->adjointmonitor[ts->numberadjointmonitors] = adjointmonitor;
1352: ts->adjointmonitordestroy[ts->numberadjointmonitors] = adjointmdestroy;
1353: ts->adjointmonitorcontext[ts->numberadjointmonitors++] = adjointmctx;
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: /*@C
1358: TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.
1360: Logically Collective
1362: Input Parameter:
1363: . ts - the `TS` context obtained from `TSCreate()`
1365: Notes:
1366: There is no way to remove a single, specific monitor.
1368: Level: intermediate
1370: .seealso: [](ch_ts), `TS`, `TSAdjointSolve()`, `TSAdjointMonitorSet()`
1371: @*/
1372: PetscErrorCode TSAdjointMonitorCancel(TS ts)
1373: {
1374: PetscInt i;
1376: PetscFunctionBegin;
1378: for (i = 0; i < ts->numberadjointmonitors; i++) {
1379: if (ts->adjointmonitordestroy[i]) PetscCall((*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]));
1380: }
1381: ts->numberadjointmonitors = 0;
1382: PetscFunctionReturn(PETSC_SUCCESS);
1383: }
1385: /*@C
1386: TSAdjointMonitorDefault - the default monitor of adjoint computations
1388: Input Parameters:
1389: + ts - the `TS` context
1390: . step - iteration number (after the final time step the monitor routine is called with a
1391: step of -1, this is at the final time which may have been interpolated to)
1392: . time - current time
1393: . v - current iterate
1394: . numcost - number of cost functions
1395: . lambda - sensitivities to initial conditions
1396: . mu - sensitivities to parameters
1397: - vf - the viewer and format
1399: Level: intermediate
1401: .seealso: [](ch_ts), `TS`, `TSAdjointSolve()`, `TSAdjointMonitorSet()`
1402: @*/
1403: PetscErrorCode TSAdjointMonitorDefault(TS ts, PetscInt step, PetscReal time, Vec v, PetscInt numcost, Vec lambda[], Vec mu[], PetscViewerAndFormat *vf)
1404: {
1405: PetscViewer viewer = vf->viewer;
1407: PetscFunctionBegin;
1408: (void)v;
1409: (void)numcost;
1410: (void)lambda;
1411: (void)mu;
1413: PetscCall(PetscViewerPushFormat(viewer, vf->format));
1414: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel));
1415: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " TS dt %g time %g%s", step, (double)ts->time_step, (double)time, ts->steprollback ? " (r)\n" : "\n"));
1416: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel));
1417: PetscCall(PetscViewerPopFormat(viewer));
1418: PetscFunctionReturn(PETSC_SUCCESS);
1419: }
1421: /*@C
1422: TSAdjointMonitorDrawSensi - Monitors progress of the adjoint `TS` solvers by calling
1423: `VecView()` for the sensitivities to initial states at each timestep
1425: Collective
1427: Input Parameters:
1428: + ts - the `TS` context
1429: . step - current time-step
1430: . ptime - current time
1431: . u - current state
1432: . numcost - number of cost functions
1433: . lambda - sensitivities to initial conditions
1434: . mu - sensitivities to parameters
1435: - dummy - either a viewer or `NULL`
1437: Level: intermediate
1439: .seealso: [](ch_ts), `TSAdjointSolve()`, `TSAdjointMonitorSet()`, `TSAdjointMonitorDefault()`, `VecView()`
1440: @*/
1441: PetscErrorCode TSAdjointMonitorDrawSensi(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscInt numcost, Vec lambda[], Vec mu[], void *dummy)
1442: {
1443: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
1444: PetscDraw draw;
1445: PetscReal xl, yl, xr, yr, h;
1446: char time[32];
1448: PetscFunctionBegin;
1449: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) PetscFunctionReturn(PETSC_SUCCESS);
1451: PetscCall(VecView(lambda[0], ictx->viewer));
1452: PetscCall(PetscViewerDrawGetDraw(ictx->viewer, 0, &draw));
1453: PetscCall(PetscSNPrintf(time, 32, "Timestep %" PetscInt_FMT " Time %g", step, (double)ptime));
1454: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
1455: h = yl + .95 * (yr - yl);
1456: PetscCall(PetscDrawStringCentered(draw, .5 * (xl + xr), h, PETSC_DRAW_BLACK, time));
1457: PetscCall(PetscDrawFlush(draw));
1458: PetscFunctionReturn(PETSC_SUCCESS);
1459: }
1461: /*@C
1462: TSAdjointSetFromOptions - Sets various `TS` adjoint parameters from options database.
1464: Collective
1466: Input Parameters:
1467: + ts - the `TS` context
1468: - PetscOptionsObject - the options context
1470: Options Database Keys:
1471: + -ts_adjoint_solve <yes,no> - After solving the ODE/DAE solve the adjoint problem (requires `-ts_save_trajectory`)
1472: . -ts_adjoint_monitor - print information at each adjoint time step
1473: - -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically
1475: Level: developer
1477: Note:
1478: This is not normally called directly by users
1480: .seealso: [](ch_ts), `TSSetSaveTrajectory()`, `TSTrajectorySetUp()`
1481: @*/
1482: PetscErrorCode TSAdjointSetFromOptions(TS ts, PetscOptionItems PetscOptionsObject)
1483: {
1484: PetscBool tflg, opt;
1486: PetscFunctionBegin;
1488: PetscOptionsHeadBegin(PetscOptionsObject, "TS Adjoint options");
1489: tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
1490: PetscCall(PetscOptionsBool("-ts_adjoint_solve", "Solve the adjoint problem immediately after solving the forward problem", "", tflg, &tflg, &opt));
1491: if (opt) {
1492: PetscCall(TSSetSaveTrajectory(ts));
1493: ts->adjoint_solve = tflg;
1494: }
1495: PetscCall(TSAdjointMonitorSetFromOptions(ts, "-ts_adjoint_monitor", "Monitor adjoint timestep size", "TSAdjointMonitorDefault", TSAdjointMonitorDefault, NULL));
1496: PetscCall(TSAdjointMonitorSetFromOptions(ts, "-ts_adjoint_monitor_sensi", "Monitor sensitivity in the adjoint computation", "TSAdjointMonitorSensi", TSAdjointMonitorSensi, NULL));
1497: opt = PETSC_FALSE;
1498: PetscCall(PetscOptionsName("-ts_adjoint_monitor_draw_sensi", "Monitor adjoint sensitivities (lambda only) graphically", "TSAdjointMonitorDrawSensi", &opt));
1499: if (opt) {
1500: TSMonitorDrawCtx ctx;
1501: PetscInt howoften = 1;
1503: PetscCall(PetscOptionsInt("-ts_adjoint_monitor_draw_sensi", "Monitor adjoint sensitivities (lambda only) graphically", "TSAdjointMonitorDrawSensi", howoften, &howoften, NULL));
1504: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
1505: PetscCall(TSAdjointMonitorSet(ts, TSAdjointMonitorDrawSensi, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
1506: }
1507: PetscFunctionReturn(PETSC_SUCCESS);
1508: }
1510: /*@
1511: TSAdjointStep - Steps one time step backward in the adjoint run
1513: Collective
1515: Input Parameter:
1516: . ts - the `TS` context obtained from `TSCreate()`
1518: Level: intermediate
1520: .seealso: [](ch_ts), `TSAdjointSetUp()`, `TSAdjointSolve()`
1521: @*/
1522: PetscErrorCode TSAdjointStep(TS ts)
1523: {
1524: DM dm;
1526: PetscFunctionBegin;
1528: PetscCall(TSGetDM(ts, &dm));
1529: PetscCall(TSAdjointSetUp(ts));
1530: ts->steps--; /* must decrease the step index before the adjoint step is taken. */
1532: ts->reason = TS_CONVERGED_ITERATING;
1533: ts->ptime_prev = ts->ptime;
1534: PetscCall(PetscLogEventBegin(TS_AdjointStep, ts, 0, 0, 0));
1535: PetscUseTypeMethod(ts, adjointstep);
1536: PetscCall(PetscLogEventEnd(TS_AdjointStep, ts, 0, 0, 0));
1537: ts->adjoint_steps++;
1539: if (ts->reason < 0) {
1540: PetscCheck(!ts->errorifstepfailed, PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSAdjointStep has failed due to %s", TSConvergedReasons[ts->reason]);
1541: } else if (!ts->reason) {
1542: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1543: }
1544: PetscFunctionReturn(PETSC_SUCCESS);
1545: }
1547: /*@
1548: TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE
1550: Collective
1551: `
1553: Input Parameter:
1554: . ts - the `TS` context obtained from `TSCreate()`
1556: Options Database Key:
1557: . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial values
1559: Level: intermediate
1561: Notes:
1562: This must be called after a call to `TSSolve()` that solves the forward problem
1564: By default this will integrate back to the initial time, one can use `TSAdjointSetSteps()` to step back to a later time
1566: .seealso: [](ch_ts), `TSCreate()`, `TSSetCostGradients()`, `TSSetSolution()`, `TSAdjointStep()`
1567: @*/
1568: PetscErrorCode TSAdjointSolve(TS ts)
1569: {
1570: static PetscBool cite = PETSC_FALSE;
1571: #if defined(TSADJOINT_STAGE)
1572: PetscLogStage adjoint_stage;
1573: #endif
1575: PetscFunctionBegin;
1577: PetscCall(PetscCitationsRegister("@article{Zhang2022tsadjoint,\n"
1578: " title = {{PETSc TSAdjoint: A Discrete Adjoint ODE Solver for First-Order and Second-Order Sensitivity Analysis}},\n"
1579: " author = {Zhang, Hong and Constantinescu, Emil M. and Smith, Barry F.},\n"
1580: " journal = {SIAM Journal on Scientific Computing},\n"
1581: " volume = {44},\n"
1582: " number = {1},\n"
1583: " pages = {C1-C24},\n"
1584: " doi = {10.1137/21M140078X},\n"
1585: " year = {2022}\n}\n",
1586: &cite));
1587: #if defined(TSADJOINT_STAGE)
1588: PetscCall(PetscLogStageRegister("TSAdjoint", &adjoint_stage));
1589: PetscCall(PetscLogStagePush(adjoint_stage));
1590: #endif
1591: PetscCall(TSAdjointSetUp(ts));
1593: /* reset time step and iteration counters */
1594: ts->adjoint_steps = 0;
1595: ts->ksp_its = 0;
1596: ts->snes_its = 0;
1597: ts->num_snes_failures = 0;
1598: ts->reject = 0;
1599: ts->reason = TS_CONVERGED_ITERATING;
1601: if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->steps;
1602: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1604: while (!ts->reason) {
1605: PetscCall(TSTrajectoryGet(ts->trajectory, ts, ts->steps, &ts->ptime));
1606: PetscCall(TSAdjointMonitor(ts, ts->steps, ts->ptime, ts->vec_sol, ts->numcost, ts->vecs_sensi, ts->vecs_sensip));
1607: PetscCall(TSAdjointEventHandler(ts));
1608: PetscCall(TSAdjointStep(ts));
1609: if ((ts->vec_costintegral || ts->quadraturets) && !ts->costintegralfwd) PetscCall(TSAdjointCostIntegral(ts));
1610: }
1611: if (!ts->steps) {
1612: PetscCall(TSTrajectoryGet(ts->trajectory, ts, ts->steps, &ts->ptime));
1613: PetscCall(TSAdjointMonitor(ts, ts->steps, ts->ptime, ts->vec_sol, ts->numcost, ts->vecs_sensi, ts->vecs_sensip));
1614: }
1615: ts->solvetime = ts->ptime;
1616: PetscCall(TSTrajectoryViewFromOptions(ts->trajectory, NULL, "-ts_trajectory_view"));
1617: PetscCall(VecViewFromOptions(ts->vecs_sensi[0], (PetscObject)ts, "-ts_adjoint_view_solution"));
1618: ts->adjoint_max_steps = 0;
1619: #if defined(TSADJOINT_STAGE)
1620: PetscCall(PetscLogStagePop());
1621: #endif
1622: PetscFunctionReturn(PETSC_SUCCESS);
1623: }
1625: /*@C
1626: TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using `TSAdjointMonitorSet()`
1628: Collective
1630: Input Parameters:
1631: + ts - time stepping context obtained from `TSCreate()`
1632: . step - step number that has just completed
1633: . ptime - model time of the state
1634: . u - state at the current model time
1635: . numcost - number of cost functions (dimension of lambda or mu)
1636: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
1637: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
1639: Level: developer
1641: Note:
1642: `TSAdjointMonitor()` is typically used automatically within the time stepping implementations.
1643: Users would almost never call this routine directly.
1645: .seealso: `TSAdjointMonitorSet()`, `TSAdjointSolve()`
1646: @*/
1647: PetscErrorCode TSAdjointMonitor(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscInt numcost, Vec lambda[], Vec mu[])
1648: {
1649: PetscInt i, n = ts->numberadjointmonitors;
1651: PetscFunctionBegin;
1654: PetscCall(VecLockReadPush(u));
1655: for (i = 0; i < n; i++) PetscCall((*ts->adjointmonitor[i])(ts, step, ptime, u, numcost, lambda, mu, ts->adjointmonitorcontext[i]));
1656: PetscCall(VecLockReadPop(u));
1657: PetscFunctionReturn(PETSC_SUCCESS);
1658: }
1660: /*@
1661: TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.
1663: Collective
1665: Input Parameter:
1666: . ts - time stepping context
1668: Level: advanced
1670: Notes:
1671: This function cannot be called until `TSAdjointStep()` has been completed.
1673: .seealso: [](ch_ts), `TSAdjointSolve()`, `TSAdjointStep()`
1674: @*/
1675: PetscErrorCode TSAdjointCostIntegral(TS ts)
1676: {
1677: PetscFunctionBegin;
1679: PetscUseTypeMethod(ts, adjointintegral);
1680: PetscFunctionReturn(PETSC_SUCCESS);
1681: }
1683: /* ------------------ Forward (tangent linear) sensitivity ------------------*/
1685: /*@
1686: TSForwardSetUp - Sets up the internal data structures for the later use
1687: of forward sensitivity analysis
1689: Collective
1691: Input Parameter:
1692: . ts - the `TS` context obtained from `TSCreate()`
1694: Level: advanced
1696: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSDestroy()`, `TSSetUp()`
1697: @*/
1698: PetscErrorCode TSForwardSetUp(TS ts)
1699: {
1700: PetscFunctionBegin;
1702: if (ts->forwardsetupcalled) PetscFunctionReturn(PETSC_SUCCESS);
1703: PetscTryTypeMethod(ts, forwardsetup);
1704: PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sensip_col));
1705: ts->forwardsetupcalled = PETSC_TRUE;
1706: PetscFunctionReturn(PETSC_SUCCESS);
1707: }
1709: /*@
1710: TSForwardReset - Reset the internal data structures used by forward sensitivity analysis
1712: Collective
1714: Input Parameter:
1715: . ts - the `TS` context obtained from `TSCreate()`
1717: Level: advanced
1719: .seealso: [](ch_ts), `TSCreate()`, `TSDestroy()`, `TSForwardSetUp()`
1720: @*/
1721: PetscErrorCode TSForwardReset(TS ts)
1722: {
1723: TS quadts = ts->quadraturets;
1725: PetscFunctionBegin;
1727: PetscTryTypeMethod(ts, forwardreset);
1728: PetscCall(MatDestroy(&ts->mat_sensip));
1729: if (quadts) PetscCall(MatDestroy(&quadts->mat_sensip));
1730: PetscCall(VecDestroy(&ts->vec_sensip_col));
1731: ts->forward_solve = PETSC_FALSE;
1732: ts->forwardsetupcalled = PETSC_FALSE;
1733: PetscFunctionReturn(PETSC_SUCCESS);
1734: }
1736: /*@
1737: TSForwardSetIntegralGradients - Set the vectors holding forward sensitivities of the integral term.
1739: Input Parameters:
1740: + ts - the `TS` context obtained from `TSCreate()`
1741: . numfwdint - number of integrals
1742: - vp - the vectors containing the gradients for each integral w.r.t. parameters
1744: Level: deprecated
1746: .seealso: [](ch_ts), `TSForwardGetSensitivities()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1747: @*/
1748: PetscErrorCode TSForwardSetIntegralGradients(TS ts, PetscInt numfwdint, Vec vp[])
1749: {
1750: PetscFunctionBegin;
1752: PetscCheck(!ts->numcost || ts->numcost == numfwdint, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "The number of cost functions (2nd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand()");
1753: if (!ts->numcost) ts->numcost = numfwdint;
1755: ts->vecs_integral_sensip = vp;
1756: PetscFunctionReturn(PETSC_SUCCESS);
1757: }
1759: /*@
1760: TSForwardGetIntegralGradients - Returns the forward sensitivities of the integral term.
1762: Input Parameter:
1763: . ts - the `TS` context obtained from `TSCreate()`
1765: Output Parameters:
1766: + numfwdint - number of integrals
1767: - vp - the vectors containing the gradients for each integral w.r.t. parameters
1769: Level: deprecated
1771: .seealso: [](ch_ts), `TSForwardSetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardStep()`
1772: @*/
1773: PetscErrorCode TSForwardGetIntegralGradients(TS ts, PetscInt *numfwdint, Vec *vp[])
1774: {
1775: PetscFunctionBegin;
1777: PetscAssertPointer(vp, 3);
1778: if (numfwdint) *numfwdint = ts->numcost;
1779: if (vp) *vp = ts->vecs_integral_sensip;
1780: PetscFunctionReturn(PETSC_SUCCESS);
1781: }
1783: /*@
1784: TSForwardStep - Compute the forward sensitivity for one time step.
1786: Collective
1788: Input Parameter:
1789: . ts - time stepping context
1791: Level: advanced
1793: Notes:
1794: This function cannot be called until `TSStep()` has been completed.
1796: .seealso: [](ch_ts), `TSForwardSetSensitivities()`, `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardSetUp()`
1797: @*/
1798: PetscErrorCode TSForwardStep(TS ts)
1799: {
1800: PetscFunctionBegin;
1802: PetscCall(PetscLogEventBegin(TS_ForwardStep, ts, 0, 0, 0));
1803: PetscUseTypeMethod(ts, forwardstep);
1804: PetscCall(PetscLogEventEnd(TS_ForwardStep, ts, 0, 0, 0));
1805: PetscCheck(ts->reason >= 0 || !ts->errorifstepfailed, PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSFowardStep has failed due to %s", TSConvergedReasons[ts->reason]);
1806: PetscFunctionReturn(PETSC_SUCCESS);
1807: }
1809: /*@
1810: TSForwardSetSensitivities - Sets the initial value of the trajectory sensitivities of solution w.r.t. the problem parameters and initial values.
1812: Logically Collective
1814: Input Parameters:
1815: + ts - the `TS` context obtained from `TSCreate()`
1816: . nump - number of parameters
1817: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1819: Level: beginner
1821: Notes:
1822: Use `PETSC_DETERMINE` to use the number of columns of `Smat` for `nump`
1824: Forward sensitivity is also called 'trajectory sensitivity' in some fields such as power systems.
1825: This function turns on a flag to trigger `TSSolve()` to compute forward sensitivities automatically.
1826: You must call this function before `TSSolve()`.
1827: The entries in the sensitivity matrix must be correctly initialized with the values S = dy/dp|startingtime.
1829: .seealso: [](ch_ts), `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1830: @*/
1831: PetscErrorCode TSForwardSetSensitivities(TS ts, PetscInt nump, Mat Smat)
1832: {
1833: PetscFunctionBegin;
1836: ts->forward_solve = PETSC_TRUE;
1837: if (nump == PETSC_DEFAULT || nump == PETSC_DETERMINE) PetscCall(MatGetSize(Smat, NULL, &ts->num_parameters));
1838: else ts->num_parameters = nump;
1839: PetscCall(PetscObjectReference((PetscObject)Smat));
1840: PetscCall(MatDestroy(&ts->mat_sensip));
1841: ts->mat_sensip = Smat;
1842: PetscFunctionReturn(PETSC_SUCCESS);
1843: }
1845: /*@
1846: TSForwardGetSensitivities - Returns the trajectory sensitivities
1848: Not Collective, but Smat returned is parallel if ts is parallel
1850: Output Parameters:
1851: + ts - the `TS` context obtained from `TSCreate()`
1852: . nump - number of parameters
1853: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1855: Level: intermediate
1857: .seealso: [](ch_ts), `TSForwardSetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1858: @*/
1859: PetscErrorCode TSForwardGetSensitivities(TS ts, PetscInt *nump, Mat *Smat)
1860: {
1861: PetscFunctionBegin;
1863: if (nump) *nump = ts->num_parameters;
1864: if (Smat) *Smat = ts->mat_sensip;
1865: PetscFunctionReturn(PETSC_SUCCESS);
1866: }
1868: /*@
1869: TSForwardCostIntegral - Evaluate the cost integral in the forward run.
1871: Collective
1873: Input Parameter:
1874: . ts - time stepping context
1876: Level: advanced
1878: Note:
1879: This function cannot be called until `TSStep()` has been completed.
1881: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSAdjointCostIntegral()`
1882: @*/
1883: PetscErrorCode TSForwardCostIntegral(TS ts)
1884: {
1885: PetscFunctionBegin;
1887: PetscUseTypeMethod(ts, forwardintegral);
1888: PetscFunctionReturn(PETSC_SUCCESS);
1889: }
1891: /*@
1892: TSForwardSetInitialSensitivities - Set initial values for tangent linear sensitivities
1894: Collective
1896: Input Parameters:
1897: + ts - the `TS` context obtained from `TSCreate()`
1898: - didp - parametric sensitivities of the initial condition
1900: Level: intermediate
1902: Notes:
1903: `TSSolve()` allows users to pass the initial solution directly to `TS`. But the tangent linear variables cannot be initialized in this way.
1904: This function is used to set initial values for tangent linear variables.
1906: .seealso: [](ch_ts), `TS`, `TSForwardSetSensitivities()`
1907: @*/
1908: PetscErrorCode TSForwardSetInitialSensitivities(TS ts, Mat didp)
1909: {
1910: PetscFunctionBegin;
1913: if (!ts->mat_sensip) PetscCall(TSForwardSetSensitivities(ts, PETSC_DETERMINE, didp));
1914: PetscFunctionReturn(PETSC_SUCCESS);
1915: }
1917: /*@
1918: TSForwardGetStages - Get the number of stages and the tangent linear sensitivities at the intermediate stages
1920: Input Parameter:
1921: . ts - the `TS` context obtained from `TSCreate()`
1923: Output Parameters:
1924: + ns - number of stages
1925: - S - tangent linear sensitivities at the intermediate stages
1927: Level: advanced
1929: .seealso: `TS`
1930: @*/
1931: PetscErrorCode TSForwardGetStages(TS ts, PetscInt *ns, Mat **S)
1932: {
1933: PetscFunctionBegin;
1936: if (!ts->ops->getstages) *S = NULL;
1937: else PetscUseTypeMethod(ts, forwardgetstages, ns, S);
1938: PetscFunctionReturn(PETSC_SUCCESS);
1939: }
1941: /*@
1942: TSCreateQuadratureTS - Create a sub-`TS` that evaluates integrals over time
1944: Input Parameters:
1945: + ts - the `TS` context obtained from `TSCreate()`
1946: - fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1948: Output Parameter:
1949: . quadts - the child `TS` context
1951: Level: intermediate
1953: .seealso: [](ch_ts), `TSGetQuadratureTS()`
1954: @*/
1955: PetscErrorCode TSCreateQuadratureTS(TS ts, PetscBool fwd, TS *quadts)
1956: {
1957: char prefix[128];
1959: PetscFunctionBegin;
1961: PetscAssertPointer(quadts, 3);
1962: PetscCall(TSDestroy(&ts->quadraturets));
1963: PetscCall(TSCreate(PetscObjectComm((PetscObject)ts), &ts->quadraturets));
1964: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->quadraturets, (PetscObject)ts, 1));
1965: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%squad_", ((PetscObject)ts)->prefix ? ((PetscObject)ts)->prefix : ""));
1966: PetscCall(TSSetOptionsPrefix(ts->quadraturets, prefix));
1967: *quadts = ts->quadraturets;
1969: if (ts->numcost) {
1970: PetscCall(VecCreateSeq(PETSC_COMM_SELF, ts->numcost, &(*quadts)->vec_sol));
1971: } else {
1972: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &(*quadts)->vec_sol));
1973: }
1974: ts->costintegralfwd = fwd;
1975: PetscFunctionReturn(PETSC_SUCCESS);
1976: }
1978: /*@
1979: TSGetQuadratureTS - Return the sub-`TS` that evaluates integrals over time
1981: Input Parameter:
1982: . ts - the `TS` context obtained from `TSCreate()`
1984: Output Parameters:
1985: + fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1986: - quadts - the child `TS` context
1988: Level: intermediate
1990: .seealso: [](ch_ts), `TSCreateQuadratureTS()`
1991: @*/
1992: PetscErrorCode TSGetQuadratureTS(TS ts, PetscBool *fwd, TS *quadts)
1993: {
1994: PetscFunctionBegin;
1996: if (fwd) *fwd = ts->costintegralfwd;
1997: if (quadts) *quadts = ts->quadraturets;
1998: PetscFunctionReturn(PETSC_SUCCESS);
1999: }
2001: /*@
2002: TSComputeSNESJacobian - Compute the Jacobian needed for the `SNESSolve()` in `TS`
2004: Collective
2006: Input Parameters:
2007: + ts - the `TS` context obtained from `TSCreate()`
2008: - x - state vector
2010: Output Parameters:
2011: + J - Jacobian matrix
2012: - Jpre - matrix used to compute the preconditioner for `J` (may be same as `J`)
2014: Level: developer
2016: Note:
2017: Uses finite differencing when `TS` Jacobian is not available.
2019: .seealso: `SNES`, `TS`, `SNESSetJacobian()`, `TSSetRHSJacobian()`, `TSSetIJacobian()`
2020: @*/
2021: PetscErrorCode TSComputeSNESJacobian(TS ts, Vec x, Mat J, Mat Jpre)
2022: {
2023: SNES snes = ts->snes;
2024: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *) = NULL;
2026: PetscFunctionBegin;
2027: /*
2028: Unlike implicit methods, explicit methods do not have SNESMatFDColoring in the snes object
2029: because SNESSolve() has not been called yet; so querying SNESMatFDColoring does not work for
2030: explicit methods. Instead, we check the Jacobian compute function directly to determine if FD
2031: coloring is used.
2032: */
2033: PetscCall(SNESGetJacobian(snes, NULL, NULL, &jac, NULL));
2034: if (jac == SNESComputeJacobianDefaultColor) {
2035: Vec f;
2036: PetscCall(SNESSetSolution(snes, x));
2037: PetscCall(SNESGetFunction(snes, &f, NULL, NULL));
2038: /* Force MatFDColoringApply to evaluate the SNES residual function for the base vector */
2039: PetscCall(SNESComputeFunction(snes, x, f));
2040: }
2041: PetscCall(SNESComputeJacobian(snes, x, J, Jpre));
2042: PetscFunctionReturn(PETSC_SUCCESS);
2043: }