1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
2 #include <petscdraw.h>
3
4 PetscLogEvent TS_AdjointStep,TS_ForwardStep,TS_JacobianPEval;
5
6 /* #define TSADJOINT_STAGE */
7
8 /* ------------------------ Sensitivity Context ---------------------------*/
9
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.
12
13 Logically Collective on TS
14
15 Input Parameters:
16 + ts - TS context obtained from TSCreate()
17 . Amat - JacobianP matrix
18 . func - function
19 - ctx - [optional] user-defined function context
20
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
27
28 Level: intermediate
29
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
32
33 .seealso: TSGetRHSJacobianP()
34 @*/
TSSetRHSJacobianP(TS ts,Mat Amat,PetscErrorCode (* func)(TS,PetscReal,Vec,Mat,void *),void * ctx)35 PetscErrorCode TSSetRHSJacobianP(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx)
36 {
37 PetscErrorCode ierr;
38
39 PetscFunctionBegin;
40 PetscValidHeaderSpecific(ts, TS_CLASSID,1);
41 PetscValidHeaderSpecific(Amat,MAT_CLASSID,2);
42
43 ts->rhsjacobianp = func;
44 ts->rhsjacobianpctx = ctx;
45 if (Amat) {
46 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr);
47 ierr = MatDestroy(&ts->Jacprhs);CHKERRQ(ierr);
48 ts->Jacprhs = Amat;
49 }
50 PetscFunctionReturn(0);
51 }
52
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.
55
56 Logically Collective on TS
57
58 Input Parameters:
59 . ts - TS context obtained from TSCreate()
60
61 Output Parameters:
62 + Amat - JacobianP matrix
63 . func - function
64 - ctx - [optional] user-defined function context
65
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
72
73 Level: intermediate
74
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
77
78 .seealso: TSSetRHSJacobianP()
79 @*/
TSGetRHSJacobianP(TS ts,Mat * Amat,PetscErrorCode (** func)(TS,PetscReal,Vec,Mat,void *),void ** ctx)80 PetscErrorCode TSGetRHSJacobianP(TS ts,Mat *Amat,PetscErrorCode (**func)(TS,PetscReal,Vec,Mat,void*),void **ctx)
81 {
82 PetscFunctionBegin;
83 if (func) *func = ts->rhsjacobianp;
84 if (ctx) *ctx = ts->rhsjacobianpctx;
85 if (Amat) *Amat = ts->Jacprhs;
86 PetscFunctionReturn(0);
87 }
88
89 /*@C
90 TSComputeRHSJacobianP - Runs the user-defined JacobianP function.
91
92 Collective on TS
93
94 Input Parameters:
95 . ts - The TS context obtained from TSCreate()
96
97 Level: developer
98
99 .seealso: TSSetRHSJacobianP()
100 @*/
TSComputeRHSJacobianP(TS ts,PetscReal t,Vec U,Mat Amat)101 PetscErrorCode TSComputeRHSJacobianP(TS ts,PetscReal t,Vec U,Mat Amat)
102 {
103 PetscErrorCode ierr;
104
105 PetscFunctionBegin;
106 if (!Amat) PetscFunctionReturn(0);
107 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
108 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
109
110 PetscStackPush("TS user JacobianP function for sensitivity analysis");
111 ierr = (*ts->rhsjacobianp)(ts,t,U,Amat,ts->rhsjacobianpctx);CHKERRQ(ierr);
112 PetscStackPop;
113 PetscFunctionReturn(0);
114 }
115
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.
118
119 Logically Collective on TS
120
121 Input Parameters:
122 + ts - TS context obtained from TSCreate()
123 . Amat - JacobianP matrix
124 . func - function
125 - ctx - [optional] user-defined function context
126
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
135
136 Level: intermediate
137
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
140
141 .seealso:
142 @*/
TSSetIJacobianP(TS ts,Mat Amat,PetscErrorCode (* func)(TS,PetscReal,Vec,Vec,PetscReal,Mat,void *),void * ctx)143 PetscErrorCode TSSetIJacobianP(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Vec,PetscReal,Mat,void*),void *ctx)
144 {
145 PetscErrorCode ierr;
146
147 PetscFunctionBegin;
148 PetscValidHeaderSpecific(ts, TS_CLASSID,1);
149 PetscValidHeaderSpecific(Amat,MAT_CLASSID,2);
150
151 ts->ijacobianp = func;
152 ts->ijacobianpctx = ctx;
153 if (Amat) {
154 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr);
155 ierr = MatDestroy(&ts->Jacp);CHKERRQ(ierr);
156 ts->Jacp = Amat;
157 }
158 PetscFunctionReturn(0);
159 }
160
161 /*@C
162 TSComputeIJacobianP - Runs the user-defined IJacobianP function.
163
164 Collective on TS
165
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
173
174 Output Parameters:
175 . A - Jacobian matrix
176
177 Level: developer
178
179 .seealso: TSSetIJacobianP()
180 @*/
TSComputeIJacobianP(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat Amat,PetscBool imex)181 PetscErrorCode TSComputeIJacobianP(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat Amat,PetscBool imex)
182 {
183 PetscErrorCode ierr;
184
185 PetscFunctionBegin;
186 if (!Amat) PetscFunctionReturn(0);
187 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
188 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
189 PetscValidHeaderSpecific(Udot,VEC_CLASSID,4);
190
191 ierr = PetscLogEventBegin(TS_JacobianPEval,ts,U,Amat,0);CHKERRQ(ierr);
192 if (ts->ijacobianp) {
193 PetscStackPush("TS user JacobianP function for sensitivity analysis");
194 ierr = (*ts->ijacobianp)(ts,t,U,Udot,shift,Amat,ts->ijacobianpctx);CHKERRQ(ierr);
195 PetscStackPop;
196 }
197 if (imex) {
198 if (!ts->ijacobianp) { /* system was written as Udot = G(t,U) */
199 PetscBool assembled;
200 ierr = MatZeroEntries(Amat);CHKERRQ(ierr);
201 ierr = MatAssembled(Amat,&assembled);CHKERRQ(ierr);
202 if (!assembled) {
203 ierr = MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
204 ierr = MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
205 }
206 }
207 } else {
208 if (ts->rhsjacobianp) {
209 ierr = TSComputeRHSJacobianP(ts,t,U,ts->Jacprhs);CHKERRQ(ierr);
210 }
211 if (ts->Jacprhs == Amat) { /* No IJacobian, so we only have the RHS matrix */
212 ierr = MatScale(Amat,-1);CHKERRQ(ierr);
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 ierr = MatZeroEntries(Amat);CHKERRQ(ierr);
217 }
218 ierr = MatAXPY(Amat,-1,ts->Jacprhs,axpy);CHKERRQ(ierr);
219 }
220 }
221 ierr = PetscLogEventEnd(TS_JacobianPEval,ts,U,Amat,0);CHKERRQ(ierr);
222 PetscFunctionReturn(0);
223 }
224
225 /*@C
226 TSSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions
227
228 Logically Collective on TS
229
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)
239
240 Calling sequence of rf:
241 $ PetscErrorCode rf(TS ts,PetscReal t,Vec U,Vec F,void *ctx);
242
243 Calling sequence of drduf:
244 $ PetscErroCode drduf(TS ts,PetscReal t,Vec U,Vec *dRdU,void *ctx);
245
246 Calling sequence of drdpf:
247 $ PetscErroCode drdpf(TS ts,PetscReal t,Vec U,Vec *dRdP,void *ctx);
248
249 Level: deprecated
250
251 Notes:
252 For optimization there is usually a single cost function (numcost = 1). For sensitivities there may be multiple cost functions
253
254 .seealso: TSSetRHSJacobianP(), TSGetCostGradients(), TSSetCostGradients()
255 @*/
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)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 {
261 PetscErrorCode ierr;
262
263 PetscFunctionBegin;
264 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
265 if (costintegral) PetscValidHeaderSpecific(costintegral,VEC_CLASSID,3);
266 if (ts->numcost && ts->numcost!=numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2rd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostGradients() or TSForwardSetIntegralGradients()");
267 if (!ts->numcost) ts->numcost=numcost;
268
269 if (costintegral) {
270 ierr = PetscObjectReference((PetscObject)costintegral);CHKERRQ(ierr);
271 ierr = VecDestroy(&ts->vec_costintegral);CHKERRQ(ierr);
272 ts->vec_costintegral = costintegral;
273 } else {
274 if (!ts->vec_costintegral) { /* Create a seq vec if user does not provide one */
275 ierr = VecCreateSeq(PETSC_COMM_SELF,numcost,&ts->vec_costintegral);CHKERRQ(ierr);
276 } else {
277 ierr = VecSet(ts->vec_costintegral,0.0);CHKERRQ(ierr);
278 }
279 }
280 if (!ts->vec_costintegrand) {
281 ierr = VecDuplicate(ts->vec_costintegral,&ts->vec_costintegrand);CHKERRQ(ierr);
282 } else {
283 ierr = VecSet(ts->vec_costintegrand,0.0);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
291 }
292
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.
296
297 Not Collective
298
299 Input Parameter:
300 . ts - the TS context obtained from TSCreate()
301
302 Output Parameter:
303 . v - the vector containing the integrals for each cost function
304
305 Level: intermediate
306
307 .seealso: TSSetCostIntegrand()
308
309 @*/
TSGetCostIntegral(TS ts,Vec * v)310 PetscErrorCode TSGetCostIntegral(TS ts,Vec *v)
311 {
312 TS quadts;
313 PetscErrorCode ierr;
314
315 PetscFunctionBegin;
316 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
317 PetscValidPointer(v,2);
318 ierr = TSGetQuadratureTS(ts,NULL,&quadts);CHKERRQ(ierr);
319 *v = quadts->vec_sol;
320 PetscFunctionReturn(0);
321 }
322
323 /*@C
324 TSComputeCostIntegrand - Evaluates the integral function in the cost functions.
325
326 Input Parameters:
327 + ts - the TS context
328 . t - current time
329 - U - state vector, i.e. current solution
330
331 Output Parameter:
332 . Q - vector of size numcost to hold the outputs
333
334 Notes:
335 Most users should not need to explicitly call this routine, as it
336 is used internally within the sensitivity analysis context.
337
338 Level: deprecated
339
340 .seealso: TSSetCostIntegrand()
341 @*/
TSComputeCostIntegrand(TS ts,PetscReal t,Vec U,Vec Q)342 PetscErrorCode TSComputeCostIntegrand(TS ts,PetscReal t,Vec U,Vec Q)
343 {
344 PetscErrorCode ierr;
345
346 PetscFunctionBegin;
347 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
348 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
349 PetscValidHeaderSpecific(Q,VEC_CLASSID,4);
350
351 ierr = PetscLogEventBegin(TS_FunctionEval,ts,U,Q,0);CHKERRQ(ierr);
352 if (ts->costintegrand) {
353 PetscStackPush("TS user integrand in the cost function");
354 ierr = (*ts->costintegrand)(ts,t,U,Q,ts->costintegrandctx);CHKERRQ(ierr);
355 PetscStackPop;
356 } else {
357 ierr = VecZeroEntries(Q);CHKERRQ(ierr);
358 }
359
360 ierr = PetscLogEventEnd(TS_FunctionEval,ts,U,Q,0);CHKERRQ(ierr);
361 PetscFunctionReturn(0);
362 }
363
364 /*@C
365 TSComputeDRDUFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobian()
366
367 Level: deprecated
368
369 @*/
TSComputeDRDUFunction(TS ts,PetscReal t,Vec U,Vec * DRDU)370 PetscErrorCode TSComputeDRDUFunction(TS ts,PetscReal t,Vec U,Vec *DRDU)
371 {
372 PetscErrorCode ierr;
373
374 PetscFunctionBegin;
375 if (!DRDU) PetscFunctionReturn(0);
376 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
377 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
378
379 PetscStackPush("TS user DRDU function for sensitivity analysis");
380 ierr = (*ts->drdufunction)(ts,t,U,DRDU,ts->costintegrandctx);CHKERRQ(ierr);
381 PetscStackPop;
382 PetscFunctionReturn(0);
383 }
384
385 /*@C
386 TSComputeDRDPFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobianP()
387
388 Level: deprecated
389
390 @*/
TSComputeDRDPFunction(TS ts,PetscReal t,Vec U,Vec * DRDP)391 PetscErrorCode TSComputeDRDPFunction(TS ts,PetscReal t,Vec U,Vec *DRDP)
392 {
393 PetscErrorCode ierr;
394
395 PetscFunctionBegin;
396 if (!DRDP) PetscFunctionReturn(0);
397 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
398 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
399
400 PetscStackPush("TS user DRDP function for sensitivity analysis");
401 ierr = (*ts->drdpfunction)(ts,t,U,DRDP,ts->costintegrandctx);CHKERRQ(ierr);
402 PetscStackPop;
403 PetscFunctionReturn(0);
404 }
405
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.
408
409 Logically Collective on TS
410
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
421
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
430
431 Level: intermediate
432
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.
443
444 .seealso:
445 @*/
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)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 {
452 PetscFunctionBegin;
453 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
454 PetscValidPointer(ihp1,2);
455
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 PetscFunctionReturn(0);
466 }
467
468 /*@C
469 TSComputeIHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Fuu.
470
471 Collective on TS
472
473 Input Parameters:
474 . ts - The TS context obtained from TSCreate()
475
476 Notes:
477 TSComputeIHessianProductFunctionUU() is typically used for sensitivity implementation,
478 so most users would not generally call this routine themselves.
479
480 Level: developer
481
482 .seealso: TSSetIHessianProduct()
483 @*/
TSComputeIHessianProductFunctionUU(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV)484 PetscErrorCode TSComputeIHessianProductFunctionUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
485 {
486 PetscErrorCode ierr;
487
488 PetscFunctionBegin;
489 if (!VHV) PetscFunctionReturn(0);
490 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
491 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
492
493 if (ts->ihessianproduct_fuu) {
494 PetscStackPush("TS user IHessianProduct function 1 for sensitivity analysis");
495 ierr = (*ts->ihessianproduct_fuu)(ts,t,U,Vl,Vr,VHV,ts->ihessianproductctx);CHKERRQ(ierr);
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 ierr = TSComputeRHSHessianProductFunctionUU(ts,t,U,Vl,Vr,VHV);CHKERRQ(ierr);
502 for (nadj=0; nadj<ts->numcost; nadj++) {
503 ierr = VecScale(VHV[nadj],-1);CHKERRQ(ierr);
504 }
505 }
506 PetscFunctionReturn(0);
507 }
508
509 /*@C
510 TSComputeIHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Fup.
511
512 Collective on TS
513
514 Input Parameters:
515 . ts - The TS context obtained from TSCreate()
516
517 Notes:
518 TSComputeIHessianProductFunctionUP() is typically used for sensitivity implementation,
519 so most users would not generally call this routine themselves.
520
521 Level: developer
522
523 .seealso: TSSetIHessianProduct()
524 @*/
TSComputeIHessianProductFunctionUP(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV)525 PetscErrorCode TSComputeIHessianProductFunctionUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
526 {
527 PetscErrorCode ierr;
528
529 PetscFunctionBegin;
530 if (!VHV) PetscFunctionReturn(0);
531 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
532 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
533
534 if (ts->ihessianproduct_fup) {
535 PetscStackPush("TS user IHessianProduct function 2 for sensitivity analysis");
536 ierr = (*ts->ihessianproduct_fup)(ts,t,U,Vl,Vr,VHV,ts->ihessianproductctx);CHKERRQ(ierr);
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 ierr = TSComputeRHSHessianProductFunctionUP(ts,t,U,Vl,Vr,VHV);CHKERRQ(ierr);
543 for (nadj=0; nadj<ts->numcost; nadj++) {
544 ierr = VecScale(VHV[nadj],-1);CHKERRQ(ierr);
545 }
546 }
547 PetscFunctionReturn(0);
548 }
549
550 /*@C
551 TSComputeIHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Fpu.
552
553 Collective on TS
554
555 Input Parameters:
556 . ts - The TS context obtained from TSCreate()
557
558 Notes:
559 TSComputeIHessianProductFunctionPU() is typically used for sensitivity implementation,
560 so most users would not generally call this routine themselves.
561
562 Level: developer
563
564 .seealso: TSSetIHessianProduct()
565 @*/
TSComputeIHessianProductFunctionPU(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV)566 PetscErrorCode TSComputeIHessianProductFunctionPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
567 {
568 PetscErrorCode ierr;
569
570 PetscFunctionBegin;
571 if (!VHV) PetscFunctionReturn(0);
572 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
573 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
574
575 if (ts->ihessianproduct_fpu) {
576 PetscStackPush("TS user IHessianProduct function 3 for sensitivity analysis");
577 ierr = (*ts->ihessianproduct_fpu)(ts,t,U,Vl,Vr,VHV,ts->ihessianproductctx);CHKERRQ(ierr);
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 ierr = TSComputeRHSHessianProductFunctionPU(ts,t,U,Vl,Vr,VHV);CHKERRQ(ierr);
584 for (nadj=0; nadj<ts->numcost; nadj++) {
585 ierr = VecScale(VHV[nadj],-1);CHKERRQ(ierr);
586 }
587 }
588 PetscFunctionReturn(0);
589 }
590
591 /*@C
592 TSComputeIHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Fpp.
593
594 Collective on TS
595
596 Input Parameters:
597 . ts - The TS context obtained from TSCreate()
598
599 Notes:
600 TSComputeIHessianProductFunctionPP() is typically used for sensitivity implementation,
601 so most users would not generally call this routine themselves.
602
603 Level: developer
604
605 .seealso: TSSetIHessianProduct()
606 @*/
TSComputeIHessianProductFunctionPP(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV)607 PetscErrorCode TSComputeIHessianProductFunctionPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
608 {
609 PetscErrorCode ierr;
610
611 PetscFunctionBegin;
612 if (!VHV) PetscFunctionReturn(0);
613 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
614 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
615
616 if (ts->ihessianproduct_fpp) {
617 PetscStackPush("TS user IHessianProduct function 3 for sensitivity analysis");
618 ierr = (*ts->ihessianproduct_fpp)(ts,t,U,Vl,Vr,VHV,ts->ihessianproductctx);CHKERRQ(ierr);
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 ierr = TSComputeRHSHessianProductFunctionPP(ts,t,U,Vl,Vr,VHV);CHKERRQ(ierr);
625 for (nadj=0; nadj<ts->numcost; nadj++) {
626 ierr = VecScale(VHV[nadj],-1);CHKERRQ(ierr);
627 }
628 }
629 PetscFunctionReturn(0);
630 }
631
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.
634
635 Logically Collective on TS
636
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
647
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
656
657 Level: intermediate
658
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.
669
670 .seealso:
671 @*/
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)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 {
678 PetscFunctionBegin;
679 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
680 PetscValidPointer(rhshp1,2);
681
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 PetscFunctionReturn(0);
692 }
693
694 /*@C
695 TSComputeRHSHessianProductFunctionUU - Runs the user-defined vector-Hessian-vector product function for Guu.
696
697 Collective on TS
698
699 Input Parameters:
700 . ts - The TS context obtained from TSCreate()
701
702 Notes:
703 TSComputeRHSHessianProductFunctionUU() is typically used for sensitivity implementation,
704 so most users would not generally call this routine themselves.
705
706 Level: developer
707
708 .seealso: TSSetRHSHessianProduct()
709 @*/
TSComputeRHSHessianProductFunctionUU(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV)710 PetscErrorCode TSComputeRHSHessianProductFunctionUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
711 {
712 PetscErrorCode ierr;
713
714 PetscFunctionBegin;
715 if (!VHV) PetscFunctionReturn(0);
716 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
717 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
718
719 PetscStackPush("TS user RHSHessianProduct function 1 for sensitivity analysis");
720 ierr = (*ts->rhshessianproduct_guu)(ts,t,U,Vl,Vr,VHV,ts->rhshessianproductctx);CHKERRQ(ierr);
721 PetscStackPop;
722 PetscFunctionReturn(0);
723 }
724
725 /*@C
726 TSComputeRHSHessianProductFunctionUP - Runs the user-defined vector-Hessian-vector product function for Gup.
727
728 Collective on TS
729
730 Input Parameters:
731 . ts - The TS context obtained from TSCreate()
732
733 Notes:
734 TSComputeRHSHessianProductFunctionUP() is typically used for sensitivity implementation,
735 so most users would not generally call this routine themselves.
736
737 Level: developer
738
739 .seealso: TSSetRHSHessianProduct()
740 @*/
TSComputeRHSHessianProductFunctionUP(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV)741 PetscErrorCode TSComputeRHSHessianProductFunctionUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
742 {
743 PetscErrorCode ierr;
744
745 PetscFunctionBegin;
746 if (!VHV) PetscFunctionReturn(0);
747 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
748 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
749
750 PetscStackPush("TS user RHSHessianProduct function 2 for sensitivity analysis");
751 ierr = (*ts->rhshessianproduct_gup)(ts,t,U,Vl,Vr,VHV,ts->rhshessianproductctx);CHKERRQ(ierr);
752 PetscStackPop;
753 PetscFunctionReturn(0);
754 }
755
756 /*@C
757 TSComputeRHSHessianProductFunctionPU - Runs the user-defined vector-Hessian-vector product function for Gpu.
758
759 Collective on TS
760
761 Input Parameters:
762 . ts - The TS context obtained from TSCreate()
763
764 Notes:
765 TSComputeRHSHessianProductFunctionPU() is typically used for sensitivity implementation,
766 so most users would not generally call this routine themselves.
767
768 Level: developer
769
770 .seealso: TSSetRHSHessianProduct()
771 @*/
TSComputeRHSHessianProductFunctionPU(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV)772 PetscErrorCode TSComputeRHSHessianProductFunctionPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
773 {
774 PetscErrorCode ierr;
775
776 PetscFunctionBegin;
777 if (!VHV) PetscFunctionReturn(0);
778 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
779 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
780
781 PetscStackPush("TS user RHSHessianProduct function 3 for sensitivity analysis");
782 ierr = (*ts->rhshessianproduct_gpu)(ts,t,U,Vl,Vr,VHV,ts->rhshessianproductctx);CHKERRQ(ierr);
783 PetscStackPop;
784 PetscFunctionReturn(0);
785 }
786
787 /*@C
788 TSComputeRHSHessianProductFunctionPP - Runs the user-defined vector-Hessian-vector product function for Gpp.
789
790 Collective on TS
791
792 Input Parameters:
793 . ts - The TS context obtained from TSCreate()
794
795 Notes:
796 TSComputeRHSHessianProductFunctionPP() is typically used for sensitivity implementation,
797 so most users would not generally call this routine themselves.
798
799 Level: developer
800
801 .seealso: TSSetRHSHessianProduct()
802 @*/
TSComputeRHSHessianProductFunctionPP(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV)803 PetscErrorCode TSComputeRHSHessianProductFunctionPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV)
804 {
805 PetscErrorCode ierr;
806
807 PetscFunctionBegin;
808 if (!VHV) PetscFunctionReturn(0);
809 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
810 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
811
812 PetscStackPush("TS user RHSHessianProduct function 3 for sensitivity analysis");
813 ierr = (*ts->rhshessianproduct_gpp)(ts,t,U,Vl,Vr,VHV,ts->rhshessianproductctx);CHKERRQ(ierr);
814 PetscStackPop;
815 PetscFunctionReturn(0);
816 }
817
818 /* --------------------------- Adjoint sensitivity ---------------------------*/
819
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.
823
824 Logically Collective on TS
825
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
831
832 Level: beginner
833
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
836
837 After TSAdjointSolve() is called the lamba and the mu contain the computed sensitivities
838
839 .seealso TSGetCostGradients()
840 @*/
TSSetCostGradients(TS ts,PetscInt numcost,Vec * lambda,Vec * mu)841 PetscErrorCode TSSetCostGradients(TS ts,PetscInt numcost,Vec *lambda,Vec *mu)
842 {
843 PetscFunctionBegin;
844 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
845 PetscValidPointer(lambda,2);
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 (2rd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand");
849 ts->numcost = numcost;
850 PetscFunctionReturn(0);
851 }
852
853 /*@
854 TSGetCostGradients - Returns the gradients from the TSAdjointSolve()
855
856 Not Collective, but Vec returned is parallel if TS is parallel
857
858 Input Parameter:
859 . ts - the TS context obtained from TSCreate()
860
861 Output Parameter:
862 + lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
863 - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
864
865 Level: intermediate
866
867 .seealso: TSSetCostGradients()
868 @*/
TSGetCostGradients(TS ts,PetscInt * numcost,Vec ** lambda,Vec ** mu)869 PetscErrorCode TSGetCostGradients(TS ts,PetscInt *numcost,Vec **lambda,Vec **mu)
870 {
871 PetscFunctionBegin;
872 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
873 if (numcost) *numcost = ts->numcost;
874 if (lambda) *lambda = ts->vecs_sensi;
875 if (mu) *mu = ts->vecs_sensip;
876 PetscFunctionReturn(0);
877 }
878
879 /*@
880 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
881 for use by the TSAdjoint routines.
882
883 Logically Collective on TS
884
885 Input Parameters:
886 + ts - the TS context obtained from TSCreate()
887 . numcost - number of cost functions
888 . 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
889 . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
890 - dir - the direction vector that are multiplied with the Hessian of the cost functions
891
892 Level: beginner
893
894 Notes: Hessian of the cost function is completely different from Hessian of the ODE/DAE system
895
896 For second-order adjoint, one needs to call this function and then TSAdjointSetForward() before TSSolve().
897
898 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.
899
900 Passing NULL for lambda2 disables the second-order calculation.
901 .seealso: TSAdjointSetForward()
902 @*/
TSSetCostHessianProducts(TS ts,PetscInt numcost,Vec * lambda2,Vec * mu2,Vec dir)903 PetscErrorCode TSSetCostHessianProducts(TS ts,PetscInt numcost,Vec *lambda2,Vec *mu2,Vec dir)
904 {
905 PetscFunctionBegin;
906 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
907 if (ts->numcost && ts->numcost!=numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2rd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand");
908 ts->numcost = numcost;
909 ts->vecs_sensi2 = lambda2;
910 ts->vecs_sensi2p = mu2;
911 ts->vec_dir = dir;
912 PetscFunctionReturn(0);
913 }
914
915 /*@
916 TSGetCostHessianProducts - Returns the gradients from the TSAdjointSolve()
917
918 Not Collective, but Vec returned is parallel if TS is parallel
919
920 Input Parameter:
921 . ts - the TS context obtained from TSCreate()
922
923 Output Parameter:
924 + numcost - number of cost functions
925 . 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
926 . mu2 - Hessian-vector product with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
927 - dir - the direction vector that are multiplied with the Hessian of the cost functions
928
929 Level: intermediate
930
931 .seealso: TSSetCostHessianProducts()
932 @*/
TSGetCostHessianProducts(TS ts,PetscInt * numcost,Vec ** lambda2,Vec ** mu2,Vec * dir)933 PetscErrorCode TSGetCostHessianProducts(TS ts,PetscInt *numcost,Vec **lambda2,Vec **mu2, Vec *dir)
934 {
935 PetscFunctionBegin;
936 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
937 if (numcost) *numcost = ts->numcost;
938 if (lambda2) *lambda2 = ts->vecs_sensi2;
939 if (mu2) *mu2 = ts->vecs_sensi2p;
940 if (dir) *dir = ts->vec_dir;
941 PetscFunctionReturn(0);
942 }
943
944 /*@
945 TSAdjointSetForward - Trigger the tangent linear solver and initialize the forward sensitivities
946
947 Logically Collective on TS
948
949 Input Parameters:
950 + ts - the TS context obtained from TSCreate()
951 - didp - the derivative of initial values w.r.t. parameters
952
953 Level: intermediate
954
955 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.
956
957 .seealso: TSSetCostHessianProducts(), TSAdjointResetForward()
958 @*/
TSAdjointSetForward(TS ts,Mat didp)959 PetscErrorCode TSAdjointSetForward(TS ts,Mat didp)
960 {
961 Mat A;
962 Vec sp;
963 PetscScalar *xarr;
964 PetscInt lsize;
965 PetscErrorCode ierr;
966
967 PetscFunctionBegin;
968 ts->forward_solve = PETSC_TRUE; /* turn on tangent linear mode */
969 if (!ts->vecs_sensi2) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetCostHessianProducts() first");
970 if (!ts->vec_dir) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Directional vector is missing. Call TSSetCostHessianProducts() to set it.");
971 /* create a single-column dense matrix */
972 ierr = VecGetLocalSize(ts->vec_sol,&lsize);CHKERRQ(ierr);
973 ierr = MatCreateDense(PetscObjectComm((PetscObject)ts),lsize,PETSC_DECIDE,PETSC_DECIDE,1,NULL,&A);CHKERRQ(ierr);
974
975 ierr = VecDuplicate(ts->vec_sol,&sp);CHKERRQ(ierr);
976 ierr = MatDenseGetColumn(A,0,&xarr);CHKERRQ(ierr);
977 ierr = VecPlaceArray(sp,xarr);CHKERRQ(ierr);
978 if (ts->vecs_sensi2p) { /* tangent linear variable initialized as 2*dIdP*dir */
979 if (didp) {
980 ierr = MatMult(didp,ts->vec_dir,sp);CHKERRQ(ierr);
981 ierr = VecScale(sp,2.);CHKERRQ(ierr);
982 } else {
983 ierr = VecZeroEntries(sp);CHKERRQ(ierr);
984 }
985 } else { /* tangent linear variable initialized as dir */
986 ierr = VecCopy(ts->vec_dir,sp);CHKERRQ(ierr);
987 }
988 ierr = VecResetArray(sp);CHKERRQ(ierr);
989 ierr = MatDenseRestoreColumn(A,&xarr);CHKERRQ(ierr);
990 ierr = VecDestroy(&sp);CHKERRQ(ierr);
991
992 ierr = TSForwardSetInitialSensitivities(ts,A);CHKERRQ(ierr); /* if didp is NULL, identity matrix is assumed */
993
994 ierr = MatDestroy(&A);CHKERRQ(ierr);
995 PetscFunctionReturn(0);
996 }
997
998 /*@
999 TSAdjointResetForward - Reset the tangent linear solver and destroy the tangent linear context
1000
1001 Logically Collective on TS
1002
1003 Input Parameters:
1004 . ts - the TS context obtained from TSCreate()
1005
1006 Level: intermediate
1007
1008 .seealso: TSAdjointSetForward()
1009 @*/
TSAdjointResetForward(TS ts)1010 PetscErrorCode TSAdjointResetForward(TS ts)
1011 {
1012 PetscErrorCode ierr;
1013
1014 PetscFunctionBegin;
1015 ts->forward_solve = PETSC_FALSE; /* turn off tangent linear mode */
1016 ierr = TSForwardReset(ts);CHKERRQ(ierr);
1017 PetscFunctionReturn(0);
1018 }
1019
1020 /*@
1021 TSAdjointSetUp - Sets up the internal data structures for the later use
1022 of an adjoint solver
1023
1024 Collective on TS
1025
1026 Input Parameter:
1027 . ts - the TS context obtained from TSCreate()
1028
1029 Level: advanced
1030
1031 .seealso: TSCreate(), TSAdjointStep(), TSSetCostGradients()
1032 @*/
TSAdjointSetUp(TS ts)1033 PetscErrorCode TSAdjointSetUp(TS ts)
1034 {
1035 TSTrajectory tj;
1036 PetscBool match;
1037 PetscErrorCode ierr;
1038
1039 PetscFunctionBegin;
1040 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1041 if (ts->adjointsetupcalled) PetscFunctionReturn(0);
1042 if (!ts->vecs_sensi) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetCostGradients() first");
1043 if (ts->vecs_sensip && !ts->Jacp && !ts->Jacprhs) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetRHSJacobianP() or TSSetIJacobianP() first");
1044 ierr = TSGetTrajectory(ts,&tj);CHKERRQ(ierr);
1045 ierr = PetscObjectTypeCompare((PetscObject)tj,TSTRAJECTORYBASIC,&match);CHKERRQ(ierr);
1046 if (match) {
1047 PetscBool solution_only;
1048 ierr = TSTrajectoryGetSolutionOnly(tj,&solution_only);CHKERRQ(ierr);
1049 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");
1050 }
1051 ierr = TSTrajectorySetUseHistory(tj,PETSC_FALSE);CHKERRQ(ierr); /* not use TSHistory */
1052
1053 if (ts->quadraturets) { /* if there is integral in the cost function */
1054 ierr = VecDuplicate(ts->vecs_sensi[0],&ts->vec_drdu_col);CHKERRQ(ierr);
1055 if (ts->vecs_sensip){
1056 ierr = VecDuplicate(ts->vecs_sensip[0],&ts->vec_drdp_col);CHKERRQ(ierr);
1057 }
1058 }
1059
1060 if (ts->ops->adjointsetup) {
1061 ierr = (*ts->ops->adjointsetup)(ts);CHKERRQ(ierr);
1062 }
1063 ts->adjointsetupcalled = PETSC_TRUE;
1064 PetscFunctionReturn(0);
1065 }
1066
1067 /*@
1068 TSAdjointReset - Resets a TSAdjoint context and removes any allocated Vecs and Mats.
1069
1070 Collective on TS
1071
1072 Input Parameter:
1073 . ts - the TS context obtained from TSCreate()
1074
1075 Level: beginner
1076
1077 .seealso: TSCreate(), TSAdjointSetUp(), TSADestroy()
1078 @*/
TSAdjointReset(TS ts)1079 PetscErrorCode TSAdjointReset(TS ts)
1080 {
1081 PetscErrorCode ierr;
1082
1083 PetscFunctionBegin;
1084 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1085 if (ts->ops->adjointreset) {
1086 ierr = (*ts->ops->adjointreset)(ts);CHKERRQ(ierr);
1087 }
1088 if (ts->quadraturets) { /* if there is integral in the cost function */
1089 ierr = VecDestroy(&ts->vec_drdu_col);CHKERRQ(ierr);
1090 if (ts->vecs_sensip){
1091 ierr = VecDestroy(&ts->vec_drdp_col);CHKERRQ(ierr);
1092 }
1093 }
1094 ts->vecs_sensi = NULL;
1095 ts->vecs_sensip = NULL;
1096 ts->vecs_sensi2 = NULL;
1097 ts->vecs_sensi2p = NULL;
1098 ts->vec_dir = NULL;
1099 ts->adjointsetupcalled = PETSC_FALSE;
1100 PetscFunctionReturn(0);
1101 }
1102
1103 /*@
1104 TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time
1105
1106 Logically Collective on TS
1107
1108 Input Parameters:
1109 + ts - the TS context obtained from TSCreate()
1110 - steps - number of steps to use
1111
1112 Level: intermediate
1113
1114 Notes:
1115 Normally one does not call this and TSAdjointSolve() integrates back to the original timestep. One can call this
1116 so as to integrate back to less than the original timestep
1117
1118 .seealso: TSSetExactFinalTime()
1119 @*/
TSAdjointSetSteps(TS ts,PetscInt steps)1120 PetscErrorCode TSAdjointSetSteps(TS ts,PetscInt steps)
1121 {
1122 PetscFunctionBegin;
1123 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1124 PetscValidLogicalCollectiveInt(ts,steps,2);
1125 if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back a negative number of steps");
1126 if (steps > ts->steps) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back more than the total number of forward steps");
1127 ts->adjoint_max_steps = steps;
1128 PetscFunctionReturn(0);
1129 }
1130
1131 /*@C
1132 TSAdjointSetRHSJacobian - Deprecated, use TSSetRHSJacobianP()
1133
1134 Level: deprecated
1135
1136 @*/
TSAdjointSetRHSJacobian(TS ts,Mat Amat,PetscErrorCode (* func)(TS,PetscReal,Vec,Mat,void *),void * ctx)1137 PetscErrorCode TSAdjointSetRHSJacobian(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx)
1138 {
1139 PetscErrorCode ierr;
1140
1141 PetscFunctionBegin;
1142 PetscValidHeaderSpecific(ts, TS_CLASSID,1);
1143 PetscValidHeaderSpecific(Amat,MAT_CLASSID,2);
1144
1145 ts->rhsjacobianp = func;
1146 ts->rhsjacobianpctx = ctx;
1147 if (Amat) {
1148 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr);
1149 ierr = MatDestroy(&ts->Jacp);CHKERRQ(ierr);
1150 ts->Jacp = Amat;
1151 }
1152 PetscFunctionReturn(0);
1153 }
1154
1155 /*@C
1156 TSAdjointComputeRHSJacobian - Deprecated, use TSComputeRHSJacobianP()
1157
1158 Level: deprecated
1159
1160 @*/
TSAdjointComputeRHSJacobian(TS ts,PetscReal t,Vec U,Mat Amat)1161 PetscErrorCode TSAdjointComputeRHSJacobian(TS ts,PetscReal t,Vec U,Mat Amat)
1162 {
1163 PetscErrorCode ierr;
1164
1165 PetscFunctionBegin;
1166 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1167 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
1168 PetscValidPointer(Amat,4);
1169
1170 PetscStackPush("TS user JacobianP function for sensitivity analysis");
1171 ierr = (*ts->rhsjacobianp)(ts,t,U,Amat,ts->rhsjacobianpctx);CHKERRQ(ierr);
1172 PetscStackPop;
1173 PetscFunctionReturn(0);
1174 }
1175
1176 /*@
1177 TSAdjointComputeDRDYFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobian()
1178
1179 Level: deprecated
1180
1181 @*/
TSAdjointComputeDRDYFunction(TS ts,PetscReal t,Vec U,Vec * DRDU)1182 PetscErrorCode TSAdjointComputeDRDYFunction(TS ts,PetscReal t,Vec U,Vec *DRDU)
1183 {
1184 PetscErrorCode ierr;
1185
1186 PetscFunctionBegin;
1187 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1188 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
1189
1190 PetscStackPush("TS user DRDY function for sensitivity analysis");
1191 ierr = (*ts->drdufunction)(ts,t,U,DRDU,ts->costintegrandctx);CHKERRQ(ierr);
1192 PetscStackPop;
1193 PetscFunctionReturn(0);
1194 }
1195
1196 /*@
1197 TSAdjointComputeDRDPFunction - Deprecated, use TSGetQuadratureTS() then TSComputeRHSJacobianP()
1198
1199 Level: deprecated
1200
1201 @*/
TSAdjointComputeDRDPFunction(TS ts,PetscReal t,Vec U,Vec * DRDP)1202 PetscErrorCode TSAdjointComputeDRDPFunction(TS ts,PetscReal t,Vec U,Vec *DRDP)
1203 {
1204 PetscErrorCode ierr;
1205
1206 PetscFunctionBegin;
1207 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1208 PetscValidHeaderSpecific(U,VEC_CLASSID,3);
1209
1210 PetscStackPush("TS user DRDP function for sensitivity analysis");
1211 ierr = (*ts->drdpfunction)(ts,t,U,DRDP,ts->costintegrandctx);CHKERRQ(ierr);
1212 PetscStackPop;
1213 PetscFunctionReturn(0);
1214 }
1215
1216 /*@C
1217 TSAdjointMonitorSensi - monitors the first lambda sensitivity
1218
1219 Level: intermediate
1220
1221 .seealso: TSAdjointMonitorSet()
1222 @*/
TSAdjointMonitorSensi(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec * lambda,Vec * mu,PetscViewerAndFormat * vf)1223 PetscErrorCode TSAdjointMonitorSensi(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
1224 {
1225 PetscErrorCode ierr;
1226 PetscViewer viewer = vf->viewer;
1227
1228 PetscFunctionBegin;
1229 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
1230 ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
1231 ierr = VecView(lambda[0],viewer);CHKERRQ(ierr);
1232 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
1233 PetscFunctionReturn(0);
1234 }
1235
1236 /*@C
1237 TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
1238
1239 Collective on TS
1240
1241 Input Parameters:
1242 + ts - TS object you wish to monitor
1243 . name - the monitor type one is seeking
1244 . help - message indicating what monitoring is done
1245 . manual - manual page for the monitor
1246 . monitor - the monitor function
1247 - 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
1248
1249 Level: developer
1250
1251 .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
1252 PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
1253 PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
1254 PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
1255 PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
1256 PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
1257 PetscOptionsFList(), PetscOptionsEList()
1258 @*/
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 *))1259 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*))
1260 {
1261 PetscErrorCode ierr;
1262 PetscViewer viewer;
1263 PetscViewerFormat format;
1264 PetscBool flg;
1265
1266 PetscFunctionBegin;
1267 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject) ts)->options,((PetscObject)ts)->prefix,name,&viewer,&format,&flg);CHKERRQ(ierr);
1268 if (flg) {
1269 PetscViewerAndFormat *vf;
1270 ierr = PetscViewerAndFormatCreate(viewer,format,&vf);CHKERRQ(ierr);
1271 ierr = PetscObjectDereference((PetscObject)viewer);CHKERRQ(ierr);
1272 if (monitorsetup) {
1273 ierr = (*monitorsetup)(ts,vf);CHKERRQ(ierr);
1274 }
1275 ierr = TSAdjointMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);CHKERRQ(ierr);
1276 }
1277 PetscFunctionReturn(0);
1278 }
1279
1280 /*@C
1281 TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
1282 timestep to display the iteration's progress.
1283
1284 Logically Collective on TS
1285
1286 Input Parameters:
1287 + ts - the TS context obtained from TSCreate()
1288 . adjointmonitor - monitoring routine
1289 . adjointmctx - [optional] user-defined context for private data for the
1290 monitor routine (use NULL if no context is desired)
1291 - adjointmonitordestroy - [optional] routine that frees monitor context
1292 (may be NULL)
1293
1294 Calling sequence of monitor:
1295 $ int adjointmonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *adjointmctx)
1296
1297 + ts - the TS context
1298 . 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
1299 been interpolated to)
1300 . time - current time
1301 . u - current iterate
1302 . numcost - number of cost functionos
1303 . lambda - sensitivities to initial conditions
1304 . mu - sensitivities to parameters
1305 - adjointmctx - [optional] adjoint monitoring context
1306
1307 Notes:
1308 This routine adds an additional monitor to the list of monitors that
1309 already has been loaded.
1310
1311 Fortran Notes:
1312 Only a single monitor function can be set for each TS object
1313
1314 Level: intermediate
1315
1316 .seealso: TSAdjointMonitorCancel()
1317 @*/
TSAdjointMonitorSet(TS ts,PetscErrorCode (* adjointmonitor)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec *,Vec *,void *),void * adjointmctx,PetscErrorCode (* adjointmdestroy)(void **))1318 PetscErrorCode TSAdjointMonitorSet(TS ts,PetscErrorCode (*adjointmonitor)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*),void *adjointmctx,PetscErrorCode (*adjointmdestroy)(void**))
1319 {
1320 PetscErrorCode ierr;
1321 PetscInt i;
1322 PetscBool identical;
1323
1324 PetscFunctionBegin;
1325 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1326 for (i=0; i<ts->numbermonitors;i++) {
1327 ierr = PetscMonitorCompare((PetscErrorCode (*)(void))adjointmonitor,adjointmctx,adjointmdestroy,(PetscErrorCode (*)(void))ts->adjointmonitor[i],ts->adjointmonitorcontext[i],ts->adjointmonitordestroy[i],&identical);CHKERRQ(ierr);
1328 if (identical) PetscFunctionReturn(0);
1329 }
1330 if (ts->numberadjointmonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many adjoint monitors set");
1331 ts->adjointmonitor[ts->numberadjointmonitors] = adjointmonitor;
1332 ts->adjointmonitordestroy[ts->numberadjointmonitors] = adjointmdestroy;
1333 ts->adjointmonitorcontext[ts->numberadjointmonitors++] = (void*)adjointmctx;
1334 PetscFunctionReturn(0);
1335 }
1336
1337 /*@C
1338 TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.
1339
1340 Logically Collective on TS
1341
1342 Input Parameters:
1343 . ts - the TS context obtained from TSCreate()
1344
1345 Notes:
1346 There is no way to remove a single, specific monitor.
1347
1348 Level: intermediate
1349
1350 .seealso: TSAdjointMonitorSet()
1351 @*/
TSAdjointMonitorCancel(TS ts)1352 PetscErrorCode TSAdjointMonitorCancel(TS ts)
1353 {
1354 PetscErrorCode ierr;
1355 PetscInt i;
1356
1357 PetscFunctionBegin;
1358 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1359 for (i=0; i<ts->numberadjointmonitors; i++) {
1360 if (ts->adjointmonitordestroy[i]) {
1361 ierr = (*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]);CHKERRQ(ierr);
1362 }
1363 }
1364 ts->numberadjointmonitors = 0;
1365 PetscFunctionReturn(0);
1366 }
1367
1368 /*@C
1369 TSAdjointMonitorDefault - the default monitor of adjoint computations
1370
1371 Level: intermediate
1372
1373 .seealso: TSAdjointMonitorSet()
1374 @*/
TSAdjointMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec * lambda,Vec * mu,PetscViewerAndFormat * vf)1375 PetscErrorCode TSAdjointMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
1376 {
1377 PetscErrorCode ierr;
1378 PetscViewer viewer = vf->viewer;
1379
1380 PetscFunctionBegin;
1381 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
1382 ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
1383 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr);
1384 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");CHKERRQ(ierr);
1385 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr);
1386 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
1387 PetscFunctionReturn(0);
1388 }
1389
1390 /*@C
1391 TSAdjointMonitorDrawSensi - Monitors progress of the adjoint TS solvers by calling
1392 VecView() for the sensitivities to initial states at each timestep
1393
1394 Collective on TS
1395
1396 Input Parameters:
1397 + ts - the TS context
1398 . step - current time-step
1399 . ptime - current time
1400 . u - current state
1401 . numcost - number of cost functions
1402 . lambda - sensitivities to initial conditions
1403 . mu - sensitivities to parameters
1404 - dummy - either a viewer or NULL
1405
1406 Level: intermediate
1407
1408 .seealso: TSAdjointMonitorSet(), TSAdjointMonitorDefault(), VecView()
1409 @*/
TSAdjointMonitorDrawSensi(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec * lambda,Vec * mu,void * dummy)1410 PetscErrorCode TSAdjointMonitorDrawSensi(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda,Vec *mu,void *dummy)
1411 {
1412 PetscErrorCode ierr;
1413 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
1414 PetscDraw draw;
1415 PetscReal xl,yl,xr,yr,h;
1416 char time[32];
1417
1418 PetscFunctionBegin;
1419 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) PetscFunctionReturn(0);
1420
1421 ierr = VecView(lambda[0],ictx->viewer);CHKERRQ(ierr);
1422 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr);
1423 ierr = PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);CHKERRQ(ierr);
1424 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr);
1425 h = yl + .95*(yr - yl);
1426 ierr = PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr);
1427 ierr = PetscDrawFlush(draw);CHKERRQ(ierr);
1428 PetscFunctionReturn(0);
1429 }
1430
1431 /*
1432 TSAdjointSetFromOptions - Sets various TSAdjoint parameters from user options.
1433
1434 Collective on TSAdjoint
1435
1436 Input Parameter:
1437 . ts - the TS context
1438
1439 Options Database Keys:
1440 + -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
1441 . -ts_adjoint_monitor - print information at each adjoint time step
1442 - -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically
1443
1444 Level: developer
1445
1446 Notes:
1447 This is not normally called directly by users
1448
1449 .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
1450 */
TSAdjointSetFromOptions(PetscOptionItems * PetscOptionsObject,TS ts)1451 PetscErrorCode TSAdjointSetFromOptions(PetscOptionItems *PetscOptionsObject,TS ts)
1452 {
1453 PetscBool tflg,opt;
1454 PetscErrorCode ierr;
1455
1456 PetscFunctionBegin;
1457 PetscValidHeaderSpecific(ts,TS_CLASSID,2);
1458 ierr = PetscOptionsHead(PetscOptionsObject,"TS Adjoint options");CHKERRQ(ierr);
1459 tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
1460 ierr = PetscOptionsBool("-ts_adjoint_solve","Solve the adjoint problem immediately after solving the forward problem","",tflg,&tflg,&opt);CHKERRQ(ierr);
1461 if (opt) {
1462 ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr);
1463 ts->adjoint_solve = tflg;
1464 }
1465 ierr = TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor","Monitor adjoint timestep size","TSAdjointMonitorDefault",TSAdjointMonitorDefault,NULL);CHKERRQ(ierr);
1466 ierr = TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor_sensi","Monitor sensitivity in the adjoint computation","TSAdjointMonitorSensi",TSAdjointMonitorSensi,NULL);CHKERRQ(ierr);
1467 opt = PETSC_FALSE;
1468 ierr = PetscOptionsName("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",&opt);CHKERRQ(ierr);
1469 if (opt) {
1470 TSMonitorDrawCtx ctx;
1471 PetscInt howoften = 1;
1472
1473 ierr = PetscOptionsInt("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",howoften,&howoften,NULL);CHKERRQ(ierr);
1474 ierr = TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);CHKERRQ(ierr);
1475 ierr = TSAdjointMonitorSet(ts,TSAdjointMonitorDrawSensi,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);CHKERRQ(ierr);
1476 }
1477 PetscFunctionReturn(0);
1478 }
1479
1480 /*@
1481 TSAdjointStep - Steps one time step backward in the adjoint run
1482
1483 Collective on TS
1484
1485 Input Parameter:
1486 . ts - the TS context obtained from TSCreate()
1487
1488 Level: intermediate
1489
1490 .seealso: TSAdjointSetUp(), TSAdjointSolve()
1491 @*/
TSAdjointStep(TS ts)1492 PetscErrorCode TSAdjointStep(TS ts)
1493 {
1494 DM dm;
1495 PetscErrorCode ierr;
1496
1497 PetscFunctionBegin;
1498 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1499 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
1500 ierr = TSAdjointSetUp(ts);CHKERRQ(ierr);
1501 ts->steps--; /* must decrease the step index before the adjoint step is taken. */
1502
1503 ts->reason = TS_CONVERGED_ITERATING;
1504 ts->ptime_prev = ts->ptime;
1505 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);
1506 ierr = PetscLogEventBegin(TS_AdjointStep,ts,0,0,0);CHKERRQ(ierr);
1507 ierr = (*ts->ops->adjointstep)(ts);CHKERRQ(ierr);
1508 ierr = PetscLogEventEnd(TS_AdjointStep,ts,0,0,0);CHKERRQ(ierr);
1509 ts->adjoint_steps++;
1510
1511 if (ts->reason < 0) {
1512 if (ts->errorifstepfailed) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSAdjointStep has failed due to %s",TSConvergedReasons[ts->reason]);
1513 } else if (!ts->reason) {
1514 if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1515 }
1516 PetscFunctionReturn(0);
1517 }
1518
1519 /*@
1520 TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE
1521
1522 Collective on TS
1523
1524 Input Parameter:
1525 . ts - the TS context obtained from TSCreate()
1526
1527 Options Database:
1528 . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial values
1529
1530 Level: intermediate
1531
1532 Notes:
1533 This must be called after a call to TSSolve() that solves the forward problem
1534
1535 By default this will integrate back to the initial time, one can use TSAdjointSetSteps() to step back to a later time
1536
1537 .seealso: TSCreate(), TSSetCostGradients(), TSSetSolution(), TSAdjointStep()
1538 @*/
TSAdjointSolve(TS ts)1539 PetscErrorCode TSAdjointSolve(TS ts)
1540 {
1541 static PetscBool cite = PETSC_FALSE;
1542 #if defined(TSADJOINT_STAGE)
1543 PetscLogStage adjoint_stage;
1544 #endif
1545 PetscErrorCode ierr;
1546
1547 PetscFunctionBegin;
1548 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1549 ierr = PetscCitationsRegister("@article{tsadjointpaper,\n"
1550 " title = {{PETSc TSAdjoint: a discrete adjoint ODE solver for first-order and second-order sensitivity analysis}},\n"
1551 " author = {Zhang, Hong and Constantinescu, Emil M. and Smith, Barry F.},\n"
1552 " journal = {arXiv e-preprints},\n"
1553 " eprint = {1912.07696},\n"
1554 " archivePrefix = {arXiv},\n"
1555 " year = {2019}\n}\n",&cite);CHKERRQ(ierr);
1556 #if defined(TSADJOINT_STAGE)
1557 ierr = PetscLogStageRegister("TSAdjoint",&adjoint_stage);CHKERRQ(ierr);
1558 ierr = PetscLogStagePush(adjoint_stage);CHKERRQ(ierr);
1559 #endif
1560 ierr = TSAdjointSetUp(ts);CHKERRQ(ierr);
1561
1562 /* reset time step and iteration counters */
1563 ts->adjoint_steps = 0;
1564 ts->ksp_its = 0;
1565 ts->snes_its = 0;
1566 ts->num_snes_failures = 0;
1567 ts->reject = 0;
1568 ts->reason = TS_CONVERGED_ITERATING;
1569
1570 if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->steps;
1571 if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
1572
1573 while (!ts->reason) {
1574 ierr = TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);CHKERRQ(ierr);
1575 ierr = TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);CHKERRQ(ierr);
1576 ierr = TSAdjointEventHandler(ts);CHKERRQ(ierr);
1577 ierr = TSAdjointStep(ts);CHKERRQ(ierr);
1578 if ((ts->vec_costintegral || ts->quadraturets) && !ts->costintegralfwd) {
1579 ierr = TSAdjointCostIntegral(ts);CHKERRQ(ierr);
1580 }
1581 }
1582 ierr = TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);CHKERRQ(ierr);
1583 ierr = TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);CHKERRQ(ierr);
1584 ts->solvetime = ts->ptime;
1585 ierr = TSTrajectoryViewFromOptions(ts->trajectory,NULL,"-ts_trajectory_view");CHKERRQ(ierr);
1586 ierr = VecViewFromOptions(ts->vecs_sensi[0],(PetscObject) ts, "-ts_adjoint_view_solution");CHKERRQ(ierr);
1587 ts->adjoint_max_steps = 0;
1588 #if defined(TSADJOINT_STAGE)
1589 ierr = PetscLogStagePop();CHKERRQ(ierr);
1590 #endif
1591 PetscFunctionReturn(0);
1592 }
1593
1594 /*@C
1595 TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using TSAdjointMonitorSet()
1596
1597 Collective on TS
1598
1599 Input Parameters:
1600 + ts - time stepping context obtained from TSCreate()
1601 . step - step number that has just completed
1602 . ptime - model time of the state
1603 . u - state at the current model time
1604 . numcost - number of cost functions (dimension of lambda or mu)
1605 . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
1606 - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
1607
1608 Notes:
1609 TSAdjointMonitor() is typically used automatically within the time stepping implementations.
1610 Users would almost never call this routine directly.
1611
1612 Level: developer
1613
1614 @*/
TSAdjointMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec * lambda,Vec * mu)1615 PetscErrorCode TSAdjointMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda, Vec *mu)
1616 {
1617 PetscErrorCode ierr;
1618 PetscInt i,n = ts->numberadjointmonitors;
1619
1620 PetscFunctionBegin;
1621 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1622 PetscValidHeaderSpecific(u,VEC_CLASSID,4);
1623 ierr = VecLockReadPush(u);CHKERRQ(ierr);
1624 for (i=0; i<n; i++) {
1625 ierr = (*ts->adjointmonitor[i])(ts,step,ptime,u,numcost,lambda,mu,ts->adjointmonitorcontext[i]);CHKERRQ(ierr);
1626 }
1627 ierr = VecLockReadPop(u);CHKERRQ(ierr);
1628 PetscFunctionReturn(0);
1629 }
1630
1631 /*@
1632 TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.
1633
1634 Collective on TS
1635
1636 Input Arguments:
1637 . ts - time stepping context
1638
1639 Level: advanced
1640
1641 Notes:
1642 This function cannot be called until TSAdjointStep() has been completed.
1643
1644 .seealso: TSAdjointSolve(), TSAdjointStep
1645 @*/
TSAdjointCostIntegral(TS ts)1646 PetscErrorCode TSAdjointCostIntegral(TS ts)
1647 {
1648 PetscErrorCode ierr;
1649 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1650 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);
1651 ierr = (*ts->ops->adjointintegral)(ts);CHKERRQ(ierr);
1652 PetscFunctionReturn(0);
1653 }
1654
1655 /* ------------------ Forward (tangent linear) sensitivity ------------------*/
1656
1657 /*@
1658 TSForwardSetUp - Sets up the internal data structures for the later use
1659 of forward sensitivity analysis
1660
1661 Collective on TS
1662
1663 Input Parameter:
1664 . ts - the TS context obtained from TSCreate()
1665
1666 Level: advanced
1667
1668 .seealso: TSCreate(), TSDestroy(), TSSetUp()
1669 @*/
TSForwardSetUp(TS ts)1670 PetscErrorCode TSForwardSetUp(TS ts)
1671 {
1672 PetscErrorCode ierr;
1673
1674 PetscFunctionBegin;
1675 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1676 if (ts->forwardsetupcalled) PetscFunctionReturn(0);
1677 if (ts->ops->forwardsetup) {
1678 ierr = (*ts->ops->forwardsetup)(ts);CHKERRQ(ierr);
1679 }
1680 ierr = VecDuplicate(ts->vec_sol,&ts->vec_sensip_col);CHKERRQ(ierr);
1681 ts->forwardsetupcalled = PETSC_TRUE;
1682 PetscFunctionReturn(0);
1683 }
1684
1685 /*@
1686 TSForwardReset - Reset the internal data structures used by forward sensitivity analysis
1687
1688 Collective on TS
1689
1690 Input Parameter:
1691 . ts - the TS context obtained from TSCreate()
1692
1693 Level: advanced
1694
1695 .seealso: TSCreate(), TSDestroy(), TSForwardSetUp()
1696 @*/
TSForwardReset(TS ts)1697 PetscErrorCode TSForwardReset(TS ts)
1698 {
1699 TS quadts = ts->quadraturets;
1700 PetscErrorCode ierr;
1701
1702 PetscFunctionBegin;
1703 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1704 if (ts->ops->forwardreset) {
1705 ierr = (*ts->ops->forwardreset)(ts);CHKERRQ(ierr);
1706 }
1707 ierr = MatDestroy(&ts->mat_sensip);CHKERRQ(ierr);
1708 if (quadts) {
1709 ierr = MatDestroy(&quadts->mat_sensip);CHKERRQ(ierr);
1710 }
1711 ierr = VecDestroy(&ts->vec_sensip_col);CHKERRQ(ierr);
1712 ts->forward_solve = PETSC_FALSE;
1713 ts->forwardsetupcalled = PETSC_FALSE;
1714 PetscFunctionReturn(0);
1715 }
1716
1717 /*@
1718 TSForwardSetIntegralGradients - Set the vectors holding forward sensitivities of the integral term.
1719
1720 Input Parameter:
1721 + ts- the TS context obtained from TSCreate()
1722 . numfwdint- number of integrals
1723 - vp = the vectors containing the gradients for each integral w.r.t. parameters
1724
1725 Level: deprecated
1726
1727 .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1728 @*/
TSForwardSetIntegralGradients(TS ts,PetscInt numfwdint,Vec * vp)1729 PetscErrorCode TSForwardSetIntegralGradients(TS ts,PetscInt numfwdint,Vec *vp)
1730 {
1731 PetscFunctionBegin;
1732 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1733 if (ts->numcost && ts->numcost!=numfwdint) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2rd parameter of TSSetCostIntegrand()) is inconsistent with the one set by TSSetCostIntegrand()");
1734 if (!ts->numcost) ts->numcost = numfwdint;
1735
1736 ts->vecs_integral_sensip = vp;
1737 PetscFunctionReturn(0);
1738 }
1739
1740 /*@
1741 TSForwardGetIntegralGradients - Returns the forward sensitivities ofthe integral term.
1742
1743 Input Parameter:
1744 . ts- the TS context obtained from TSCreate()
1745
1746 Output Parameter:
1747 . vp = the vectors containing the gradients for each integral w.r.t. parameters
1748
1749 Level: deprecated
1750
1751 .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1752 @*/
TSForwardGetIntegralGradients(TS ts,PetscInt * numfwdint,Vec ** vp)1753 PetscErrorCode TSForwardGetIntegralGradients(TS ts,PetscInt *numfwdint,Vec **vp)
1754 {
1755 PetscFunctionBegin;
1756 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1757 PetscValidPointer(vp,3);
1758 if (numfwdint) *numfwdint = ts->numcost;
1759 if (vp) *vp = ts->vecs_integral_sensip;
1760 PetscFunctionReturn(0);
1761 }
1762
1763 /*@
1764 TSForwardStep - Compute the forward sensitivity for one time step.
1765
1766 Collective on TS
1767
1768 Input Arguments:
1769 . ts - time stepping context
1770
1771 Level: advanced
1772
1773 Notes:
1774 This function cannot be called until TSStep() has been completed.
1775
1776 .seealso: TSForwardSetSensitivities(), TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardSetUp()
1777 @*/
TSForwardStep(TS ts)1778 PetscErrorCode TSForwardStep(TS ts)
1779 {
1780 PetscErrorCode ierr;
1781 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1782 if (!ts->ops->forwardstep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide forward sensitivity analysis",((PetscObject)ts)->type_name);
1783 ierr = PetscLogEventBegin(TS_ForwardStep,ts,0,0,0);CHKERRQ(ierr);
1784 ierr = (*ts->ops->forwardstep)(ts);CHKERRQ(ierr);
1785 ierr = PetscLogEventEnd(TS_ForwardStep,ts,0,0,0);CHKERRQ(ierr);
1786 if (ts->reason < 0 && ts->errorifstepfailed) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSFowardStep has failed due to %s",TSConvergedReasons[ts->reason]);
1787 PetscFunctionReturn(0);
1788 }
1789
1790 /*@
1791 TSForwardSetSensitivities - Sets the initial value of the trajectory sensitivities of solution w.r.t. the problem parameters and initial values.
1792
1793 Logically Collective on TS
1794
1795 Input Parameters:
1796 + ts - the TS context obtained from TSCreate()
1797 . nump - number of parameters
1798 - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1799
1800 Level: beginner
1801
1802 Notes:
1803 Forward sensitivity is also called 'trajectory sensitivity' in some fields such as power systems.
1804 This function turns on a flag to trigger TSSolve() to compute forward sensitivities automatically.
1805 You must call this function before TSSolve().
1806 The entries in the sensitivity matrix must be correctly initialized with the values S = dy/dp|startingtime.
1807
1808 .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1809 @*/
TSForwardSetSensitivities(TS ts,PetscInt nump,Mat Smat)1810 PetscErrorCode TSForwardSetSensitivities(TS ts,PetscInt nump,Mat Smat)
1811 {
1812 PetscErrorCode ierr;
1813
1814 PetscFunctionBegin;
1815 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1816 PetscValidHeaderSpecific(Smat,MAT_CLASSID,3);
1817 ts->forward_solve = PETSC_TRUE;
1818 if (nump == PETSC_DEFAULT) {
1819 ierr = MatGetSize(Smat,NULL,&ts->num_parameters);CHKERRQ(ierr);
1820 } else ts->num_parameters = nump;
1821 ierr = PetscObjectReference((PetscObject)Smat);CHKERRQ(ierr);
1822 ierr = MatDestroy(&ts->mat_sensip);CHKERRQ(ierr);
1823 ts->mat_sensip = Smat;
1824 PetscFunctionReturn(0);
1825 }
1826
1827 /*@
1828 TSForwardGetSensitivities - Returns the trajectory sensitivities
1829
1830 Not Collective, but Vec returned is parallel if TS is parallel
1831
1832 Output Parameter:
1833 + ts - the TS context obtained from TSCreate()
1834 . nump - number of parameters
1835 - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1836
1837 Level: intermediate
1838
1839 .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1840 @*/
TSForwardGetSensitivities(TS ts,PetscInt * nump,Mat * Smat)1841 PetscErrorCode TSForwardGetSensitivities(TS ts,PetscInt *nump,Mat *Smat)
1842 {
1843 PetscFunctionBegin;
1844 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1845 if (nump) *nump = ts->num_parameters;
1846 if (Smat) *Smat = ts->mat_sensip;
1847 PetscFunctionReturn(0);
1848 }
1849
1850 /*@
1851 TSForwardCostIntegral - Evaluate the cost integral in the forward run.
1852
1853 Collective on TS
1854
1855 Input Arguments:
1856 . ts - time stepping context
1857
1858 Level: advanced
1859
1860 Notes:
1861 This function cannot be called until TSStep() has been completed.
1862
1863 .seealso: TSSolve(), TSAdjointCostIntegral()
1864 @*/
TSForwardCostIntegral(TS ts)1865 PetscErrorCode TSForwardCostIntegral(TS ts)
1866 {
1867 PetscErrorCode ierr;
1868 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1869 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);
1870 ierr = (*ts->ops->forwardintegral)(ts);CHKERRQ(ierr);
1871 PetscFunctionReturn(0);
1872 }
1873
1874 /*@
1875 TSForwardSetInitialSensitivities - Set initial values for tangent linear sensitivities
1876
1877 Collective on TS
1878
1879 Input Parameter:
1880 + ts - the TS context obtained from TSCreate()
1881 - didp - parametric sensitivities of the initial condition
1882
1883 Level: intermediate
1884
1885 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.
1886
1887 .seealso: TSForwardSetSensitivities()
1888 @*/
TSForwardSetInitialSensitivities(TS ts,Mat didp)1889 PetscErrorCode TSForwardSetInitialSensitivities(TS ts,Mat didp)
1890 {
1891 PetscErrorCode ierr;
1892
1893 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1894 PetscValidHeaderSpecific(didp,MAT_CLASSID,2);
1895 if (!ts->mat_sensip) {
1896 ierr = TSForwardSetSensitivities(ts,PETSC_DEFAULT,didp);CHKERRQ(ierr);
1897 }
1898 PetscFunctionReturn(0);
1899 }
1900
1901 /*@
1902 TSForwardGetStages - Get the number of stages and the tangent linear sensitivities at the intermediate stages
1903
1904 Input Parameter:
1905 . ts - the TS context obtained from TSCreate()
1906
1907 Output Parameters:
1908 + ns - number of stages
1909 - S - tangent linear sensitivities at the intermediate stages
1910
1911 Level: advanced
1912
1913 @*/
TSForwardGetStages(TS ts,PetscInt * ns,Mat ** S)1914 PetscErrorCode TSForwardGetStages(TS ts,PetscInt *ns,Mat **S)
1915 {
1916 PetscErrorCode ierr;
1917
1918 PetscFunctionBegin;
1919 PetscValidHeaderSpecific(ts, TS_CLASSID,1);
1920
1921 if (!ts->ops->getstages) *S=NULL;
1922 else {
1923 ierr = (*ts->ops->forwardgetstages)(ts,ns,S);CHKERRQ(ierr);
1924 }
1925 PetscFunctionReturn(0);
1926 }
1927
1928 /*@
1929 TSCreateQuadratureTS - Create a sub-TS that evaluates integrals over time
1930
1931 Input Parameter:
1932 + ts - the TS context obtained from TSCreate()
1933 - fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1934
1935 Output Parameters:
1936 . quadts - the child TS context
1937
1938 Level: intermediate
1939
1940 .seealso: TSGetQuadratureTS()
1941 @*/
TSCreateQuadratureTS(TS ts,PetscBool fwd,TS * quadts)1942 PetscErrorCode TSCreateQuadratureTS(TS ts,PetscBool fwd,TS *quadts)
1943 {
1944 char prefix[128];
1945 PetscErrorCode ierr;
1946
1947 PetscFunctionBegin;
1948 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1949 PetscValidPointer(quadts,2);
1950 ierr = TSDestroy(&ts->quadraturets);CHKERRQ(ierr);
1951 ierr = TSCreate(PetscObjectComm((PetscObject)ts),&ts->quadraturets);CHKERRQ(ierr);
1952 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->quadraturets,(PetscObject)ts,1);CHKERRQ(ierr);
1953 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->quadraturets);CHKERRQ(ierr);
1954 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%squad_",((PetscObject)ts)->prefix ? ((PetscObject)ts)->prefix : "");CHKERRQ(ierr);
1955 ierr = TSSetOptionsPrefix(ts->quadraturets,prefix);CHKERRQ(ierr);
1956 *quadts = ts->quadraturets;
1957
1958 if (ts->numcost) {
1959 ierr = VecCreateSeq(PETSC_COMM_SELF,ts->numcost,&(*quadts)->vec_sol);CHKERRQ(ierr);
1960 } else {
1961 ierr = VecCreateSeq(PETSC_COMM_SELF,1,&(*quadts)->vec_sol);CHKERRQ(ierr);
1962 }
1963 ts->costintegralfwd = fwd;
1964 PetscFunctionReturn(0);
1965 }
1966
1967 /*@
1968 TSGetQuadratureTS - Return the sub-TS that evaluates integrals over time
1969
1970 Input Parameter:
1971 . ts - the TS context obtained from TSCreate()
1972
1973 Output Parameters:
1974 + fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
1975 - quadts - the child TS context
1976
1977 Level: intermediate
1978
1979 .seealso: TSCreateQuadratureTS()
1980 @*/
TSGetQuadratureTS(TS ts,PetscBool * fwd,TS * quadts)1981 PetscErrorCode TSGetQuadratureTS(TS ts,PetscBool *fwd,TS *quadts)
1982 {
1983 PetscFunctionBegin;
1984 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1985 if (fwd) *fwd = ts->costintegralfwd;
1986 if (quadts) *quadts = ts->quadraturets;
1987 PetscFunctionReturn(0);
1988 }
1989
1990 /*@
1991 TSComputeSNESJacobian - Compute the SNESJacobian
1992
1993 Input Parameters:
1994 + ts - the TS context obtained from TSCreate()
1995 - x - state vector
1996
1997 Output Parameters:
1998 + J - Jacobian matrix
1999 - Jpre - preconditioning matrix for J (may be same as J)
2000
2001 Level: developer
2002
2003 Notes:
2004 Using SNES to compute the Jacobian enables finite differencing when TS Jacobian is not available.
2005 @*/
TSComputeSNESJacobian(TS ts,Vec x,Mat J,Mat Jpre)2006 PetscErrorCode TSComputeSNESJacobian(TS ts,Vec x,Mat J,Mat Jpre)
2007 {
2008 SNES snes = ts->snes;
2009 PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*) = NULL;
2010 PetscErrorCode ierr;
2011
2012 PetscFunctionBegin;
2013 /*
2014 Unlike implicit methods, explicit methods do not have SNESMatFDColoring in the snes object
2015 because SNESSolve() has not been called yet; so querying SNESMatFDColoring does not work for
2016 explicit methods. Instead, we check the Jacobian compute function directly to determin if FD
2017 coloring is used.
2018 */
2019 ierr = SNESGetJacobian(snes,NULL,NULL,&jac,NULL);CHKERRQ(ierr);
2020 if (jac == SNESComputeJacobianDefaultColor) {
2021 Vec f;
2022 ierr = SNESSetSolution(snes,x);CHKERRQ(ierr);
2023 ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr);
2024 /* Force MatFDColoringApply to evaluate the SNES residual function for the base vector */
2025 ierr = SNESComputeFunction(snes,x,f);CHKERRQ(ierr);
2026 }
2027 ierr = SNESComputeJacobian(snes,x,J,Jpre);CHKERRQ(ierr);
2028 PetscFunctionReturn(0);
2029 }
2030