Actual source code: tssen.c
petsc-3.11.4 2019-09-28
1: #include <petsc/private/tsimpl.h>
2: #include <petscdraw.h>
4: PetscLogEvent TS_AdjointStep, TS_ForwardStep;
6: /* ------------------------ Sensitivity Context ---------------------------*/
8: /*@C
9: TSSetRHSJacobianP - Sets the function that computes the Jacobian of G w.r.t. the parameters p where y_t = G(y,p,t), as well as the location to store the matrix.
11: Logically Collective on TS
13: Input Parameters:
14: + ts - The TS context obtained from TSCreate()
15: - func - The function
17: Calling sequence of func:
18: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
19: + t - current timestep
20: . y - input vector (current ODE solution)
21: . A - output matrix
22: - ctx - [optional] user-defined function context
24: Level: intermediate
26: Notes:
27: 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
29: .keywords: TS, sensitivity
30: .seealso:
31: @*/
32: PetscErrorCode TSSetRHSJacobianP(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx)
33: {
40: ts->rhsjacobianp = func;
41: ts->rhsjacobianpctx = ctx;
42: if(Amat) {
43: PetscObjectReference((PetscObject)Amat);
44: MatDestroy(&ts->Jacp);
45: ts->Jacp = Amat;
46: }
47: return(0);
48: }
50: /*@C
51: TSComputeRHSJacobianP - Runs the user-defined JacobianP function.
53: Collective on TS
55: Input Parameters:
56: . ts - The TS context obtained from TSCreate()
58: Level: developer
60: .keywords: TS, sensitivity
61: .seealso: TSSetRHSJacobianP()
62: @*/
63: PetscErrorCode TSComputeRHSJacobianP(TS ts,PetscReal t,Vec X,Mat Amat)
64: {
72: PetscStackPush("TS user JacobianP function for sensitivity analysis");
73: (*ts->rhsjacobianp)(ts,t,X,Amat,ts->rhsjacobianpctx);
74: PetscStackPop;
75: return(0);
76: }
78: /*@C
79: TSSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions
81: Logically Collective on TS
83: Input Parameters:
84: + ts - the TS context obtained from TSCreate()
85: . numcost - number of gradients to be computed, this is the number of cost functions
86: . costintegral - vector that stores the integral values
87: . rf - routine for evaluating the integrand function
88: . drdyf - function that computes the gradients of the r's with respect to y
89: . 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)
90: . fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
91: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
93: Calling sequence of rf:
94: $ PetscErrorCode rf(TS ts,PetscReal t,Vec y,Vec f,void *ctx);
96: Calling sequence of drdyf:
97: $ PetscErroCode drdyf(TS ts,PetscReal t,Vec y,Vec *drdy,void *ctx);
99: Calling sequence of drdpf:
100: $ PetscErroCode drdpf(TS ts,PetscReal t,Vec y,Vec *drdp,void *ctx);
102: Level: intermediate
104: Notes:
105: For optimization there is usually a single cost function (numcost = 1). For sensitivities there may be multiple cost functions
107: .keywords: TS, sensitivity analysis, timestep, set, quadrature, function
109: .seealso: TSSetRHSJacobianP(), TSGetCostGradients(), TSSetCostGradients()
110: @*/
111: PetscErrorCode TSSetCostIntegrand(TS ts,PetscInt numcost,Vec costintegral,PetscErrorCode (*rf)(TS,PetscReal,Vec,Vec,void*),
112: PetscErrorCode (*drdyf)(TS,PetscReal,Vec,Vec*,void*),
113: PetscErrorCode (*drdpf)(TS,PetscReal,Vec,Vec*,void*),
114: PetscBool fwd,void *ctx)
115: {
121: 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()");
122: if (!ts->numcost) ts->numcost=numcost;
124: if (costintegral) {
125: PetscObjectReference((PetscObject)costintegral);
126: VecDestroy(&ts->vec_costintegral);
127: ts->vec_costintegral = costintegral;
128: } else {
129: if (!ts->vec_costintegral) { /* Create a seq vec if user does not provide one */
130: VecCreateSeq(PETSC_COMM_SELF,numcost,&ts->vec_costintegral);
131: } else {
132: VecSet(ts->vec_costintegral,0.0);
133: }
134: }
135: if (!ts->vec_costintegrand) {
136: VecDuplicate(ts->vec_costintegral,&ts->vec_costintegrand);
137: } else {
138: VecSet(ts->vec_costintegrand,0.0);
139: }
140: ts->costintegralfwd = fwd; /* Evaluate the cost integral in forward run if fwd is true */
141: ts->costintegrand = rf;
142: ts->costintegrandctx = ctx;
143: ts->drdyfunction = drdyf;
144: ts->drdpfunction = drdpf;
145: return(0);
146: }
148: /*@
149: TSGetCostIntegral - Returns the values of the integral term in the cost functions.
150: It is valid to call the routine after a backward run.
152: Not Collective
154: Input Parameter:
155: . ts - the TS context obtained from TSCreate()
157: Output Parameter:
158: . v - the vector containing the integrals for each cost function
160: Level: intermediate
162: .seealso: TSSetCostIntegrand()
164: .keywords: TS, sensitivity analysis
165: @*/
166: PetscErrorCode TSGetCostIntegral(TS ts,Vec *v)
167: {
171: *v = ts->vec_costintegral;
172: return(0);
173: }
175: /*@
176: TSComputeCostIntegrand - Evaluates the integral function in the cost functions.
178: Input Parameters:
179: + ts - the TS context
180: . t - current time
181: - y - state vector, i.e. current solution
183: Output Parameter:
184: . q - vector of size numcost to hold the outputs
186: Note:
187: Most users should not need to explicitly call this routine, as it
188: is used internally within the sensitivity analysis context.
190: Level: developer
192: .keywords: TS, compute
194: .seealso: TSSetCostIntegrand()
195: @*/
196: PetscErrorCode TSComputeCostIntegrand(TS ts,PetscReal t,Vec y,Vec q)
197: {
205: PetscLogEventBegin(TS_FunctionEval,ts,y,q,0);
206: if (ts->costintegrand) {
207: PetscStackPush("TS user integrand in the cost function");
208: (*ts->costintegrand)(ts,t,y,q,ts->costintegrandctx);
209: PetscStackPop;
210: } else {
211: VecZeroEntries(q);
212: }
214: PetscLogEventEnd(TS_FunctionEval,ts,y,q,0);
215: return(0);
216: }
218: /*@
219: TSComputeDRDYFunction - Runs the user-defined DRDY function.
221: Collective on TS
223: Input Parameters:
224: . ts - The TS context obtained from TSCreate()
226: Notes:
227: TSAdjointComputeDRDYFunction() is typically used for sensitivity implementation,
228: so most users would not generally call this routine themselves.
230: Level: developer
232: .keywords: TS, sensitivity
233: .seealso: TSSetCostIntegrand()
234: @*/
235: PetscErrorCode TSComputeDRDYFunction(TS ts,PetscReal t,Vec y,Vec *drdy)
236: {
243: PetscStackPush("TS user DRDY function for sensitivity analysis");
244: (*ts->drdyfunction)(ts,t,y,drdy,ts->costintegrandctx);
245: PetscStackPop;
246: return(0);
247: }
249: /*@
250: TSComputeDRDPFunction - Runs the user-defined DRDP function.
252: Collective on TS
254: Input Parameters:
255: . ts - The TS context obtained from TSCreate()
257: Notes:
258: TSDRDPFunction() is typically used for sensitivity implementation,
259: so most users would not generally call this routine themselves.
261: Level: developer
263: .keywords: TS, sensitivity
264: .seealso: TSSetCostIntegrand()
265: @*/
266: PetscErrorCode TSComputeDRDPFunction(TS ts,PetscReal t,Vec y,Vec *drdp)
267: {
274: PetscStackPush("TS user DRDP function for sensitivity analysis");
275: (*ts->drdpfunction)(ts,t,y,drdp,ts->costintegrandctx);
276: PetscStackPop;
277: return(0);
278: }
280: /* --------------------------- Adjoint sensitivity ---------------------------*/
282: /*@
283: TSSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial values and w.r.t. the problem parameters
284: for use by the TSAdjoint routines.
286: Logically Collective on TS and Vec
288: Input Parameters:
289: + ts - the TS context obtained from TSCreate()
290: . 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
291: - mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
293: Level: beginner
295: Notes:
296: the entries in these vectors must be correctly initialized with the values lamda_i = df/dy|finaltime mu_i = df/dp|finaltime
298: After TSAdjointSolve() is called the lamba and the mu contain the computed sensitivities
300: .keywords: TS, timestep, set, sensitivity, initial values
301: @*/
302: PetscErrorCode TSSetCostGradients(TS ts,PetscInt numcost,Vec *lambda,Vec *mu)
303: {
307: ts->vecs_sensi = lambda;
308: ts->vecs_sensip = mu;
309: 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");
310: ts->numcost = numcost;
311: return(0);
312: }
314: /*@
315: TSGetCostGradients - Returns the gradients from the TSAdjointSolve()
317: Not Collective, but Vec returned is parallel if TS is parallel
319: Input Parameter:
320: . ts - the TS context obtained from TSCreate()
322: Output Parameter:
323: + lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
324: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
326: Level: intermediate
328: .seealso: TSGetTimeStep()
330: .keywords: TS, timestep, get, sensitivity
331: @*/
332: PetscErrorCode TSGetCostGradients(TS ts,PetscInt *numcost,Vec **lambda,Vec **mu)
333: {
336: if (numcost) *numcost = ts->numcost;
337: if (lambda) *lambda = ts->vecs_sensi;
338: if (mu) *mu = ts->vecs_sensip;
339: return(0);
340: }
342: /*@
343: TSAdjointSetUp - Sets up the internal data structures for the later use
344: of an adjoint solver
346: Collective on TS
348: Input Parameter:
349: . ts - the TS context obtained from TSCreate()
351: Level: advanced
353: .keywords: TS, timestep, setup
355: .seealso: TSCreate(), TSAdjointStep(), TSSetCostGradients()
356: @*/
357: PetscErrorCode TSAdjointSetUp(TS ts)
358: {
363: if (ts->adjointsetupcalled) return(0);
364: if (!ts->vecs_sensi) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetCostGradients() first");
365: if (ts->vecs_sensip && !ts->Jacp) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSAdjointSetRHSJacobian() first");
367: if (ts->vec_costintegral) { /* if there is integral in the cost function */
368: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&ts->vecs_drdy);
369: if (ts->vecs_sensip){
370: VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&ts->vecs_drdp);
371: }
372: }
374: if (ts->ops->adjointsetup) {
375: (*ts->ops->adjointsetup)(ts);
376: }
377: ts->adjointsetupcalled = PETSC_TRUE;
378: return(0);
379: }
381: /*@
382: TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time
384: Logically Collective on TS
386: Input Parameters:
387: + ts - the TS context obtained from TSCreate()
388: . steps - number of steps to use
390: Level: intermediate
392: Notes:
393: Normally one does not call this and TSAdjointSolve() integrates back to the original timestep. One can call this
394: so as to integrate back to less than the original timestep
396: .keywords: TS, timestep, set, maximum, iterations
398: .seealso: TSSetExactFinalTime()
399: @*/
400: PetscErrorCode TSAdjointSetSteps(TS ts,PetscInt steps)
401: {
405: if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back a negative number of steps");
406: if (steps > ts->steps) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back more than the total number of forward steps");
407: ts->adjoint_max_steps = steps;
408: return(0);
409: }
411: /*@C
412: TSAdjointSetRHSJacobian - Deprecated, use TSSetRHSJacobianP()
414: Level: deprecated
416: @*/
417: PetscErrorCode TSAdjointSetRHSJacobian(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx)
418: {
425: ts->rhsjacobianp = func;
426: ts->rhsjacobianpctx = ctx;
427: if(Amat) {
428: PetscObjectReference((PetscObject)Amat);
429: MatDestroy(&ts->Jacp);
430: ts->Jacp = Amat;
431: }
432: return(0);
433: }
435: /*@C
436: TSAdjointComputeRHSJacobian - Deprecated, use TSComputeRHSJacobianP()
438: Level: deprecated
440: @*/
441: PetscErrorCode TSAdjointComputeRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat)
442: {
450: PetscStackPush("TS user JacobianP function for sensitivity analysis");
451: (*ts->rhsjacobianp)(ts,t,X,Amat,ts->rhsjacobianpctx);
452: PetscStackPop;
453: return(0);
454: }
456: /*@
457: TSAdjointComputeDRDYFunction - Deprecated, use TSComputeDRDYFunction()
459: Level: deprecated
461: @*/
462: PetscErrorCode TSAdjointComputeDRDYFunction(TS ts,PetscReal t,Vec y,Vec *drdy)
463: {
470: PetscStackPush("TS user DRDY function for sensitivity analysis");
471: (*ts->drdyfunction)(ts,t,y,drdy,ts->costintegrandctx);
472: PetscStackPop;
473: return(0);
474: }
476: /*@
477: TSAdjointComputeDRDPFunction - Deprecated, use TSComputeDRDPFunction()
479: Level: deprecated
481: @*/
482: PetscErrorCode TSAdjointComputeDRDPFunction(TS ts,PetscReal t,Vec y,Vec *drdp)
483: {
490: PetscStackPush("TS user DRDP function for sensitivity analysis");
491: (*ts->drdpfunction)(ts,t,y,drdp,ts->costintegrandctx);
492: PetscStackPop;
493: return(0);
494: }
496: /*@C
497: TSAdjointMonitorSensi - monitors the first lambda sensitivity
499: Level: intermediate
501: .keywords: TS, set, monitor
503: .seealso: TSAdjointMonitorSet()
504: @*/
505: PetscErrorCode TSAdjointMonitorSensi(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
506: {
508: PetscViewer viewer = vf->viewer;
512: PetscViewerPushFormat(viewer,vf->format);
513: VecView(lambda[0],viewer);
514: PetscViewerPopFormat(viewer);
515: return(0);
516: }
518: /*@C
519: TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
521: Collective on TS
523: Input Parameters:
524: + ts - TS object you wish to monitor
525: . name - the monitor type one is seeking
526: . help - message indicating what monitoring is done
527: . manual - manual page for the monitor
528: . monitor - the monitor function
529: - 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
531: Level: developer
533: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
534: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
535: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
536: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
537: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
538: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
539: PetscOptionsFList(), PetscOptionsEList()
540: @*/
541: 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*))
542: {
543: PetscErrorCode ierr;
544: PetscViewer viewer;
545: PetscViewerFormat format;
546: PetscBool flg;
549: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject) ts)->options,((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
550: if (flg) {
551: PetscViewerAndFormat *vf;
552: PetscViewerAndFormatCreate(viewer,format,&vf);
553: PetscObjectDereference((PetscObject)viewer);
554: if (monitorsetup) {
555: (*monitorsetup)(ts,vf);
556: }
557: TSAdjointMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
558: }
559: return(0);
560: }
562: /*@C
563: TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
564: timestep to display the iteration's progress.
566: Logically Collective on TS
568: Input Parameters:
569: + ts - the TS context obtained from TSCreate()
570: . adjointmonitor - monitoring routine
571: . adjointmctx - [optional] user-defined context for private data for the
572: monitor routine (use NULL if no context is desired)
573: - adjointmonitordestroy - [optional] routine that frees monitor context
574: (may be NULL)
576: Calling sequence of monitor:
577: $ int adjointmonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *adjointmctx)
579: + ts - the TS context
580: . 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
581: been interpolated to)
582: . time - current time
583: . u - current iterate
584: . numcost - number of cost functionos
585: . lambda - sensitivities to initial conditions
586: . mu - sensitivities to parameters
587: - adjointmctx - [optional] adjoint monitoring context
589: Notes:
590: This routine adds an additional monitor to the list of monitors that
591: already has been loaded.
593: Fortran Notes:
594: Only a single monitor function can be set for each TS object
596: Level: intermediate
598: .keywords: TS, timestep, set, adjoint, monitor
600: .seealso: TSAdjointMonitorCancel()
601: @*/
602: PetscErrorCode TSAdjointMonitorSet(TS ts,PetscErrorCode (*adjointmonitor)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*),void *adjointmctx,PetscErrorCode (*adjointmdestroy)(void**))
603: {
605: PetscInt i;
606: PetscBool identical;
610: for (i=0; i<ts->numbermonitors;i++) {
611: PetscMonitorCompare((PetscErrorCode (*)(void))adjointmonitor,adjointmctx,adjointmdestroy,(PetscErrorCode (*)(void))ts->adjointmonitor[i],ts->adjointmonitorcontext[i],ts->adjointmonitordestroy[i],&identical);
612: if (identical) return(0);
613: }
614: if (ts->numberadjointmonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many adjoint monitors set");
615: ts->adjointmonitor[ts->numberadjointmonitors] = adjointmonitor;
616: ts->adjointmonitordestroy[ts->numberadjointmonitors] = adjointmdestroy;
617: ts->adjointmonitorcontext[ts->numberadjointmonitors++] = (void*)adjointmctx;
618: return(0);
619: }
621: /*@C
622: TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.
624: Logically Collective on TS
626: Input Parameters:
627: . ts - the TS context obtained from TSCreate()
629: Notes:
630: There is no way to remove a single, specific monitor.
632: Level: intermediate
634: .keywords: TS, timestep, set, adjoint, monitor
636: .seealso: TSAdjointMonitorSet()
637: @*/
638: PetscErrorCode TSAdjointMonitorCancel(TS ts)
639: {
641: PetscInt i;
645: for (i=0; i<ts->numberadjointmonitors; i++) {
646: if (ts->adjointmonitordestroy[i]) {
647: (*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]);
648: }
649: }
650: ts->numberadjointmonitors = 0;
651: return(0);
652: }
654: /*@C
655: TSAdjointMonitorDefault - the default monitor of adjoint computations
657: Level: intermediate
659: .keywords: TS, set, monitor
661: .seealso: TSAdjointMonitorSet()
662: @*/
663: PetscErrorCode TSAdjointMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
664: {
666: PetscViewer viewer = vf->viewer;
670: PetscViewerPushFormat(viewer,vf->format);
671: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
672: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
673: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
674: PetscViewerPopFormat(viewer);
675: return(0);
676: }
678: /*@C
679: TSAdjointMonitorDrawSensi - Monitors progress of the adjoint TS solvers by calling
680: VecView() for the sensitivities to initial states at each timestep
682: Collective on TS
684: Input Parameters:
685: + ts - the TS context
686: . step - current time-step
687: . ptime - current time
688: . u - current state
689: . numcost - number of cost functions
690: . lambda - sensitivities to initial conditions
691: . mu - sensitivities to parameters
692: - dummy - either a viewer or NULL
694: Level: intermediate
696: .keywords: TS, vector, adjoint, monitor, view
698: .seealso: TSAdjointMonitorSet(), TSAdjointMonitorDefault(), VecView()
699: @*/
700: PetscErrorCode TSAdjointMonitorDrawSensi(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda,Vec *mu,void *dummy)
701: {
702: PetscErrorCode ierr;
703: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
704: PetscDraw draw;
705: PetscReal xl,yl,xr,yr,h;
706: char time[32];
709: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);
711: VecView(lambda[0],ictx->viewer);
712: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
713: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
714: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
715: h = yl + .95*(yr - yl);
716: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
717: PetscDrawFlush(draw);
718: return(0);
719: }
721: /*
722: TSAdjointSetFromOptions - Sets various TSAdjoint parameters from user options.
724: Collective on TSAdjoint
726: Input Parameter:
727: . ts - the TS context
729: Options Database Keys:
730: + -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
731: . -ts_adjoint_monitor - print information at each adjoint time step
732: - -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically
734: Level: developer
736: Notes:
737: This is not normally called directly by users
739: .keywords: TS, trajectory, timestep, set, options, database
741: .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
742: */
743: PetscErrorCode TSAdjointSetFromOptions(PetscOptionItems *PetscOptionsObject,TS ts)
744: {
745: PetscBool tflg,opt;
750: PetscOptionsHead(PetscOptionsObject,"TS Adjoint options");
751: tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
752: PetscOptionsBool("-ts_adjoint_solve","Solve the adjoint problem immediately after solving the forward problem","",tflg,&tflg,&opt);
753: if (opt) {
754: TSSetSaveTrajectory(ts);
755: ts->adjoint_solve = tflg;
756: }
757: TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor","Monitor adjoint timestep size","TSAdjointMonitorDefault",TSAdjointMonitorDefault,NULL);
758: TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor_sensi","Monitor sensitivity in the adjoint computation","TSAdjointMonitorSensi",TSAdjointMonitorSensi,NULL);
759: opt = PETSC_FALSE;
760: PetscOptionsName("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",&opt);
761: if (opt) {
762: TSMonitorDrawCtx ctx;
763: PetscInt howoften = 1;
765: PetscOptionsInt("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",howoften,&howoften,NULL);
766: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
767: TSAdjointMonitorSet(ts,TSAdjointMonitorDrawSensi,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
768: }
769: return(0);
770: }
772: /*@
773: TSAdjointStep - Steps one time step backward in the adjoint run
775: Collective on TS
777: Input Parameter:
778: . ts - the TS context obtained from TSCreate()
780: Level: intermediate
782: .keywords: TS, adjoint, step
784: .seealso: TSAdjointSetUp(), TSAdjointSolve()
785: @*/
786: PetscErrorCode TSAdjointStep(TS ts)
787: {
788: DM dm;
789: PetscErrorCode ierr;
793: TSGetDM(ts,&dm);
794: TSAdjointSetUp(ts);
796: VecViewFromOptions(ts->vec_sol,(PetscObject)ts,"-ts_view_solution");
798: ts->reason = TS_CONVERGED_ITERATING;
799: ts->ptime_prev = ts->ptime;
800: 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);
801: PetscLogEventBegin(TS_AdjointStep,ts,0,0,0);
802: (*ts->ops->adjointstep)(ts);
803: PetscLogEventEnd(TS_AdjointStep,ts,0,0,0);
804: ts->adjoint_steps++; ts->steps--;
806: if (ts->reason < 0) {
807: if (ts->errorifstepfailed) {
808: if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s, increase -ts_max_snes_failures or make negative to attempt recovery",TSConvergedReasons[ts->reason]);
809: else if (ts->reason == TS_DIVERGED_STEP_REJECTED) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s, increase -ts_max_reject or make negative to attempt recovery",TSConvergedReasons[ts->reason]);
810: else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
811: }
812: } else if (!ts->reason) {
813: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
814: }
815: return(0);
816: }
818: /*@
819: TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE
821: Collective on TS
823: Input Parameter:
824: . ts - the TS context obtained from TSCreate()
826: Options Database:
827: . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial values
829: Level: intermediate
831: Notes:
832: This must be called after a call to TSSolve() that solves the forward problem
834: By default this will integrate back to the initial time, one can use TSAdjointSetSteps() to step back to a later time
836: .keywords: TS, timestep, solve
838: .seealso: TSCreate(), TSSetCostGradients(), TSSetSolution(), TSAdjointStep()
839: @*/
840: PetscErrorCode TSAdjointSolve(TS ts)
841: {
842: PetscErrorCode ierr;
846: TSAdjointSetUp(ts);
848: /* reset time step and iteration counters */
849: ts->adjoint_steps = 0;
850: ts->ksp_its = 0;
851: ts->snes_its = 0;
852: ts->num_snes_failures = 0;
853: ts->reject = 0;
854: ts->reason = TS_CONVERGED_ITERATING;
856: if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->steps;
857: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
859: while (!ts->reason) {
860: TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);
861: TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
862: TSAdjointEventHandler(ts);
863: TSAdjointStep(ts);
864: if (ts->vec_costintegral && !ts->costintegralfwd) {
865: TSAdjointCostIntegral(ts);
866: }
867: }
868: TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);
869: TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
870: ts->solvetime = ts->ptime;
871: TSTrajectoryViewFromOptions(ts->trajectory,NULL,"-ts_trajectory_view");
872: VecViewFromOptions(ts->vecs_sensi[0],(PetscObject) ts, "-ts_adjoint_view_solution");
873: ts->adjoint_max_steps = 0;
874: return(0);
875: }
877: /*@C
878: TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using TSAdjointMonitorSet()
880: Collective on TS
882: Input Parameters:
883: + ts - time stepping context obtained from TSCreate()
884: . step - step number that has just completed
885: . ptime - model time of the state
886: . u - state at the current model time
887: . numcost - number of cost functions (dimension of lambda or mu)
888: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
889: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
891: Notes:
892: TSAdjointMonitor() is typically used automatically within the time stepping implementations.
893: Users would almost never call this routine directly.
895: Level: developer
897: .keywords: TS, timestep
898: @*/
899: PetscErrorCode TSAdjointMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda, Vec *mu)
900: {
902: PetscInt i,n = ts->numberadjointmonitors;
907: VecLockReadPush(u);
908: for (i=0; i<n; i++) {
909: (*ts->adjointmonitor[i])(ts,step,ptime,u,numcost,lambda,mu,ts->adjointmonitorcontext[i]);
910: }
911: VecLockReadPop(u);
912: return(0);
913: }
915: /*@
916: TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.
918: Collective on TS
920: Input Arguments:
921: . ts - time stepping context
923: Level: advanced
925: Notes:
926: This function cannot be called until TSAdjointStep() has been completed.
928: .seealso: TSAdjointSolve(), TSAdjointStep
929: @*/
930: PetscErrorCode TSAdjointCostIntegral(TS ts)
931: {
934: 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);
935: (*ts->ops->adjointintegral)(ts);
936: return(0);
937: }
939: /* ------------------ Forward (tangent linear) sensitivity ------------------*/
941: /*@
942: TSForwardSetUp - Sets up the internal data structures for the later use
943: of forward sensitivity analysis
945: Collective on TS
947: Input Parameter:
948: . ts - the TS context obtained from TSCreate()
950: Level: advanced
952: .keywords: TS, forward sensitivity, setup
954: .seealso: TSCreate(), TSDestroy(), TSSetUp()
955: @*/
956: PetscErrorCode TSForwardSetUp(TS ts)
957: {
962: if (ts->forwardsetupcalled) return(0);
963: if (ts->vec_costintegral && !ts->vecs_integral_sensip ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSForwardSetIntegralGradients() before TSSetCostIntegrand()");
964: if (ts->vecs_integral_sensip) {
965: VecDuplicateVecs(ts->vec_sol,ts->numcost,&ts->vecs_drdy);
966: VecDuplicateVecs(ts->vecs_integral_sensip[0],ts->numcost,&ts->vecs_drdp);
967: }
969: if (ts->ops->forwardsetup) {
970: (*ts->ops->forwardsetup)(ts);
971: }
972: ts->forwardsetupcalled = PETSC_TRUE;
973: return(0);
974: }
976: /*@
977: TSForwardSetIntegralGradients - Set the vectors holding forward sensitivities of the integral term.
979: Input Parameter:
980: . ts- the TS context obtained from TSCreate()
981: . numfwdint- number of integrals
982: . vp = the vectors containing the gradients for each integral w.r.t. parameters
984: Level: intermediate
986: .keywords: TS, forward sensitivity
988: .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
989: @*/
990: PetscErrorCode TSForwardSetIntegralGradients(TS ts,PetscInt numfwdint,Vec *vp)
991: {
994: 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()");
995: if (!ts->numcost) ts->numcost = numfwdint;
997: ts->vecs_integral_sensip = vp;
998: return(0);
999: }
1001: /*@
1002: TSForwardGetIntegralGradients - Returns the forward sensitivities ofthe integral term.
1004: Input Parameter:
1005: . ts- the TS context obtained from TSCreate()
1007: Output Parameter:
1008: . vp = the vectors containing the gradients for each integral w.r.t. parameters
1010: Level: intermediate
1012: .keywords: TS, forward sensitivity
1014: .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1015: @*/
1016: PetscErrorCode TSForwardGetIntegralGradients(TS ts,PetscInt *numfwdint,Vec **vp)
1017: {
1021: if (numfwdint) *numfwdint = ts->numcost;
1022: if (vp) *vp = ts->vecs_integral_sensip;
1023: return(0);
1024: }
1026: /*@
1027: TSForwardStep - Compute the forward sensitivity for one time step.
1029: Collective on TS
1031: Input Arguments:
1032: . ts - time stepping context
1034: Level: advanced
1036: Notes:
1037: This function cannot be called until TSStep() has been completed.
1039: .keywords: TS, forward sensitivity
1041: .seealso: TSForwardSetSensitivities(), TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardSetUp()
1042: @*/
1043: PetscErrorCode TSForwardStep(TS ts)
1044: {
1047: if (!ts->ops->forwardstep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide forward sensitivity analysis",((PetscObject)ts)->type_name);
1048: PetscLogEventBegin(TS_ForwardStep,ts,0,0,0);
1049: (*ts->ops->forwardstep)(ts);
1050: PetscLogEventEnd(TS_ForwardStep,ts,0,0,0);
1051: return(0);
1052: }
1054: /*@
1055: TSForwardSetSensitivities - Sets the initial value of the trajectory sensitivities of solution w.r.t. the problem parameters and initial values.
1057: Logically Collective on TS and Vec
1059: Input Parameters:
1060: + ts - the TS context obtained from TSCreate()
1061: . nump - number of parameters
1062: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1064: Level: beginner
1066: Notes:
1067: Forward sensitivity is also called 'trajectory sensitivity' in some fields such as power systems.
1068: This function turns on a flag to trigger TSSolve() to compute forward sensitivities automatically.
1069: You must call this function before TSSolve().
1070: The entries in the sensitivity matrix must be correctly initialized with the values S = dy/dp|startingtime.
1072: .keywords: TS, timestep, set, forward sensitivity, initial values
1074: .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1075: @*/
1076: PetscErrorCode TSForwardSetSensitivities(TS ts,PetscInt nump,Mat Smat)
1077: {
1083: ts->forward_solve = PETSC_TRUE;
1084: ts->num_parameters = nump;
1085: PetscObjectReference((PetscObject)Smat);
1086: MatDestroy(&ts->mat_sensip);
1087: ts->mat_sensip = Smat;
1088: return(0);
1089: }
1091: /*@
1092: TSForwardGetSensitivities - Returns the trajectory sensitivities
1094: Not Collective, but Vec returned is parallel if TS is parallel
1096: Output Parameter:
1097: + ts - the TS context obtained from TSCreate()
1098: . nump - number of parameters
1099: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1101: Level: intermediate
1103: .keywords: TS, forward sensitivity
1105: .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1106: @*/
1107: PetscErrorCode TSForwardGetSensitivities(TS ts,PetscInt *nump,Mat *Smat)
1108: {
1111: if (nump) *nump = ts->num_parameters;
1112: if (Smat) *Smat = ts->mat_sensip;
1113: return(0);
1114: }
1116: /*@
1117: TSForwardCostIntegral - Evaluate the cost integral in the forward run.
1119: Collective on TS
1121: Input Arguments:
1122: . ts - time stepping context
1124: Level: advanced
1126: Notes:
1127: This function cannot be called until TSStep() has been completed.
1129: .seealso: TSSolve(), TSAdjointCostIntegral()
1130: @*/
1131: PetscErrorCode TSForwardCostIntegral(TS ts)
1132: {
1135: 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);
1136: (*ts->ops->forwardintegral)(ts);
1137: return(0);
1138: }