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: {
43: ts->rhsjacobianp = func;
44: ts->rhsjacobianpctx = ctx;
45: if (Amat) {
46: PetscObjectReference((PetscObject)Amat);
47: MatDestroy(&ts->Jacprhs);
48: ts->Jacprhs = Amat;
49: }
50: return(0);
51: }
53: /*@C
54: 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.
56: Logically Collective on TS
58: Input Parameter:
59: . ts - TS context obtained from TSCreate()
61: Output Parameters:
62: + Amat - JacobianP matrix
63: . func - function
64: - ctx - [optional] user-defined function context
66: Calling sequence of func:
67: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
68: + t - current timestep
69: . U - input vector (current ODE solution)
70: . A - output matrix
71: - ctx - [optional] user-defined function context
73: Level: intermediate
75: Notes:
76: 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
78: .seealso: TSSetRHSJacobianP()
79: @*/
80: PetscErrorCode TSGetRHSJacobianP(TS ts,Mat *Amat,PetscErrorCode (**func)(TS,PetscReal,Vec,Mat,void*),void **ctx)
81: {
83: if (func) *func = ts->rhsjacobianp;
84: if (ctx) *ctx = ts->rhsjacobianpctx;
85: if (Amat) *Amat = ts->Jacprhs;
86: return(0);
87: }
89: /*@C
90: TSComputeRHSJacobianP - Runs the user-defined JacobianP function.
92: Collective on TS
94: Input Parameters:
95: . ts - The TS context obtained from TSCreate()
97: Level: developer
99: .seealso: TSSetRHSJacobianP()
100: @*/
101: PetscErrorCode TSComputeRHSJacobianP(TS ts,PetscReal t,Vec U,Mat Amat)
102: {
106: if (!Amat) return(0);
110: PetscStackPush("TS user JacobianP function for sensitivity analysis");
111: (*ts->rhsjacobianp)(ts,t,U,Amat,ts->rhsjacobianpctx);
112: PetscStackPop;
113: return(0);
114: }
116: /*@C
117: 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.
119: Logically Collective on TS
121: Input Parameters:
122: + ts - TS context obtained from TSCreate()
123: . Amat - JacobianP matrix
124: . func - function
125: - ctx - [optional] user-defined function context
127: Calling sequence of func:
128: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
129: + t - current timestep
130: . U - input vector (current ODE solution)
131: . Udot - time derivative of state vector
132: . shift - shift to apply, see note below
133: . A - output matrix
134: - ctx - [optional] user-defined function context
136: Level: intermediate
138: Notes:
139: 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
141: .seealso:
142: @*/
143: PetscErrorCode TSSetIJacobianP(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Vec,PetscReal,Mat,void*),void *ctx)
144: {
151: ts->ijacobianp = func;
152: ts->ijacobianpctx = ctx;
153: if (Amat) {
154: PetscObjectReference((PetscObject)Amat);
155: MatDestroy(&ts->Jacp);
156: ts->Jacp = Amat;
157: }
158: return(0);
159: }
161: /*@C
162: TSComputeIJacobianP - Runs the user-defined IJacobianP function.
164: Collective on TS
166: Input Parameters:
167: + ts - the TS context
168: . t - current timestep
169: . U - state vector
170: . Udot - time derivative of state vector
171: . shift - shift to apply, see note below
172: - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate
174: Output Parameters:
175: . A - Jacobian matrix
177: Level: developer
179: .seealso: TSSetIJacobianP()
180: @*/
181: PetscErrorCode TSComputeIJacobianP(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat Amat,PetscBool imex)
182: {
186: if (!Amat) return(0);
191: PetscLogEventBegin(TS_JacobianPEval,ts,U,Amat,0);
192: if (ts->ijacobianp) {
193: PetscStackPush("TS user JacobianP function for sensitivity analysis");
194: (*ts->ijacobianp)(ts,t,U,Udot,shift,Amat,ts->ijacobianpctx);
195: PetscStackPop;
196: }
197: if (imex) {
198: if (!ts->ijacobianp) { /* system was written as Udot = G(t,U) */
199: PetscBool assembled;
200: MatZeroEntries(Amat);
201: MatAssembled(Amat,&assembled);
202: if (!assembled) {
203: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
204: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
205: }
206: }
207: } else {
208: if (ts->rhsjacobianp) {
209: TSComputeRHSJacobianP(ts,t,U,ts->Jacprhs);
210: }
211: if (ts->Jacprhs == Amat) { /* No IJacobian, so we only have the RHS matrix */
212: MatScale(Amat,-1);
213: } else if (ts->Jacprhs) { /* Both IJacobian and RHSJacobian */
214: MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
215: if (!ts->ijacobianp) { /* No IJacobianp provided, but we have a separate RHS matrix */
216: MatZeroEntries(Amat);
217: }
218: MatAXPY(Amat,-1,ts->Jacprhs,axpy);
219: }
220: }
221: PetscLogEventEnd(TS_JacobianPEval,ts,U,Amat,0);
222: return(0);
223: }
225: /*@C
226: TSSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions
228: Logically Collective on TS
230: Input Parameters:
231: + ts - the TS context obtained from TSCreate()
232: . numcost - number of gradients to be computed, this is the number of cost functions
233: . costintegral - vector that stores the integral values
234: . rf - routine for evaluating the integrand function
235: . drduf - function that computes the gradients of the r's with respect to u
236: . 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)
237: . fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
238: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
240: Calling sequence of rf:
241: $ PetscErrorCode rf(TS ts,PetscReal t,Vec U,Vec F,void *ctx);
243: Calling sequence of drduf:
244: $ PetscErroCode drduf(TS ts,PetscReal t,Vec U,Vec *dRdU,void *ctx);
246: Calling sequence of drdpf:
247: $ PetscErroCode drdpf(TS ts,PetscReal t,Vec U,Vec *dRdP,void *ctx);
249: Level: deprecated
251: Notes:
252: For optimization there is usually a single cost function (numcost = 1). For sensitivities there may be multiple cost functions
254: .seealso: TSSetRHSJacobianP(), TSGetCostGradients(), TSSetCostGradients()
255: @*/
256: PetscErrorCode TSSetCostIntegrand(TS ts,PetscInt numcost,Vec costintegral,PetscErrorCode (*rf)(TS,PetscReal,Vec,Vec,void*),
257: PetscErrorCode (*drduf)(TS,PetscReal,Vec,Vec*,void*),
258: PetscErrorCode (*drdpf)(TS,PetscReal,Vec,Vec*,void*),
259: PetscBool fwd,void *ctx)
260: {
266: if (ts->numcost && ts->numcost!=numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2nd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostGradients() or TSForwardSetIntegralGradients()");
267: if (!ts->numcost) ts->numcost=numcost;
269: if (costintegral) {
270: PetscObjectReference((PetscObject)costintegral);
271: VecDestroy(&ts->vec_costintegral);
272: ts->vec_costintegral = costintegral;
273: } else {
274: if (!ts->vec_costintegral) { /* Create a seq vec if user does not provide one */
275: VecCreateSeq(PETSC_COMM_SELF,numcost,&ts->vec_costintegral);
276: } else {
277: VecSet(ts->vec_costintegral,0.0);
278: }
279: }
280: if (!ts->vec_costintegrand) {
281: VecDuplicate(ts->vec_costintegral,&ts->vec_costintegrand);
282: } else {
283: VecSet(ts->vec_costintegrand,0.0);
284: }
285: ts->costintegralfwd = fwd; /* Evaluate the cost integral in forward run if fwd is true */
286: ts->costintegrand = rf;
287: ts->costintegrandctx = ctx;
288: ts->drdufunction = drduf;
289: ts->drdpfunction = drdpf;
290: return(0);
291: }
293: /*@C
294: TSGetCostIntegral - Returns the values of the integral term in the cost functions.
295: It is valid to call the routine after a backward run.
297: Not Collective
299: Input Parameter:
300: . ts - the TS context obtained from TSCreate()
302: Output Parameter:
303: . v - the vector containing the integrals for each cost function
305: Level: intermediate
307: .seealso: TSSetCostIntegrand()
309: @*/
310: PetscErrorCode TSGetCostIntegral(TS ts,Vec *v)
311: {
312: TS quadts;
318: TSGetQuadratureTS(ts,NULL,&quadts);
319: *v = quadts->vec_sol;
320: return(0);
321: }
323: /*@C
324: TSComputeCostIntegrand - Evaluates the integral function in the cost functions.
326: Input Parameters:
327: + ts - the TS context
328: . t - current time
329: - U - state vector, i.e. current solution
331: Output Parameter:
332: . Q - vector of size numcost to hold the outputs
334: Notes:
335: Most users should not need to explicitly call this routine, as it
336: is used internally within the sensitivity analysis context.
338: Level: deprecated
340: .seealso: TSSetCostIntegrand()
341: @*/
342: PetscErrorCode TSComputeCostIntegrand(TS ts,PetscReal t,Vec U,Vec Q)
343: {
351: PetscLogEventBegin(TS_FunctionEval,ts,U,Q,0);
352: if (ts->costintegrand) {
353: PetscStackPush("TS user integrand in the cost function");
354: (*ts->costintegrand)(ts,t,U,Q,ts->costintegrandctx);
355: PetscStackPop;
356: } else {
357: VecZeroEntries(Q);
358: }
360: PetscLogEventEnd(TS_FunctionEval,ts,U,Q,0);
361: return(0);
362: }
364: /*@C
365: TSComputeDRDUFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobian()
367: Level: deprecated
369: @*/
370: PetscErrorCode TSComputeDRDUFunction(TS ts,PetscReal t,Vec U,Vec *DRDU)
371: {
375: if (!DRDU) return(0);
379: PetscStackPush("TS user DRDU function for sensitivity analysis");
380: (*ts->drdufunction)(ts,t,U,DRDU,ts->costintegrandctx);
381: PetscStackPop;
382: return(0);
383: }
385: /*@C
386: TSComputeDRDPFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobianP()
388: Level: deprecated
390: @*/
391: PetscErrorCode TSComputeDRDPFunction(TS ts,PetscReal t,Vec U,Vec *DRDP)
392: {
396: if (!DRDP) return(0);
400: PetscStackPush("TS user DRDP function for sensitivity analysis");
401: (*ts->drdpfunction)(ts,t,U,DRDP,ts->costintegrandctx);
402: PetscStackPop;
403: return(0);
404: }
406: /*@C
407: 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.
409: Logically Collective on TS
411: Input Parameters:
412: + ts - TS context obtained from TSCreate()
413: . ihp1 - an array of vectors storing the result of vector-Hessian-vector product for F_UU
414: . hessianproductfunc1 - vector-Hessian-vector product function for F_UU
415: . ihp2 - an array of vectors storing the result of vector-Hessian-vector product for F_UP
416: . hessianproductfunc2 - vector-Hessian-vector product function for F_UP
417: . ihp3 - an array of vectors storing the result of vector-Hessian-vector product for F_PU
418: . hessianproductfunc3 - vector-Hessian-vector product function for F_PU
419: . ihp4 - an array of vectors storing the result of vector-Hessian-vector product for F_PP
420: - hessianproductfunc4 - vector-Hessian-vector product function for F_PP
422: Calling sequence of ihessianproductfunc:
423: $ ihessianproductfunc (TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx);
424: + t - current timestep
425: . U - input vector (current ODE solution)
426: . Vl - an array of input vectors to be left-multiplied with the Hessian
427: . Vr - input vector to be right-multiplied with the Hessian
428: . VHV - an array of output vectors for vector-Hessian-vector product
429: - ctx - [optional] user-defined function context
431: Level: intermediate
433: Notes:
434: The first Hessian function and the working array are required.
435: As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
436: $ Vl_n^T*F_UP*Vr
437: 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.
438: Each entry of F_UP corresponds to the derivative
439: $ F_UP[i][j][k] = \frac{\partial^2 F[i]}{\partial U[j] \partial P[k]}.
440: 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
441: $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * F_UP[i][j][k] * Vr[k]}
442: If the cost function is a scalar, there will be only one vector in Vl and VHV.
444: .seealso:
445: @*/
446: PetscErrorCode TSSetIHessianProduct(TS ts,Vec *ihp1,PetscErrorCode (*ihessianproductfunc1)(TS,PetscReal,Vec,Vec*,Vec,Vec*,void*),
447: Vec *ihp2,PetscErrorCode (*ihessianproductfunc2)(TS,PetscReal,Vec,Vec*,Vec,Vec*,void*),
448: Vec *ihp3,PetscErrorCode (*ihessianproductfunc3)(TS,PetscReal,Vec,Vec*,Vec,Vec*,void*),
449: Vec *ihp4,PetscErrorCode (*ihessianproductfunc4)(TS,PetscReal,Vec,Vec*,Vec,Vec*,void*),
450: void *ctx)
451: {
456: ts->ihessianproductctx = ctx;
457: if (ihp1) ts->vecs_fuu = ihp1;
458: if (ihp2) ts->vecs_fup = ihp2;
459: if (ihp3) ts->vecs_fpu = ihp3;
460: if (ihp4) ts->vecs_fpp = ihp4;
461: ts->ihessianproduct_fuu = ihessianproductfunc1;
462: ts->ihessianproduct_fup = ihessianproductfunc2;
463: ts->ihessianproduct_fpu = ihessianproductfunc3;
464: ts->ihessianproduct_fpp = ihessianproductfunc4;
465: return(0);
466: }
468: /*@C
469: TSComputeIHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Fuu.
471: Collective on TS
473: Input Parameters:
474: . ts - The TS context obtained from TSCreate()
476: Notes:
477: TSComputeIHessianProductFunctionUU() is typically used for sensitivity implementation,
478: so most users would not generally call this routine themselves.
480: Level: developer
482: .seealso: TSSetIHessianProduct()
483: @*/
484: PetscErrorCode TSComputeIHessianProductFunctionUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
485: {
489: if (!VHV) return(0);
493: if (ts->ihessianproduct_fuu) {
494: PetscStackPush("TS user IHessianProduct function 1 for sensitivity analysis");
495: (*ts->ihessianproduct_fuu)(ts,t,U,Vl,Vr,VHV,ts->ihessianproductctx);
496: PetscStackPop;
497: }
498: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
499: if (ts->rhshessianproduct_guu) {
500: PetscInt nadj;
501: TSComputeRHSHessianProductFunctionUU(ts,t,U,Vl,Vr,VHV);
502: for (nadj=0; nadj<ts->numcost; nadj++) {
503: VecScale(VHV[nadj],-1);
504: }
505: }
506: return(0);
507: }
509: /*@C
510: TSComputeIHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Fup.
512: Collective on TS
514: Input Parameters:
515: . ts - The TS context obtained from TSCreate()
517: Notes:
518: TSComputeIHessianProductFunctionUP() is typically used for sensitivity implementation,
519: so most users would not generally call this routine themselves.
521: Level: developer
523: .seealso: TSSetIHessianProduct()
524: @*/
525: PetscErrorCode TSComputeIHessianProductFunctionUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
526: {
530: if (!VHV) return(0);
534: if (ts->ihessianproduct_fup) {
535: PetscStackPush("TS user IHessianProduct function 2 for sensitivity analysis");
536: (*ts->ihessianproduct_fup)(ts,t,U,Vl,Vr,VHV,ts->ihessianproductctx);
537: PetscStackPop;
538: }
539: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
540: if (ts->rhshessianproduct_gup) {
541: PetscInt nadj;
542: TSComputeRHSHessianProductFunctionUP(ts,t,U,Vl,Vr,VHV);
543: for (nadj=0; nadj<ts->numcost; nadj++) {
544: VecScale(VHV[nadj],-1);
545: }
546: }
547: return(0);
548: }
550: /*@C
551: TSComputeIHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Fpu.
553: Collective on TS
555: Input Parameters:
556: . ts - The TS context obtained from TSCreate()
558: Notes:
559: TSComputeIHessianProductFunctionPU() is typically used for sensitivity implementation,
560: so most users would not generally call this routine themselves.
562: Level: developer
564: .seealso: TSSetIHessianProduct()
565: @*/
566: PetscErrorCode TSComputeIHessianProductFunctionPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
567: {
571: if (!VHV) return(0);
575: if (ts->ihessianproduct_fpu) {
576: PetscStackPush("TS user IHessianProduct function 3 for sensitivity analysis");
577: (*ts->ihessianproduct_fpu)(ts,t,U,Vl,Vr,VHV,ts->ihessianproductctx);
578: PetscStackPop;
579: }
580: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
581: if (ts->rhshessianproduct_gpu) {
582: PetscInt nadj;
583: TSComputeRHSHessianProductFunctionPU(ts,t,U,Vl,Vr,VHV);
584: for (nadj=0; nadj<ts->numcost; nadj++) {
585: VecScale(VHV[nadj],-1);
586: }
587: }
588: return(0);
589: }
591: /*@C
592: TSComputeIHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Fpp.
594: Collective on TS
596: Input Parameters:
597: . ts - The TS context obtained from TSCreate()
599: Notes:
600: TSComputeIHessianProductFunctionPP() is typically used for sensitivity implementation,
601: so most users would not generally call this routine themselves.
603: Level: developer
605: .seealso: TSSetIHessianProduct()
606: @*/
607: PetscErrorCode TSComputeIHessianProductFunctionPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
608: {
612: if (!VHV) return(0);
616: if (ts->ihessianproduct_fpp) {
617: PetscStackPush("TS user IHessianProduct function 3 for sensitivity analysis");
618: (*ts->ihessianproduct_fpp)(ts,t,U,Vl,Vr,VHV,ts->ihessianproductctx);
619: PetscStackPop;
620: }
621: /* does not consider IMEX for now, so either IHessian or RHSHessian will be calculated, using the same output VHV */
622: if (ts->rhshessianproduct_gpp) {
623: PetscInt nadj;
624: TSComputeRHSHessianProductFunctionPP(ts,t,U,Vl,Vr,VHV);
625: for (nadj=0; nadj<ts->numcost; nadj++) {
626: VecScale(VHV[nadj],-1);
627: }
628: }
629: return(0);
630: }
632: /*@C
633: 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.
635: Logically Collective on TS
637: Input Parameters:
638: + ts - TS context obtained from TSCreate()
639: . rhshp1 - an array of vectors storing the result of vector-Hessian-vector product for G_UU
640: . hessianproductfunc1 - vector-Hessian-vector product function for G_UU
641: . rhshp2 - an array of vectors storing the result of vector-Hessian-vector product for G_UP
642: . hessianproductfunc2 - vector-Hessian-vector product function for G_UP
643: . rhshp3 - an array of vectors storing the result of vector-Hessian-vector product for G_PU
644: . hessianproductfunc3 - vector-Hessian-vector product function for G_PU
645: . rhshp4 - an array of vectors storing the result of vector-Hessian-vector product for G_PP
646: - hessianproductfunc4 - vector-Hessian-vector product function for G_PP
648: Calling sequence of ihessianproductfunc:
649: $ rhshessianproductfunc (TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx);
650: + t - current timestep
651: . U - input vector (current ODE solution)
652: . Vl - an array of input vectors to be left-multiplied with the Hessian
653: . Vr - input vector to be right-multiplied with the Hessian
654: . VHV - an array of output vectors for vector-Hessian-vector product
655: - ctx - [optional] user-defined function context
657: Level: intermediate
659: Notes:
660: The first Hessian function and the working array are required.
661: As an example to implement the callback functions, the second callback function calculates the vector-Hessian-vector product
662: $ Vl_n^T*G_UP*Vr
663: 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.
664: Each entry of G_UP corresponds to the derivative
665: $ G_UP[i][j][k] = \frac{\partial^2 G[i]}{\partial U[j] \partial P[k]}.
666: The result of the vector-Hessian-vector product for Vl_n needs to be stored in vector VHV_n with j-th entry being
667: $ VHV_n[j] = \sum_i \sum_k {Vl_n[i] * G_UP[i][j][k] * Vr[k]}
668: If the cost function is a scalar, there will be only one vector in Vl and VHV.
670: .seealso:
671: @*/
672: PetscErrorCode TSSetRHSHessianProduct(TS ts,Vec *rhshp1,PetscErrorCode (*rhshessianproductfunc1)(TS,PetscReal,Vec,Vec*,Vec,Vec*,void*),
673: Vec *rhshp2,PetscErrorCode (*rhshessianproductfunc2)(TS,PetscReal,Vec,Vec*,Vec,Vec*,void*),
674: Vec *rhshp3,PetscErrorCode (*rhshessianproductfunc3)(TS,PetscReal,Vec,Vec*,Vec,Vec*,void*),
675: Vec *rhshp4,PetscErrorCode (*rhshessianproductfunc4)(TS,PetscReal,Vec,Vec*,Vec,Vec*,void*),
676: void *ctx)
677: {
682: ts->rhshessianproductctx = ctx;
683: if (rhshp1) ts->vecs_guu = rhshp1;
684: if (rhshp2) ts->vecs_gup = rhshp2;
685: if (rhshp3) ts->vecs_gpu = rhshp3;
686: if (rhshp4) ts->vecs_gpp = rhshp4;
687: ts->rhshessianproduct_guu = rhshessianproductfunc1;
688: ts->rhshessianproduct_gup = rhshessianproductfunc2;
689: ts->rhshessianproduct_gpu = rhshessianproductfunc3;
690: ts->rhshessianproduct_gpp = rhshessianproductfunc4;
691: return(0);
692: }
694: /*@C
695: TSComputeRHSHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Guu.
697: Collective on TS
699: Input Parameters:
700: . ts - The TS context obtained from TSCreate()
702: Notes:
703: TSComputeRHSHessianProductFunctionUU() is typically used for sensitivity implementation,
704: so most users would not generally call this routine themselves.
706: Level: developer
708: .seealso: TSSetRHSHessianProduct()
709: @*/
710: PetscErrorCode TSComputeRHSHessianProductFunctionUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
711: {
715: if (!VHV) return(0);
719: PetscStackPush("TS user RHSHessianProduct function 1 for sensitivity analysis");
720: (*ts->rhshessianproduct_guu)(ts,t,U,Vl,Vr,VHV,ts->rhshessianproductctx);
721: PetscStackPop;
722: return(0);
723: }
725: /*@C
726: TSComputeRHSHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Gup.
728: Collective on TS
730: Input Parameters:
731: . ts - The TS context obtained from TSCreate()
733: Notes:
734: TSComputeRHSHessianProductFunctionUP() is typically used for sensitivity implementation,
735: so most users would not generally call this routine themselves.
737: Level: developer
739: .seealso: TSSetRHSHessianProduct()
740: @*/
741: PetscErrorCode TSComputeRHSHessianProductFunctionUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
742: {
746: if (!VHV) return(0);
750: PetscStackPush("TS user RHSHessianProduct function 2 for sensitivity analysis");
751: (*ts->rhshessianproduct_gup)(ts,t,U,Vl,Vr,VHV,ts->rhshessianproductctx);
752: PetscStackPop;
753: return(0);
754: }
756: /*@C
757: TSComputeRHSHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Gpu.
759: Collective on TS
761: Input Parameters:
762: . ts - The TS context obtained from TSCreate()
764: Notes:
765: TSComputeRHSHessianProductFunctionPU() is typically used for sensitivity implementation,
766: so most users would not generally call this routine themselves.
768: Level: developer
770: .seealso: TSSetRHSHessianProduct()
771: @*/
772: PetscErrorCode TSComputeRHSHessianProductFunctionPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
773: {
777: if (!VHV) return(0);
781: PetscStackPush("TS user RHSHessianProduct function 3 for sensitivity analysis");
782: (*ts->rhshessianproduct_gpu)(ts,t,U,Vl,Vr,VHV,ts->rhshessianproductctx);
783: PetscStackPop;
784: return(0);
785: }
787: /*@C
788: TSComputeRHSHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Gpp.
790: Collective on TS
792: Input Parameters:
793: . ts - The TS context obtained from TSCreate()
795: Notes:
796: TSComputeRHSHessianProductFunctionPP() is typically used for sensitivity implementation,
797: so most users would not generally call this routine themselves.
799: Level: developer
801: .seealso: TSSetRHSHessianProduct()
802: @*/
803: PetscErrorCode TSComputeRHSHessianProductFunctionPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
804: {
808: if (!VHV) return(0);
812: PetscStackPush("TS user RHSHessianProduct function 3 for sensitivity analysis");
813: (*ts->rhshessianproduct_gpp)(ts,t,U,Vl,Vr,VHV,ts->rhshessianproductctx);
814: PetscStackPop;
815: return(0);
816: }
818: /* --------------------------- Adjoint sensitivity ---------------------------*/
820: /*@
821: TSSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial values and w.r.t. the problem parameters
822: for use by the TSAdjoint routines.
824: Logically Collective on TS
826: Input Parameters:
827: + ts - the TS context obtained from TSCreate()
828: . numcost - number of gradients to be computed, this is the number of cost functions
829: . 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
830: - mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
832: Level: beginner
834: Notes:
835: the entries in these vectors must be correctly initialized with the values lamda_i = df/dy|finaltime mu_i = df/dp|finaltime
837: After TSAdjointSolve() is called the lamba and the mu contain the computed sensitivities
839: .seealso TSGetCostGradients()
840: @*/
841: PetscErrorCode TSSetCostGradients(TS ts,PetscInt numcost,Vec *lambda,Vec *mu)
842: {
846: ts->vecs_sensi = lambda;
847: ts->vecs_sensip = mu;
848: if (ts->numcost && ts->numcost!=numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2nd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand");
849: ts->numcost = numcost;
850: return(0);
851: }
853: /*@
854: TSGetCostGradients - Returns the gradients from the TSAdjointSolve()
856: Not Collective, but Vec returned is parallel if TS is parallel
858: Input Parameter:
859: . ts - the TS context obtained from TSCreate()
861: Output Parameters:
862: + numcost - size of returned arrays
863: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
864: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
866: Level: intermediate
868: .seealso: TSSetCostGradients()
869: @*/
870: PetscErrorCode TSGetCostGradients(TS ts,PetscInt *numcost,Vec **lambda,Vec **mu)
871: {
874: if (numcost) *numcost = ts->numcost;
875: if (lambda) *lambda = ts->vecs_sensi;
876: if (mu) *mu = ts->vecs_sensip;
877: return(0);
878: }
880: /*@
881: 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
882: for use by the TSAdjoint routines.
884: Logically Collective on TS
886: Input Parameters:
887: + ts - the TS context obtained from TSCreate()
888: . numcost - number of cost functions
889: . 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
890: . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
891: - dir - the direction vector that are multiplied with the Hessian of the cost functions
893: Level: beginner
895: Notes: Hessian of the cost function is completely different from Hessian of the ODE/DAE system
897: For second-order adjoint, one needs to call this function and then TSAdjointSetForward() before TSSolve().
899: 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.
901: Passing NULL for lambda2 disables the second-order calculation.
902: .seealso: TSAdjointSetForward()
903: @*/
904: PetscErrorCode TSSetCostHessianProducts(TS ts,PetscInt numcost,Vec *lambda2,Vec *mu2,Vec dir)
905: {
908: if (ts->numcost && ts->numcost!=numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2nd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand");
909: ts->numcost = numcost;
910: ts->vecs_sensi2 = lambda2;
911: ts->vecs_sensi2p = mu2;
912: ts->vec_dir = dir;
913: return(0);
914: }
916: /*@
917: TSGetCostHessianProducts - Returns the gradients from the TSAdjointSolve()
919: Not Collective, but Vec returned is parallel if TS is parallel
921: Input Parameter:
922: . ts - the TS context obtained from TSCreate()
924: Output Parameters:
925: + numcost - number of cost functions
926: . lambda2 - Hessian-vector product with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector
927: . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
928: - dir - the direction vector that are multiplied with the Hessian of the cost functions
930: Level: intermediate
932: .seealso: TSSetCostHessianProducts()
933: @*/
934: PetscErrorCode TSGetCostHessianProducts(TS ts,PetscInt *numcost,Vec **lambda2,Vec **mu2, Vec *dir)
935: {
938: if (numcost) *numcost = ts->numcost;
939: if (lambda2) *lambda2 = ts->vecs_sensi2;
940: if (mu2) *mu2 = ts->vecs_sensi2p;
941: if (dir) *dir = ts->vec_dir;
942: return(0);
943: }
945: /*@
946: TSAdjointSetForward - Trigger the tangent linear solver and initialize the forward sensitivities
948: Logically Collective on TS
950: Input Parameters:
951: + ts - the TS context obtained from TSCreate()
952: - didp - the derivative of initial values w.r.t. parameters
954: Level: intermediate
956: 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.
958: .seealso: TSSetCostHessianProducts(), TSAdjointResetForward()
959: @*/
960: PetscErrorCode TSAdjointSetForward(TS ts,Mat didp)
961: {
962: Mat A;
963: Vec sp;
964: PetscScalar *xarr;
965: PetscInt lsize;
969: ts->forward_solve = PETSC_TRUE; /* turn on tangent linear mode */
970: if (!ts->vecs_sensi2) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetCostHessianProducts() first");
971: if (!ts->vec_dir) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Directional vector is missing. Call TSSetCostHessianProducts() to set it.");
972: /* create a single-column dense matrix */
973: VecGetLocalSize(ts->vec_sol,&lsize);
974: MatCreateDense(PetscObjectComm((PetscObject)ts),lsize,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&A);
976: VecDuplicate(ts->vec_sol,&sp);
977: MatDenseGetColumn(A,0,&xarr);
978: VecPlaceArray(sp,xarr);
979: if (ts->vecs_sensi2p) { /* tangent linear variable initialized as 2*dIdP*dir */
980: if (didp) {
981: MatMult(didp,ts->vec_dir,sp);
982: VecScale(sp,2.);
983: } else {
984: VecZeroEntries(sp);
985: }
986: } else { /* tangent linear variable initialized as dir */
987: VecCopy(ts->vec_dir,sp);
988: }
989: VecResetArray(sp);
990: MatDenseRestoreColumn(A,&xarr);
991: VecDestroy(&sp);
993: TSForwardSetInitialSensitivities(ts,A); /* if didp is NULL, identity matrix is assumed */
995: MatDestroy(&A);
996: return(0);
997: }
999: /*@
1000: TSAdjointResetForward - Reset the tangent linear solver and destroy the tangent linear context
1002: Logically Collective on TS
1004: Input Parameters:
1005: . ts - the TS context obtained from TSCreate()
1007: Level: intermediate
1009: .seealso: TSAdjointSetForward()
1010: @*/
1011: PetscErrorCode TSAdjointResetForward(TS ts)
1012: {
1016: ts->forward_solve = PETSC_FALSE; /* turn off tangent linear mode */
1017: TSForwardReset(ts);
1018: return(0);
1019: }
1021: /*@
1022: TSAdjointSetUp - Sets up the internal data structures for the later use
1023: of an adjoint solver
1025: Collective on TS
1027: Input Parameter:
1028: . ts - the TS context obtained from TSCreate()
1030: Level: advanced
1032: .seealso: TSCreate(), TSAdjointStep(), TSSetCostGradients()
1033: @*/
1034: PetscErrorCode TSAdjointSetUp(TS ts)
1035: {
1036: TSTrajectory tj;
1037: PetscBool match;
1038: PetscErrorCode ierr;
1042: if (ts->adjointsetupcalled) return(0);
1043: if (!ts->vecs_sensi) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetCostGradients() first");
1044: if (ts->vecs_sensip && !ts->Jacp && !ts->Jacprhs) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetRHSJacobianP() or TSSetIJacobianP() first");
1045: TSGetTrajectory(ts,&tj);
1046: PetscObjectTypeCompare((PetscObject)tj,TSTRAJECTORYBASIC,&match);
1047: if (match) {
1048: PetscBool solution_only;
1049: TSTrajectoryGetSolutionOnly(tj,&solution_only);
1050: if (solution_only) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"TSAdjoint cannot use the solution-only mode when choosing the Basic TSTrajectory type. Turn it off with -ts_trajectory_solution_only 0");
1051: }
1052: TSTrajectorySetUseHistory(tj,PETSC_FALSE); /* not use TSHistory */
1054: if (ts->quadraturets) { /* if there is integral in the cost function */
1055: VecDuplicate(ts->vecs_sensi[0],&ts->vec_drdu_col);
1056: if (ts->vecs_sensip) {
1057: VecDuplicate(ts->vecs_sensip[0],&ts->vec_drdp_col);
1058: }
1059: }
1061: if (ts->ops->adjointsetup) {
1062: (*ts->ops->adjointsetup)(ts);
1063: }
1064: ts->adjointsetupcalled = PETSC_TRUE;
1065: return(0);
1066: }
1068: /*@
1069: TSAdjointReset - Resets a TSAdjoint context and removes any allocated Vecs and Mats.
1071: Collective on TS
1073: Input Parameter:
1074: . ts - the TS context obtained from TSCreate()
1076: Level: beginner
1078: .seealso: TSCreate(), TSAdjointSetUp(), TSADestroy()
1079: @*/
1080: PetscErrorCode TSAdjointReset(TS ts)
1081: {
1086: if (ts->ops->adjointreset) {
1087: (*ts->ops->adjointreset)(ts);
1088: }
1089: if (ts->quadraturets) { /* if there is integral in the cost function */
1090: VecDestroy(&ts->vec_drdu_col);
1091: if (ts->vecs_sensip) {
1092: VecDestroy(&ts->vec_drdp_col);
1093: }
1094: }
1095: ts->vecs_sensi = NULL;
1096: ts->vecs_sensip = NULL;
1097: ts->vecs_sensi2 = NULL;
1098: ts->vecs_sensi2p = NULL;
1099: ts->vec_dir = NULL;
1100: ts->adjointsetupcalled = PETSC_FALSE;
1101: return(0);
1102: }
1104: /*@
1105: TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time
1107: Logically Collective on TS
1109: Input Parameters:
1110: + ts - the TS context obtained from TSCreate()
1111: - steps - number of steps to use
1113: Level: intermediate
1115: Notes:
1116: Normally one does not call this and TSAdjointSolve() integrates back to the original timestep. One can call this
1117: so as to integrate back to less than the original timestep
1119: .seealso: TSSetExactFinalTime()
1120: @*/
1121: PetscErrorCode TSAdjointSetSteps(TS ts,PetscInt steps)
1122: {
1126: if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back a negative number of steps");
1127: if (steps > ts->steps) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back more than the total number of forward steps");
1128: ts->adjoint_max_steps = steps;
1129: return(0);
1130: }
1132: /*@C
1133: TSAdjointSetRHSJacobian - Deprecated, use TSSetRHSJacobianP()
1135: Level: deprecated
1137: @*/
1138: PetscErrorCode TSAdjointSetRHSJacobian(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx)
1139: {
1146: ts->rhsjacobianp = func;
1147: ts->rhsjacobianpctx = ctx;
1148: if (Amat) {
1149: PetscObjectReference((PetscObject)Amat);
1150: MatDestroy(&ts->Jacp);
1151: ts->Jacp = Amat;
1152: }
1153: return(0);
1154: }
1156: /*@C
1157: TSAdjointComputeRHSJacobian - Deprecated, use TSComputeRHSJacobianP()
1159: Level: deprecated
1161: @*/
1162: PetscErrorCode TSAdjointComputeRHSJacobian(TS ts,PetscReal t,Vec U,Mat Amat)
1163: {
1171: PetscStackPush("TS user JacobianP function for sensitivity analysis");
1172: (*ts->rhsjacobianp)(ts,t,U,Amat,ts->rhsjacobianpctx);
1173: PetscStackPop;
1174: return(0);
1175: }
1177: /*@
1178: TSAdjointComputeDRDYFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobian()
1180: Level: deprecated
1182: @*/
1183: PetscErrorCode TSAdjointComputeDRDYFunction(TS ts,PetscReal t,Vec U,Vec *DRDU)
1184: {
1191: PetscStackPush("TS user DRDY function for sensitivity analysis");
1192: (*ts->drdufunction)(ts,t,U,DRDU,ts->costintegrandctx);
1193: PetscStackPop;
1194: return(0);
1195: }
1197: /*@
1198: TSAdjointComputeDRDPFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobianP()
1200: Level: deprecated
1202: @*/
1203: PetscErrorCode TSAdjointComputeDRDPFunction(TS ts,PetscReal t,Vec U,Vec *DRDP)
1204: {
1211: PetscStackPush("TS user DRDP function for sensitivity analysis");
1212: (*ts->drdpfunction)(ts,t,U,DRDP,ts->costintegrandctx);
1213: PetscStackPop;
1214: return(0);
1215: }
1217: /*@C
1218: TSAdjointMonitorSensi - monitors the first lambda sensitivity
1220: Level: intermediate
1222: .seealso: TSAdjointMonitorSet()
1223: @*/
1224: PetscErrorCode TSAdjointMonitorSensi(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
1225: {
1227: PetscViewer viewer = vf->viewer;
1231: PetscViewerPushFormat(viewer,vf->format);
1232: VecView(lambda[0],viewer);
1233: PetscViewerPopFormat(viewer);
1234: return(0);
1235: }
1237: /*@C
1238: TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
1240: Collective on TS
1242: Input Parameters:
1243: + ts - TS object you wish to monitor
1244: . name - the monitor type one is seeking
1245: . help - message indicating what monitoring is done
1246: . manual - manual page for the monitor
1247: . monitor - the monitor function
1248: - 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
1250: Level: developer
1252: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
1253: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
1254: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
1255: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
1256: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
1257: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
1258: PetscOptionsFList(), PetscOptionsEList()
1259: @*/
1260: 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*))
1261: {
1262: PetscErrorCode ierr;
1263: PetscViewer viewer;
1264: PetscViewerFormat format;
1265: PetscBool flg;
1268: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject) ts)->options,((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
1269: if (flg) {
1270: PetscViewerAndFormat *vf;
1271: PetscViewerAndFormatCreate(viewer,format,&vf);
1272: PetscObjectDereference((PetscObject)viewer);
1273: if (monitorsetup) {
1274: (*monitorsetup)(ts,vf);
1275: }
1276: TSAdjointMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
1277: }
1278: return(0);
1279: }
1281: /*@C
1282: TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
1283: timestep to display the iteration's progress.
1285: Logically Collective on TS
1287: Input Parameters:
1288: + ts - the TS context obtained from TSCreate()
1289: . adjointmonitor - monitoring routine
1290: . adjointmctx - [optional] user-defined context for private data for the
1291: monitor routine (use NULL if no context is desired)
1292: - adjointmonitordestroy - [optional] routine that frees monitor context
1293: (may be NULL)
1295: Calling sequence of monitor:
1296: $ int adjointmonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *adjointmctx)
1298: + ts - the TS context
1299: . 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
1300: been interpolated to)
1301: . time - current time
1302: . u - current iterate
1303: . numcost - number of cost functionos
1304: . lambda - sensitivities to initial conditions
1305: . mu - sensitivities to parameters
1306: - adjointmctx - [optional] adjoint monitoring context
1308: Notes:
1309: This routine adds an additional monitor to the list of monitors that
1310: already has been loaded.
1312: Fortran Notes:
1313: Only a single monitor function can be set for each TS object
1315: Level: intermediate
1317: .seealso: TSAdjointMonitorCancel()
1318: @*/
1319: PetscErrorCode TSAdjointMonitorSet(TS ts,PetscErrorCode (*adjointmonitor)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*),void *adjointmctx,PetscErrorCode (*adjointmdestroy)(void**))
1320: {
1322: PetscInt i;
1323: PetscBool identical;
1327: for (i=0; i<ts->numbermonitors;i++) {
1328: PetscMonitorCompare((PetscErrorCode (*)(void))adjointmonitor,adjointmctx,adjointmdestroy,(PetscErrorCode (*)(void))ts->adjointmonitor[i],ts->adjointmonitorcontext[i],ts->adjointmonitordestroy[i],&identical);
1329: if (identical) return(0);
1330: }
1331: if (ts->numberadjointmonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many adjoint monitors set");
1332: ts->adjointmonitor[ts->numberadjointmonitors] = adjointmonitor;
1333: ts->adjointmonitordestroy[ts->numberadjointmonitors] = adjointmdestroy;
1334: ts->adjointmonitorcontext[ts->numberadjointmonitors++] = (void*)adjointmctx;
1335: return(0);
1336: }
1338: /*@C
1339: TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.
1341: Logically Collective on TS
1343: Input Parameters:
1344: . ts - the TS context obtained from TSCreate()
1346: Notes:
1347: There is no way to remove a single, specific monitor.
1349: Level: intermediate
1351: .seealso: TSAdjointMonitorSet()
1352: @*/
1353: PetscErrorCode TSAdjointMonitorCancel(TS ts)
1354: {
1356: PetscInt i;
1360: for (i=0; i<ts->numberadjointmonitors; i++) {
1361: if (ts->adjointmonitordestroy[i]) {
1362: (*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]);
1363: }
1364: }
1365: ts->numberadjointmonitors = 0;
1366: return(0);
1367: }
1369: /*@C
1370: TSAdjointMonitorDefault - the default monitor of adjoint computations
1372: Level: intermediate
1374: .seealso: TSAdjointMonitorSet()
1375: @*/
1376: PetscErrorCode TSAdjointMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
1377: {
1379: PetscViewer viewer = vf->viewer;
1383: PetscViewerPushFormat(viewer,vf->format);
1384: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
1385: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
1386: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
1387: PetscViewerPopFormat(viewer);
1388: return(0);
1389: }
1391: /*@C
1392: TSAdjointMonitorDrawSensi - Monitors progress of the adjoint TS solvers by calling
1393: VecView() for the sensitivities to initial states at each timestep
1395: Collective on TS
1397: Input Parameters:
1398: + ts - the TS context
1399: . step - current time-step
1400: . ptime - current time
1401: . u - current state
1402: . numcost - number of cost functions
1403: . lambda - sensitivities to initial conditions
1404: . mu - sensitivities to parameters
1405: - dummy - either a viewer or NULL
1407: Level: intermediate
1409: .seealso: TSAdjointMonitorSet(), TSAdjointMonitorDefault(), VecView()
1410: @*/
1411: PetscErrorCode TSAdjointMonitorDrawSensi(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda,Vec *mu,void *dummy)
1412: {
1413: PetscErrorCode ierr;
1414: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
1415: PetscDraw draw;
1416: PetscReal xl,yl,xr,yr,h;
1417: char time[32];
1420: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);
1422: VecView(lambda[0],ictx->viewer);
1423: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
1424: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
1425: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1426: h = yl + .95*(yr - yl);
1427: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
1428: PetscDrawFlush(draw);
1429: return(0);
1430: }
1432: /*
1433: TSAdjointSetFromOptions - Sets various TSAdjoint parameters from user options.
1435: Collective on TSAdjoint
1437: Input Parameter:
1438: . ts - the TS context
1440: Options Database Keys:
1441: + -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
1442: . -ts_adjoint_monitor - print information at each adjoint time step
1443: - -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically
1445: Level: developer
1447: Notes:
1448: This is not normally called directly by users
1450: .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
1451: */
1452: PetscErrorCode TSAdjointSetFromOptions(PetscOptionItems *PetscOptionsObject,TS ts)
1453: {
1454: PetscBool tflg,opt;
1459: PetscOptionsHead(PetscOptionsObject,"TS Adjoint options");
1460: tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
1461: PetscOptionsBool("-ts_adjoint_solve","Solve the adjoint problem immediately after solving the forward problem","",tflg,&tflg,&opt);
1462: if (opt) {
1463: TSSetSaveTrajectory(ts);
1464: ts->adjoint_solve = tflg;
1465: }
1466: TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor","Monitor adjoint timestep size","TSAdjointMonitorDefault",TSAdjointMonitorDefault,NULL);
1467: TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor_sensi","Monitor sensitivity in the adjoint computation","TSAdjointMonitorSensi",TSAdjointMonitorSensi,NULL);
1468: opt = PETSC_FALSE;
1469: PetscOptionsName("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",&opt);
1470: if (opt) {
1471: TSMonitorDrawCtx ctx;
1472: PetscInt howoften = 1;
1474: PetscOptionsInt("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",howoften,&howoften,NULL);
1475: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
1476: TSAdjointMonitorSet(ts,TSAdjointMonitorDrawSensi,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
1477: }
1478: return(0);
1479: }
1481: /*@
1482: TSAdjointStep - Steps one time step backward in the adjoint run
1484: Collective on TS
1486: Input Parameter:
1487: . ts - the TS context obtained from TSCreate()
1489: Level: intermediate
1491: .seealso: TSAdjointSetUp(), TSAdjointSolve()
1492: @*/
1493: PetscErrorCode TSAdjointStep(TS ts)
1494: {
1495: DM dm;
1496: PetscErrorCode ierr;
1500: TSGetDM(ts,&dm);
1501: TSAdjointSetUp(ts);
1502: ts->steps--; /* must decrease the step index before the adjoint step is taken. */
1504: ts->reason = TS_CONVERGED_ITERATING;
1505: ts->ptime_prev = ts->ptime;
1506: if (!ts->ops->adjointstep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed because the adjoint of %s has not been implemented, try other time stepping methods for adjoint sensitivity analysis",((PetscObject)ts)->type_name);
1507: PetscLogEventBegin(TS_AdjointStep,ts,0,0,0);
1508: (*ts->ops->adjointstep)(ts);
1509: PetscLogEventEnd(TS_AdjointStep,ts,0,0,0);
1510: ts->adjoint_steps++;
1512: if (ts->reason < 0) {
1513: if (ts->errorifstepfailed) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSAdjointStep has failed due to %s",TSConvergedReasons[ts->reason]);
1514: } else if (!ts->reason) {
1515: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1516: }
1517: return(0);
1518: }
1520: /*@
1521: TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE
1523: Collective on TS
1525: Input Parameter:
1526: . ts - the TS context obtained from TSCreate()
1528: Options Database:
1529: . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial values
1531: Level: intermediate
1533: Notes:
1534: This must be called after a call to TSSolve() that solves the forward problem
1536: By default this will integrate back to the initial time, one can use TSAdjointSetSteps() to step back to a later time
1538: .seealso: TSCreate(), TSSetCostGradients(), TSSetSolution(), TSAdjointStep()
1539: @*/
1540: PetscErrorCode TSAdjointSolve(TS ts)
1541: {
1542: static PetscBool cite = PETSC_FALSE;
1543: #if defined(TSADJOINT_STAGE)
1544: PetscLogStage adjoint_stage;
1545: #endif
1550: PetscCitationsRegister("@article{tsadjointpaper,\n"
1551: " title = {{PETSc TSAdjoint: a discrete adjoint ODE solver for first-order and second-order sensitivity analysis}},\n"
1552: " author = {Zhang, Hong and Constantinescu, Emil M. and Smith, Barry F.},\n"
1553: " journal = {arXiv e-preprints},\n"
1554: " eprint = {1912.07696},\n"
1555: " archivePrefix = {arXiv},\n"
1556: " year = {2019}\n}\n",&cite);
1557: #if defined(TSADJOINT_STAGE)
1558: PetscLogStageRegister("TSAdjoint",&adjoint_stage);
1559: PetscLogStagePush(adjoint_stage);
1560: #endif
1561: TSAdjointSetUp(ts);
1563: /* reset time step and iteration counters */
1564: ts->adjoint_steps = 0;
1565: ts->ksp_its = 0;
1566: ts->snes_its = 0;
1567: ts->num_snes_failures = 0;
1568: ts->reject = 0;
1569: ts->reason = TS_CONVERGED_ITERATING;
1571: if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->steps;
1572: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1574: while (!ts->reason) {
1575: TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);
1576: TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
1577: TSAdjointEventHandler(ts);
1578: TSAdjointStep(ts);
1579: if ((ts->vec_costintegral || ts->quadraturets) && !ts->costintegralfwd) {
1580: TSAdjointCostIntegral(ts);
1581: }
1582: }
1583: if (!ts->steps) {
1584: TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);
1585: TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
1586: }
1587: ts->solvetime = ts->ptime;
1588: TSTrajectoryViewFromOptions(ts->trajectory,NULL,"-ts_trajectory_view");
1589: VecViewFromOptions(ts->vecs_sensi[0],(PetscObject) ts, "-ts_adjoint_view_solution");
1590: ts->adjoint_max_steps = 0;
1591: #if defined(TSADJOINT_STAGE)
1592: PetscLogStagePop();
1593: #endif
1594: return(0);
1595: }
1597: /*@C
1598: TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using TSAdjointMonitorSet()
1600: Collective on TS
1602: Input Parameters:
1603: + ts - time stepping context obtained from TSCreate()
1604: . step - step number that has just completed
1605: . ptime - model time of the state
1606: . u - state at the current model time
1607: . numcost - number of cost functions (dimension of lambda or mu)
1608: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
1609: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
1611: Notes:
1612: TSAdjointMonitor() is typically used automatically within the time stepping implementations.
1613: Users would almost never call this routine directly.
1615: Level: developer
1617: @*/
1618: PetscErrorCode TSAdjointMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda, Vec *mu)
1619: {
1621: PetscInt i,n = ts->numberadjointmonitors;
1626: VecLockReadPush(u);
1627: for (i=0; i<n; i++) {
1628: (*ts->adjointmonitor[i])(ts,step,ptime,u,numcost,lambda,mu,ts->adjointmonitorcontext[i]);
1629: }
1630: VecLockReadPop(u);
1631: return(0);
1632: }
1634: /*@
1635: TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.
1637: Collective on TS
1639: Input Parameter:
1640: . ts - time stepping context
1642: Level: advanced
1644: Notes:
1645: This function cannot be called until TSAdjointStep() has been completed.
1647: .seealso: TSAdjointSolve(), TSAdjointStep
1648: @*/
1649: PetscErrorCode TSAdjointCostIntegral(TS ts)
1650: {
1654: if (!ts->ops->adjointintegral) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide integral evaluation in the adjoint run",((PetscObject)ts)->type_name);
1655: (*ts->ops->adjointintegral)(ts);
1656: return(0);
1657: }
1659: /* ------------------ Forward (tangent linear) sensitivity ------------------*/
1661: /*@
1662: TSForwardSetUp - Sets up the internal data structures for the later use
1663: of forward sensitivity analysis
1665: Collective on TS
1667: Input Parameter:
1668: . ts - the TS context obtained from TSCreate()
1670: Level: advanced
1672: .seealso: TSCreate(), TSDestroy(), TSSetUp()
1673: @*/
1674: PetscErrorCode TSForwardSetUp(TS ts)
1675: {
1680: if (ts->forwardsetupcalled) return(0);
1681: if (ts->ops->forwardsetup) {
1682: (*ts->ops->forwardsetup)(ts);
1683: }
1684: VecDuplicate(ts->vec_sol,&ts->vec_sensip_col);
1685: ts->forwardsetupcalled = PETSC_TRUE;
1686: return(0);
1687: }
1689: /*@
1690: TSForwardReset - Reset the internal data structures used by forward sensitivity analysis
1692: Collective on TS
1694: Input Parameter:
1695: . ts - the TS context obtained from TSCreate()
1697: Level: advanced
1699: .seealso: TSCreate(), TSDestroy(), TSForwardSetUp()
1700: @*/
1701: PetscErrorCode TSForwardReset(TS ts)
1702: {
1703: TS quadts = ts->quadraturets;
1708: if (ts->ops->forwardreset) {
1709: (*ts->ops->forwardreset)(ts);
1710: }
1711: MatDestroy(&ts->mat_sensip);
1712: if (quadts) {
1713: MatDestroy(&quadts->mat_sensip);
1714: }
1715: VecDestroy(&ts->vec_sensip_col);
1716: ts->forward_solve = PETSC_FALSE;
1717: ts->forwardsetupcalled = PETSC_FALSE;
1718: return(0);
1719: }
1721: /*@
1722: TSForwardSetIntegralGradients - Set the vectors holding forward sensitivities of the integral term.
1724: Input Parameters:
1725: + ts- the TS context obtained from TSCreate()
1726: . numfwdint- number of integrals
1727: - vp = the vectors containing the gradients for each integral w.r.t. parameters
1729: Level: deprecated
1731: .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1732: @*/
1733: PetscErrorCode TSForwardSetIntegralGradients(TS ts,PetscInt numfwdint,Vec *vp)
1734: {
1737: if (ts->numcost && ts->numcost!=numfwdint) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2nd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand()");
1738: if (!ts->numcost) ts->numcost = numfwdint;
1740: ts->vecs_integral_sensip = vp;
1741: return(0);
1742: }
1744: /*@
1745: TSForwardGetIntegralGradients - Returns the forward sensitivities ofthe integral term.
1747: Input Parameter:
1748: . ts- the TS context obtained from TSCreate()
1750: Output Parameter:
1751: . vp = the vectors containing the gradients for each integral w.r.t. parameters
1753: Level: deprecated
1755: .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1756: @*/
1757: PetscErrorCode TSForwardGetIntegralGradients(TS ts,PetscInt *numfwdint,Vec **vp)
1758: {
1762: if (numfwdint) *numfwdint = ts->numcost;
1763: if (vp) *vp = ts->vecs_integral_sensip;
1764: return(0);
1765: }
1767: /*@
1768: TSForwardStep - Compute the forward sensitivity for one time step.
1770: Collective on TS
1772: Input Parameter:
1773: . ts - time stepping context
1775: Level: advanced
1777: Notes:
1778: This function cannot be called until TSStep() has been completed.
1780: .seealso: TSForwardSetSensitivities(), TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardSetUp()
1781: @*/
1782: PetscErrorCode TSForwardStep(TS ts)
1783: {
1787: if (!ts->ops->forwardstep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide forward sensitivity analysis",((PetscObject)ts)->type_name);
1788: PetscLogEventBegin(TS_ForwardStep,ts,0,0,0);
1789: (*ts->ops->forwardstep)(ts);
1790: PetscLogEventEnd(TS_ForwardStep,ts,0,0,0);
1791: if (ts->reason < 0 && ts->errorifstepfailed) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSFowardStep has failed due to %s",TSConvergedReasons[ts->reason]);
1792: return(0);
1793: }
1795: /*@
1796: TSForwardSetSensitivities - Sets the initial value of the trajectory sensitivities of solution w.r.t. the problem parameters and initial values.
1798: Logically Collective on TS
1800: Input Parameters:
1801: + ts - the TS context obtained from TSCreate()
1802: . nump - number of parameters
1803: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1805: Level: beginner
1807: Notes:
1808: Forward sensitivity is also called 'trajectory sensitivity' in some fields such as power systems.
1809: This function turns on a flag to trigger TSSolve() to compute forward sensitivities automatically.
1810: You must call this function before TSSolve().
1811: The entries in the sensitivity matrix must be correctly initialized with the values S = dy/dp|startingtime.
1813: .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1814: @*/
1815: PetscErrorCode TSForwardSetSensitivities(TS ts,PetscInt nump,Mat Smat)
1816: {
1822: ts->forward_solve = PETSC_TRUE;
1823: if (nump == PETSC_DEFAULT) {
1824: MatGetSize(Smat,NULL,&ts->num_parameters);
1825: } else ts->num_parameters = nump;
1826: PetscObjectReference((PetscObject)Smat);
1827: MatDestroy(&ts->mat_sensip);
1828: ts->mat_sensip = Smat;
1829: return(0);
1830: }
1832: /*@
1833: TSForwardGetSensitivities - Returns the trajectory sensitivities
1835: Not Collective, but Vec returned is parallel if TS is parallel
1837: Output Parameters:
1838: + ts - the TS context obtained from TSCreate()
1839: . nump - number of parameters
1840: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1842: Level: intermediate
1844: .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1845: @*/
1846: PetscErrorCode TSForwardGetSensitivities(TS ts,PetscInt *nump,Mat *Smat)
1847: {
1850: if (nump) *nump = ts->num_parameters;
1851: if (Smat) *Smat = ts->mat_sensip;
1852: return(0);
1853: }
1855: /*@
1856: TSForwardCostIntegral - Evaluate the cost integral in the forward run.
1858: Collective on TS
1860: Input Parameter:
1861: . ts - time stepping context
1863: Level: advanced
1865: Notes:
1866: This function cannot be called until TSStep() has been completed.
1868: .seealso: TSSolve(), TSAdjointCostIntegral()
1869: @*/
1870: PetscErrorCode TSForwardCostIntegral(TS ts)
1871: {
1876: if (!ts->ops->forwardintegral) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide integral evaluation in the forward run",((PetscObject)ts)->type_name);
1877: (*ts->ops->forwardintegral)(ts);
1878: return(0);
1879: }
1881: /*@
1882: TSForwardSetInitialSensitivities - Set initial values for tangent linear sensitivities
1884: Collective on TS
1886: Input Parameters:
1887: + ts - the TS context obtained from TSCreate()
1888: - didp - parametric sensitivities of the initial condition
1890: Level: intermediate
1892: 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.
1894: .seealso: TSForwardSetSensitivities()
1895: @*/
1896: PetscErrorCode TSForwardSetInitialSensitivities(TS ts,Mat didp)
1897: {
1903: if (!ts->mat_sensip) {
1904: TSForwardSetSensitivities(ts,PETSC_DEFAULT,didp);
1905: }
1906: return(0);
1907: }
1909: /*@
1910: TSForwardGetStages - Get the number of stages and the tangent linear sensitivities at the intermediate stages
1912: Input Parameter:
1913: . ts - the TS context obtained from TSCreate()
1915: Output Parameters:
1916: + ns - number of stages
1917: - S - tangent linear sensitivities at the intermediate stages
1919: Level: advanced
1921: @*/
1922: PetscErrorCode TSForwardGetStages(TS ts,PetscInt *ns,Mat **S)
1923: {
1929: if (!ts->ops->getstages) *S=NULL;
1930: else {
1931: (*ts->ops->forwardgetstages)(ts,ns,S);
1932: }
1933: return(0);
1934: }
1936: /*@
1937: TSCreateQuadratureTS - Create a sub-TS that evaluates integrals over time
1939: Input Parameters:
1940: + ts - the TS context obtained from TSCreate()
1941: - fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1943: Output Parameters:
1944: . quadts - the child TS context
1946: Level: intermediate
1948: .seealso: TSGetQuadratureTS()
1949: @*/
1950: PetscErrorCode TSCreateQuadratureTS(TS ts,PetscBool fwd,TS *quadts)
1951: {
1952: char prefix[128];
1958: TSDestroy(&ts->quadraturets);
1959: TSCreate(PetscObjectComm((PetscObject)ts),&ts->quadraturets);
1960: PetscObjectIncrementTabLevel((PetscObject)ts->quadraturets,(PetscObject)ts,1);
1961: PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->quadraturets);
1962: PetscSNPrintf(prefix,sizeof(prefix),"%squad_",((PetscObject)ts)->prefix ? ((PetscObject)ts)->prefix : "");
1963: TSSetOptionsPrefix(ts->quadraturets,prefix);
1964: *quadts = ts->quadraturets;
1966: if (ts->numcost) {
1967: VecCreateSeq(PETSC_COMM_SELF,ts->numcost,&(*quadts)->vec_sol);
1968: } else {
1969: VecCreateSeq(PETSC_COMM_SELF,1,&(*quadts)->vec_sol);
1970: }
1971: ts->costintegralfwd = fwd;
1972: return(0);
1973: }
1975: /*@
1976: TSGetQuadratureTS - Return the sub-TS that evaluates integrals over time
1978: Input Parameter:
1979: . ts - the TS context obtained from TSCreate()
1981: Output Parameters:
1982: + fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1983: - quadts - the child TS context
1985: Level: intermediate
1987: .seealso: TSCreateQuadratureTS()
1988: @*/
1989: PetscErrorCode TSGetQuadratureTS(TS ts,PetscBool *fwd,TS *quadts)
1990: {
1993: if (fwd) *fwd = ts->costintegralfwd;
1994: if (quadts) *quadts = ts->quadraturets;
1995: return(0);
1996: }
1998: /*@
1999: TSComputeSNESJacobian - Compute the SNESJacobian
2001: Input Parameters:
2002: + ts - the TS context obtained from TSCreate()
2003: - x - state vector
2005: Output Parameters:
2006: + J - Jacobian matrix
2007: - Jpre - preconditioning matrix for J (may be same as J)
2009: Level: developer
2011: Notes:
2012: Using SNES to compute the Jacobian enables finite differencing when TS Jacobian is not available.
2013: @*/
2014: PetscErrorCode TSComputeSNESJacobian(TS ts,Vec x,Mat J,Mat Jpre)
2015: {
2016: SNES snes = ts->snes;
2017: PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*) = NULL;
2021: /*
2022: Unlike implicit methods, explicit methods do not have SNESMatFDColoring in the snes object
2023: because SNESSolve() has not been called yet; so querying SNESMatFDColoring does not work for
2024: explicit methods. Instead, we check the Jacobian compute function directly to determin if FD
2025: coloring is used.
2026: */
2027: SNESGetJacobian(snes,NULL,NULL,&jac,NULL);
2028: if (jac == SNESComputeJacobianDefaultColor) {
2029: Vec f;
2030: SNESSetSolution(snes,x);
2031: SNESGetFunction(snes,&f,NULL,NULL);
2032: /* Force MatFDColoringApply to evaluate the SNES residual function for the base vector */
2033: SNESComputeFunction(snes,x,f);
2034: }
2035: SNESComputeJacobian(snes,x,J,Jpre);
2036: return(0);
2037: }