Actual source code: tssen.c
1: #include <petsc/private/tsimpl.h>
2: #include <petscdraw.h>
4: PetscLogEvent TS_AdjointStep, TS_ForwardStep, TS_JacobianPEval;
6: /* #define TSADJOINT_STAGE */
8: /* ------------------------ Sensitivity Context ---------------------------*/
10: /*@C
11: TSSetRHSJacobianP - Sets the function that computes the Jacobian of G w.r.t. the parameters P where U_t = G(U,P,t), as well as the location to store the matrix.
13: Logically Collective
15: Input Parameters:
16: + ts - `TS` context obtained from `TSCreate()`
17: . Amat - JacobianP matrix
18: . func - function
19: - ctx - [optional] user-defined function context
21: Calling sequence of func:
22: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
23: + t - current timestep
24: . U - input vector (current ODE solution)
25: . A - output matrix
26: - ctx - [optional] user-defined function context
28: Level: intermediate
30: Note:
31: Amat has the same number of rows and the same row parallel layout as u, Amat has the same number of columns and parallel layout as p
33: .seealso: [](chapter_ts), `TS`, `TSGetRHSJacobianP()`
34: @*/
35: PetscErrorCode TSSetRHSJacobianP(TS ts, Mat Amat, PetscErrorCode (*func)(TS, PetscReal, Vec, Mat, void *), void *ctx)
36: {
40: ts->rhsjacobianp = func;
41: ts->rhsjacobianpctx = ctx;
42: if (Amat) {
43: PetscObjectReference((PetscObject)Amat);
44: MatDestroy(&ts->Jacprhs);
45: ts->Jacprhs = Amat;
46: }
47: return 0;
48: }
50: /*@C
51: TSGetRHSJacobianP - Gets the function that computes the Jacobian of G w.r.t. the parameters P where U_t = G(U,P,t), as well as the location to store the matrix.
53: Logically Collective
55: Input Parameter:
56: . ts - `TS` context obtained from `TSCreate()`
58: Output Parameters:
59: + Amat - JacobianP matrix
60: . func - function
61: - ctx - [optional] user-defined function context
63: Calling sequence of func:
64: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
65: + t - current timestep
66: . U - input vector (current ODE solution)
67: . A - output matrix
68: - ctx - [optional] user-defined function context
70: Level: intermediate
72: Note:
73: Amat has the same number of rows and the same row parallel layout as u, Amat has the same number of columns and parallel layout as p
75: .seealso: [](chapter_ts), `TSSetRHSJacobianP()`, `TS`, `TSGetRHSJacobianP()`
76: @*/
77: PetscErrorCode TSGetRHSJacobianP(TS ts, Mat *Amat, PetscErrorCode (**func)(TS, PetscReal, Vec, Mat, void *), void **ctx)
78: {
79: if (func) *func = ts->rhsjacobianp;
80: if (ctx) *ctx = ts->rhsjacobianpctx;
81: if (Amat) *Amat = ts->Jacprhs;
82: return 0;
83: }
85: /*@C
86: TSComputeRHSJacobianP - Runs the user-defined JacobianP function.
88: Collective
90: Input Parameters:
91: . ts - The `TS` context obtained from `TSCreate()`
93: Level: developer
95: .seealso: [](chapter_ts), `TSSetRHSJacobianP()`, `TS`
96: @*/
97: PetscErrorCode TSComputeRHSJacobianP(TS ts, PetscReal t, Vec U, Mat Amat)
98: {
99: if (!Amat) return 0;
103: PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->rhsjacobianp)(ts, t, U, Amat, ts->rhsjacobianpctx));
104: return 0;
105: }
107: /*@C
108: TSSetIJacobianP - Sets the function that computes the Jacobian of F w.r.t. the parameters P where F(Udot,U,t) = G(U,P,t), as well as the location to store the matrix.
110: Logically Collective
112: Input Parameters:
113: + ts - `TS` context obtained from `TSCreate()`
114: . Amat - JacobianP matrix
115: . func - function
116: - ctx - [optional] user-defined function context
118: Calling sequence of func:
119: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
120: + t - current timestep
121: . U - input vector (current ODE solution)
122: . Udot - time derivative of state vector
123: . shift - shift to apply, see note below
124: . A - output matrix
125: - ctx - [optional] user-defined function context
127: Level: intermediate
129: Note:
130: Amat has the same number of rows and the same row parallel layout as u, Amat has the same number of columns and parallel layout as p
132: .seealso: [](chapter_ts), `TSSetRHSJacobianP()`, `TS`
133: @*/
134: PetscErrorCode TSSetIJacobianP(TS ts, Mat Amat, PetscErrorCode (*func)(TS, PetscReal, Vec, Vec, PetscReal, Mat, void *), void *ctx)
135: {
139: ts->ijacobianp = func;
140: ts->ijacobianpctx = ctx;
141: if (Amat) {
142: PetscObjectReference((PetscObject)Amat);
143: MatDestroy(&ts->Jacp);
144: ts->Jacp = Amat;
145: }
146: return 0;
147: }
149: /*@C
150: TSComputeIJacobianP - Runs the user-defined IJacobianP function.
152: Collective
154: Input Parameters:
155: + ts - the `TS` context
156: . t - current timestep
157: . U - state vector
158: . Udot - time derivative of state vector
159: . shift - shift to apply, see note below
160: - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate
162: Output Parameters:
163: . A - Jacobian matrix
165: Level: developer
167: .seealso: [](chapter_ts), `TS`, `TSSetIJacobianP()`
168: @*/
169: PetscErrorCode TSComputeIJacobianP(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat Amat, PetscBool imex)
170: {
171: if (!Amat) return 0;
176: PetscLogEventBegin(TS_JacobianPEval, ts, U, Amat, 0);
177: if (ts->ijacobianp) PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->ijacobianp)(ts, t, U, Udot, shift, Amat, ts->ijacobianpctx));
178: if (imex) {
179: if (!ts->ijacobianp) { /* system was written as Udot = G(t,U) */
180: PetscBool assembled;
181: MatZeroEntries(Amat);
182: MatAssembled(Amat, &assembled);
183: if (!assembled) {
184: MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
185: MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
186: }
187: }
188: } else {
189: if (ts->rhsjacobianp) TSComputeRHSJacobianP(ts, t, U, ts->Jacprhs);
190: if (ts->Jacprhs == Amat) { /* No IJacobian, so we only have the RHS matrix */
191: MatScale(Amat, -1);
192: } else if (ts->Jacprhs) { /* Both IJacobian and RHSJacobian */
193: MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
194: if (!ts->ijacobianp) { /* No IJacobianp provided, but we have a separate RHS matrix */
195: MatZeroEntries(Amat);
196: }
197: MatAXPY(Amat, -1, ts->Jacprhs, axpy);
198: }
199: }
200: PetscLogEventEnd(TS_JacobianPEval, ts, U, Amat, 0);
201: return 0;
202: }
204: /*@C
205: TSSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions
207: Logically Collective
209: Input Parameters:
210: + ts - the `TS` context obtained from `TSCreate()`
211: . numcost - number of gradients to be computed, this is the number of cost functions
212: . costintegral - vector that stores the integral values
213: . rf - routine for evaluating the integrand function
214: . drduf - function that computes the gradients of the r's with respect to u
215: . drdpf - function that computes the gradients of the r's with respect to p, can be NULL if parametric sensitivity is not desired (mu=NULL)
216: . fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
217: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
219: Calling sequence of rf:
220: $ PetscErrorCode rf(TS ts,PetscReal t,Vec U,Vec F,void *ctx);
222: Calling sequence of drduf:
223: $ PetscErroCode drduf(TS ts,PetscReal t,Vec U,Vec *dRdU,void *ctx);
225: Calling sequence of drdpf:
226: $ PetscErroCode drdpf(TS ts,PetscReal t,Vec U,Vec *dRdP,void *ctx);
228: Level: deprecated
230: Note:
231: For optimization there is usually a single cost function (numcost = 1). For sensitivities there may be multiple cost functions
233: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobianP()`, `TSGetCostGradients()`, `TSSetCostGradients()`
234: @*/
235: PetscErrorCode TSSetCostIntegrand(TS ts, PetscInt numcost, Vec costintegral, PetscErrorCode (*rf)(TS, PetscReal, Vec, Vec, void *), PetscErrorCode (*drduf)(TS, PetscReal, Vec, Vec *, void *), PetscErrorCode (*drdpf)(TS, PetscReal, Vec, Vec *, void *), PetscBool fwd, void *ctx)
236: {
240: if (!ts->numcost) ts->numcost = numcost;
242: if (costintegral) {
243: PetscObjectReference((PetscObject)costintegral);
244: VecDestroy(&ts->vec_costintegral);
245: ts->vec_costintegral = costintegral;
246: } else {
247: if (!ts->vec_costintegral) { /* Create a seq vec if user does not provide one */
248: VecCreateSeq(PETSC_COMM_SELF, numcost, &ts->vec_costintegral);
249: } else {
250: VecSet(ts->vec_costintegral, 0.0);
251: }
252: }
253: if (!ts->vec_costintegrand) {
254: VecDuplicate(ts->vec_costintegral, &ts->vec_costintegrand);
255: } else {
256: VecSet(ts->vec_costintegrand, 0.0);
257: }
258: ts->costintegralfwd = fwd; /* Evaluate the cost integral in forward run if fwd is true */
259: ts->costintegrand = rf;
260: ts->costintegrandctx = ctx;
261: ts->drdufunction = drduf;
262: ts->drdpfunction = drdpf;
263: return 0;
264: }
266: /*@C
267: TSGetCostIntegral - Returns the values of the integral term in the cost functions.
268: It is valid to call the routine after a backward run.
270: Not Collective
272: Input Parameter:
273: . ts - the `TS` context obtained from `TSCreate()`
275: Output Parameter:
276: . v - the vector containing the integrals for each cost function
278: Level: intermediate
280: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, ``TSSetCostIntegrand()`
281: @*/
282: PetscErrorCode TSGetCostIntegral(TS ts, Vec *v)
283: {
284: TS quadts;
288: TSGetQuadratureTS(ts, NULL, &quadts);
289: *v = quadts->vec_sol;
290: return 0;
291: }
293: /*@C
294: TSComputeCostIntegrand - Evaluates the integral function in the cost functions.
296: Input Parameters:
297: + ts - the `TS` context
298: . t - current time
299: - U - state vector, i.e. current solution
301: Output Parameter:
302: . Q - vector of size numcost to hold the outputs
304: Level: deprecated
306: Note:
307: Most users should not need to explicitly call this routine, as it
308: is used internally within the sensitivity analysis context.
310: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSSetCostIntegrand()`
311: @*/
312: PetscErrorCode TSComputeCostIntegrand(TS ts, PetscReal t, Vec U, Vec Q)
313: {
318: PetscLogEventBegin(TS_FunctionEval, ts, U, Q, 0);
319: if (ts->costintegrand) PetscCallBack("TS callback integrand in the cost function", (*ts->costintegrand)(ts, t, U, Q, ts->costintegrandctx));
320: else VecZeroEntries(Q);
321: PetscLogEventEnd(TS_FunctionEval, ts, U, Q, 0);
322: return 0;
323: }
325: /*@C
326: TSComputeDRDUFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobian()`
328: Level: deprecated
330: @*/
331: PetscErrorCode TSComputeDRDUFunction(TS ts, PetscReal t, Vec U, Vec *DRDU)
332: {
333: if (!DRDU) return 0;
337: PetscCallBack("TS callback DRDU for sensitivity analysis", (*ts->drdufunction)(ts, t, U, DRDU, ts->costintegrandctx));
338: return 0;
339: }
341: /*@C
342: TSComputeDRDPFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobianP()`
344: Level: deprecated
346: @*/
347: PetscErrorCode TSComputeDRDPFunction(TS ts, PetscReal t, Vec U, Vec *DRDP)
348: {
349: if (!DRDP) return 0;
353: PetscCallBack("TS callback DRDP for sensitivity analysis", (*ts->drdpfunction)(ts, t, U, DRDP, ts->costintegrandctx));
354: return 0;
355: }
357: /*@C
358: TSSetIHessianProduct - Sets the function that computes the vector-Hessian-vector product. The Hessian is the second-order derivative of F (IFunction) w.r.t. the state variable.
360: Logically Collective
362: Input Parameters:
363: + ts - `TS` context obtained from `TSCreate()`
364: . ihp1 - an array of vectors storing the result of vector-Hessian-vector product for F_UU
365: . hessianproductfunc1 - vector-Hessian-vector product function for F_UU
366: . ihp2 - an array of vectors storing the result of vector-Hessian-vector product for F_UP
367: . hessianproductfunc2 - vector-Hessian-vector product function for F_UP
368: . ihp3 - an array of vectors storing the result of vector-Hessian-vector product for F_PU
369: . hessianproductfunc3 - vector-Hessian-vector product function for F_PU
370: . ihp4 - an array of vectors storing the result of vector-Hessian-vector product for F_PP
371: - hessianproductfunc4 - vector-Hessian-vector product function for F_PP
373: Calling sequence of ihessianproductfunc:
374: $ ihessianproductfunc (TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx);
375: + t - current timestep
376: . U - input vector (current ODE solution)
377: . Vl - an array of input vectors to be left-multiplied with the Hessian
378: . Vr - input vector to be right-multiplied with the Hessian
379: . VHV - an array of output vectors for vector-Hessian-vector product
380: - ctx - [optional] user-defined function context
382: Level: intermediate
384: Notes:
385: The first Hessian function and the working array are required.
386: As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
387: $ Vl_n^T*F_UP*Vr
388: where the vector Vl_n (n-th element in the array Vl) and Vr are of size N and M respectively, and the Hessian F_UP is of size N x N x M.
389: Each entry of F_UP corresponds to the derivative
390: $ F_UP[i][j][k] = \frac{\partial^2 F[i]}{\partial U[j] \partial P[k]}.
391: The result of the vector-Hessian-vector product for Vl_n needs to be stored in vector VHV_n with the j-th entry being
392: $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * F_UP[i][j][k] * Vr[k]}
393: If the cost function is a scalar, there will be only one vector in Vl and VHV.
395: .seealso: [](chapter_ts), `TS`
396: @*/
397: PetscErrorCode TSSetIHessianProduct(TS ts, Vec *ihp1, PetscErrorCode (*ihessianproductfunc1)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *ihp2, PetscErrorCode (*ihessianproductfunc2)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *ihp3, PetscErrorCode (*ihessianproductfunc3)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *ihp4, PetscErrorCode (*ihessianproductfunc4)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), void *ctx)
398: {
402: ts->ihessianproductctx = ctx;
403: if (ihp1) ts->vecs_fuu = ihp1;
404: if (ihp2) ts->vecs_fup = ihp2;
405: if (ihp3) ts->vecs_fpu = ihp3;
406: if (ihp4) ts->vecs_fpp = ihp4;
407: ts->ihessianproduct_fuu = ihessianproductfunc1;
408: ts->ihessianproduct_fup = ihessianproductfunc2;
409: ts->ihessianproduct_fpu = ihessianproductfunc3;
410: ts->ihessianproduct_fpp = ihessianproductfunc4;
411: return 0;
412: }
414: /*@C
415: TSComputeIHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Fuu.
417: Collective
419: Input Parameters:
420: . ts - The `TS` context obtained from `TSCreate()`
422: Level: developer
424: Note:
425: `TSComputeIHessianProductFunctionUU()` is typically used for sensitivity implementation,
426: so most users would not generally call this routine themselves.
428: .seealso: [](chapter_ts), `TSSetIHessianProduct()`
429: @*/
430: PetscErrorCode TSComputeIHessianProductFunctionUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
431: {
432: if (!VHV) return 0;
436: if (ts->ihessianproduct_fuu) PetscCallBack("TS callback IHessianProduct 1 for sensitivity analysis", (*ts->ihessianproduct_fuu)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
438: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
439: if (ts->rhshessianproduct_guu) {
440: PetscInt nadj;
441: TSComputeRHSHessianProductFunctionUU(ts, t, U, Vl, Vr, VHV);
442: for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
443: }
444: return 0;
445: }
447: /*@C
448: TSComputeIHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Fup.
450: Collective
452: Input Parameters:
453: . ts - The `TS` context obtained from `TSCreate()`
455: Level: developer
457: Note:
458: `TSComputeIHessianProductFunctionUP()` is typically used for sensitivity implementation,
459: so most users would not generally call this routine themselves.
461: .seealso: [](chapter_ts), `TSSetIHessianProduct()`
462: @*/
463: PetscErrorCode TSComputeIHessianProductFunctionUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
464: {
465: if (!VHV) return 0;
469: if (ts->ihessianproduct_fup) PetscCallBack("TS callback IHessianProduct 2 for sensitivity analysis", (*ts->ihessianproduct_fup)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
471: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
472: if (ts->rhshessianproduct_gup) {
473: PetscInt nadj;
474: TSComputeRHSHessianProductFunctionUP(ts, t, U, Vl, Vr, VHV);
475: for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
476: }
477: return 0;
478: }
480: /*@C
481: TSComputeIHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Fpu.
483: Collective
485: Input Parameters:
486: . ts - The `TS` context obtained from `TSCreate()`
488: Level: developer
490: Note:
491: `TSComputeIHessianProductFunctionPU()` is typically used for sensitivity implementation,
492: so most users would not generally call this routine themselves.
494: .seealso: [](chapter_ts), `TSSetIHessianProduct()`
495: @*/
496: PetscErrorCode TSComputeIHessianProductFunctionPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
497: {
498: if (!VHV) return 0;
502: if (ts->ihessianproduct_fpu) PetscCallBack("TS callback IHessianProduct 3 for sensitivity analysis", (*ts->ihessianproduct_fpu)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
504: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
505: if (ts->rhshessianproduct_gpu) {
506: PetscInt nadj;
507: TSComputeRHSHessianProductFunctionPU(ts, t, U, Vl, Vr, VHV);
508: for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
509: }
510: return 0;
511: }
513: /*@C
514: TSComputeIHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Fpp.
516: Collective
518: Input Parameters:
519: . ts - The `TS` context obtained from `TSCreate()`
521: Level: developer
523: Note:
524: `TSComputeIHessianProductFunctionPP()` is typically used for sensitivity implementation,
525: so most users would not generally call this routine themselves.
527: .seealso: [](chapter_ts), `TSSetIHessianProduct()`
528: @*/
529: PetscErrorCode TSComputeIHessianProductFunctionPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
530: {
531: if (!VHV) return 0;
535: if (ts->ihessianproduct_fpp) PetscCallBack("TS callback IHessianProduct 3 for sensitivity analysis", (*ts->ihessianproduct_fpp)(ts, t, U, Vl, Vr, VHV, ts->ihessianproductctx));
537: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
538: if (ts->rhshessianproduct_gpp) {
539: PetscInt nadj;
540: TSComputeRHSHessianProductFunctionPP(ts, t, U, Vl, Vr, VHV);
541: for (nadj = 0; nadj < ts->numcost; nadj++) VecScale(VHV[nadj], -1);
542: }
543: return 0;
544: }
546: /*@C
547: TSSetRHSHessianProduct - Sets the function that computes the vector-Hessian-vector product. The Hessian is the second-order derivative of G (RHSFunction) w.r.t. the state variable.
549: Logically Collective
551: Input Parameters:
552: + ts - `TS` context obtained from `TSCreate()`
553: . rhshp1 - an array of vectors storing the result of vector-Hessian-vector product for G_UU
554: . hessianproductfunc1 - vector-Hessian-vector product function for G_UU
555: . rhshp2 - an array of vectors storing the result of vector-Hessian-vector product for G_UP
556: . hessianproductfunc2 - vector-Hessian-vector product function for G_UP
557: . rhshp3 - an array of vectors storing the result of vector-Hessian-vector product for G_PU
558: . hessianproductfunc3 - vector-Hessian-vector product function for G_PU
559: . rhshp4 - an array of vectors storing the result of vector-Hessian-vector product for G_PP
560: - hessianproductfunc4 - vector-Hessian-vector product function for G_PP
562: Calling sequence of ihessianproductfunc:
563: $ rhshessianproductfunc (TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx);
564: + t - current timestep
565: . U - input vector (current ODE solution)
566: . Vl - an array of input vectors to be left-multiplied with the Hessian
567: . Vr - input vector to be right-multiplied with the Hessian
568: . VHV - an array of output vectors for vector-Hessian-vector product
569: - ctx - [optional] user-defined function context
571: Level: intermediate
573: Notes:
574: The first Hessian function and the working array are required.
575: As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
576: $ Vl_n^T*G_UP*Vr
577: where the vector Vl_n (n-th element in the array Vl) and Vr are of size N and M respectively, and the Hessian G_UP is of size N x N x M.
578: Each entry of G_UP corresponds to the derivative
579: $ G_UP[i][j][k] = \frac{\partial^2 G[i]}{\partial U[j] \partial P[k]}.
580: The result of the vector-Hessian-vector product for Vl_n needs to be stored in vector VHV_n with j-th entry being
581: $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * G_UP[i][j][k] * Vr[k]}
582: If the cost function is a scalar, there will be only one vector in Vl and VHV.
584: .seealso: `TS`
585: @*/
586: PetscErrorCode TSSetRHSHessianProduct(TS ts, Vec *rhshp1, PetscErrorCode (*rhshessianproductfunc1)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *rhshp2, PetscErrorCode (*rhshessianproductfunc2)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *rhshp3, PetscErrorCode (*rhshessianproductfunc3)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), Vec *rhshp4, PetscErrorCode (*rhshessianproductfunc4)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *), void *ctx)
587: {
591: ts->rhshessianproductctx = ctx;
592: if (rhshp1) ts->vecs_guu = rhshp1;
593: if (rhshp2) ts->vecs_gup = rhshp2;
594: if (rhshp3) ts->vecs_gpu = rhshp3;
595: if (rhshp4) ts->vecs_gpp = rhshp4;
596: ts->rhshessianproduct_guu = rhshessianproductfunc1;
597: ts->rhshessianproduct_gup = rhshessianproductfunc2;
598: ts->rhshessianproduct_gpu = rhshessianproductfunc3;
599: ts->rhshessianproduct_gpp = rhshessianproductfunc4;
600: return 0;
601: }
603: /*@C
604: TSComputeRHSHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Guu.
606: Collective
608: Input Parameters:
609: . ts - The `TS` context obtained from `TSCreate()`
611: Level: developer
613: Note:
614: `TSComputeRHSHessianProductFunctionUU()` is typically used for sensitivity implementation,
615: so most users would not generally call this routine themselves.
617: .seealso: [](chapter_ts), `TS`, `TSSetRHSHessianProduct()`
618: @*/
619: PetscErrorCode TSComputeRHSHessianProductFunctionUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
620: {
621: if (!VHV) return 0;
625: PetscCallBack("TS callback RHSHessianProduct 1 for sensitivity analysis", (*ts->rhshessianproduct_guu)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
626: return 0;
627: }
629: /*@C
630: TSComputeRHSHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Gup.
632: Collective
634: Input Parameters:
635: . ts - The `TS` context obtained from `TSCreate()`
637: Level: developer
639: Note:
640: `TSComputeRHSHessianProductFunctionUP()` is typically used for sensitivity implementation,
641: so most users would not generally call this routine themselves.
643: .seealso: [](chapter_ts), `TS`, `TSSetRHSHessianProduct()`
644: @*/
645: PetscErrorCode TSComputeRHSHessianProductFunctionUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
646: {
647: if (!VHV) return 0;
651: PetscCallBack("TS callback RHSHessianProduct 2 for sensitivity analysis", (*ts->rhshessianproduct_gup)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
652: return 0;
653: }
655: /*@C
656: TSComputeRHSHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Gpu.
658: Collective
660: Input Parameters:
661: . ts - The `TS` context obtained from `TSCreate()`
663: Level: developer
665: Note:
666: `TSComputeRHSHessianProductFunctionPU()` is typically used for sensitivity implementation,
667: so most users would not generally call this routine themselves.
669: .seealso: [](chapter_ts), `TSSetRHSHessianProduct()`
670: @*/
671: PetscErrorCode TSComputeRHSHessianProductFunctionPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
672: {
673: if (!VHV) return 0;
677: PetscCallBack("TS callback RHSHessianProduct 3 for sensitivity analysis", (*ts->rhshessianproduct_gpu)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
678: return 0;
679: }
681: /*@C
682: TSComputeRHSHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Gpp.
684: Collective
686: Input Parameters:
687: . ts - The `TS` context obtained from `TSCreate()`
689: Level: developer
691: Note:
692: `TSComputeRHSHessianProductFunctionPP()` is typically used for sensitivity implementation,
693: so most users would not generally call this routine themselves.
695: .seealso: [](chapter_ts), `TSSetRHSHessianProduct()`
696: @*/
697: PetscErrorCode TSComputeRHSHessianProductFunctionPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV)
698: {
699: if (!VHV) return 0;
703: PetscCallBack("TS callback RHSHessianProduct 3 for sensitivity analysis", (*ts->rhshessianproduct_gpp)(ts, t, U, Vl, Vr, VHV, ts->rhshessianproductctx));
704: return 0;
705: }
707: /* --------------------------- Adjoint sensitivity ---------------------------*/
709: /*@
710: TSSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial values and w.r.t. the problem parameters
711: for use by the `TS` adjoint routines.
713: Logically Collective
715: Input Parameters:
716: + ts - the `TS` context obtained from `TSCreate()`
717: . numcost - number of gradients to be computed, this is the number of cost functions
718: . lambda - gradients with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
719: - mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
721: Level: beginner
723: Notes:
724: the entries in these vectors must be correctly initialized with the values lambda_i = df/dy|finaltime mu_i = df/dp|finaltime
726: After `TSAdjointSolve()` is called the lambda and the mu contain the computed sensitivities
728: .seealso: `TS`, `TSAdjointSolve()`, `TSGetCostGradients()`
729: @*/
730: PetscErrorCode TSSetCostGradients(TS ts, PetscInt numcost, Vec *lambda, Vec *mu)
731: {
734: ts->vecs_sensi = lambda;
735: ts->vecs_sensip = mu;
737: ts->numcost = numcost;
738: return 0;
739: }
741: /*@
742: TSGetCostGradients - Returns the gradients from the `TSAdjointSolve()`
744: Not Collective, but the vectors returned are parallel if `TS` is parallel
746: Input Parameter:
747: . ts - the `TS` context obtained from `TSCreate()`
749: Output Parameters:
750: + numcost - size of returned arrays
751: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
752: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
754: Level: intermediate
756: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSSetCostGradients()`
757: @*/
758: PetscErrorCode TSGetCostGradients(TS ts, PetscInt *numcost, Vec **lambda, Vec **mu)
759: {
761: if (numcost) *numcost = ts->numcost;
762: if (lambda) *lambda = ts->vecs_sensi;
763: if (mu) *mu = ts->vecs_sensip;
764: return 0;
765: }
767: /*@
768: TSSetCostHessianProducts - Sets the initial value of the Hessian-vector products of the cost function w.r.t. initial values and w.r.t. the problem parameters
769: for use by the `TS` adjoint routines.
771: Logically Collective
773: Input Parameters:
774: + ts - the `TS` context obtained from `TSCreate()`
775: . numcost - number of cost functions
776: . lambda2 - Hessian-vector product with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
777: . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
778: - dir - the direction vector that are multiplied with the Hessian of the cost functions
780: Level: beginner
782: Notes:
783: Hessian of the cost function is completely different from Hessian of the ODE/DAE system
785: For second-order adjoint, one needs to call this function and then `TSAdjointSetForward()` before `TSSolve()`.
787: After `TSAdjointSolve()` is called, the lambda2 and the mu2 will contain the computed second-order adjoint sensitivities, and can be used to produce Hessian-vector product (not the full Hessian matrix). Users must provide a direction vector; it is usually generated by an optimization solver.
789: Passing NULL for lambda2 disables the second-order calculation.
791: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSAdjointSetForward()`
792: @*/
793: PetscErrorCode TSSetCostHessianProducts(TS ts, PetscInt numcost, Vec *lambda2, Vec *mu2, Vec dir)
794: {
797: ts->numcost = numcost;
798: ts->vecs_sensi2 = lambda2;
799: ts->vecs_sensi2p = mu2;
800: ts->vec_dir = dir;
801: return 0;
802: }
804: /*@
805: TSGetCostHessianProducts - Returns the gradients from the `TSAdjointSolve()`
807: Not Collective, but vectors returned are parallel if `TS` is parallel
809: Input Parameter:
810: . ts - the `TS` context obtained from `TSCreate()`
812: Output Parameters:
813: + numcost - number of cost functions
814: . lambda2 - Hessian-vector product with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
815: . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
816: - dir - the direction vector that are multiplied with the Hessian of the cost functions
818: Level: intermediate
820: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSSetCostHessianProducts()`
821: @*/
822: PetscErrorCode TSGetCostHessianProducts(TS ts, PetscInt *numcost, Vec **lambda2, Vec **mu2, Vec *dir)
823: {
825: if (numcost) *numcost = ts->numcost;
826: if (lambda2) *lambda2 = ts->vecs_sensi2;
827: if (mu2) *mu2 = ts->vecs_sensi2p;
828: if (dir) *dir = ts->vec_dir;
829: return 0;
830: }
832: /*@
833: TSAdjointSetForward - Trigger the tangent linear solver and initialize the forward sensitivities
835: Logically Collective
837: Input Parameters:
838: + ts - the `TS` context obtained from `TSCreate()`
839: - didp - the derivative of initial values w.r.t. parameters
841: Level: intermediate
843: Notes:
844: When computing sensitivies w.r.t. initial condition, set didp to NULL so that the solver will take it as an identity matrix mathematically.
845: `TS` adjoint does not reset the tangent linear solver automatically, `TSAdjointResetForward()` should be called to reset the tangent linear solver.
847: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSSetCostHessianProducts()`, `TSAdjointResetForward()`
848: @*/
849: PetscErrorCode TSAdjointSetForward(TS ts, Mat didp)
850: {
851: Mat A;
852: Vec sp;
853: PetscScalar *xarr;
854: PetscInt lsize;
856: ts->forward_solve = PETSC_TRUE; /* turn on tangent linear mode */
859: /* create a single-column dense matrix */
860: VecGetLocalSize(ts->vec_sol, &lsize);
861: MatCreateDense(PetscObjectComm((PetscObject)ts), lsize, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, &A);
863: VecDuplicate(ts->vec_sol, &sp);
864: MatDenseGetColumn(A, 0, &xarr);
865: VecPlaceArray(sp, xarr);
866: if (ts->vecs_sensi2p) { /* tangent linear variable initialized as 2*dIdP*dir */
867: if (didp) {
868: MatMult(didp, ts->vec_dir, sp);
869: VecScale(sp, 2.);
870: } else {
871: VecZeroEntries(sp);
872: }
873: } else { /* tangent linear variable initialized as dir */
874: VecCopy(ts->vec_dir, sp);
875: }
876: VecResetArray(sp);
877: MatDenseRestoreColumn(A, &xarr);
878: VecDestroy(&sp);
880: TSForwardSetInitialSensitivities(ts, A); /* if didp is NULL, identity matrix is assumed */
882: MatDestroy(&A);
883: return 0;
884: }
886: /*@
887: TSAdjointResetForward - Reset the tangent linear solver and destroy the tangent linear context
889: Logically Collective
891: Input Parameters:
892: . ts - the `TS` context obtained from `TSCreate()`
894: Level: intermediate
896: .seealso: [](chapter_ts), `TSAdjointSetForward()`
897: @*/
898: PetscErrorCode TSAdjointResetForward(TS ts)
899: {
900: ts->forward_solve = PETSC_FALSE; /* turn off tangent linear mode */
901: TSForwardReset(ts);
902: return 0;
903: }
905: /*@
906: TSAdjointSetUp - Sets up the internal data structures for the later use
907: of an adjoint solver
909: Collective
911: Input Parameter:
912: . ts - the `TS` context obtained from `TSCreate()`
914: Level: advanced
916: .seealso: [](chapter_ts), `TSCreate()`, `TSAdjointStep()`, `TSSetCostGradients()`
917: @*/
918: PetscErrorCode TSAdjointSetUp(TS ts)
919: {
920: TSTrajectory tj;
921: PetscBool match;
924: if (ts->adjointsetupcalled) return 0;
927: TSGetTrajectory(ts, &tj);
928: PetscObjectTypeCompare((PetscObject)tj, TSTRAJECTORYBASIC, &match);
929: if (match) {
930: PetscBool solution_only;
931: TSTrajectoryGetSolutionOnly(tj, &solution_only);
933: }
934: TSTrajectorySetUseHistory(tj, PETSC_FALSE); /* not use TSHistory */
936: if (ts->quadraturets) { /* if there is integral in the cost function */
937: VecDuplicate(ts->vecs_sensi[0], &ts->vec_drdu_col);
938: if (ts->vecs_sensip) VecDuplicate(ts->vecs_sensip[0], &ts->vec_drdp_col);
939: }
941: PetscTryTypeMethod(ts, adjointsetup);
942: ts->adjointsetupcalled = PETSC_TRUE;
943: return 0;
944: }
946: /*@
947: TSAdjointReset - Resets a `TS` adjoint context and removes any allocated `Vec`s and `Mat`s.
949: Collective
951: Input Parameter:
952: . ts - the `TS` context obtained from `TSCreate()`
954: Level: beginner
956: .seealso: [](chapter_ts), `TSCreate()`, `TSAdjointSetUp()`, `TSADestroy()`
957: @*/
958: PetscErrorCode TSAdjointReset(TS ts)
959: {
961: PetscTryTypeMethod(ts, adjointreset);
962: if (ts->quadraturets) { /* if there is integral in the cost function */
963: VecDestroy(&ts->vec_drdu_col);
964: if (ts->vecs_sensip) VecDestroy(&ts->vec_drdp_col);
965: }
966: ts->vecs_sensi = NULL;
967: ts->vecs_sensip = NULL;
968: ts->vecs_sensi2 = NULL;
969: ts->vecs_sensi2p = NULL;
970: ts->vec_dir = NULL;
971: ts->adjointsetupcalled = PETSC_FALSE;
972: return 0;
973: }
975: /*@
976: TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time
978: Logically Collective
980: Input Parameters:
981: + ts - the `TS` context obtained from `TSCreate()`
982: - steps - number of steps to use
984: Level: intermediate
986: Notes:
987: Normally one does not call this and `TSAdjointSolve()` integrates back to the original timestep. One can call this
988: so as to integrate back to less than the original timestep
990: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TS`, `TSSetExactFinalTime()`
991: @*/
992: PetscErrorCode TSAdjointSetSteps(TS ts, PetscInt steps)
993: {
998: ts->adjoint_max_steps = steps;
999: return 0;
1000: }
1002: /*@C
1003: TSAdjointSetRHSJacobian - Deprecated, use `TSSetRHSJacobianP()`
1005: Level: deprecated
1007: @*/
1008: PetscErrorCode TSAdjointSetRHSJacobian(TS ts, Mat Amat, PetscErrorCode (*func)(TS, PetscReal, Vec, Mat, void *), void *ctx)
1009: {
1013: ts->rhsjacobianp = func;
1014: ts->rhsjacobianpctx = ctx;
1015: if (Amat) {
1016: PetscObjectReference((PetscObject)Amat);
1017: MatDestroy(&ts->Jacp);
1018: ts->Jacp = Amat;
1019: }
1020: return 0;
1021: }
1023: /*@C
1024: TSAdjointComputeRHSJacobian - Deprecated, use `TSComputeRHSJacobianP()`
1026: Level: deprecated
1028: @*/
1029: PetscErrorCode TSAdjointComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat Amat)
1030: {
1035: PetscCallBack("TS callback JacobianP for sensitivity analysis", (*ts->rhsjacobianp)(ts, t, U, Amat, ts->rhsjacobianpctx));
1036: return 0;
1037: }
1039: /*@
1040: TSAdjointComputeDRDYFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobian()`
1042: Level: deprecated
1044: @*/
1045: PetscErrorCode TSAdjointComputeDRDYFunction(TS ts, PetscReal t, Vec U, Vec *DRDU)
1046: {
1050: PetscCallBack("TS callback DRDY for sensitivity analysis", (*ts->drdufunction)(ts, t, U, DRDU, ts->costintegrandctx));
1051: return 0;
1052: }
1054: /*@
1055: TSAdjointComputeDRDPFunction - Deprecated, use `TSGetQuadratureTS()` then `TSComputeRHSJacobianP()`
1057: Level: deprecated
1059: @*/
1060: PetscErrorCode TSAdjointComputeDRDPFunction(TS ts, PetscReal t, Vec U, Vec *DRDP)
1061: {
1065: PetscCallBack("TS callback DRDP for sensitivity analysis", (*ts->drdpfunction)(ts, t, U, DRDP, ts->costintegrandctx));
1066: return 0;
1067: }
1069: /*@C
1070: TSAdjointMonitorSensi - monitors the first lambda sensitivity
1072: Level: intermediate
1074: .seealso: [](chapter_ts), `TSAdjointMonitorSet()`
1075: @*/
1076: PetscErrorCode TSAdjointMonitorSensi(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscInt numcost, Vec *lambda, Vec *mu, PetscViewerAndFormat *vf)
1077: {
1078: PetscViewer viewer = vf->viewer;
1081: PetscViewerPushFormat(viewer, vf->format);
1082: VecView(lambda[0], viewer);
1083: PetscViewerPopFormat(viewer);
1084: return 0;
1085: }
1087: /*@C
1088: TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
1090: Collective
1092: Input Parameters:
1093: + ts - `TS` object you wish to monitor
1094: . name - the monitor type one is seeking
1095: . help - message indicating what monitoring is done
1096: . manual - manual page for the monitor
1097: . monitor - the monitor function
1098: - monitorsetup - a function that is called once ONLY if the user selected this monitor that may set additional features of the `TS` or `PetscViewer` objects
1100: Level: developer
1102: .seealso: [](chapter_ts), `PetscOptionsGetViewer()`, `PetscOptionsGetReal()`, `PetscOptionsHasName()`, `PetscOptionsGetString()`,
1103: `PetscOptionsGetIntArray()`, `PetscOptionsGetRealArray()`, `PetscOptionsBool()`
1104: `PetscOptionsInt()`, `PetscOptionsString()`, `PetscOptionsReal()`, `PetscOptionsBool()`,
1105: `PetscOptionsName()`, `PetscOptionsBegin()`, `PetscOptionsEnd()`, `PetscOptionsHeadBegin()`,
1106: `PetscOptionsStringArray()`, `PetscOptionsRealArray()`, `PetscOptionsScalar()`,
1107: `PetscOptionsBoolGroupBegin()`, `PetscOptionsBoolGroup()`, `PetscOptionsBoolGroupEnd()`,
1108: `PetscOptionsFList()`, `PetscOptionsEList()`
1109: @*/
1110: PetscErrorCode TSAdjointMonitorSetFromOptions(TS ts, const char name[], const char help[], const char manual[], PetscErrorCode (*monitor)(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, PetscViewerAndFormat *), PetscErrorCode (*monitorsetup)(TS, PetscViewerAndFormat *))
1111: {
1112: PetscViewer viewer;
1113: PetscViewerFormat format;
1114: PetscBool flg;
1116: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, name, &viewer, &format, &flg);
1117: if (flg) {
1118: PetscViewerAndFormat *vf;
1119: PetscViewerAndFormatCreate(viewer, format, &vf);
1120: PetscObjectDereference((PetscObject)viewer);
1121: if (monitorsetup) (*monitorsetup)(ts, vf);
1122: TSAdjointMonitorSet(ts, (PetscErrorCode(*)(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, void *))monitor, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1123: }
1124: return 0;
1125: }
1127: /*@C
1128: TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
1129: timestep to display the iteration's progress.
1131: Logically Collective
1133: Input Parameters:
1134: + ts - the `TS` context obtained from `TSCreate()`
1135: . adjointmonitor - monitoring routine
1136: . adjointmctx - [optional] user-defined context for private data for the
1137: monitor routine (use NULL if no context is desired)
1138: - adjointmonitordestroy - [optional] routine that frees monitor context
1139: (may be NULL)
1141: Calling sequence of monitor:
1142: $ int adjointmonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *adjointmctx)
1144: + ts - the `TS` context
1145: . steps - iteration number (after the final time step the monitor routine is called with a step of -1, this is at the final time which may have
1146: been interpolated to)
1147: . time - current time
1148: . u - current iterate
1149: . numcost - number of cost functionos
1150: . lambda - sensitivities to initial conditions
1151: . mu - sensitivities to parameters
1152: - adjointmctx - [optional] adjoint monitoring context
1154: Level: intermediate
1156: Note:
1157: This routine adds an additional monitor to the list of monitors that
1158: already has been loaded.
1160: Fortran Note:
1161: Only a single monitor function can be set for each TS object
1163: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSAdjointMonitorCancel()`
1164: @*/
1165: PetscErrorCode TSAdjointMonitorSet(TS ts, PetscErrorCode (*adjointmonitor)(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, void *), void *adjointmctx, PetscErrorCode (*adjointmdestroy)(void **))
1166: {
1167: PetscInt i;
1168: PetscBool identical;
1171: for (i = 0; i < ts->numbermonitors; i++) {
1172: PetscMonitorCompare((PetscErrorCode(*)(void))adjointmonitor, adjointmctx, adjointmdestroy, (PetscErrorCode(*)(void))ts->adjointmonitor[i], ts->adjointmonitorcontext[i], ts->adjointmonitordestroy[i], &identical);
1173: if (identical) return 0;
1174: }
1176: ts->adjointmonitor[ts->numberadjointmonitors] = adjointmonitor;
1177: ts->adjointmonitordestroy[ts->numberadjointmonitors] = adjointmdestroy;
1178: ts->adjointmonitorcontext[ts->numberadjointmonitors++] = (void *)adjointmctx;
1179: return 0;
1180: }
1182: /*@C
1183: TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.
1185: Logically Collective
1187: Input Parameters:
1188: . ts - the `TS` context obtained from `TSCreate()`
1190: Notes:
1191: There is no way to remove a single, specific monitor.
1193: Level: intermediate
1195: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSAdjointMonitorSet()`
1196: @*/
1197: PetscErrorCode TSAdjointMonitorCancel(TS ts)
1198: {
1199: PetscInt i;
1202: for (i = 0; i < ts->numberadjointmonitors; i++) {
1203: if (ts->adjointmonitordestroy[i]) (*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]);
1204: }
1205: ts->numberadjointmonitors = 0;
1206: return 0;
1207: }
1209: /*@C
1210: TSAdjointMonitorDefault - the default monitor of adjoint computations
1212: Level: intermediate
1214: .seealso: [](chapter_ts), `TS`, `TSAdjointSolve()`, `TSAdjointMonitorSet()`
1215: @*/
1216: PetscErrorCode TSAdjointMonitorDefault(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscInt numcost, Vec *lambda, Vec *mu, PetscViewerAndFormat *vf)
1217: {
1218: PetscViewer viewer = vf->viewer;
1221: PetscViewerPushFormat(viewer, vf->format);
1222: PetscViewerASCIIAddTab(viewer, ((PetscObject)ts)->tablevel);
1223: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " TS dt %g time %g%s", step, (double)ts->time_step, (double)ptime, ts->steprollback ? " (r)\n" : "\n");
1224: PetscViewerASCIISubtractTab(viewer, ((PetscObject)ts)->tablevel);
1225: PetscViewerPopFormat(viewer);
1226: return 0;
1227: }
1229: /*@C
1230: TSAdjointMonitorDrawSensi - Monitors progress of the adjoint `TS` solvers by calling
1231: `VecView()` for the sensitivities to initial states at each timestep
1233: Collective
1235: Input Parameters:
1236: + ts - the `TS` context
1237: . step - current time-step
1238: . ptime - current time
1239: . u - current state
1240: . numcost - number of cost functions
1241: . lambda - sensitivities to initial conditions
1242: . mu - sensitivities to parameters
1243: - dummy - either a viewer or NULL
1245: Level: intermediate
1247: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSAdjointMonitorSet()`, `TSAdjointMonitorDefault()`, `VecView()`
1248: @*/
1249: PetscErrorCode TSAdjointMonitorDrawSensi(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscInt numcost, Vec *lambda, Vec *mu, void *dummy)
1250: {
1251: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
1252: PetscDraw draw;
1253: PetscReal xl, yl, xr, yr, h;
1254: char time[32];
1256: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return 0;
1258: VecView(lambda[0], ictx->viewer);
1259: PetscViewerDrawGetDraw(ictx->viewer, 0, &draw);
1260: PetscSNPrintf(time, 32, "Timestep %d Time %g", (int)step, (double)ptime);
1261: PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr);
1262: h = yl + .95 * (yr - yl);
1263: PetscDrawStringCentered(draw, .5 * (xl + xr), h, PETSC_DRAW_BLACK, time);
1264: PetscDrawFlush(draw);
1265: return 0;
1266: }
1268: /*
1269: TSAdjointSetFromOptions - Sets various `TS` adjoint parameters from user options.
1271: Collective
1273: Input Parameter:
1274: . ts - the `TS` context
1276: Options Database Keys:
1277: + -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
1278: . -ts_adjoint_monitor - print information at each adjoint time step
1279: - -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically
1281: Level: developer
1283: Note:
1284: This is not normally called directly by users
1286: .seealso: [](chapter_ts), `TSSetSaveTrajectory()`, `TSTrajectorySetUp()`
1287: */
1288: PetscErrorCode TSAdjointSetFromOptions(TS ts, PetscOptionItems *PetscOptionsObject)
1289: {
1290: PetscBool tflg, opt;
1293: PetscOptionsHeadBegin(PetscOptionsObject, "TS Adjoint options");
1294: tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
1295: PetscOptionsBool("-ts_adjoint_solve", "Solve the adjoint problem immediately after solving the forward problem", "", tflg, &tflg, &opt);
1296: if (opt) {
1297: TSSetSaveTrajectory(ts);
1298: ts->adjoint_solve = tflg;
1299: }
1300: TSAdjointMonitorSetFromOptions(ts, "-ts_adjoint_monitor", "Monitor adjoint timestep size", "TSAdjointMonitorDefault", TSAdjointMonitorDefault, NULL);
1301: TSAdjointMonitorSetFromOptions(ts, "-ts_adjoint_monitor_sensi", "Monitor sensitivity in the adjoint computation", "TSAdjointMonitorSensi", TSAdjointMonitorSensi, NULL);
1302: opt = PETSC_FALSE;
1303: PetscOptionsName("-ts_adjoint_monitor_draw_sensi", "Monitor adjoint sensitivities (lambda only) graphically", "TSAdjointMonitorDrawSensi", &opt);
1304: if (opt) {
1305: TSMonitorDrawCtx ctx;
1306: PetscInt howoften = 1;
1308: PetscOptionsInt("-ts_adjoint_monitor_draw_sensi", "Monitor adjoint sensitivities (lambda only) graphically", "TSAdjointMonitorDrawSensi", howoften, &howoften, NULL);
1309: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
1310: TSAdjointMonitorSet(ts, TSAdjointMonitorDrawSensi, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
1311: }
1312: return 0;
1313: }
1315: /*@
1316: TSAdjointStep - Steps one time step backward in the adjoint run
1318: Collective
1320: Input Parameter:
1321: . ts - the `TS` context obtained from `TSCreate()`
1323: Level: intermediate
1325: .seealso: [](chapter_ts), `TSAdjointSetUp()`, `TSAdjointSolve()`
1326: @*/
1327: PetscErrorCode TSAdjointStep(TS ts)
1328: {
1329: DM dm;
1332: TSGetDM(ts, &dm);
1333: TSAdjointSetUp(ts);
1334: ts->steps--; /* must decrease the step index before the adjoint step is taken. */
1336: ts->reason = TS_CONVERGED_ITERATING;
1337: ts->ptime_prev = ts->ptime;
1338: PetscLogEventBegin(TS_AdjointStep, ts, 0, 0, 0);
1339: PetscUseTypeMethod(ts, adjointstep);
1340: PetscLogEventEnd(TS_AdjointStep, ts, 0, 0, 0);
1341: ts->adjoint_steps++;
1343: if (ts->reason < 0) {
1345: } else if (!ts->reason) {
1346: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1347: }
1348: return 0;
1349: }
1351: /*@
1352: TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE
1354: Collective
1355: `
1356: Input Parameter:
1357: . ts - the `TS` context obtained from `TSCreate()`
1359: Options Database Key:
1360: . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial values
1362: Level: intermediate
1364: Notes:
1365: This must be called after a call to `TSSolve()` that solves the forward problem
1367: By default this will integrate back to the initial time, one can use `TSAdjointSetSteps()` to step back to a later time
1369: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSCreate()`, `TSSetCostGradients()`, `TSSetSolution()`, `TSAdjointStep()`
1370: @*/
1371: PetscErrorCode TSAdjointSolve(TS ts)
1372: {
1373: static PetscBool cite = PETSC_FALSE;
1374: #if defined(TSADJOINT_STAGE)
1375: PetscLogStage adjoint_stage;
1376: #endif
1379: PetscCall(PetscCitationsRegister("@article{Zhang2022tsadjoint,\n"
1380: " title = {{PETSc TSAdjoint: A Discrete Adjoint ODE Solver for First-Order and Second-Order Sensitivity Analysis}},\n"
1381: " author = {Zhang, Hong and Constantinescu, Emil M. and Smith, Barry F.},\n"
1382: " journal = {SIAM Journal on Scientific Computing},\n"
1383: " volume = {44},\n"
1384: " number = {1},\n"
1385: " pages = {C1-C24},\n"
1386: " doi = {10.1137/21M140078X},\n"
1387: " year = {2022}\n}\n",
1388: &cite));
1389: #if defined(TSADJOINT_STAGE)
1390: PetscLogStageRegister("TSAdjoint", &adjoint_stage);
1391: PetscLogStagePush(adjoint_stage);
1392: #endif
1393: TSAdjointSetUp(ts);
1395: /* reset time step and iteration counters */
1396: ts->adjoint_steps = 0;
1397: ts->ksp_its = 0;
1398: ts->snes_its = 0;
1399: ts->num_snes_failures = 0;
1400: ts->reject = 0;
1401: ts->reason = TS_CONVERGED_ITERATING;
1403: if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->steps;
1404: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1406: while (!ts->reason) {
1407: TSTrajectoryGet(ts->trajectory, ts, ts->steps, &ts->ptime);
1408: TSAdjointMonitor(ts, ts->steps, ts->ptime, ts->vec_sol, ts->numcost, ts->vecs_sensi, ts->vecs_sensip);
1409: TSAdjointEventHandler(ts);
1410: TSAdjointStep(ts);
1411: if ((ts->vec_costintegral || ts->quadraturets) && !ts->costintegralfwd) TSAdjointCostIntegral(ts);
1412: }
1413: if (!ts->steps) {
1414: TSTrajectoryGet(ts->trajectory, ts, ts->steps, &ts->ptime);
1415: TSAdjointMonitor(ts, ts->steps, ts->ptime, ts->vec_sol, ts->numcost, ts->vecs_sensi, ts->vecs_sensip);
1416: }
1417: ts->solvetime = ts->ptime;
1418: TSTrajectoryViewFromOptions(ts->trajectory, NULL, "-ts_trajectory_view");
1419: VecViewFromOptions(ts->vecs_sensi[0], (PetscObject)ts, "-ts_adjoint_view_solution");
1420: ts->adjoint_max_steps = 0;
1421: #if defined(TSADJOINT_STAGE)
1422: PetscLogStagePop();
1423: #endif
1424: return 0;
1425: }
1427: /*@C
1428: TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using `TSAdjointMonitorSet()`
1430: Collective
1432: Input Parameters:
1433: + ts - time stepping context obtained from `TSCreate()`
1434: . step - step number that has just completed
1435: . ptime - model time of the state
1436: . u - state at the current model time
1437: . numcost - number of cost functions (dimension of lambda or mu)
1438: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
1439: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
1441: Level: developer
1443: Note:
1444: `TSAdjointMonitor()` is typically used automatically within the time stepping implementations.
1445: Users would almost never call this routine directly.
1447: .seealso: `TSAdjointMonitorSet()`, `TSAdjointSolve()`
1448: @*/
1449: PetscErrorCode TSAdjointMonitor(TS ts, PetscInt step, PetscReal ptime, Vec u, PetscInt numcost, Vec *lambda, Vec *mu)
1450: {
1451: PetscInt i, n = ts->numberadjointmonitors;
1455: VecLockReadPush(u);
1456: for (i = 0; i < n; i++) (*ts->adjointmonitor[i])(ts, step, ptime, u, numcost, lambda, mu, ts->adjointmonitorcontext[i]);
1457: VecLockReadPop(u);
1458: return 0;
1459: }
1461: /*@
1462: TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.
1464: Collective
1466: Input Parameter:
1467: . ts - time stepping context
1469: Level: advanced
1471: Notes:
1472: This function cannot be called until `TSAdjointStep()` has been completed.
1474: .seealso: [](chapter_ts), `TSAdjointSolve()`, `TSAdjointStep()`
1475: @*/
1476: PetscErrorCode TSAdjointCostIntegral(TS ts)
1477: {
1479: PetscUseTypeMethod(ts, adjointintegral);
1480: return 0;
1481: }
1483: /* ------------------ Forward (tangent linear) sensitivity ------------------*/
1485: /*@
1486: TSForwardSetUp - Sets up the internal data structures for the later use
1487: of forward sensitivity analysis
1489: Collective
1491: Input Parameter:
1492: . ts - the `TS` context obtained from `TSCreate()`
1494: Level: advanced
1496: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSDestroy()`, `TSSetUp()`
1497: @*/
1498: PetscErrorCode TSForwardSetUp(TS ts)
1499: {
1501: if (ts->forwardsetupcalled) return 0;
1502: PetscTryTypeMethod(ts, forwardsetup);
1503: VecDuplicate(ts->vec_sol, &ts->vec_sensip_col);
1504: ts->forwardsetupcalled = PETSC_TRUE;
1505: return 0;
1506: }
1508: /*@
1509: TSForwardReset - Reset the internal data structures used by forward sensitivity analysis
1511: Collective
1513: Input Parameter:
1514: . ts - the `TS` context obtained from `TSCreate()`
1516: Level: advanced
1518: .seealso: [](chapter_ts), `TSCreate()`, `TSDestroy()`, `TSForwardSetUp()`
1519: @*/
1520: PetscErrorCode TSForwardReset(TS ts)
1521: {
1522: TS quadts = ts->quadraturets;
1525: PetscTryTypeMethod(ts, forwardreset);
1526: MatDestroy(&ts->mat_sensip);
1527: if (quadts) MatDestroy(&quadts->mat_sensip);
1528: VecDestroy(&ts->vec_sensip_col);
1529: ts->forward_solve = PETSC_FALSE;
1530: ts->forwardsetupcalled = PETSC_FALSE;
1531: return 0;
1532: }
1534: /*@
1535: TSForwardSetIntegralGradients - Set the vectors holding forward sensitivities of the integral term.
1537: Input Parameters:
1538: + ts - the `TS` context obtained from `TSCreate()`
1539: . numfwdint - number of integrals
1540: - vp - the vectors containing the gradients for each integral w.r.t. parameters
1542: Level: deprecated
1544: .seealso: [](chapter_ts), `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1545: @*/
1546: PetscErrorCode TSForwardSetIntegralGradients(TS ts, PetscInt numfwdint, Vec *vp)
1547: {
1550: if (!ts->numcost) ts->numcost = numfwdint;
1552: ts->vecs_integral_sensip = vp;
1553: return 0;
1554: }
1556: /*@
1557: TSForwardGetIntegralGradients - Returns the forward sensitivities ofthe integral term.
1559: Input Parameter:
1560: . ts - the `TS` context obtained from `TSCreate()`
1562: Output Parameter:
1563: . vp - the vectors containing the gradients for each integral w.r.t. parameters
1565: Level: deprecated
1567: .seealso: [](chapter_ts), `TSForwardSetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1568: @*/
1569: PetscErrorCode TSForwardGetIntegralGradients(TS ts, PetscInt *numfwdint, Vec **vp)
1570: {
1573: if (numfwdint) *numfwdint = ts->numcost;
1574: if (vp) *vp = ts->vecs_integral_sensip;
1575: return 0;
1576: }
1578: /*@
1579: TSForwardStep - Compute the forward sensitivity for one time step.
1581: Collective
1583: Input Parameter:
1584: . ts - time stepping context
1586: Level: advanced
1588: Notes:
1589: This function cannot be called until `TSStep()` has been completed.
1591: .seealso: [](chapter_ts), `TSForwardSetSensitivities()`, `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardSetUp()`
1592: @*/
1593: PetscErrorCode TSForwardStep(TS ts)
1594: {
1596: PetscLogEventBegin(TS_ForwardStep, ts, 0, 0, 0);
1597: PetscUseTypeMethod(ts, forwardstep);
1598: PetscLogEventEnd(TS_ForwardStep, ts, 0, 0, 0);
1600: return 0;
1601: }
1603: /*@
1604: TSForwardSetSensitivities - Sets the initial value of the trajectory sensitivities of solution w.r.t. the problem parameters and initial values.
1606: Logically Collective
1608: Input Parameters:
1609: + ts - the `TS` context obtained from `TSCreate()`
1610: . nump - number of parameters
1611: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1613: Level: beginner
1615: Notes:
1616: Forward sensitivity is also called 'trajectory sensitivity' in some fields such as power systems.
1617: This function turns on a flag to trigger `TSSolve()` to compute forward sensitivities automatically.
1618: You must call this function before `TSSolve()`.
1619: The entries in the sensitivity matrix must be correctly initialized with the values S = dy/dp|startingtime.
1621: .seealso: [](chapter_ts), `TSForwardGetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1622: @*/
1623: PetscErrorCode TSForwardSetSensitivities(TS ts, PetscInt nump, Mat Smat)
1624: {
1627: ts->forward_solve = PETSC_TRUE;
1628: if (nump == PETSC_DEFAULT) {
1629: MatGetSize(Smat, NULL, &ts->num_parameters);
1630: } else ts->num_parameters = nump;
1631: PetscObjectReference((PetscObject)Smat);
1632: MatDestroy(&ts->mat_sensip);
1633: ts->mat_sensip = Smat;
1634: return 0;
1635: }
1637: /*@
1638: TSForwardGetSensitivities - Returns the trajectory sensitivities
1640: Not Collective, but Smat returned is parallel if ts is parallel
1642: Output Parameters:
1643: + ts - the `TS` context obtained from `TSCreate()`
1644: . nump - number of parameters
1645: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1647: Level: intermediate
1649: .seealso: [](chapter_ts), `TSForwardSetSensitivities()`, `TSForwardSetIntegralGradients()`, `TSForwardGetIntegralGradients()`, `TSForwardStep()`
1650: @*/
1651: PetscErrorCode TSForwardGetSensitivities(TS ts, PetscInt *nump, Mat *Smat)
1652: {
1654: if (nump) *nump = ts->num_parameters;
1655: if (Smat) *Smat = ts->mat_sensip;
1656: return 0;
1657: }
1659: /*@
1660: TSForwardCostIntegral - Evaluate the cost integral in the forward run.
1662: Collective
1664: Input Parameter:
1665: . ts - time stepping context
1667: Level: advanced
1669: Note:
1670: This function cannot be called until `TSStep()` has been completed.
1672: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSAdjointCostIntegral()`
1673: @*/
1674: PetscErrorCode TSForwardCostIntegral(TS ts)
1675: {
1677: PetscUseTypeMethod(ts, forwardintegral);
1678: return 0;
1679: }
1681: /*@
1682: TSForwardSetInitialSensitivities - Set initial values for tangent linear sensitivities
1684: Collective
1686: Input Parameters:
1687: + ts - the `TS` context obtained from `TSCreate()`
1688: - didp - parametric sensitivities of the initial condition
1690: Level: intermediate
1692: Notes:
1693: `TSSolve()` allows users to pass the initial solution directly to `TS`. But the tangent linear variables cannot be initialized in this way.
1694: This function is used to set initial values for tangent linear variables.
1696: .seealso: [](chapter_ts), `TS`, `TSForwardSetSensitivities()`
1697: @*/
1698: PetscErrorCode TSForwardSetInitialSensitivities(TS ts, Mat didp)
1699: {
1702: if (!ts->mat_sensip) TSForwardSetSensitivities(ts, PETSC_DEFAULT, didp);
1703: return 0;
1704: }
1706: /*@
1707: TSForwardGetStages - Get the number of stages and the tangent linear sensitivities at the intermediate stages
1709: Input Parameter:
1710: . ts - the `TS` context obtained from `TSCreate()`
1712: Output Parameters:
1713: + ns - number of stages
1714: - S - tangent linear sensitivities at the intermediate stages
1716: Level: advanced
1718: .seealso: `TS`
1719: @*/
1720: PetscErrorCode TSForwardGetStages(TS ts, PetscInt *ns, Mat **S)
1721: {
1724: if (!ts->ops->getstages) *S = NULL;
1725: else PetscUseTypeMethod(ts, forwardgetstages, ns, S);
1726: return 0;
1727: }
1729: /*@
1730: TSCreateQuadratureTS - Create a sub-`TS` that evaluates integrals over time
1732: Input Parameters:
1733: + ts - the `TS` context obtained from `TSCreate()`
1734: - fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1736: Output Parameters:
1737: . quadts - the child `TS` context
1739: Level: intermediate
1741: .seealso: [](chapter_ts), `TSGetQuadratureTS()`
1742: @*/
1743: PetscErrorCode TSCreateQuadratureTS(TS ts, PetscBool fwd, TS *quadts)
1744: {
1745: char prefix[128];
1749: TSDestroy(&ts->quadraturets);
1750: TSCreate(PetscObjectComm((PetscObject)ts), &ts->quadraturets);
1751: PetscObjectIncrementTabLevel((PetscObject)ts->quadraturets, (PetscObject)ts, 1);
1752: PetscSNPrintf(prefix, sizeof(prefix), "%squad_", ((PetscObject)ts)->prefix ? ((PetscObject)ts)->prefix : "");
1753: TSSetOptionsPrefix(ts->quadraturets, prefix);
1754: *quadts = ts->quadraturets;
1756: if (ts->numcost) {
1757: VecCreateSeq(PETSC_COMM_SELF, ts->numcost, &(*quadts)->vec_sol);
1758: } else {
1759: VecCreateSeq(PETSC_COMM_SELF, 1, &(*quadts)->vec_sol);
1760: }
1761: ts->costintegralfwd = fwd;
1762: return 0;
1763: }
1765: /*@
1766: TSGetQuadratureTS - Return the sub-`TS` that evaluates integrals over time
1768: Input Parameter:
1769: . ts - the `TS` context obtained from `TSCreate()`
1771: Output Parameters:
1772: + fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1773: - quadts - the child `TS` context
1775: Level: intermediate
1777: .seealso: [](chapter_ts), `TSCreateQuadratureTS()`
1778: @*/
1779: PetscErrorCode TSGetQuadratureTS(TS ts, PetscBool *fwd, TS *quadts)
1780: {
1782: if (fwd) *fwd = ts->costintegralfwd;
1783: if (quadts) *quadts = ts->quadraturets;
1784: return 0;
1785: }
1787: /*@
1788: TSComputeSNESJacobian - Compute the Jacobian needed for the `SNESSolve()` in `TS`
1790: Collective
1792: Input Parameters:
1793: + ts - the `TS` context obtained from `TSCreate()`
1794: - x - state vector
1796: Output Parameters:
1797: + J - Jacobian matrix
1798: - Jpre - preconditioning matrix for J (may be same as J)
1800: Level: developer
1802: Note:
1803: Uses finite differencing when `TS` Jacobian is not available.
1805: .seealso: `SNES`, `TS`, `SNESSetJacobian()`, TSSetRHSJacobian()`, `TSSetIJacobian()`
1806: @*/
1807: PetscErrorCode TSComputeSNESJacobian(TS ts, Vec x, Mat J, Mat Jpre)
1808: {
1809: SNES snes = ts->snes;
1810: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *) = NULL;
1812: /*
1813: Unlike implicit methods, explicit methods do not have SNESMatFDColoring in the snes object
1814: because SNESSolve() has not been called yet; so querying SNESMatFDColoring does not work for
1815: explicit methods. Instead, we check the Jacobian compute function directly to determine if FD
1816: coloring is used.
1817: */
1818: SNESGetJacobian(snes, NULL, NULL, &jac, NULL);
1819: if (jac == SNESComputeJacobianDefaultColor) {
1820: Vec f;
1821: SNESSetSolution(snes, x);
1822: SNESGetFunction(snes, &f, NULL, NULL);
1823: /* Force MatFDColoringApply to evaluate the SNES residual function for the base vector */
1824: SNESComputeFunction(snes, x, f);
1825: }
1826: SNESComputeJacobian(snes, x, J, Jpre);
1827: return 0;
1828: }