Actual source code: tssen.c

petsc-3.11.2 2019-05-18
Report Typos and Errors
  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: }