Actual source code: tssen.c

petsc-3.9.4 2018-09-11
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: 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: }