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 on TS
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: Notes:
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: `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 on TS
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: Notes:
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: `TSSetRHSJacobianP()`
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 on TS
90: Input Parameters:
91: . ts - The TS context obtained from TSCreate()
93: Level: developer
95: .seealso: `TSSetRHSJacobianP()`
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 on TS
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: Notes:
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:
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 on TS
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: `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 on TS
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: Notes:
231: For optimization there is usually a single cost function (numcost = 1). For sensitivities there may be multiple cost functions
233: .seealso: `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: `TSSetCostIntegrand()`
282: @*/
283: PetscErrorCode TSGetCostIntegral(TS ts, Vec *v)
284: {
285: TS quadts;
289: TSGetQuadratureTS(ts, NULL, &quadts);
290: *v = quadts->vec_sol;
291: return 0;
292: }
294: /*@C
295: TSComputeCostIntegrand - Evaluates the integral function in the cost functions.
297: Input Parameters:
298: + ts - the TS context
299: . t - current time
300: - U - state vector, i.e. current solution
302: Output Parameter:
303: . Q - vector of size numcost to hold the outputs
305: Notes:
306: Most users should not need to explicitly call this routine, as it
307: is used internally within the sensitivity analysis context.
309: Level: deprecated
311: .seealso: `TSSetCostIntegrand()`
312: @*/
313: PetscErrorCode TSComputeCostIntegrand(TS ts, PetscReal t, Vec U, Vec Q)
314: {
319: PetscLogEventBegin(TS_FunctionEval, ts, U, Q, 0);
320: if (ts->costintegrand) PetscCallBack("TS callback integrand in the cost function", (*ts->costintegrand)(ts, t, U, Q, ts->costintegrandctx));
321: else VecZeroEntries(Q);
322: PetscLogEventEnd(TS_FunctionEval, ts, U, Q, 0);
323: return 0;
324: }
326: /*@C
327: TSComputeDRDUFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobian()
329: Level: deprecated
331: @*/
332: PetscErrorCode TSComputeDRDUFunction(TS ts, PetscReal t, Vec U, Vec *DRDU)
333: {
334: if (!DRDU) return 0;
338: PetscCallBack("TS callback DRDU for sensitivity analysis", (*ts->drdufunction)(ts, t, U, DRDU, ts->costintegrandctx));
339: return 0;
340: }
342: /*@C
343: TSComputeDRDPFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobianP()
345: Level: deprecated
347: @*/
348: PetscErrorCode TSComputeDRDPFunction(TS ts, PetscReal t, Vec U, Vec *DRDP)
349: {
350: if (!DRDP) return 0;
354: PetscCallBack("TS callback DRDP for sensitivity analysis", (*ts->drdpfunction)(ts, t, U, DRDP, ts->costintegrandctx));
355: return 0;
356: }
358: /*@C
359: 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.
361: Logically Collective on TS
363: Input Parameters:
364: + ts - TS context obtained from TSCreate()
365: . ihp1 - an array of vectors storing the result of vector-Hessian-vector product for F_UU
366: . hessianproductfunc1 - vector-Hessian-vector product function for F_UU
367: . ihp2 - an array of vectors storing the result of vector-Hessian-vector product for F_UP
368: . hessianproductfunc2 - vector-Hessian-vector product function for F_UP
369: . ihp3 - an array of vectors storing the result of vector-Hessian-vector product for F_PU
370: . hessianproductfunc3 - vector-Hessian-vector product function for F_PU
371: . ihp4 - an array of vectors storing the result of vector-Hessian-vector product for F_PP
372: - hessianproductfunc4 - vector-Hessian-vector product function for F_PP
374: Calling sequence of ihessianproductfunc:
375: $ ihessianproductfunc (TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx);
376: + t - current timestep
377: . U - input vector (current ODE solution)
378: . Vl - an array of input vectors to be left-multiplied with the Hessian
379: . Vr - input vector to be right-multiplied with the Hessian
380: . VHV - an array of output vectors for vector-Hessian-vector product
381: - ctx - [optional] user-defined function context
383: Level: intermediate
385: Notes:
386: The first Hessian function and the working array are required.
387: As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
388: $ Vl_n^T*F_UP*Vr
389: 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.
390: Each entry of F_UP corresponds to the derivative
391: $ F_UP[i][j][k] = \frac{\partial^2 F[i]}{\partial U[j] \partial P[k]}.
392: 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
393: $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * F_UP[i][j][k] * Vr[k]}
394: If the cost function is a scalar, there will be only one vector in Vl and VHV.
396: .seealso:
397: @*/
398: 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)
399: {
403: ts->ihessianproductctx = ctx;
404: if (ihp1) ts->vecs_fuu = ihp1;
405: if (ihp2) ts->vecs_fup = ihp2;
406: if (ihp3) ts->vecs_fpu = ihp3;
407: if (ihp4) ts->vecs_fpp = ihp4;
408: ts->ihessianproduct_fuu = ihessianproductfunc1;
409: ts->ihessianproduct_fup = ihessianproductfunc2;
410: ts->ihessianproduct_fpu = ihessianproductfunc3;
411: ts->ihessianproduct_fpp = ihessianproductfunc4;
412: return 0;
413: }
415: /*@C
416: TSComputeIHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Fuu.
418: Collective on TS
420: Input Parameters:
421: . ts - The TS context obtained from TSCreate()
423: Notes:
424: TSComputeIHessianProductFunctionUU() is typically used for sensitivity implementation,
425: so most users would not generally call this routine themselves.
427: Level: developer
429: .seealso: `TSSetIHessianProduct()`
430: @*/
431: PetscErrorCode TSComputeIHessianProductFunctionUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
432: {
433: if (!VHV) return 0;
437: if (ts->ihessianproduct_fuu) PetscCallBack("TS callback IHessianProduct 1 for sensitivity analysis", (*ts->ihessianproduct_fuu)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
439: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
440: if (ts->rhshessianproduct_guu) {
441: PetscInt nadj;
442: TSComputeRHSHessianProductFunctionUU(ts, t, U, Vl, Vr, VHV);
443: for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
444: }
445: return 0;
446: }
448: /*@C
449: TSComputeIHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Fup.
451: Collective on TS
453: Input Parameters:
454: . ts - The TS context obtained from TSCreate()
456: Notes:
457: TSComputeIHessianProductFunctionUP() is typically used for sensitivity implementation,
458: so most users would not generally call this routine themselves.
460: Level: developer
462: .seealso: `TSSetIHessianProduct()`
463: @*/
464: PetscErrorCode TSComputeIHessianProductFunctionUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
465: {
466: if (!VHV) return 0;
470: if (ts->ihessianproduct_fup) PetscCallBack("TS callback IHessianProduct 2 for sensitivity analysis", (*ts->ihessianproduct_fup)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
472: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
473: if (ts->rhshessianproduct_gup) {
474: PetscInt nadj;
475: TSComputeRHSHessianProductFunctionUP(ts, t, U, Vl, Vr, VHV);
476: for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
477: }
478: return 0;
479: }
481: /*@C
482: TSComputeIHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Fpu.
484: Collective on TS
486: Input Parameters:
487: . ts - The TS context obtained from TSCreate()
489: Notes:
490: TSComputeIHessianProductFunctionPU() is typically used for sensitivity implementation,
491: so most users would not generally call this routine themselves.
493: Level: developer
495: .seealso: `TSSetIHessianProduct()`
496: @*/
497: PetscErrorCode TSComputeIHessianProductFunctionPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
498: {
499: if (!VHV) return 0;
503: if (ts->ihessianproduct_fpu) PetscCallBack("TS callback IHessianProduct 3 for sensitivity analysis", (*ts->ihessianproduct_fpu)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
505: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
506: if (ts->rhshessianproduct_gpu) {
507: PetscInt nadj;
508: TSComputeRHSHessianProductFunctionPU(ts, t, U, Vl, Vr, VHV);
509: for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
510: }
511: return 0;
512: }
514: /*@C
515: TSComputeIHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Fpp.
517: Collective on TS
519: Input Parameters:
520: . ts - The TS context obtained from TSCreate()
522: Notes:
523: TSComputeIHessianProductFunctionPP() is typically used for sensitivity implementation,
524: so most users would not generally call this routine themselves.
526: Level: developer
528: .seealso: `TSSetIHessianProduct()`
529: @*/
530: PetscErrorCode TSComputeIHessianProductFunctionPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
531: {
532: if (!VHV) return 0;
536: if (ts->ihessianproduct_fpp) PetscCallBack("TS callback IHessianProduct 3 for sensitivity analysis", (*ts->ihessianproduct_fpp)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
538: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
539: if (ts->rhshessianproduct_gpp) {
540: PetscInt nadj;
541: TSComputeRHSHessianProductFunctionPP(ts, t, U, Vl, Vr, VHV);
542: for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
543: }
544: return 0;
545: }
547: /*@C
548: 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.
550: Logically Collective on TS
552: Input Parameters:
553: + ts - TS context obtained from TSCreate()
554: . rhshp1 - an array of vectors storing the result of vector-Hessian-vector product for G_UU
555: . hessianproductfunc1 - vector-Hessian-vector product function for G_UU
556: . rhshp2 - an array of vectors storing the result of vector-Hessian-vector product for G_UP
557: . hessianproductfunc2 - vector-Hessian-vector product function for G_UP
558: . rhshp3 - an array of vectors storing the result of vector-Hessian-vector product for G_PU
559: . hessianproductfunc3 - vector-Hessian-vector product function for G_PU
560: . rhshp4 - an array of vectors storing the result of vector-Hessian-vector product for G_PP
561: - hessianproductfunc4 - vector-Hessian-vector product function for G_PP
563: Calling sequence of ihessianproductfunc:
564: $ rhshessianproductfunc (TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx);
565: + t - current timestep
566: . U - input vector (current ODE solution)
567: . Vl - an array of input vectors to be left-multiplied with the Hessian
568: . Vr - input vector to be right-multiplied with the Hessian
569: . VHV - an array of output vectors for vector-Hessian-vector product
570: - ctx - [optional] user-defined function context
572: Level: intermediate
574: Notes:
575: The first Hessian function and the working array are required.
576: As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
577: $ Vl_n^T*G_UP*Vr
578: 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.
579: Each entry of G_UP corresponds to the derivative
580: $ G_UP[i][j][k] = \frac{\partial^2 G[i]}{\partial U[j] \partial P[k]}.
581: The result of the vector-Hessian-vector product for Vl_n needs to be stored in vector VHV_n with j-th entry being
582: $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * G_UP[i][j][k] * Vr[k]}
583: If the cost function is a scalar, there will be only one vector in Vl and VHV.
585: .seealso:
586: @*/
587: 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)
588: {
592: ts->rhshessianproductctx = ctx;
593: if (rhshp1) ts->vecs_guu = rhshp1;
594: if (rhshp2) ts->vecs_gup = rhshp2;
595: if (rhshp3) ts->vecs_gpu = rhshp3;
596: if (rhshp4) ts->vecs_gpp = rhshp4;
597: ts->rhshessianproduct_guu = rhshessianproductfunc1;
598: ts->rhshessianproduct_gup = rhshessianproductfunc2;
599: ts->rhshessianproduct_gpu = rhshessianproductfunc3;
600: ts->rhshessianproduct_gpp = rhshessianproductfunc4;
601: return 0;
602: }
604: /*@C
605: TSComputeRHSHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Guu.
607: Collective on TS
609: Input Parameters:
610: . ts - The TS context obtained from TSCreate()
612: Notes:
613: TSComputeRHSHessianProductFunctionUU() is typically used for sensitivity implementation,
614: so most users would not generally call this routine themselves.
616: Level: developer
618: .seealso: `TSSetRHSHessianProduct()`
619: @*/
620: PetscErrorCode TSComputeRHSHessianProductFunctionUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
621: {
622: if (!VHV) return 0;
626: PetscCallBack("TS callback RHSHessianProduct 1 for sensitivity analysis", (*ts->rhshessianproduct_guu)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
627: return 0;
628: }
630: /*@C
631: TSComputeRHSHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Gup.
633: Collective on TS
635: Input Parameters:
636: . ts - The TS context obtained from TSCreate()
638: Notes:
639: TSComputeRHSHessianProductFunctionUP() is typically used for sensitivity implementation,
640: so most users would not generally call this routine themselves.
642: Level: developer
644: .seealso: `TSSetRHSHessianProduct()`
645: @*/
646: PetscErrorCode TSComputeRHSHessianProductFunctionUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
647: {
648: if (!VHV) return 0;
652: PetscCallBack("TS callback RHSHessianProduct 2 for sensitivity analysis", (*ts->rhshessianproduct_gup)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
653: return 0;
654: }
656: /*@C
657: TSComputeRHSHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Gpu.
659: Collective on TS
661: Input Parameters:
662: . ts - The TS context obtained from TSCreate()
664: Notes:
665: TSComputeRHSHessianProductFunctionPU() is typically used for sensitivity implementation,
666: so most users would not generally call this routine themselves.
668: Level: developer
670: .seealso: `TSSetRHSHessianProduct()`
671: @*/
672: PetscErrorCode TSComputeRHSHessianProductFunctionPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
673: {
674: if (!VHV) return 0;
678: PetscCallBack("TS callback RHSHessianProduct 3 for sensitivity analysis", (*ts->rhshessianproduct_gpu)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
679: return 0;
680: }
682: /*@C
683: TSComputeRHSHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Gpp.
685: Collective on TS
687: Input Parameters:
688: . ts - The TS context obtained from TSCreate()
690: Notes:
691: TSComputeRHSHessianProductFunctionPP() is typically used for sensitivity implementation,
692: so most users would not generally call this routine themselves.
694: Level: developer
696: .seealso: `TSSetRHSHessianProduct()`
697: @*/
698: PetscErrorCode TSComputeRHSHessianProductFunctionPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
699: {
700: if (!VHV) return 0;
704: PetscCallBack("TS callback RHSHessianProduct 3 for sensitivity analysis", (*ts->rhshessianproduct_gpp)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
705: return 0;
706: }
708: /* --------------------------- Adjoint sensitivity ---------------------------*/
710: /*@
711: TSSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial values and w.r.t. the problem parameters
712: for use by the TSAdjoint routines.
714: Logically Collective on TS
716: Input Parameters:
717: + ts - the TS context obtained from TSCreate()
718: . numcost - number of gradients to be computed, this is the number of cost functions
719: . 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
720: - mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
722: Level: beginner
724: Notes:
725: the entries in these vectors must be correctly initialized with the values lamda_i = df/dy|finaltime mu_i = df/dp|finaltime
727: After TSAdjointSolve() is called the lamba and the mu contain the computed sensitivities
729: .seealso `TSGetCostGradients()`
730: @*/
731: PetscErrorCode TSSetCostGradients(TS ts, PetscInt numcost, Vec *lambda, Vec *mu)
732: {
735: ts->vecs_sensi = lambda;
736: ts->vecs_sensip = mu;
738: ts->numcost = numcost;
739: return 0;
740: }
742: /*@
743: TSGetCostGradients - Returns the gradients from the TSAdjointSolve()
745: Not Collective, but Vec returned is parallel if TS is parallel
747: Input Parameter:
748: . ts - the TS context obtained from TSCreate()
750: Output Parameters:
751: + numcost - size of returned arrays
752: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
753: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
755: Level: intermediate
757: .seealso: `TSSetCostGradients()`
758: @*/
759: PetscErrorCode TSGetCostGradients(TS ts, PetscInt *numcost, Vec **lambda, Vec **mu)
760: {
762: if (numcost) *numcost = ts->numcost;
763: if (lambda) *lambda = ts->vecs_sensi;
764: if (mu) *mu = ts->vecs_sensip;
765: return 0;
766: }
768: /*@
769: 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
770: for use by the TSAdjoint routines.
772: Logically Collective on TS
774: Input Parameters:
775: + ts - the TS context obtained from TSCreate()
776: . numcost - number of cost functions
777: . 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
778: . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
779: - dir - the direction vector that are multiplied with the Hessian of the cost functions
781: Level: beginner
783: Notes: 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 lamba2 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.
790: .seealso: `TSAdjointSetForward()`
791: @*/
792: PetscErrorCode TSSetCostHessianProducts(TS ts, PetscInt numcost, Vec *lambda2, Vec *mu2, Vec dir)
793: {
796: ts->numcost = numcost;
797: ts->vecs_sensi2 = lambda2;
798: ts->vecs_sensi2p = mu2;
799: ts->vec_dir = dir;
800: return 0;
801: }
803: /*@
804: TSGetCostHessianProducts - Returns the gradients from the TSAdjointSolve()
806: Not Collective, but Vec returned is parallel if TS is parallel
808: Input Parameter:
809: . ts - the TS context obtained from TSCreate()
811: Output Parameters:
812: + numcost - number of cost functions
813: . 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
814: . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
815: - dir - the direction vector that are multiplied with the Hessian of the cost functions
817: Level: intermediate
819: .seealso: `TSSetCostHessianProducts()`
820: @*/
821: PetscErrorCode TSGetCostHessianProducts(TS ts, PetscInt *numcost, Vec **lambda2, Vec **mu2, Vec *dir)
822: {
824: if (numcost) *numcost = ts->numcost;
825: if (lambda2) *lambda2 = ts->vecs_sensi2;
826: if (mu2) *mu2 = ts->vecs_sensi2p;
827: if (dir) *dir = ts->vec_dir;
828: return 0;
829: }
831: /*@
832: TSAdjointSetForward - Trigger the tangent linear solver and initialize the forward sensitivities
834: Logically Collective on TS
836: Input Parameters:
837: + ts - the TS context obtained from TSCreate()
838: - didp - the derivative of initial values w.r.t. parameters
840: Level: intermediate
842: Notes: When computing sensitivies w.r.t. initial condition, set didp to NULL so that the solver will take it as an identity matrix mathematically. TSAdjoint does not reset the tangent linear solver automatically, TSAdjointResetForward() should be called to reset the tangent linear solver.
844: .seealso: `TSSetCostHessianProducts()`, `TSAdjointResetForward()`
845: @*/
846: PetscErrorCode TSAdjointSetForward(TS ts, Mat didp)
847: {
848: Mat A;
849: Vec sp;
850: PetscScalar *xarr;
851: PetscInt lsize;
853: ts->forward_solve = PETSC_TRUE; /* turn on tangent linear mode */
856: /* create a single-column dense matrix */
857: VecGetLocalSize(ts->vec_sol, &lsize);
858: MatCreateDense(PetscObjectComm((PetscObject)ts), lsize, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, &A);
860: VecDuplicate(ts->vec_sol, &sp);
861: MatDenseGetColumn(A, 0, &xarr);
862: VecPlaceArray(sp, xarr);
863: if (ts->vecs_sensi2p) { /* tangent linear variable initialized as 2*dIdP*dir */
864: if (didp) {
865: MatMult(didp, ts->vec_dir, sp);
866: VecScale(sp, 2.);
867: } else {
868: VecZeroEntries(sp);
869: }
870: } else { /* tangent linear variable initialized as dir */
871: VecCopy(ts->vec_dir, sp);
872: }
873: VecResetArray(sp);
874: MatDenseRestoreColumn(A, &xarr);
875: VecDestroy(&sp);
877: TSForwardSetInitialSensitivities(ts, A); /* if didp is NULL, identity matrix is assumed */
879: MatDestroy(&A);
880: return 0;
881: }
883: /*@
884: TSAdjointResetForward - Reset the tangent linear solver and destroy the tangent linear context
886: Logically Collective on TS
888: Input Parameters:
889: . ts - the TS context obtained from TSCreate()
891: Level: intermediate
893: .seealso: `TSAdjointSetForward()`
894: @*/
895: PetscErrorCode TSAdjointResetForward(TS ts)
896: {
897: ts->forward_solve = PETSC_FALSE; /* turn off tangent linear mode */
898: TSForwardReset(ts);
899: return 0;
900: }
902: /*@
903: TSAdjointSetUp - Sets up the internal data structures for the later use
904: of an adjoint solver
906: Collective on TS
908: Input Parameter:
909: . ts - the TS context obtained from TSCreate()
911: Level: advanced
913: .seealso: `TSCreate()`, `TSAdjointStep()`, `TSSetCostGradients()`
914: @*/
915: PetscErrorCode TSAdjointSetUp(TS ts)
916: {
917: TSTrajectory tj;
918: PetscBool match;
921: if (ts->adjointsetupcalled) return 0;
924: TSGetTrajectory(ts, &tj);
925: PetscObjectTypeCompare((PetscObject)tj, TSTRAJECTORYBASIC, &match);
926: if (match) {
927: PetscBool solution_only;
928: TSTrajectoryGetSolutionOnly(tj, &solution_only);
930: }
931: TSTrajectorySetUseHistory(tj, PETSC_FALSE); /* not use TSHistory */
933: if (ts->quadraturets) { /* if there is integral in the cost function */
934: VecDuplicate(ts->vecs_sensi[0], &ts->vec_drdu_col);
935: if (ts->vecs_sensip) VecDuplicate(ts->vecs_sensip[0], &ts->vec_drdp_col);
936: }
938: PetscTryTypeMethod(ts, adjointsetup);
939: ts->adjointsetupcalled = PETSC_TRUE;
940: return 0;
941: }
943: /*@
944: TSAdjointReset - Resets a TSAdjoint context and removes any allocated Vecs and Mats.
946: Collective on TS
948: Input Parameter:
949: . ts - the TS context obtained from TSCreate()
951: Level: beginner
953: .seealso: `TSCreate()`, `TSAdjointSetUp()`, `TSADestroy()`
954: @*/
955: PetscErrorCode TSAdjointReset(TS ts)
956: {
958: PetscTryTypeMethod(ts, adjointreset);
959: if (ts->quadraturets) { /* if there is integral in the cost function */
960: VecDestroy(&ts->vec_drdu_col);
961: if (ts->vecs_sensip) VecDestroy(&ts->vec_drdp_col);
962: }
963: ts->vecs_sensi = NULL;
964: ts->vecs_sensip = NULL;
965: ts->vecs_sensi2 = NULL;
966: ts->vecs_sensi2p = NULL;
967: ts->vec_dir = NULL;
968: ts->adjointsetupcalled = PETSC_FALSE;
969: return 0;
970: }
972: /*@
973: TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time
975: Logically Collective on TS
977: Input Parameters:
978: + ts - the TS context obtained from TSCreate()
979: - steps - number of steps to use
981: Level: intermediate
983: Notes:
984: Normally one does not call this and TSAdjointSolve() integrates back to the original timestep. One can call this
985: so as to integrate back to less than the original timestep
987: .seealso: `TSSetExactFinalTime()`
988: @*/
989: PetscErrorCode TSAdjointSetSteps(TS ts, PetscInt steps)
990: {
995: ts->adjoint_max_steps = steps;
996: return 0;
997: }
999: /*@C
1000: TSAdjointSetRHSJacobian - Deprecated, use TSSetRHSJacobianP()
1002: Level: deprecated
1004: @*/
1005: PetscErrorCode TSAdjointSetRHSJacobian(TS ts, Mat Amat, PetscErrorCode (*func)(TS, PetscReal, Vec, Mat, void *), void *ctx)
1006: {
1010: ts->rhsjacobianp = func;
1011: ts->rhsjacobianpctx = ctx;
1012: if (Amat) {
1013: PetscObjectReference((PetscObject)Amat);
1014: MatDestroy(&ts->Jacp);
1015: ts->Jacp = Amat;
1016: }
1017: return 0;
1018: }
1020: /*@C
1021: TSAdjointComputeRHSJacobian - Deprecated, use TSComputeRHSJacobianP()
1023: Level: deprecated
1025: @*/
1026: PetscErrorCode TSAdjointComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat Amat)
1027: {
1032: PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->rhsjacobianp)(ts, t, U, Amat, ts->rhsjacobianpctx));
1033: return 0;
1034: }
1036: /*@
1037: TSAdjointComputeDRDYFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobian()
1039: Level: deprecated
1041: @*/
1042: PetscErrorCode TSAdjointComputeDRDYFunction(TS ts, PetscReal t, Vec U, Vec *DRDU)
1043: {
1047: PetscCallBack("TS callback DRDY for sensitivity analysis", (*ts->drdufunction)(ts, t, U, DRDU, ts->costintegrandctx));
1048: return 0;
1049: }
1051: /*@
1052: TSAdjointComputeDRDPFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobianP()
1054: Level: deprecated
1056: @*/
1057: PetscErrorCode TSAdjointComputeDRDPFunction(TS ts, PetscReal t, Vec U, Vec *DRDP)
1058: {
1062: PetscCallBack("TS callback DRDP for sensitivity analysis", (*ts->drdpfunction)(ts, t, U, DRDP, ts->costintegrandctx));
1063: return 0;
1064: }
1066: /*@C
1067: TSAdjointMonitorSensi - monitors the first lambda sensitivity
1069: Level: intermediate
1071: .seealso: `TSAdjointMonitorSet()`
1072: @*/
1073: PetscErrorCode TSAdjointMonitorSensi(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscInt numcost, Vec *lambda, Vec *mu, PetscViewerAndFormat *vf)
1074: {
1075: PetscViewer viewer = vf->viewer;
1078: PetscViewerPushFormat(viewer, vf->format);
1079: VecView(lambda[0], viewer);
1080: PetscViewerPopFormat(viewer);
1081: return 0;
1082: }
1084: /*@C
1085: TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
1087: Collective on TS
1089: Input Parameters:
1090: + ts - TS object you wish to monitor
1091: . name - the monitor type one is seeking
1092: . help - message indicating what monitoring is done
1093: . manual - manual page for the monitor
1094: . monitor - the monitor function
1095: - 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
1097: Level: developer
1099: .seealso: `PetscOptionsGetViewer()`, `PetscOptionsGetReal()`, `PetscOptionsHasName()`, `PetscOptionsGetString()`,
1100: `PetscOptionsGetIntArray()`, `PetscOptionsGetRealArray()`, `PetscOptionsBool()`
1101: `PetscOptionsInt()`, `PetscOptionsString()`, `PetscOptionsReal()`, `PetscOptionsBool()`,
1102: `PetscOptionsName()`, `PetscOptionsBegin()`, `PetscOptionsEnd()`, `PetscOptionsHeadBegin()`,
1103: `PetscOptionsStringArray()`, `PetscOptionsRealArray()`, `PetscOptionsScalar()`,
1104: `PetscOptionsBoolGroupBegin()`, `PetscOptionsBoolGroup()`, `PetscOptionsBoolGroupEnd()`,
1105: `PetscOptionsFList()`, `PetscOptionsEList()`
1106: @*/
1107: 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 *))
1108: {
1109: PetscViewer viewer;
1110: PetscViewerFormat format;
1111: PetscBool flg;
1113: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, name, &viewer, &format, &flg);
1114: if (flg) {
1115: PetscViewerAndFormat *vf;
1116: PetscViewerAndFormatCreate(viewer, format, &vf);
1117: PetscObjectDereference((PetscObject)viewer);
1118: if (monitorsetup) (*monitorsetup)(ts, vf);
1119: TSAdjointMonitorSet(ts, (PetscErrorCode(*)(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, void *))monitor, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1120: }
1121: return 0;
1122: }
1124: /*@C
1125: TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
1126: timestep to display the iteration's progress.
1128: Logically Collective on TS
1130: Input Parameters:
1131: + ts - the TS context obtained from TSCreate()
1132: . adjointmonitor - monitoring routine
1133: . adjointmctx - [optional] user-defined context for private data for the
1134: monitor routine (use NULL if no context is desired)
1135: - adjointmonitordestroy - [optional] routine that frees monitor context
1136: (may be NULL)
1138: Calling sequence of monitor:
1139: $ int adjointmonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *adjointmctx)
1141: + ts - the TS context
1142: . 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
1143: been interpolated to)
1144: . time - current time
1145: . u - current iterate
1146: . numcost - number of cost functionos
1147: . lambda - sensitivities to initial conditions
1148: . mu - sensitivities to parameters
1149: - adjointmctx - [optional] adjoint monitoring context
1151: Notes:
1152: This routine adds an additional monitor to the list of monitors that
1153: already has been loaded.
1155: Fortran Notes:
1156: Only a single monitor function can be set for each TS object
1158: Level: intermediate
1160: .seealso: `TSAdjointMonitorCancel()`
1161: @*/
1162: PetscErrorCode TSAdjointMonitorSet(TS ts, PetscErrorCode (*adjointmonitor)(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, void *), void *adjointmctx, PetscErrorCode (*adjointmdestroy)(void **))
1163: {
1164: PetscInt i;
1165: PetscBool identical;
1168: for (i = 0; i < ts->numbermonitors; i++) {
1169: PetscMonitorCompare((PetscErrorCode(*)(void))adjointmonitor, adjointmctx, adjointmdestroy, (PetscErrorCode(*)(void))ts->adjointmonitor[i], ts->adjointmonitorcontext[i], ts->adjointmonitordestroy[i], &identical);
1170: if (identical) return 0;
1171: }
1173: ts->adjointmonitor[ts->numberadjointmonitors] = adjointmonitor;
1174: ts->adjointmonitordestroy[ts->numberadjointmonitors] = adjointmdestroy;
1175: ts->adjointmonitorcontext[ts->numberadjointmonitors++] = (void *)adjointmctx;
1176: return 0;
1177: }
1179: /*@C
1180: TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.
1182: Logically Collective on TS
1184: Input Parameters:
1185: . ts - the TS context obtained from TSCreate()
1187: Notes:
1188: There is no way to remove a single, specific monitor.
1190: Level: intermediate
1192: .seealso: `TSAdjointMonitorSet()`
1193: @*/
1194: PetscErrorCode TSAdjointMonitorCancel(TS ts)
1195: {
1196: PetscInt i;
1199: for (i = 0; i < ts->numberadjointmonitors; i++) {
1200: if (ts->adjointmonitordestroy[i]) (*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]);
1201: }
1202: ts->numberadjointmonitors = 0;
1203: return 0;
1204: }
1206: /*@C
1207: TSAdjointMonitorDefault - the default monitor of adjoint computations
1209: Level: intermediate
1211: .seealso: `TSAdjointMonitorSet()`
1212: @*/
1213: PetscErrorCode TSAdjointMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscInt numcost, Vec *lambda, Vec *mu, PetscViewerAndFormat *vf)
1214: {
1215: PetscViewer viewer = vf->viewer;
1218: PetscViewerPushFormat(viewer, vf->format);
1219: PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel);
1220: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " TS dt %g time %g%s", step, (double)ts->time_step, (double)ptime, ts->steprollback ? " (r)\n" : "\n");
1221: PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel);
1222: PetscViewerPopFormat(viewer);
1223: return 0;
1224: }
1226: /*@C
1227: TSAdjointMonitorDrawSensi - Monitors progress of the adjoint TS solvers by calling
1228: VecView() for the sensitivities to initial states at each timestep
1230: Collective on TS
1232: Input Parameters:
1233: + ts - the TS context
1234: . step - current time-step
1235: . ptime - current time
1236: . u - current state
1237: . numcost - number of cost functions
1238: . lambda - sensitivities to initial conditions
1239: . mu - sensitivities to parameters
1240: - dummy - either a viewer or NULL
1242: Level: intermediate
1244: .seealso: `TSAdjointMonitorSet()`, `TSAdjointMonitorDefault()`, `VecView()`
1245: @*/
1246: PetscErrorCode TSAdjointMonitorDrawSensi(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscInt numcost, Vec *lambda, Vec *mu, void *dummy)
1247: {
1248: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
1249: PetscDraw draw;
1250: PetscReal xl, yl, xr, yr, h;
1251: char time[32];
1253: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return 0;
1255: VecView(lambda[0], ictx->viewer);
1256: PetscViewerDrawGetDraw(ictx->viewer, 0, &draw);
1257: PetscSNPrintf(time, 32, "Timestep %d Time %g", (int)step, (double)ptime);
1258: PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr);
1259: h = yl + .95 * (yr - yl);
1260: PetscDrawStringCentered(draw, .5 * (xl + xr), h, PETSC_DRAW_BLACK, time);
1261: PetscDrawFlush(draw);
1262: return 0;
1263: }
1265: /*
1266: TSAdjointSetFromOptions - Sets various TSAdjoint parameters from user options.
1268: Collective on TSAdjoint
1270: Input Parameter:
1271: . ts - the TS context
1273: Options Database Keys:
1274: + -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
1275: . -ts_adjoint_monitor - print information at each adjoint time step
1276: - -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically
1278: Level: developer
1280: Notes:
1281: This is not normally called directly by users
1283: .seealso: `TSSetSaveTrajectory()`, `TSTrajectorySetUp()`
1284: */
1285: PetscErrorCode TSAdjointSetFromOptions(TS ts, PetscOptionItems *PetscOptionsObject)
1286: {
1287: PetscBool tflg, opt;
1290: PetscOptionsHeadBegin(PetscOptionsObject, "TS Adjoint options");
1291: tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
1292: PetscOptionsBool("-ts_adjoint_solve", "Solve the adjoint problem immediately after solving the forward problem", "", tflg, &tflg, &opt);
1293: if (opt) {
1294: TSSetSaveTrajectory(ts);
1295: ts->adjoint_solve = tflg;
1296: }
1297: TSAdjointMonitorSetFromOptions(ts, "-ts_adjoint_monitor", "Monitor adjoint timestep size", "TSAdjointMonitorDefault", TSAdjointMonitorDefault, NULL);
1298: TSAdjointMonitorSetFromOptions(ts, "-ts_adjoint_monitor_sensi", "Monitor sensitivity in the adjoint computation", "TSAdjointMonitorSensi", TSAdjointMonitorSensi, NULL);
1299: opt = PETSC_FALSE;
1300: PetscOptionsName("-ts_adjoint_monitor_draw_sensi", "Monitor adjoint sensitivities (lambda only) graphically", "TSAdjointMonitorDrawSensi", &opt);
1301: if (opt) {
1302: TSMonitorDrawCtx ctx;
1303: PetscInt howoften = 1;
1305: PetscOptionsInt("-ts_adjoint_monitor_draw_sensi", "Monitor adjoint sensitivities (lambda only) graphically", "TSAdjointMonitorDrawSensi", howoften, &howoften, NULL);
1306: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
1307: TSAdjointMonitorSet(ts, TSAdjointMonitorDrawSensi, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
1308: }
1309: return 0;
1310: }
1312: /*@
1313: TSAdjointStep - Steps one time step backward in the adjoint run
1315: Collective on TS
1317: Input Parameter:
1318: . ts - the TS context obtained from TSCreate()
1320: Level: intermediate
1322: .seealso: `TSAdjointSetUp()`, `TSAdjointSolve()`
1323: @*/
1324: PetscErrorCode TSAdjointStep(TS ts)
1325: {
1326: DM dm;
1329: TSGetDM(ts, &dm);
1330: TSAdjointSetUp(ts);
1331: ts->steps--; /* must decrease the step index before the adjoint step is taken. */
1333: ts->reason = TS_CONVERGED_ITERATING;
1334: ts->ptime_prev = ts->ptime;
1335: PetscLogEventBegin(TS_AdjointStep, ts, 0, 0, 0);
1336: PetscUseTypeMethod(ts, adjointstep);
1337: PetscLogEventEnd(TS_AdjointStep, ts, 0, 0, 0);
1338: ts->adjoint_steps++;
1340: if (ts->reason < 0) {
1342: } else if (!ts->reason) {
1343: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1344: }
1345: return 0;
1346: }
1348: /*@
1349: TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE
1351: Collective on TS
1353: Input Parameter:
1354: . ts - the TS context obtained from TSCreate()
1356: Options Database:
1357: . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial values
1359: Level: intermediate
1361: Notes:
1362: This must be called after a call to TSSolve() that solves the forward problem
1364: By default this will integrate back to the initial time, one can use TSAdjointSetSteps() to step back to a later time
1366: .seealso: `TSCreate()`, `TSSetCostGradients()`, `TSSetSolution()`, `TSAdjointStep()`
1367: @*/
1368: PetscErrorCode TSAdjointSolve(TS ts)
1369: {
1370: static PetscBool cite = PETSC_FALSE;
1371: #if defined(TSADJOINT_STAGE)
1372: PetscLogStage adjoint_stage;
1373: #endif
1376: PetscCall(PetscCitationsRegister("@article{Zhang2022tsadjoint,\n"
1377: " title = {{PETSc TSAdjoint: A Discrete Adjoint ODE Solver for First-Order and Second-Order Sensitivity Analysis}},\n"
1378: " author = {Zhang, Hong and Constantinescu, Emil M. and Smith, Barry F.},\n"
1379: " journal = {SIAM Journal on Scientific Computing},\n"
1380: " volume = {44},\n"
1381: " number = {1},\n"
1382: " pages = {C1-C24},\n"
1383: " doi = {10.1137/21M140078X},\n"
1384: " year = {2022}\n}\n",
1385: &cite));
1386: #if defined(TSADJOINT_STAGE)
1387: PetscLogStageRegister("TSAdjoint", &adjoint_stage);
1388: PetscLogStagePush(adjoint_stage);
1389: #endif
1390: TSAdjointSetUp(ts);
1392: /* reset time step and iteration counters */
1393: ts->adjoint_steps = 0;
1394: ts->ksp_its = 0;
1395: ts->snes_its = 0;
1396: ts->num_snes_failures = 0;
1397: ts->reject = 0;
1398: ts->reason = TS_CONVERGED_ITERATING;
1400: if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->steps;
1401: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1403: while (!ts->reason) {
1404: TSTrajectoryGet(ts->trajectory, ts, ts->steps, &ts->ptime);
1405: TSAdjointMonitor(ts, ts->steps, ts->ptime, ts->vec_sol, ts->numcost, ts->vecs_sensi, ts->vecs_sensip);
1406: TSAdjointEventHandler(ts);
1407: TSAdjointStep(ts);
1408: if ((ts->vec_costintegral || ts->quadraturets) && !ts->costintegralfwd) TSAdjointCostIntegral(ts);
1409: }
1410: if (!ts->steps) {
1411: TSTrajectoryGet(ts->trajectory, ts, ts->steps, &ts->ptime);
1412: TSAdjointMonitor(ts, ts->steps, ts->ptime, ts->vec_sol, ts->numcost, ts->vecs_sensi, ts->vecs_sensip);
1413: }
1414: ts->solvetime = ts->ptime;
1415: TSTrajectoryViewFromOptions(ts->trajectory, NULL, "-ts_trajectory_view");
1416: VecViewFromOptions(ts->vecs_sensi[0], (PetscObject)ts, "-ts_adjoint_view_solution");
1417: ts->adjoint_max_steps = 0;
1418: #if defined(TSADJOINT_STAGE)
1419: PetscLogStagePop();
1420: #endif
1421: return 0;
1422: }
1424: /*@C
1425: TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using TSAdjointMonitorSet()
1427: Collective on TS
1429: Input Parameters:
1430: + ts - time stepping context obtained from TSCreate()
1431: . step - step number that has just completed
1432: . ptime - model time of the state
1433: . u - state at the current model time
1434: . numcost - number of cost functions (dimension of lambda or mu)
1435: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
1436: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
1438: Notes:
1439: TSAdjointMonitor() is typically used automatically within the time stepping implementations.
1440: Users would almost never call this routine directly.
1442: Level: developer
1444: @*/
1445: PetscErrorCode TSAdjointMonitor(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscInt numcost, Vec *lambda, Vec *mu)
1446: {
1447: PetscInt i, n = ts->numberadjointmonitors;
1451: VecLockReadPush(u);
1452: for (i = 0; i < n; i++) (*ts->adjointmonitor[i])(ts, step, ptime, u, numcost, lambda, mu, ts->adjointmonitorcontext[i]);
1453: VecLockReadPop(u);
1454: return 0;
1455: }
1457: /*@
1458: TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.
1460: Collective on TS
1462: Input Parameter:
1463: . ts - time stepping context
1465: Level: advanced
1467: Notes:
1468: This function cannot be called until TSAdjointStep() has been completed.
1470: .seealso: `TSAdjointSolve()`, `TSAdjointStep`
1471: @*/
1472: PetscErrorCode TSAdjointCostIntegral(TS ts)
1473: {
1475: PetscUseTypeMethod(ts, adjointintegral);
1476: return 0;
1477: }
1479: /* ------------------ Forward (tangent linear) sensitivity ------------------*/
1481: /*@
1482: TSForwardSetUp - Sets up the internal data structures for the later use
1483: of forward sensitivity analysis
1485: Collective on TS
1487: Input Parameter:
1488: . ts - the TS context obtained from TSCreate()
1490: Level: advanced
1492: .seealso: `TSCreate()`, `TSDestroy()`, `TSSetUp()`
1493: @*/
1494: PetscErrorCode TSForwardSetUp(TS ts)
1495: {
1497: if (ts->forwardsetupcalled) return 0;
1498: PetscTryTypeMethod(ts, forwardsetup);
1499: VecDuplicate(ts->vec_sol, &ts->vec_sensip_col);
1500: ts->forwardsetupcalled = PETSC_TRUE;
1501: return 0;
1502: }
1504: /*@
1505: TSForwardReset - Reset the internal data structures used by forward sensitivity analysis
1507: Collective on TS
1509: Input Parameter:
1510: . ts - the TS context obtained from TSCreate()
1512: Level: advanced
1514: .seealso: `TSCreate()`, `TSDestroy()`, `TSForwardSetUp()`
1515: @*/
1516: PetscErrorCode TSForwardReset(TS ts)
1517: {
1518: TS quadts = ts->quadraturets;
1521: PetscTryTypeMethod(ts, forwardreset);
1522: MatDestroy(&ts->mat_sensip);
1523: if (quadts) MatDestroy(&quadts->mat_sensip);
1524: VecDestroy(&ts->vec_sensip_col);
1525: ts->forward_solve = PETSC_FALSE;
1526: ts->forwardsetupcalled = PETSC_FALSE;
1527: return 0;
1528: }
1530: /*@
1531: TSForwardSetIntegralGradients - Set the vectors holding forward sensitivities of the integral term.
1533: Input Parameters:
1534: + ts - the TS context obtained from TSCreate()
1535: . numfwdint - number of integrals
1536: - vp - the vectors containing the gradients for each integral w.r.t. parameters
1538: Level: deprecated
1540: .seealso: `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1541: @*/
1542: PetscErrorCode TSForwardSetIntegralGradients(TS ts, PetscInt numfwdint, Vec *vp)
1543: {
1546: if (!ts->numcost) ts->numcost = numfwdint;
1548: ts->vecs_integral_sensip = vp;
1549: return 0;
1550: }
1552: /*@
1553: TSForwardGetIntegralGradients - Returns the forward sensitivities ofthe integral term.
1555: Input Parameter:
1556: . ts - the TS context obtained from TSCreate()
1558: Output Parameter:
1559: . vp - the vectors containing the gradients for each integral w.r.t. parameters
1561: Level: deprecated
1563: .seealso: `TSForwardSetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1564: @*/
1565: PetscErrorCode TSForwardGetIntegralGradients(TS ts, PetscInt *numfwdint, Vec **vp)
1566: {
1569: if (numfwdint) *numfwdint = ts->numcost;
1570: if (vp) *vp = ts->vecs_integral_sensip;
1571: return 0;
1572: }
1574: /*@
1575: TSForwardStep - Compute the forward sensitivity for one time step.
1577: Collective on TS
1579: Input Parameter:
1580: . ts - time stepping context
1582: Level: advanced
1584: Notes:
1585: This function cannot be called until TSStep() has been completed.
1587: .seealso: `TSForwardSetSensitivities()`, `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardSetUp()`
1588: @*/
1589: PetscErrorCode TSForwardStep(TS ts)
1590: {
1592: PetscLogEventBegin(TS_ForwardStep, ts, 0, 0, 0);
1593: PetscUseTypeMethod(ts, forwardstep);
1594: PetscLogEventEnd(TS_ForwardStep, ts, 0, 0, 0);
1596: return 0;
1597: }
1599: /*@
1600: TSForwardSetSensitivities - Sets the initial value of the trajectory sensitivities of solution w.r.t. the problem parameters and initial values.
1602: Logically Collective on TS
1604: Input Parameters:
1605: + ts - the TS context obtained from TSCreate()
1606: . nump - number of parameters
1607: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1609: Level: beginner
1611: Notes:
1612: Forward sensitivity is also called 'trajectory sensitivity' in some fields such as power systems.
1613: This function turns on a flag to trigger TSSolve() to compute forward sensitivities automatically.
1614: You must call this function before TSSolve().
1615: The entries in the sensitivity matrix must be correctly initialized with the values S = dy/dp|startingtime.
1617: .seealso: `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1618: @*/
1619: PetscErrorCode TSForwardSetSensitivities(TS ts, PetscInt nump, Mat Smat)
1620: {
1623: ts->forward_solve = PETSC_TRUE;
1624: if (nump == PETSC_DEFAULT) {
1625: MatGetSize(Smat, NULL, &ts->num_parameters);
1626: } else ts->num_parameters = nump;
1627: PetscObjectReference((PetscObject)Smat);
1628: MatDestroy(&ts->mat_sensip);
1629: ts->mat_sensip = Smat;
1630: return 0;
1631: }
1633: /*@
1634: TSForwardGetSensitivities - Returns the trajectory sensitivities
1636: Not Collective, but Vec returned is parallel if TS is parallel
1638: Output Parameters:
1639: + ts - the TS context obtained from TSCreate()
1640: . nump - number of parameters
1641: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1643: Level: intermediate
1645: .seealso: `TSForwardSetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1646: @*/
1647: PetscErrorCode TSForwardGetSensitivities(TS ts, PetscInt *nump, Mat *Smat)
1648: {
1650: if (nump) *nump = ts->num_parameters;
1651: if (Smat) *Smat = ts->mat_sensip;
1652: return 0;
1653: }
1655: /*@
1656: TSForwardCostIntegral - Evaluate the cost integral in the forward run.
1658: Collective on TS
1660: Input Parameter:
1661: . ts - time stepping context
1663: Level: advanced
1665: Notes:
1666: This function cannot be called until TSStep() has been completed.
1668: .seealso: `TSSolve()`, `TSAdjointCostIntegral()`
1669: @*/
1670: PetscErrorCode TSForwardCostIntegral(TS ts)
1671: {
1673: PetscUseTypeMethod(ts, forwardintegral);
1674: return 0;
1675: }
1677: /*@
1678: TSForwardSetInitialSensitivities - Set initial values for tangent linear sensitivities
1680: Collective on TS
1682: Input Parameters:
1683: + ts - the TS context obtained from TSCreate()
1684: - didp - parametric sensitivities of the initial condition
1686: Level: intermediate
1688: Notes: TSSolve() allows users to pass the initial solution directly to TS. But the tangent linear variables cannot be initialized in this way. This function is used to set initial values for tangent linear variables.
1690: .seealso: `TSForwardSetSensitivities()`
1691: @*/
1692: PetscErrorCode TSForwardSetInitialSensitivities(TS ts, Mat didp)
1693: {
1696: if (!ts->mat_sensip) TSForwardSetSensitivities(ts, PETSC_DEFAULT, didp);
1697: return 0;
1698: }
1700: /*@
1701: TSForwardGetStages - Get the number of stages and the tangent linear sensitivities at the intermediate stages
1703: Input Parameter:
1704: . ts - the TS context obtained from TSCreate()
1706: Output Parameters:
1707: + ns - number of stages
1708: - S - tangent linear sensitivities at the intermediate stages
1710: Level: advanced
1712: @*/
1713: PetscErrorCode TSForwardGetStages(TS ts, PetscInt *ns, Mat **S)
1714: {
1717: if (!ts->ops->getstages) *S = NULL;
1718: else PetscUseTypeMethod(ts, forwardgetstages, ns, S);
1719: return 0;
1720: }
1722: /*@
1723: TSCreateQuadratureTS - Create a sub-TS that evaluates integrals over time
1725: Input Parameters:
1726: + ts - the TS context obtained from TSCreate()
1727: - fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1729: Output Parameters:
1730: . quadts - the child TS context
1732: Level: intermediate
1734: .seealso: `TSGetQuadratureTS()`
1735: @*/
1736: PetscErrorCode TSCreateQuadratureTS(TS ts, PetscBool fwd, TS *quadts)
1737: {
1738: char prefix[128];
1742: TSDestroy(&ts->quadraturets);
1743: TSCreate(PetscObjectComm((PetscObject)ts), &ts->quadraturets);
1744: PetscObjectIncrementTabLevel((PetscObject)ts->quadraturets, (PetscObject)ts, 1);
1745: PetscSNPrintf(prefix, sizeof(prefix), "%squad_", ((PetscObject)ts)->prefix ? ((PetscObject)ts)->prefix : "");
1746: TSSetOptionsPrefix(ts->quadraturets, prefix);
1747: *quadts = ts->quadraturets;
1749: if (ts->numcost) {
1750: VecCreateSeq(PETSC_COMM_SELF, ts->numcost, &(*quadts)->vec_sol);
1751: } else {
1752: VecCreateSeq(PETSC_COMM_SELF, 1, &(*quadts)->vec_sol);
1753: }
1754: ts->costintegralfwd = fwd;
1755: return 0;
1756: }
1758: /*@
1759: TSGetQuadratureTS - Return the sub-TS that evaluates integrals over time
1761: Input Parameter:
1762: . ts - the TS context obtained from TSCreate()
1764: Output Parameters:
1765: + fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1766: - quadts - the child TS context
1768: Level: intermediate
1770: .seealso: `TSCreateQuadratureTS()`
1771: @*/
1772: PetscErrorCode TSGetQuadratureTS(TS ts, PetscBool *fwd, TS *quadts)
1773: {
1775: if (fwd) *fwd = ts->costintegralfwd;
1776: if (quadts) *quadts = ts->quadraturets;
1777: return 0;
1778: }
1780: /*@
1781: TSComputeSNESJacobian - Compute the SNESJacobian
1783: Input Parameters:
1784: + ts - the TS context obtained from TSCreate()
1785: - x - state vector
1787: Output Parameters:
1788: + J - Jacobian matrix
1789: - Jpre - preconditioning matrix for J (may be same as J)
1791: Level: developer
1793: Notes:
1794: Using SNES to compute the Jacobian enables finite differencing when TS Jacobian is not available.
1795: @*/
1796: PetscErrorCode TSComputeSNESJacobian(TS ts, Vec x, Mat J, Mat Jpre)
1797: {
1798: SNES snes = ts->snes;
1799: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *) = NULL;
1801: /*
1802: Unlike implicit methods, explicit methods do not have SNESMatFDColoring in the snes object
1803: because SNESSolve() has not been called yet; so querying SNESMatFDColoring does not work for
1804: explicit methods. Instead, we check the Jacobian compute function directly to determin if FD
1805: coloring is used.
1806: */
1807: SNESGetJacobian(snes, NULL, NULL, &jac, NULL);
1808: if (jac == SNESComputeJacobianDefaultColor) {
1809: Vec f;
1810: SNESSetSolution(snes, x);
1811: SNESGetFunction(snes, &f, NULL, NULL);
1812: /* Force MatFDColoringApply to evaluate the SNES residual function for the base vector */
1813: SNESComputeFunction(snes, x, f);
1814: }
1815: SNESComputeJacobian(snes, x, J, Jpre);
1816: return 0;
1817: }