Actual source code: tssen.c
petsc-3.9.4 2018-09-11
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: 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
28: .keywords: TS, sensitivity
29: .seealso:
30: @*/
31: PetscErrorCode TSSetRHSJacobianP(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx)
32: {
39: ts->rhsjacobianp = func;
40: ts->rhsjacobianpctx = ctx;
41: if(Amat) {
42: PetscObjectReference((PetscObject)Amat);
43: MatDestroy(&ts->Jacp);
44: ts->Jacp = Amat;
45: }
46: return(0);
47: }
49: /*@C
50: TSComputeRHSJacobianP - Runs the user-defined JacobianP function.
52: Collective on TS
54: Input Parameters:
55: . ts - The TS context obtained from TSCreate()
57: Level: developer
59: .keywords: TS, sensitivity
60: .seealso: TSSetRHSJacobianP()
61: @*/
62: PetscErrorCode TSComputeRHSJacobianP(TS ts,PetscReal t,Vec X,Mat Amat)
63: {
71: PetscStackPush("TS user JacobianP function for sensitivity analysis");
72: (*ts->rhsjacobianp)(ts,t,X,Amat,ts->rhsjacobianpctx);
73: PetscStackPop;
74: return(0);
75: }
77: /*@C
78: TSSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions
80: Logically Collective on TS
82: Input Parameters:
83: + ts - the TS context obtained from TSCreate()
84: . numcost - number of gradients to be computed, this is the number of cost functions
85: . costintegral - vector that stores the integral values
86: . rf - routine for evaluating the integrand function
87: . drdyf - function that computes the gradients of the r's with respect to y
88: . 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)
89: . fwd - flag indicating whether to evaluate cost integral in the forward run or the adjoint run
90: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
92: Calling sequence of rf:
93: $ PetscErrorCode rf(TS ts,PetscReal t,Vec y,Vec f,void *ctx);
95: Calling sequence of drdyf:
96: $ PetscErroCode drdyf(TS ts,PetscReal t,Vec y,Vec *drdy,void *ctx);
98: Calling sequence of drdpf:
99: $ PetscErroCode drdpf(TS ts,PetscReal t,Vec y,Vec *drdp,void *ctx);
101: Level: intermediate
103: Notes: For optimization there is usually a single cost function (numcost = 1). For sensitivities there may be multiple cost functions
105: .keywords: TS, sensitivity analysis, timestep, set, quadrature, function
107: .seealso: TSSetRHSJacobianP(), TSGetCostGradients(), TSSetCostGradients()
108: @*/
109: PetscErrorCode TSSetCostIntegrand(TS ts,PetscInt numcost,Vec costintegral,PetscErrorCode (*rf)(TS,PetscReal,Vec,Vec,void*),
110: PetscErrorCode (*drdyf)(TS,PetscReal,Vec,Vec*,void*),
111: PetscErrorCode (*drdpf)(TS,PetscReal,Vec,Vec*,void*),
112: PetscBool fwd,void *ctx)
113: {
119: 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()");
120: if (!ts->numcost) ts->numcost=numcost;
122: if (costintegral) {
123: PetscObjectReference((PetscObject)costintegral);
124: VecDestroy(&ts->vec_costintegral);
125: ts->vec_costintegral = costintegral;
126: } else {
127: if (!ts->vec_costintegral) { /* Create a seq vec if user does not provide one */
128: VecCreateSeq(PETSC_COMM_SELF,numcost,&ts->vec_costintegral);
129: } else {
130: VecSet(ts->vec_costintegral,0.0);
131: }
132: }
133: if (!ts->vec_costintegrand) {
134: VecDuplicate(ts->vec_costintegral,&ts->vec_costintegrand);
135: } else {
136: VecSet(ts->vec_costintegrand,0.0);
137: }
138: ts->costintegralfwd = fwd; /* Evaluate the cost integral in forward run if fwd is true */
139: ts->costintegrand = rf;
140: ts->costintegrandctx = ctx;
141: ts->drdyfunction = drdyf;
142: ts->drdpfunction = drdpf;
143: return(0);
144: }
146: /*@
147: TSGetCostIntegral - Returns the values of the integral term in the cost functions.
148: It is valid to call the routine after a backward run.
150: Not Collective
152: Input Parameter:
153: . ts - the TS context obtained from TSCreate()
155: Output Parameter:
156: . v - the vector containing the integrals for each cost function
158: Level: intermediate
160: .seealso: TSSetCostIntegrand()
162: .keywords: TS, sensitivity analysis
163: @*/
164: PetscErrorCode TSGetCostIntegral(TS ts,Vec *v)
165: {
169: *v = ts->vec_costintegral;
170: return(0);
171: }
173: /*@
174: TSComputeCostIntegrand - Evaluates the integral function in the cost functions.
176: Input Parameters:
177: + ts - the TS context
178: . t - current time
179: - y - state vector, i.e. current solution
181: Output Parameter:
182: . q - vector of size numcost to hold the outputs
184: Note:
185: Most users should not need to explicitly call this routine, as it
186: is used internally within the sensitivity analysis context.
188: Level: developer
190: .keywords: TS, compute
192: .seealso: TSSetCostIntegrand()
193: @*/
194: PetscErrorCode TSComputeCostIntegrand(TS ts,PetscReal t,Vec y,Vec q)
195: {
203: PetscLogEventBegin(TS_FunctionEval,ts,y,q,0);
204: if (ts->costintegrand) {
205: PetscStackPush("TS user integrand in the cost function");
206: (*ts->costintegrand)(ts,t,y,q,ts->costintegrandctx);
207: PetscStackPop;
208: } else {
209: VecZeroEntries(q);
210: }
212: PetscLogEventEnd(TS_FunctionEval,ts,y,q,0);
213: return(0);
214: }
216: /*@
217: TSComputeDRDYFunction - Runs the user-defined DRDY function.
219: Collective on TS
221: Input Parameters:
222: . ts - The TS context obtained from TSCreate()
224: Notes:
225: TSAdjointComputeDRDYFunction() is typically used for sensitivity implementation,
226: so most users would not generally call this routine themselves.
228: Level: developer
230: .keywords: TS, sensitivity
231: .seealso: TSSetCostIntegrand()
232: @*/
233: PetscErrorCode TSComputeDRDYFunction(TS ts,PetscReal t,Vec y,Vec *drdy)
234: {
241: PetscStackPush("TS user DRDY function for sensitivity analysis");
242: (*ts->drdyfunction)(ts,t,y,drdy,ts->costintegrandctx);
243: PetscStackPop;
244: return(0);
245: }
247: /*@
248: TSComputeDRDPFunction - Runs the user-defined DRDP function.
250: Collective on TS
252: Input Parameters:
253: . ts - The TS context obtained from TSCreate()
255: Notes:
256: TSDRDPFunction() is typically used for sensitivity implementation,
257: so most users would not generally call this routine themselves.
259: Level: developer
261: .keywords: TS, sensitivity
262: .seealso: TSSetCostIntegrand()
263: @*/
264: PetscErrorCode TSComputeDRDPFunction(TS ts,PetscReal t,Vec y,Vec *drdp)
265: {
272: PetscStackPush("TS user DRDP function for sensitivity analysis");
273: (*ts->drdpfunction)(ts,t,y,drdp,ts->costintegrandctx);
274: PetscStackPop;
275: return(0);
276: }
278: /* --------------------------- Adjoint sensitivity ---------------------------*/
280: /*@
281: TSSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial values and w.r.t. the problem parameters
282: for use by the TSAdjoint routines.
284: Logically Collective on TS and Vec
286: Input Parameters:
287: + ts - the TS context obtained from TSCreate()
288: . 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
289: - mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
291: Level: beginner
293: Notes: the entries in these vectors must be correctly initialized with the values lamda_i = df/dy|finaltime mu_i = df/dp|finaltime
295: After TSAdjointSolve() is called the lamba and the mu contain the computed sensitivities
297: .keywords: TS, timestep, set, sensitivity, initial values
298: @*/
299: PetscErrorCode TSSetCostGradients(TS ts,PetscInt numcost,Vec *lambda,Vec *mu)
300: {
304: ts->vecs_sensi = lambda;
305: ts->vecs_sensip = mu;
306: 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");
307: ts->numcost = numcost;
308: return(0);
309: }
311: /*@
312: TSGetCostGradients - Returns the gradients from the TSAdjointSolve()
314: Not Collective, but Vec returned is parallel if TS is parallel
316: Input Parameter:
317: . ts - the TS context obtained from TSCreate()
319: Output Parameter:
320: + lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
321: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
323: Level: intermediate
325: .seealso: TSGetTimeStep()
327: .keywords: TS, timestep, get, sensitivity
328: @*/
329: PetscErrorCode TSGetCostGradients(TS ts,PetscInt *numcost,Vec **lambda,Vec **mu)
330: {
333: if (numcost) *numcost = ts->numcost;
334: if (lambda) *lambda = ts->vecs_sensi;
335: if (mu) *mu = ts->vecs_sensip;
336: return(0);
337: }
339: /*@
340: TSAdjointSetUp - Sets up the internal data structures for the later use
341: of an adjoint solver
343: Collective on TS
345: Input Parameter:
346: . ts - the TS context obtained from TSCreate()
348: Level: advanced
350: .keywords: TS, timestep, setup
352: .seealso: TSCreate(), TSAdjointStep(), TSSetCostGradients()
353: @*/
354: PetscErrorCode TSAdjointSetUp(TS ts)
355: {
360: if (ts->adjointsetupcalled) return(0);
361: if (!ts->vecs_sensi) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetCostGradients() first");
362: if (ts->vecs_sensip && !ts->Jacp) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must call TSAdjointSetRHSJacobian() first");
364: if (ts->vec_costintegral) { /* if there is integral in the cost function */
365: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&ts->vecs_drdy);
366: if (ts->vecs_sensip){
367: VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&ts->vecs_drdp);
368: }
369: }
371: if (ts->ops->adjointsetup) {
372: (*ts->ops->adjointsetup)(ts);
373: }
374: ts->adjointsetupcalled = PETSC_TRUE;
375: return(0);
376: }
378: /*@
379: TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time
381: Logically Collective on TS
383: Input Parameters:
384: + ts - the TS context obtained from TSCreate()
385: . steps - number of steps to use
387: Level: intermediate
389: Notes: Normally one does not call this and TSAdjointSolve() integrates back to the original timestep. One can call this
390: so as to integrate back to less than the original timestep
392: .keywords: TS, timestep, set, maximum, iterations
394: .seealso: TSSetExactFinalTime()
395: @*/
396: PetscErrorCode TSAdjointSetSteps(TS ts,PetscInt steps)
397: {
401: if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back a negative number of steps");
402: if (steps > ts->steps) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back more than the total number of forward steps");
403: ts->adjoint_max_steps = steps;
404: return(0);
405: }
407: /*@C
408: TSAdjointSetRHSJacobian - Deprecated, use TSSetRHSJacobianP()
410: Level: deprecated
412: @*/
413: PetscErrorCode TSAdjointSetRHSJacobian(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx)
414: {
421: ts->rhsjacobianp = func;
422: ts->rhsjacobianpctx = ctx;
423: if(Amat) {
424: PetscObjectReference((PetscObject)Amat);
425: MatDestroy(&ts->Jacp);
426: ts->Jacp = Amat;
427: }
428: return(0);
429: }
431: /*@C
432: TSAdjointComputeRHSJacobian - Deprecated, use TSComputeRHSJacobianP()
434: Level: deprecated
436: @*/
437: PetscErrorCode TSAdjointComputeRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat)
438: {
446: PetscStackPush("TS user JacobianP function for sensitivity analysis");
447: (*ts->rhsjacobianp)(ts,t,X,Amat,ts->rhsjacobianpctx);
448: PetscStackPop;
449: return(0);
450: }
452: /*@
453: TSAdjointComputeDRDYFunction - Deprecated, use TSComputeDRDYFunction()
455: Level: deprecated
457: @*/
458: PetscErrorCode TSAdjointComputeDRDYFunction(TS ts,PetscReal t,Vec y,Vec *drdy)
459: {
466: PetscStackPush("TS user DRDY function for sensitivity analysis");
467: (*ts->drdyfunction)(ts,t,y,drdy,ts->costintegrandctx);
468: PetscStackPop;
469: return(0);
470: }
472: /*@
473: TSAdjointComputeDRDPFunction - Deprecated, use TSComputeDRDPFunction()
475: Level: deprecated
477: @*/
478: PetscErrorCode TSAdjointComputeDRDPFunction(TS ts,PetscReal t,Vec y,Vec *drdp)
479: {
486: PetscStackPush("TS user DRDP function for sensitivity analysis");
487: (*ts->drdpfunction)(ts,t,y,drdp,ts->costintegrandctx);
488: PetscStackPop;
489: return(0);
490: }
492: /*@C
493: TSAdjointMonitorSensi - monitors the first lambda sensitivity
495: Level: intermediate
497: .keywords: TS, set, monitor
499: .seealso: TSAdjointMonitorSet()
500: @*/
501: PetscErrorCode TSAdjointMonitorSensi(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
502: {
504: PetscViewer viewer = vf->viewer;
508: PetscViewerPushFormat(viewer,vf->format);
509: VecView(lambda[0],viewer);
510: PetscViewerPopFormat(viewer);
511: return(0);
512: }
514: /*@C
515: TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user
517: Collective on TS
519: Input Parameters:
520: + ts - TS object you wish to monitor
521: . name - the monitor type one is seeking
522: . help - message indicating what monitoring is done
523: . manual - manual page for the monitor
524: . monitor - the monitor function
525: - 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
527: Level: developer
529: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
530: PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
531: PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
532: PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
533: PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
534: PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
535: PetscOptionsFList(), PetscOptionsEList()
536: @*/
537: 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*))
538: {
539: PetscErrorCode ierr;
540: PetscViewer viewer;
541: PetscViewerFormat format;
542: PetscBool flg;
545: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
546: if (flg) {
547: PetscViewerAndFormat *vf;
548: PetscViewerAndFormatCreate(viewer,format,&vf);
549: PetscObjectDereference((PetscObject)viewer);
550: if (monitorsetup) {
551: (*monitorsetup)(ts,vf);
552: }
553: TSAdjointMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
554: }
555: return(0);
556: }
558: /*@C
559: TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
560: timestep to display the iteration's progress.
562: Logically Collective on TS
564: Input Parameters:
565: + ts - the TS context obtained from TSCreate()
566: . adjointmonitor - monitoring routine
567: . adjointmctx - [optional] user-defined context for private data for the
568: monitor routine (use NULL if no context is desired)
569: - adjointmonitordestroy - [optional] routine that frees monitor context
570: (may be NULL)
572: Calling sequence of monitor:
573: $ int adjointmonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *adjointmctx)
575: + ts - the TS context
576: . 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
577: been interpolated to)
578: . time - current time
579: . u - current iterate
580: . numcost - number of cost functionos
581: . lambda - sensitivities to initial conditions
582: . mu - sensitivities to parameters
583: - adjointmctx - [optional] adjoint monitoring context
585: Notes:
586: This routine adds an additional monitor to the list of monitors that
587: already has been loaded.
589: Fortran notes: Only a single monitor function can be set for each TS object
591: Level: intermediate
593: .keywords: TS, timestep, set, adjoint, monitor
595: .seealso: TSAdjointMonitorCancel()
596: @*/
597: PetscErrorCode TSAdjointMonitorSet(TS ts,PetscErrorCode (*adjointmonitor)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*),void *adjointmctx,PetscErrorCode (*adjointmdestroy)(void**))
598: {
600: PetscInt i;
601: PetscBool identical;
605: for (i=0; i<ts->numbermonitors;i++) {
606: PetscMonitorCompare((PetscErrorCode (*)(void))adjointmonitor,adjointmctx,adjointmdestroy,(PetscErrorCode (*)(void))ts->adjointmonitor[i],ts->adjointmonitorcontext[i],ts->adjointmonitordestroy[i],&identical);
607: if (identical) return(0);
608: }
609: if (ts->numberadjointmonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many adjoint monitors set");
610: ts->adjointmonitor[ts->numberadjointmonitors] = adjointmonitor;
611: ts->adjointmonitordestroy[ts->numberadjointmonitors] = adjointmdestroy;
612: ts->adjointmonitorcontext[ts->numberadjointmonitors++] = (void*)adjointmctx;
613: return(0);
614: }
616: /*@C
617: TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.
619: Logically Collective on TS
621: Input Parameters:
622: . ts - the TS context obtained from TSCreate()
624: Notes:
625: There is no way to remove a single, specific monitor.
627: Level: intermediate
629: .keywords: TS, timestep, set, adjoint, monitor
631: .seealso: TSAdjointMonitorSet()
632: @*/
633: PetscErrorCode TSAdjointMonitorCancel(TS ts)
634: {
636: PetscInt i;
640: for (i=0; i<ts->numberadjointmonitors; i++) {
641: if (ts->adjointmonitordestroy[i]) {
642: (*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]);
643: }
644: }
645: ts->numberadjointmonitors = 0;
646: return(0);
647: }
649: /*@C
650: TSAdjointMonitorDefault - the default monitor of adjoint computations
652: Level: intermediate
654: .keywords: TS, set, monitor
656: .seealso: TSAdjointMonitorSet()
657: @*/
658: PetscErrorCode TSAdjointMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
659: {
661: PetscViewer viewer = vf->viewer;
665: PetscViewerPushFormat(viewer,vf->format);
666: PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
667: PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
668: PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
669: PetscViewerPopFormat(viewer);
670: return(0);
671: }
673: /*@C
674: TSAdjointMonitorDrawSensi - Monitors progress of the adjoint TS solvers by calling
675: VecView() for the sensitivities to initial states at each timestep
677: Collective on TS
679: Input Parameters:
680: + ts - the TS context
681: . step - current time-step
682: . ptime - current time
683: . u - current state
684: . numcost - number of cost functions
685: . lambda - sensitivities to initial conditions
686: . mu - sensitivities to parameters
687: - dummy - either a viewer or NULL
689: Level: intermediate
691: .keywords: TS, vector, adjoint, monitor, view
693: .seealso: TSAdjointMonitorSet(), TSAdjointMonitorDefault(), VecView()
694: @*/
695: PetscErrorCode TSAdjointMonitorDrawSensi(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda,Vec *mu,void *dummy)
696: {
697: PetscErrorCode ierr;
698: TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
699: PetscDraw draw;
700: PetscReal xl,yl,xr,yr,h;
701: char time[32];
704: if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);
706: VecView(lambda[0],ictx->viewer);
707: PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
708: PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
709: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
710: h = yl + .95*(yr - yl);
711: PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
712: PetscDrawFlush(draw);
713: return(0);
714: }
716: /*
717: TSAdjointSetFromOptions - Sets various TSAdjoint parameters from user options.
719: Collective on TSAdjoint
721: Input Parameter:
722: . ts - the TS context
724: Options Database Keys:
725: + -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
726: . -ts_adjoint_monitor - print information at each adjoint time step
727: - -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically
729: Level: developer
731: Notes: This is not normally called directly by users
733: .keywords: TS, trajectory, timestep, set, options, database
735: .seealso: TSSetSaveTrajectory(), TSTrajectorySetUp()
736: */
737: PetscErrorCode TSAdjointSetFromOptions(PetscOptionItems *PetscOptionsObject,TS ts)
738: {
739: PetscBool tflg,opt;
744: PetscOptionsHead(PetscOptionsObject,"TS Adjoint options");
745: tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
746: PetscOptionsBool("-ts_adjoint_solve","Solve the adjoint problem immediately after solving the forward problem","",tflg,&tflg,&opt);
747: if (opt) {
748: TSSetSaveTrajectory(ts);
749: ts->adjoint_solve = tflg;
750: }
751: TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor","Monitor adjoint timestep size","TSAdjointMonitorDefault",TSAdjointMonitorDefault,NULL);
752: TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor_sensi","Monitor sensitivity in the adjoint computation","TSAdjointMonitorSensi",TSAdjointMonitorSensi,NULL);
753: opt = PETSC_FALSE;
754: PetscOptionsName("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",&opt);
755: if (opt) {
756: TSMonitorDrawCtx ctx;
757: PetscInt howoften = 1;
759: PetscOptionsInt("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",howoften,&howoften,NULL);
760: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
761: TSAdjointMonitorSet(ts,TSAdjointMonitorDrawSensi,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
762: }
763: return(0);
764: }
766: /*@
767: TSAdjointStep - Steps one time step backward in the adjoint run
769: Collective on TS
771: Input Parameter:
772: . ts - the TS context obtained from TSCreate()
774: Level: intermediate
776: .keywords: TS, adjoint, step
778: .seealso: TSAdjointSetUp(), TSAdjointSolve()
779: @*/
780: PetscErrorCode TSAdjointStep(TS ts)
781: {
782: DM dm;
783: PetscErrorCode ierr;
787: TSGetDM(ts,&dm);
788: TSAdjointSetUp(ts);
790: VecViewFromOptions(ts->vec_sol,(PetscObject)ts,"-ts_view_solution");
792: ts->reason = TS_CONVERGED_ITERATING;
793: ts->ptime_prev = ts->ptime;
794: 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);
795: PetscLogEventBegin(TS_AdjointStep,ts,0,0,0);
796: (*ts->ops->adjointstep)(ts);
797: PetscLogEventEnd(TS_AdjointStep,ts,0,0,0);
798: ts->adjoint_steps++; ts->steps--;
800: if (ts->reason < 0) {
801: if (ts->errorifstepfailed) {
802: 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]);
803: 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]);
804: else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
805: }
806: } else if (!ts->reason) {
807: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
808: }
809: return(0);
810: }
812: /*@
813: TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE
815: Collective on TS
817: Input Parameter:
818: . ts - the TS context obtained from TSCreate()
820: Options Database:
821: . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial values
823: Level: intermediate
825: Notes:
826: This must be called after a call to TSSolve() that solves the forward problem
828: By default this will integrate back to the initial time, one can use TSAdjointSetSteps() to step back to a later time
830: .keywords: TS, timestep, solve
832: .seealso: TSCreate(), TSSetCostGradients(), TSSetSolution(), TSAdjointStep()
833: @*/
834: PetscErrorCode TSAdjointSolve(TS ts)
835: {
836: PetscErrorCode ierr;
840: TSAdjointSetUp(ts);
842: /* reset time step and iteration counters */
843: ts->adjoint_steps = 0;
844: ts->ksp_its = 0;
845: ts->snes_its = 0;
846: ts->num_snes_failures = 0;
847: ts->reject = 0;
848: ts->reason = TS_CONVERGED_ITERATING;
850: if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->steps;
851: if (ts->adjoint_steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
853: while (!ts->reason) {
854: TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);
855: TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
856: TSAdjointEventHandler(ts);
857: TSAdjointStep(ts);
858: if (ts->vec_costintegral && !ts->costintegralfwd) {
859: TSAdjointCostIntegral(ts);
860: }
861: }
862: TSTrajectoryGet(ts->trajectory,ts,ts->steps,&ts->ptime);
863: TSAdjointMonitor(ts,ts->steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
864: ts->solvetime = ts->ptime;
865: TSTrajectoryViewFromOptions(ts->trajectory,NULL,"-ts_trajectory_view");
866: VecViewFromOptions(ts->vecs_sensi[0],(PetscObject) ts, "-ts_adjoint_view_solution");
867: ts->adjoint_max_steps = 0;
868: return(0);
869: }
871: /*@C
872: TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using TSAdjointMonitorSet()
874: Collective on TS
876: Input Parameters:
877: + ts - time stepping context obtained from TSCreate()
878: . step - step number that has just completed
879: . ptime - model time of the state
880: . u - state at the current model time
881: . numcost - number of cost functions (dimension of lambda or mu)
882: . lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
883: - mu - vectors containing the gradients of the cost functions with respect to the problem parameters
885: Notes:
886: TSAdjointMonitor() is typically used automatically within the time stepping implementations.
887: Users would almost never call this routine directly.
889: Level: developer
891: .keywords: TS, timestep
892: @*/
893: PetscErrorCode TSAdjointMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda, Vec *mu)
894: {
896: PetscInt i,n = ts->numberadjointmonitors;
901: VecLockPush(u);
902: for (i=0; i<n; i++) {
903: (*ts->adjointmonitor[i])(ts,step,ptime,u,numcost,lambda,mu,ts->adjointmonitorcontext[i]);
904: }
905: VecLockPop(u);
906: return(0);
907: }
909: /*@
910: TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.
912: Collective on TS
914: Input Arguments:
915: . ts - time stepping context
917: Level: advanced
919: Notes:
920: This function cannot be called until TSAdjointStep() has been completed.
922: .seealso: TSAdjointSolve(), TSAdjointStep
923: @*/
924: PetscErrorCode TSAdjointCostIntegral(TS ts)
925: {
928: 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);
929: (*ts->ops->adjointintegral)(ts);
930: return(0);
931: }
933: /* ------------------ Forward (tangent linear) sensitivity ------------------*/
935: /*@
936: TSForwardSetUp - Sets up the internal data structures for the later use
937: of forward sensitivity analysis
939: Collective on TS
941: Input Parameter:
942: . ts - the TS context obtained from TSCreate()
944: Level: advanced
946: .keywords: TS, forward sensitivity, setup
948: .seealso: TSCreate(), TSDestroy(), TSSetUp()
949: @*/
950: PetscErrorCode TSForwardSetUp(TS ts)
951: {
956: if (ts->forwardsetupcalled) return(0);
957: if (ts->vec_costintegral && !ts->vecs_integral_sensip ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSForwardSetIntegralGradients() before TSSetCostIntegrand()");
958: if (ts->vecs_integral_sensip) {
959: VecDuplicateVecs(ts->vec_sol,ts->numcost,&ts->vecs_drdy);
960: VecDuplicateVecs(ts->vecs_integral_sensip[0],ts->numcost,&ts->vecs_drdp);
961: }
963: if (ts->ops->forwardsetup) {
964: (*ts->ops->forwardsetup)(ts);
965: }
966: ts->forwardsetupcalled = PETSC_TRUE;
967: return(0);
968: }
970: /*@
971: TSForwardSetIntegralGradients - Set the vectors holding forward sensitivities of the integral term.
973: Input Parameter:
974: . ts- the TS context obtained from TSCreate()
975: . numfwdint- number of integrals
976: . vp = the vectors containing the gradients for each integral w.r.t. parameters
978: Level: intermediate
980: .keywords: TS, forward sensitivity
982: .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
983: @*/
984: PetscErrorCode TSForwardSetIntegralGradients(TS ts,PetscInt numfwdint,Vec *vp)
985: {
988: 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()");
989: if (!ts->numcost) ts->numcost = numfwdint;
991: ts->vecs_integral_sensip = vp;
992: return(0);
993: }
995: /*@
996: TSForwardGetIntegralGradients - Returns the forward sensitivities ofthe integral term.
998: Input Parameter:
999: . ts- the TS context obtained from TSCreate()
1001: Output Parameter:
1002: . vp = the vectors containing the gradients for each integral w.r.t. parameters
1004: Level: intermediate
1006: .keywords: TS, forward sensitivity
1008: .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1009: @*/
1010: PetscErrorCode TSForwardGetIntegralGradients(TS ts,PetscInt *numfwdint,Vec **vp)
1011: {
1015: if (numfwdint) *numfwdint = ts->numcost;
1016: if (vp) *vp = ts->vecs_integral_sensip;
1017: return(0);
1018: }
1020: /*@
1021: TSForwardStep - Compute the forward sensitivity for one time step.
1023: Collective on TS
1025: Input Arguments:
1026: . ts - time stepping context
1028: Level: advanced
1030: Notes:
1031: This function cannot be called until TSStep() has been completed.
1033: .keywords: TS, forward sensitivity
1035: .seealso: TSForwardSetSensitivities(), TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardSetUp()
1036: @*/
1037: PetscErrorCode TSForwardStep(TS ts)
1038: {
1041: if (!ts->ops->forwardstep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide forward sensitivity analysis",((PetscObject)ts)->type_name);
1042: PetscLogEventBegin(TS_ForwardStep,ts,0,0,0);
1043: (*ts->ops->forwardstep)(ts);
1044: PetscLogEventEnd(TS_ForwardStep,ts,0,0,0);
1045: return(0);
1046: }
1048: /*@
1049: TSForwardSetSensitivities - Sets the initial value of the trajectory sensitivities of solution w.r.t. the problem parameters and initial values.
1051: Logically Collective on TS and Vec
1053: Input Parameters:
1054: + ts - the TS context obtained from TSCreate()
1055: . nump - number of parameters
1056: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1058: Level: beginner
1060: Notes:
1061: Forward sensitivity is also called 'trajectory sensitivity' in some fields such as power systems.
1062: This function turns on a flag to trigger TSSolve() to compute forward sensitivities automatically.
1063: You must call this function before TSSolve().
1064: The entries in the sensitivity matrix must be correctly initialized with the values S = dy/dp|startingtime.
1066: .keywords: TS, timestep, set, forward sensitivity, initial values
1068: .seealso: TSForwardGetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1069: @*/
1070: PetscErrorCode TSForwardSetSensitivities(TS ts,PetscInt nump,Mat Smat)
1071: {
1077: ts->forward_solve = PETSC_TRUE;
1078: ts->num_parameters = nump;
1079: PetscObjectReference((PetscObject)Smat);
1080: MatDestroy(&ts->mat_sensip);
1081: ts->mat_sensip = Smat;
1082: return(0);
1083: }
1085: /*@
1086: TSForwardGetSensitivities - Returns the trajectory sensitivities
1088: Not Collective, but Vec returned is parallel if TS is parallel
1090: Output Parameter:
1091: + ts - the TS context obtained from TSCreate()
1092: . nump - number of parameters
1093: - Smat - sensitivities with respect to the parameters, the number of entries in these vectors is the same as the number of parameters
1095: Level: intermediate
1097: .keywords: TS, forward sensitivity
1099: .seealso: TSForwardSetSensitivities(), TSForwardSetIntegralGradients(), TSForwardGetIntegralGradients(), TSForwardStep()
1100: @*/
1101: PetscErrorCode TSForwardGetSensitivities(TS ts,PetscInt *nump,Mat *Smat)
1102: {
1105: if (nump) *nump = ts->num_parameters;
1106: if (Smat) *Smat = ts->mat_sensip;
1107: return(0);
1108: }
1110: /*@
1111: TSForwardCostIntegral - Evaluate the cost integral in the forward run.
1113: Collective on TS
1115: Input Arguments:
1116: . ts - time stepping context
1118: Level: advanced
1120: Notes:
1121: This function cannot be called until TSStep() has been completed.
1123: .seealso: TSSolve(), TSAdjointCostIntegral()
1124: @*/
1125: PetscErrorCode TSForwardCostIntegral(TS ts)
1126: {
1129: 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);
1130: (*ts->ops->forwardintegral)(ts);
1131: return(0);
1132: }