Actual source code: ts.c

petsc-3.7.6 2017-04-24
Report Typos and Errors
  2: #include <petsc/private/tsimpl.h>        /*I "petscts.h"  I*/
  3: #include <petscdmshell.h>
  4: #include <petscdmda.h>
  5: #include <petscviewer.h>
  6: #include <petscdraw.h>

  8: /* Logging support */
  9: PetscClassId  TS_CLASSID, DMTS_CLASSID;
 10: PetscLogEvent TS_AdjointStep, TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;

 12: const char *const TSExactFinalTimeOptions[] = {"UNSPECIFIED","STEPOVER","INTERPOLATE","MATCHSTEP","TSExactFinalTimeOption","TS_EXACTFINALTIME_",0};

 14: struct _n_TSMonitorDrawCtx {
 15:   PetscViewer   viewer;
 16:   Vec           initialsolution;
 17:   PetscBool     showinitial;
 18:   PetscInt      howoften;  /* when > 0 uses step % howoften, when negative only final solution plotted */
 19:   PetscBool     showtimestepandtime;
 20: };

 24: /*@C
 25:    TSMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user

 27:    Collective on TS

 29:    Input Parameters:
 30: +  ts - TS object you wish to monitor
 31: .  name - the monitor type one is seeking
 32: .  help - message indicating what monitoring is done
 33: .  manual - manual page for the monitor
 34: .  monitor - the monitor function
 35: -  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

 37:    Level: developer

 39: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
 40:           PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
 41:           PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
 42:           PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
 43:           PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
 44:           PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
 45:           PetscOptionsFList(), PetscOptionsEList()
 46: @*/
 47: PetscErrorCode  TSMonitorSetFromOptions(TS ts,const char name[],const char help[], const char manual[],PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,PetscViewerAndFormat*),PetscErrorCode (*monitorsetup)(TS,PetscViewerAndFormat*))
 48: {
 49:   PetscErrorCode    ierr;
 50:   PetscViewer       viewer;
 51:   PetscViewerFormat format;
 52:   PetscBool         flg;

 55:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
 56:   if (flg) {
 57:     PetscViewerAndFormat *vf;
 58:     PetscViewerAndFormatCreate(viewer,format,&vf);
 59:     PetscObjectDereference((PetscObject)viewer);
 60:     if (monitorsetup) {
 61:       (*monitorsetup)(ts,vf);
 62:     }
 63:     TSMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
 64:   }
 65:   return(0);
 66: }

 70: /*@C
 71:    TSAdjointMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user

 73:    Collective on TS

 75:    Input Parameters:
 76: +  ts - TS object you wish to monitor
 77: .  name - the monitor type one is seeking
 78: .  help - message indicating what monitoring is done
 79: .  manual - manual page for the monitor
 80: .  monitor - the monitor function
 81: -  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

 83:    Level: developer

 85: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
 86:           PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
 87:           PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
 88:           PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
 89:           PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
 90:           PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
 91:           PetscOptionsFList(), PetscOptionsEList()
 92: @*/
 93: 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*))
 94: {
 95:   PetscErrorCode    ierr;
 96:   PetscViewer       viewer;
 97:   PetscViewerFormat format;
 98:   PetscBool         flg;

101:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
102:   if (flg) {
103:     PetscViewerAndFormat *vf;
104:     PetscViewerAndFormatCreate(viewer,format,&vf);
105:     PetscObjectDereference((PetscObject)viewer);
106:     if (monitorsetup) {
107:       (*monitorsetup)(ts,vf);
108:     }
109:     TSAdjointMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
110:   }
111:   return(0);
112: }

116: /*@
117:    TSSetFromOptions - Sets various TS parameters from user options.

119:    Collective on TS

121:    Input Parameter:
122: .  ts - the TS context obtained from TSCreate()

124:    Options Database Keys:
125: +  -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSALPHA, TSGL, TSSSP
126: .  -ts_save_trajectory - checkpoint the solution at each time-step
127: .  -ts_max_steps <maxsteps> - maximum number of time-steps to take
128: .  -ts_final_time <time> - maximum time to compute to
129: .  -ts_dt <dt> - initial time step
130: .  -ts_exact_final_time <stepover,interpolate,matchstep> whether to stop at the exact given final time and how to compute the solution at that ti,e
131: .  -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
132: .  -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
133: .  -ts_error_if_step_fails <true,false> - Error if no step succeeds
134: .  -ts_rtol <rtol> - relative tolerance for local truncation error
135: .  -ts_atol <atol> Absolute tolerance for local truncation error
136: .  -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
137: .  -ts_fd_color - Use finite differences with coloring to compute IJacobian
138: .  -ts_monitor - print information at each timestep
139: .  -ts_monitor_lg_solution - Monitor solution graphically
140: .  -ts_monitor_lg_error - Monitor error graphically
141: .  -ts_monitor_lg_timestep - Monitor timestep size graphically
142: .  -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
143: .  -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
144: .  -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
145: .  -ts_monitor_draw_solution - Monitor solution graphically
146: .  -ts_monitor_draw_solution_phase  <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
147: .  -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
148: .  -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
149: .  -ts_monitor_solution_vtk <filename.vts> - Save each time step to a binary file, use filename-%%03D.vts
150: .  -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
151: .  -ts_adjoint_monitor - print information at each adjoint time step
152: -  -ts_adjoint_monitor_draw_sensi - monitor the sensitivity of the first cost function wrt initial conditions (lambda[0]) graphically

154:    Developer Note: We should unify all the -ts_monitor options in the way that -xxx_view has been unified

156:    Level: beginner

158: .keywords: TS, timestep, set, options, database

160: .seealso: TSGetType()
161: @*/
162: PetscErrorCode  TSSetFromOptions(TS ts)
163: {
164:   PetscBool              opt,flg,tflg;
165:   PetscErrorCode         ierr;
166:   char                   monfilename[PETSC_MAX_PATH_LEN];
167:   PetscReal              time_step;
168:   TSExactFinalTimeOption eftopt;
169:   char                   dir[16];
170:   TSIFunction            ifun;
171:   const char             *defaultType;
172:   char                   typeName[256];


177:   TSRegisterAll();
178:   TSGetIFunction(ts,NULL,&ifun,NULL);

180:   PetscObjectOptionsBegin((PetscObject)ts);
181:   if (((PetscObject)ts)->type_name)
182:     defaultType = ((PetscObject)ts)->type_name;
183:   else
184:     defaultType = ifun ? TSBEULER : TSEULER;
185:   PetscOptionsFList("-ts_type","TS method","TSSetType",TSList,defaultType,typeName,256,&opt);
186:   if (opt) {
187:     TSSetType(ts,typeName);
188:   } else {
189:     TSSetType(ts,defaultType);
190:   }

192:   /* Handle generic TS options */
193:   PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetDuration",ts->max_steps,&ts->max_steps,NULL);
194:   PetscOptionsReal("-ts_final_time","Time to run to","TSSetDuration",ts->max_time,&ts->max_time,NULL);
195:   PetscOptionsReal("-ts_init_time","Initial time","TSSetTime",ts->ptime,&ts->ptime,NULL);
196:   PetscOptionsReal("-ts_dt","Initial time step","TSSetTimeStep",ts->time_step,&time_step,&flg);
197:   if (flg) {TSSetTimeStep(ts,time_step);}
198:   PetscOptionsEnum("-ts_exact_final_time","Option for handling of final time step","TSSetExactFinalTime",TSExactFinalTimeOptions,(PetscEnum)ts->exact_final_time,(PetscEnum*)&eftopt,&flg);
199:   if (flg) {TSSetExactFinalTime(ts,eftopt);}
200:   PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","TSSetMaxSNESFailures",ts->max_snes_failures,&ts->max_snes_failures,NULL);
201:   PetscOptionsInt("-ts_max_reject","Maximum number of step rejections before step fails","TSSetMaxStepRejections",ts->max_reject,&ts->max_reject,NULL);
202:   PetscOptionsBool("-ts_error_if_step_fails","Error if no step succeeds","TSSetErrorIfStepFails",ts->errorifstepfailed,&ts->errorifstepfailed,NULL);
203:   PetscOptionsReal("-ts_rtol","Relative tolerance for local truncation error","TSSetTolerances",ts->rtol,&ts->rtol,NULL);
204:   PetscOptionsReal("-ts_atol","Absolute tolerance for local truncation error","TSSetTolerances",ts->atol,&ts->atol,NULL);

206: #if defined(PETSC_HAVE_SAWS)
207:   {
208:   PetscBool set;
209:   flg  = PETSC_FALSE;
210:   PetscOptionsBool("-ts_saws_block","Block for SAWs memory snooper at end of TSSolve","PetscObjectSAWsBlock",((PetscObject)ts)->amspublishblock,&flg,&set);
211:   if (set) {
212:     PetscObjectSAWsSetBlock((PetscObject)ts,flg);
213:   }
214:   }
215: #endif

217:   /* Monitor options */
218:   TSMonitorSetFromOptions(ts,"-ts_monitor","Monitor time and timestep size","TSMonitorDefault",TSMonitorDefault,NULL);
219:   TSMonitorSetFromOptions(ts,"-ts_monitor_solution","View the solution at each timestep","TSMonitorSolution",TSMonitorSolution,NULL);
220:   TSAdjointMonitorSetFromOptions(ts,"-ts_adjoint_monitor","Monitor adjoint timestep size","TSAdjointMonitorDefault",TSAdjointMonitorDefault,NULL);

222:   PetscOptionsString("-ts_monitor_python","Use Python function","TSMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
223:   if (flg) {PetscPythonMonitorSet((PetscObject)ts,monfilename);}

225:   PetscOptionsName("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",&opt);
226:   if (opt) {
227:     TSMonitorLGCtx ctx;
228:     PetscInt       howoften = 1;

230:     PetscOptionsInt("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",howoften,&howoften,NULL);
231:     TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
232:     TSMonitorSet(ts,TSMonitorLGSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
233:   }

235:   PetscOptionsName("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",&opt);
236:   if (opt) {
237:     TSMonitorLGCtx ctx;
238:     PetscInt       howoften = 1;

240:     PetscOptionsInt("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",howoften,&howoften,NULL);
241:     TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
242:     TSMonitorSet(ts,TSMonitorLGError,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
243:   }

245:   PetscOptionsName("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",&opt);
246:   if (opt) {
247:     TSMonitorLGCtx ctx;
248:     PetscInt       howoften = 1;

250:     PetscOptionsInt("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);
251:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
252:     TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
253:   }
254:   PetscOptionsName("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",&opt);
255:   if (opt) {
256:     TSMonitorLGCtx ctx;
257:     PetscInt       howoften = 1;

259:     PetscOptionsInt("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",howoften,&howoften,NULL);
260:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
261:     TSMonitorSet(ts,TSMonitorLGSNESIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
262:   }
263:   PetscOptionsName("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",&opt);
264:   if (opt) {
265:     TSMonitorLGCtx ctx;
266:     PetscInt       howoften = 1;

268:     PetscOptionsInt("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",howoften,&howoften,NULL);
269:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
270:     TSMonitorSet(ts,TSMonitorLGKSPIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
271:   }
272:   PetscOptionsName("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",&opt);
273:   if (opt) {
274:     TSMonitorSPEigCtx ctx;
275:     PetscInt          howoften = 1;

277:     PetscOptionsInt("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",howoften,&howoften,NULL);
278:     TSMonitorSPEigCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
279:     TSMonitorSet(ts,TSMonitorSPEig,ctx,(PetscErrorCode (*)(void**))TSMonitorSPEigCtxDestroy);
280:   }
281:   opt  = PETSC_FALSE;
282:   PetscOptionsName("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",&opt);
283:   if (opt) {
284:     TSMonitorDrawCtx ctx;
285:     PetscInt         howoften = 1;

287:     PetscOptionsInt("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",howoften,&howoften,NULL);
288:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
289:     TSMonitorSet(ts,TSMonitorDrawSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
290:   }
291:   opt  = PETSC_FALSE;
292:   PetscOptionsName("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",&opt);
293:   if (opt) {
294:     TSMonitorDrawCtx ctx;
295:     PetscInt         howoften = 1;

297:     PetscOptionsInt("-ts_adjoint_monitor_draw_sensi","Monitor adjoint sensitivities (lambda only) graphically","TSAdjointMonitorDrawSensi",howoften,&howoften,NULL);
298:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
299:     TSAdjointMonitorSet(ts,TSAdjointMonitorDrawSensi,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
300:   }
301:   opt  = PETSC_FALSE;
302:   PetscOptionsName("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",&opt);
303:   if (opt) {
304:     TSMonitorDrawCtx ctx;
305:     PetscReal        bounds[4];
306:     PetscInt         n = 4;
307:     PetscDraw        draw;
308:     PetscDrawAxis    axis;

310:     PetscOptionsRealArray("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",bounds,&n,NULL);
311:     if (n != 4) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Must provide bounding box of phase field");
312:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,1,&ctx);
313:     PetscViewerDrawGetDraw(ctx->viewer,0,&draw);
314:     PetscViewerDrawGetDrawAxis(ctx->viewer,0,&axis);
315:     PetscDrawAxisSetLimits(axis,bounds[0],bounds[2],bounds[1],bounds[3]);
316:     PetscDrawAxisSetLabels(axis,"Phase Diagram","Variable 1","Variable 2");
317:     TSMonitorSet(ts,TSMonitorDrawSolutionPhase,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
318:   }
319:   opt  = PETSC_FALSE;
320:   PetscOptionsName("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",&opt);
321:   if (opt) {
322:     TSMonitorDrawCtx ctx;
323:     PetscInt         howoften = 1;

325:     PetscOptionsInt("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",howoften,&howoften,NULL);
326:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
327:     TSMonitorSet(ts,TSMonitorDrawError,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
328:   }

330:   opt  = PETSC_FALSE;
331:   PetscOptionsString("-ts_monitor_solution_vtk","Save each time step to a binary file, use filename-%%03D.vts","TSMonitorSolutionVTK",0,monfilename,PETSC_MAX_PATH_LEN,&flg);
332:   if (flg) {
333:     const char *ptr,*ptr2;
334:     char       *filetemplate;
335:     if (!monfilename[0]) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
336:     /* Do some cursory validation of the input. */
337:     PetscStrstr(monfilename,"%",(char**)&ptr);
338:     if (!ptr) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
339:     for (ptr++; ptr && *ptr; ptr++) {
340:       PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
341:       if (!ptr2 && (*ptr < '0' || '9' < *ptr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Invalid file template argument to -ts_monitor_solution_vtk, should look like filename-%%03D.vts");
342:       if (ptr2) break;
343:     }
344:     PetscStrallocpy(monfilename,&filetemplate);
345:     TSMonitorSet(ts,TSMonitorSolutionVTK,filetemplate,(PetscErrorCode (*)(void**))TSMonitorSolutionVTKDestroy);
346:   }

348:   PetscOptionsString("-ts_monitor_dmda_ray","Display a ray of the solution","None","y=0",dir,16,&flg);
349:   if (flg) {
350:     TSMonitorDMDARayCtx *rayctx;
351:     int                  ray = 0;
352:     DMDADirection        ddir;
353:     DM                   da;
354:     PetscMPIInt          rank;

356:     if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
357:     if (dir[0] == 'x') ddir = DMDA_X;
358:     else if (dir[0] == 'y') ddir = DMDA_Y;
359:     else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
360:     sscanf(dir+2,"%d",&ray);

362:     PetscInfo2(((PetscObject)ts),"Displaying DMDA ray %c = %D\n",dir[0],ray);
363:     PetscNew(&rayctx);
364:     TSGetDM(ts,&da);
365:     DMDAGetRay(da,ddir,ray,&rayctx->ray,&rayctx->scatter);
366:     MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);
367:     if (!rank) {
368:       PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,600,300,&rayctx->viewer);
369:     }
370:     rayctx->lgctx = NULL;
371:     TSMonitorSet(ts,TSMonitorDMDARay,rayctx,TSMonitorDMDARayDestroy);
372:   }
373:   PetscOptionsString("-ts_monitor_lg_dmda_ray","Display a ray of the solution","None","x=0",dir,16,&flg);
374:   if (flg) {
375:     TSMonitorDMDARayCtx *rayctx;
376:     int                 ray = 0;
377:     DMDADirection       ddir;
378:     DM                  da;
379:     PetscInt            howoften = 1;

381:     if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
382:     if      (dir[0] == 'x') ddir = DMDA_X;
383:     else if (dir[0] == 'y') ddir = DMDA_Y;
384:     else SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
385:     sscanf(dir+2, "%d", &ray);

387:     PetscInfo2(((PetscObject) ts),"Displaying LG DMDA ray %c = %D\n", dir[0], ray);
388:     PetscNew(&rayctx);
389:     TSGetDM(ts, &da);
390:     DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);
391:     TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&rayctx->lgctx);
392:     TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy);
393:   }

395:   PetscOptionsName("-ts_monitor_envelope","Monitor maximum and minimum value of each component of the solution","TSMonitorEnvelope",&opt);
396:   if (opt) {
397:     TSMonitorEnvelopeCtx ctx;

399:     TSMonitorEnvelopeCtxCreate(ts,&ctx);
400:     TSMonitorSet(ts,TSMonitorEnvelope,ctx,(PetscErrorCode (*)(void**))TSMonitorEnvelopeCtxDestroy);
401:   }

403:   flg  = PETSC_FALSE;
404:   PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeJacobianDefaultColor", flg, &flg, NULL);
405:   if (flg) {
406:     DM   dm;
407:     DMTS tdm;

409:     TSGetDM(ts, &dm);
410:     DMGetDMTS(dm, &tdm);
411:     tdm->ijacobianctx = NULL;
412:     TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, 0);
413:     PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n");
414:   }

416:   if (ts->adapt) {
417:     TSAdaptSetFromOptions(PetscOptionsObject,ts->adapt);
418:   }

420:   /* Handle specific TS options */
421:   if (ts->ops->setfromoptions) {
422:     (*ts->ops->setfromoptions)(PetscOptionsObject,ts);
423:   }

425:   /* TS trajectory must be set after TS, since it may use some TS options above */
426:   tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
427:   PetscOptionsBool("-ts_save_trajectory","Save the solution at each timestep","TSSetSaveTrajectory",tflg,&tflg,NULL);
428:   if (tflg) {
429:     TSSetSaveTrajectory(ts);
430:   }
431:   tflg = ts->adjoint_solve ? PETSC_TRUE : PETSC_FALSE;
432:   PetscOptionsBool("-ts_adjoint_solve","Solve the adjoint problem immediately after solving the forward problem","",tflg,&tflg,&flg);
433:   if (flg) {
434:     TSSetSaveTrajectory(ts);
435:     ts->adjoint_solve = tflg;
436:   }

438:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
439:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)ts);
440:   PetscOptionsEnd();

442:   if (ts->trajectory) {
443:     TSTrajectorySetFromOptions(ts->trajectory,ts);
444:   }

446:   TSGetSNES(ts,&ts->snes);
447:   if (ts->problem_type == TS_LINEAR) {SNESSetType(ts->snes,SNESKSPONLY);}
448:   SNESSetFromOptions(ts->snes);
449:   return(0);
450: }

454: /*@
455:    TSSetSaveTrajectory - Causes the TS to save its solutions as it iterates forward in time in a TSTrajectory object

457:    Collective on TS

459:    Input Parameters:
460: .  ts - the TS context obtained from TSCreate()

462: Note: This routine should be called after all TS options have been set

464:    Level: intermediate

466: .seealso: TSGetTrajectory(), TSAdjointSolve()

468: .keywords: TS, set, checkpoint,
469: @*/
470: PetscErrorCode  TSSetSaveTrajectory(TS ts)
471: {

476:   if (!ts->trajectory) {
477:     TSTrajectoryCreate(PetscObjectComm((PetscObject)ts),&ts->trajectory);
478:     TSTrajectorySetFromOptions(ts->trajectory,ts);
479:   }
480:   return(0);
481: }

485: /*@
486:    TSComputeRHSJacobian - Computes the Jacobian matrix that has been
487:       set with TSSetRHSJacobian().

489:    Collective on TS and Vec

491:    Input Parameters:
492: +  ts - the TS context
493: .  t - current timestep
494: -  U - input vector

496:    Output Parameters:
497: +  A - Jacobian matrix
498: .  B - optional preconditioning matrix
499: -  flag - flag indicating matrix structure

501:    Notes:
502:    Most users should not need to explicitly call this routine, as it
503:    is used internally within the nonlinear solvers.

505:    See KSPSetOperators() for important information about setting the
506:    flag parameter.

508:    Level: developer

510: .keywords: SNES, compute, Jacobian, matrix

512: .seealso:  TSSetRHSJacobian(), KSPSetOperators()
513: @*/
514: PetscErrorCode  TSComputeRHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B)
515: {
517:   PetscObjectState Ustate;
518:   DM             dm;
519:   DMTS           tsdm;
520:   TSRHSJacobian  rhsjacobianfunc;
521:   void           *ctx;
522:   TSIJacobian    ijacobianfunc;
523:   TSRHSFunction  rhsfunction;

529:   TSGetDM(ts,&dm);
530:   DMGetDMTS(dm,&tsdm);
531:   DMTSGetRHSJacobian(dm,&rhsjacobianfunc,&ctx);
532:   DMTSGetIJacobian(dm,&ijacobianfunc,NULL);
533:   DMTSGetRHSFunction(dm,&rhsfunction,&ctx);
534:   PetscObjectStateGet((PetscObject)U,&Ustate);
535:   if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.X == U && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) {
536:     return(0);
537:   }

539:   if (!rhsjacobianfunc && !ijacobianfunc) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");

541:   if (ts->rhsjacobian.reuse) {
542:     MatShift(A,-ts->rhsjacobian.shift);
543:     MatScale(A,1./ts->rhsjacobian.scale);
544:     if (A != B) {
545:       MatShift(B,-ts->rhsjacobian.shift);
546:       MatScale(B,1./ts->rhsjacobian.scale);
547:     }
548:     ts->rhsjacobian.shift = 0;
549:     ts->rhsjacobian.scale = 1.;
550:   }

552:   if (rhsjacobianfunc) {
553:     PetscBool missing;
554:     PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);
555:     PetscStackPush("TS user Jacobian function");
556:     (*rhsjacobianfunc)(ts,t,U,A,B,ctx);
557:     PetscStackPop;
558:     PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);
559:     if (A) {
560:       MatMissingDiagonal(A,&missing,NULL);
561:       if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Amat passed to TSSetRHSJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
562:     }
563:     if (B && B != A) {
564:       MatMissingDiagonal(B,&missing,NULL);
565:       if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Bmat passed to TSSetRHSJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
566:     }
567:   } else {
568:     MatZeroEntries(A);
569:     if (A != B) {MatZeroEntries(B);}
570:   }
571:   ts->rhsjacobian.time       = t;
572:   ts->rhsjacobian.X          = U;
573:   PetscObjectStateGet((PetscObject)U,&ts->rhsjacobian.Xstate);
574:   return(0);
575: }

579: /*@
580:    TSComputeRHSFunction - Evaluates the right-hand-side function.

582:    Collective on TS and Vec

584:    Input Parameters:
585: +  ts - the TS context
586: .  t - current time
587: -  U - state vector

589:    Output Parameter:
590: .  y - right hand side

592:    Note:
593:    Most users should not need to explicitly call this routine, as it
594:    is used internally within the nonlinear solvers.

596:    Level: developer

598: .keywords: TS, compute

600: .seealso: TSSetRHSFunction(), TSComputeIFunction()
601: @*/
602: PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec U,Vec y)
603: {
605:   TSRHSFunction  rhsfunction;
606:   TSIFunction    ifunction;
607:   void           *ctx;
608:   DM             dm;

614:   TSGetDM(ts,&dm);
615:   DMTSGetRHSFunction(dm,&rhsfunction,&ctx);
616:   DMTSGetIFunction(dm,&ifunction,NULL);

618:   if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");

620:   PetscLogEventBegin(TS_FunctionEval,ts,U,y,0);
621:   if (rhsfunction) {
622:     PetscStackPush("TS user right-hand-side function");
623:     (*rhsfunction)(ts,t,U,y,ctx);
624:     PetscStackPop;
625:   } else {
626:     VecZeroEntries(y);
627:   }

629:   PetscLogEventEnd(TS_FunctionEval,ts,U,y,0);
630:   return(0);
631: }

635: /*@
636:    TSComputeSolutionFunction - Evaluates the solution function.

638:    Collective on TS and Vec

640:    Input Parameters:
641: +  ts - the TS context
642: -  t - current time

644:    Output Parameter:
645: .  U - the solution

647:    Note:
648:    Most users should not need to explicitly call this routine, as it
649:    is used internally within the nonlinear solvers.

651:    Level: developer

653: .keywords: TS, compute

655: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
656: @*/
657: PetscErrorCode TSComputeSolutionFunction(TS ts,PetscReal t,Vec U)
658: {
659:   PetscErrorCode     ierr;
660:   TSSolutionFunction solutionfunction;
661:   void               *ctx;
662:   DM                 dm;

667:   TSGetDM(ts,&dm);
668:   DMTSGetSolutionFunction(dm,&solutionfunction,&ctx);

670:   if (solutionfunction) {
671:     PetscStackPush("TS user solution function");
672:     (*solutionfunction)(ts,t,U,ctx);
673:     PetscStackPop;
674:   }
675:   return(0);
676: }
679: /*@
680:    TSComputeForcingFunction - Evaluates the forcing function.

682:    Collective on TS and Vec

684:    Input Parameters:
685: +  ts - the TS context
686: -  t - current time

688:    Output Parameter:
689: .  U - the function value

691:    Note:
692:    Most users should not need to explicitly call this routine, as it
693:    is used internally within the nonlinear solvers.

695:    Level: developer

697: .keywords: TS, compute

699: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
700: @*/
701: PetscErrorCode TSComputeForcingFunction(TS ts,PetscReal t,Vec U)
702: {
703:   PetscErrorCode     ierr, (*forcing)(TS,PetscReal,Vec,void*);
704:   void               *ctx;
705:   DM                 dm;

710:   TSGetDM(ts,&dm);
711:   DMTSGetForcingFunction(dm,&forcing,&ctx);

713:   if (forcing) {
714:     PetscStackPush("TS user forcing function");
715:     (*forcing)(ts,t,U,ctx);
716:     PetscStackPop;
717:   }
718:   return(0);
719: }

723: static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs)
724: {
725:   Vec            F;

729:   *Frhs = NULL;
730:   TSGetIFunction(ts,&F,NULL,NULL);
731:   if (!ts->Frhs) {
732:     VecDuplicate(F,&ts->Frhs);
733:   }
734:   *Frhs = ts->Frhs;
735:   return(0);
736: }

740: static PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs)
741: {
742:   Mat            A,B;

746:   if (Arhs) *Arhs = NULL;
747:   if (Brhs) *Brhs = NULL;
748:   TSGetIJacobian(ts,&A,&B,NULL,NULL);
749:   if (Arhs) {
750:     if (!ts->Arhs) {
751:       MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);
752:     }
753:     *Arhs = ts->Arhs;
754:   }
755:   if (Brhs) {
756:     if (!ts->Brhs) {
757:       if (A != B) {
758:         MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);
759:       } else {
760:         PetscObjectReference((PetscObject)ts->Arhs);
761:         ts->Brhs = ts->Arhs;
762:       }
763:     }
764:     *Brhs = ts->Brhs;
765:   }
766:   return(0);
767: }

771: /*@
772:    TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,U,Udot)=0

774:    Collective on TS and Vec

776:    Input Parameters:
777: +  ts - the TS context
778: .  t - current time
779: .  U - state vector
780: .  Udot - time derivative of state vector
781: -  imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate

783:    Output Parameter:
784: .  Y - right hand side

786:    Note:
787:    Most users should not need to explicitly call this routine, as it
788:    is used internally within the nonlinear solvers.

790:    If the user did did not write their equations in implicit form, this
791:    function recasts them in implicit form.

793:    Level: developer

795: .keywords: TS, compute

797: .seealso: TSSetIFunction(), TSComputeRHSFunction()
798: @*/
799: PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec Y,PetscBool imex)
800: {
802:   TSIFunction    ifunction;
803:   TSRHSFunction  rhsfunction;
804:   void           *ctx;
805:   DM             dm;


813:   TSGetDM(ts,&dm);
814:   DMTSGetIFunction(dm,&ifunction,&ctx);
815:   DMTSGetRHSFunction(dm,&rhsfunction,NULL);

817:   if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");

819:   PetscLogEventBegin(TS_FunctionEval,ts,U,Udot,Y);
820:   if (ifunction) {
821:     PetscStackPush("TS user implicit function");
822:     (*ifunction)(ts,t,U,Udot,Y,ctx);
823:     PetscStackPop;
824:   }
825:   if (imex) {
826:     if (!ifunction) {
827:       VecCopy(Udot,Y);
828:     }
829:   } else if (rhsfunction) {
830:     if (ifunction) {
831:       Vec Frhs;
832:       TSGetRHSVec_Private(ts,&Frhs);
833:       TSComputeRHSFunction(ts,t,U,Frhs);
834:       VecAXPY(Y,-1,Frhs);
835:     } else {
836:       TSComputeRHSFunction(ts,t,U,Y);
837:       VecAYPX(Y,-1,Udot);
838:     }
839:   }
840:   PetscLogEventEnd(TS_FunctionEval,ts,U,Udot,Y);
841:   return(0);
842: }

846: /*@
847:    TSComputeIJacobian - Evaluates the Jacobian of the DAE

849:    Collective on TS and Vec

851:    Input
852:       Input Parameters:
853: +  ts - the TS context
854: .  t - current timestep
855: .  U - state vector
856: .  Udot - time derivative of state vector
857: .  shift - shift to apply, see note below
858: -  imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate

860:    Output Parameters:
861: +  A - Jacobian matrix
862: .  B - optional preconditioning matrix
863: -  flag - flag indicating matrix structure

865:    Notes:
866:    If F(t,U,Udot)=0 is the DAE, the required Jacobian is

868:    dF/dU + shift*dF/dUdot

870:    Most users should not need to explicitly call this routine, as it
871:    is used internally within the nonlinear solvers.

873:    Level: developer

875: .keywords: TS, compute, Jacobian, matrix

877: .seealso:  TSSetIJacobian()
878: @*/
879: PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,PetscBool imex)
880: {
882:   TSIJacobian    ijacobian;
883:   TSRHSJacobian  rhsjacobian;
884:   DM             dm;
885:   void           *ctx;


896:   TSGetDM(ts,&dm);
897:   DMTSGetIJacobian(dm,&ijacobian,&ctx);
898:   DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);

900:   if (!rhsjacobian && !ijacobian) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");

902:   PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);
903:   if (ijacobian) {
904:     PetscBool missing;
905:     PetscStackPush("TS user implicit Jacobian");
906:     (*ijacobian)(ts,t,U,Udot,shift,A,B,ctx);
907:     PetscStackPop;
908:     if (A) {
909:       MatMissingDiagonal(A,&missing,NULL);
910:       if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Amat passed to TSSetIJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
911:     }
912:     if (B && B != A) {
913:       MatMissingDiagonal(B,&missing,NULL);
914:       if (missing) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Bmat passed to TSSetIJacobian() must have all diagonal entries set, if they are zero you must still set them with a zero value");
915:     }
916:   }
917:   if (imex) {
918:     if (!ijacobian) {  /* system was written as Udot = G(t,U) */
919:       PetscBool assembled;
920:       MatZeroEntries(A);
921:       MatAssembled(A,&assembled);
922:       if (!assembled) {
923:         MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
924:         MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
925:       }
926:       MatShift(A,shift);
927:       if (A != B) {
928:         MatZeroEntries(B);
929:         MatAssembled(B,&assembled);
930:         if (!assembled) {
931:           MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
932:           MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
933:         }
934:         MatShift(B,shift);
935:       }
936:     }
937:   } else {
938:     Mat Arhs = NULL,Brhs = NULL;
939:     if (rhsjacobian) {
940:       TSGetRHSMats_Private(ts,&Arhs,&Brhs);
941:       TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);
942:     }
943:     if (Arhs == A) {           /* No IJacobian, so we only have the RHS matrix */
944:       ts->rhsjacobian.scale = -1;
945:       ts->rhsjacobian.shift = shift;
946:       MatScale(A,-1);
947:       MatShift(A,shift);
948:       if (A != B) {
949:         MatScale(B,-1);
950:         MatShift(B,shift);
951:       }
952:     } else if (Arhs) {          /* Both IJacobian and RHSJacobian */
953:       MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
954:       if (!ijacobian) {         /* No IJacobian provided, but we have a separate RHS matrix */
955:         MatZeroEntries(A);
956:         MatShift(A,shift);
957:         if (A != B) {
958:           MatZeroEntries(B);
959:           MatShift(B,shift);
960:         }
961:       }
962:       MatAXPY(A,-1,Arhs,axpy);
963:       if (A != B) {
964:         MatAXPY(B,-1,Brhs,axpy);
965:       }
966:     }
967:   }
968:   PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);
969:   return(0);
970: }

974: /*@C
975:     TSSetRHSFunction - Sets the routine for evaluating the function,
976:     where U_t = G(t,u).

978:     Logically Collective on TS

980:     Input Parameters:
981: +   ts - the TS context obtained from TSCreate()
982: .   r - vector to put the computed right hand side (or NULL to have it created)
983: .   f - routine for evaluating the right-hand-side function
984: -   ctx - [optional] user-defined context for private data for the
985:           function evaluation routine (may be NULL)

987:     Calling sequence of func:
988: $     func (TS ts,PetscReal t,Vec u,Vec F,void *ctx);

990: +   t - current timestep
991: .   u - input vector
992: .   F - function vector
993: -   ctx - [optional] user-defined function context

995:     Level: beginner

997:     Notes: You must call this function or TSSetIFunction() to define your ODE. You cannot use this function when solving a DAE.

999: .keywords: TS, timestep, set, right-hand-side, function

1001: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSSetIFunction()
1002: @*/
1003: PetscErrorCode  TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx)
1004: {
1006:   SNES           snes;
1007:   Vec            ralloc = NULL;
1008:   DM             dm;


1014:   TSGetDM(ts,&dm);
1015:   DMTSSetRHSFunction(dm,f,ctx);
1016:   TSGetSNES(ts,&snes);
1017:   if (!r && !ts->dm && ts->vec_sol) {
1018:     VecDuplicate(ts->vec_sol,&ralloc);
1019:     r = ralloc;
1020:   }
1021:   SNESSetFunction(snes,r,SNESTSFormFunction,ts);
1022:   VecDestroy(&ralloc);
1023:   return(0);
1024: }

1028: /*@C
1029:     TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE

1031:     Logically Collective on TS

1033:     Input Parameters:
1034: +   ts - the TS context obtained from TSCreate()
1035: .   f - routine for evaluating the solution
1036: -   ctx - [optional] user-defined context for private data for the
1037:           function evaluation routine (may be NULL)

1039:     Calling sequence of func:
1040: $     func (TS ts,PetscReal t,Vec u,void *ctx);

1042: +   t - current timestep
1043: .   u - output vector
1044: -   ctx - [optional] user-defined function context

1046:     Notes:
1047:     This routine is used for testing accuracy of time integration schemes when you already know the solution.
1048:     If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1049:     create closed-form solutions with non-physical forcing terms.

1051:     For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history.

1053:     Level: beginner

1055: .keywords: TS, timestep, set, right-hand-side, function

1057: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetForcingFunction()
1058: @*/
1059: PetscErrorCode  TSSetSolutionFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx)
1060: {
1062:   DM             dm;

1066:   TSGetDM(ts,&dm);
1067:   DMTSSetSolutionFunction(dm,f,ctx);
1068:   return(0);
1069: }

1073: /*@C
1074:     TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE

1076:     Logically Collective on TS

1078:     Input Parameters:
1079: +   ts - the TS context obtained from TSCreate()
1080: .   f - routine for evaluating the forcing function
1081: -   ctx - [optional] user-defined context for private data for the
1082:           function evaluation routine (may be NULL)

1084:     Calling sequence of func:
1085: $     func (TS ts,PetscReal t,Vec u,void *ctx);

1087: +   t - current timestep
1088: .   u - output vector
1089: -   ctx - [optional] user-defined function context

1091:     Notes:
1092:     This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1093:     create closed-form solutions with a non-physical forcing term.

1095:     For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history.

1097:     Level: beginner

1099: .keywords: TS, timestep, set, right-hand-side, function

1101: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetSolutionFunction()
1102: @*/
1103: PetscErrorCode  TSSetForcingFunction(TS ts,TSForcingFunction f,void *ctx)
1104: {
1106:   DM             dm;

1110:   TSGetDM(ts,&dm);
1111:   DMTSSetForcingFunction(dm,f,ctx);
1112:   return(0);
1113: }

1117: /*@C
1118:    TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1119:    where U_t = G(U,t), as well as the location to store the matrix.

1121:    Logically Collective on TS

1123:    Input Parameters:
1124: +  ts  - the TS context obtained from TSCreate()
1125: .  Amat - (approximate) Jacobian matrix
1126: .  Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1127: .  f   - the Jacobian evaluation routine
1128: -  ctx - [optional] user-defined context for private data for the
1129:          Jacobian evaluation routine (may be NULL)

1131:    Calling sequence of f:
1132: $     func (TS ts,PetscReal t,Vec u,Mat A,Mat B,void *ctx);

1134: +  t - current timestep
1135: .  u - input vector
1136: .  Amat - (approximate) Jacobian matrix
1137: .  Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1138: -  ctx - [optional] user-defined context for matrix evaluation routine

1140:    Notes:
1141:    You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value

1143:    The TS solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f()
1144:    You should not assume the values are the same in the next call to f() as you set them in the previous call.

1146:    Level: beginner

1148: .keywords: TS, timestep, set, right-hand-side, Jacobian

1150: .seealso: SNESComputeJacobianDefaultColor(), TSSetRHSFunction(), TSRHSJacobianSetReuse(), TSSetIJacobian()

1152: @*/
1153: PetscErrorCode  TSSetRHSJacobian(TS ts,Mat Amat,Mat Pmat,TSRHSJacobian f,void *ctx)
1154: {
1156:   SNES           snes;
1157:   DM             dm;
1158:   TSIJacobian    ijacobian;


1167:   TSGetDM(ts,&dm);
1168:   DMTSSetRHSJacobian(dm,f,ctx);
1169:   if (f == TSComputeRHSJacobianConstant) {
1170:     /* Handle this case automatically for the user; otherwise user should call themselves. */
1171:     TSRHSJacobianSetReuse(ts,PETSC_TRUE);
1172:   }
1173:   DMTSGetIJacobian(dm,&ijacobian,NULL);
1174:   TSGetSNES(ts,&snes);
1175:   if (!ijacobian) {
1176:     SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
1177:   }
1178:   if (Amat) {
1179:     PetscObjectReference((PetscObject)Amat);
1180:     MatDestroy(&ts->Arhs);
1181:     ts->Arhs = Amat;
1182:   }
1183:   if (Pmat) {
1184:     PetscObjectReference((PetscObject)Pmat);
1185:     MatDestroy(&ts->Brhs);
1186:     ts->Brhs = Pmat;
1187:   }
1188:   return(0);
1189: }


1194: /*@C
1195:    TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.

1197:    Logically Collective on TS

1199:    Input Parameters:
1200: +  ts  - the TS context obtained from TSCreate()
1201: .  r   - vector to hold the residual (or NULL to have it created internally)
1202: .  f   - the function evaluation routine
1203: -  ctx - user-defined context for private data for the function evaluation routine (may be NULL)

1205:    Calling sequence of f:
1206: $  f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx);

1208: +  t   - time at step/stage being solved
1209: .  u   - state vector
1210: .  u_t - time derivative of state vector
1211: .  F   - function vector
1212: -  ctx - [optional] user-defined context for matrix evaluation routine

1214:    Important:
1215:    The user MUST call either this routine or TSSetRHSFunction() to define the ODE.  When solving DAEs you must use this function.

1217:    Level: beginner

1219: .keywords: TS, timestep, set, DAE, Jacobian

1221: .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian()
1222: @*/
1223: PetscErrorCode  TSSetIFunction(TS ts,Vec r,TSIFunction f,void *ctx)
1224: {
1226:   SNES           snes;
1227:   Vec            ralloc = NULL;
1228:   DM             dm;


1234:   TSGetDM(ts,&dm);
1235:   DMTSSetIFunction(dm,f,ctx);

1237:   TSGetSNES(ts,&snes);
1238:   if (!r && !ts->dm && ts->vec_sol) {
1239:     VecDuplicate(ts->vec_sol,&ralloc);
1240:     r  = ralloc;
1241:   }
1242:   SNESSetFunction(snes,r,SNESTSFormFunction,ts);
1243:   VecDestroy(&ralloc);
1244:   return(0);
1245: }

1249: /*@C
1250:    TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it.

1252:    Not Collective

1254:    Input Parameter:
1255: .  ts - the TS context

1257:    Output Parameter:
1258: +  r - vector to hold residual (or NULL)
1259: .  func - the function to compute residual (or NULL)
1260: -  ctx - the function context (or NULL)

1262:    Level: advanced

1264: .keywords: TS, nonlinear, get, function

1266: .seealso: TSSetIFunction(), SNESGetFunction()
1267: @*/
1268: PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx)
1269: {
1271:   SNES           snes;
1272:   DM             dm;

1276:   TSGetSNES(ts,&snes);
1277:   SNESGetFunction(snes,r,NULL,NULL);
1278:   TSGetDM(ts,&dm);
1279:   DMTSGetIFunction(dm,func,ctx);
1280:   return(0);
1281: }

1285: /*@C
1286:    TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it.

1288:    Not Collective

1290:    Input Parameter:
1291: .  ts - the TS context

1293:    Output Parameter:
1294: +  r - vector to hold computed right hand side (or NULL)
1295: .  func - the function to compute right hand side (or NULL)
1296: -  ctx - the function context (or NULL)

1298:    Level: advanced

1300: .keywords: TS, nonlinear, get, function

1302: .seealso: TSSetRHSFunction(), SNESGetFunction()
1303: @*/
1304: PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx)
1305: {
1307:   SNES           snes;
1308:   DM             dm;

1312:   TSGetSNES(ts,&snes);
1313:   SNESGetFunction(snes,r,NULL,NULL);
1314:   TSGetDM(ts,&dm);
1315:   DMTSGetRHSFunction(dm,func,ctx);
1316:   return(0);
1317: }

1321: /*@C
1322:    TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1323:         provided with TSSetIFunction().

1325:    Logically Collective on TS

1327:    Input Parameters:
1328: +  ts  - the TS context obtained from TSCreate()
1329: .  Amat - (approximate) Jacobian matrix
1330: .  Pmat - matrix used to compute preconditioner (usually the same as Amat)
1331: .  f   - the Jacobian evaluation routine
1332: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)

1334:    Calling sequence of f:
1335: $  f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat Amat,Mat Pmat,void *ctx);

1337: +  t    - time at step/stage being solved
1338: .  U    - state vector
1339: .  U_t  - time derivative of state vector
1340: .  a    - shift
1341: .  Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1342: .  Pmat - matrix used for constructing preconditioner, usually the same as Amat
1343: -  ctx  - [optional] user-defined context for matrix evaluation routine

1345:    Notes:
1346:    The matrices Amat and Pmat are exactly the matrices that are used by SNES for the nonlinear solve.

1348:    If you know the operator Amat has a null space you can use MatSetNullSpace() and MatSetTransposeNullSpace() to supply the null
1349:    space to Amat and the KSP solvers will automatically use that null space as needed during the solution process.

1351:    The matrix dF/dU + a*dF/dU_t you provide turns out to be
1352:    the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1353:    The time integrator internally approximates U_t by W+a*U where the positive "shift"
1354:    a and vector W depend on the integration method, step size, and past states. For example with
1355:    the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1356:    W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt

1358:    You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value

1360:    The TS solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f()
1361:    You should not assume the values are the same in the next call to f() as you set them in the previous call.

1363:    Level: beginner

1365: .keywords: TS, timestep, DAE, Jacobian

1367: .seealso: TSSetIFunction(), TSSetRHSJacobian(), SNESComputeJacobianDefaultColor(), SNESComputeJacobianDefault(), TSSetRHSFunction()

1369: @*/
1370: PetscErrorCode  TSSetIJacobian(TS ts,Mat Amat,Mat Pmat,TSIJacobian f,void *ctx)
1371: {
1373:   SNES           snes;
1374:   DM             dm;


1383:   TSGetDM(ts,&dm);
1384:   DMTSSetIJacobian(dm,f,ctx);

1386:   TSGetSNES(ts,&snes);
1387:   SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
1388:   return(0);
1389: }

1393: /*@
1394:    TSRHSJacobianSetReuse - restore RHS Jacobian before re-evaluating.  Without this flag, TS will change the sign and
1395:    shift the RHS Jacobian for a finite-time-step implicit solve, in which case the user function will need to recompute
1396:    the entire Jacobian.  The reuse flag must be set if the evaluation function will assume that the matrix entries have
1397:    not been changed by the TS.

1399:    Logically Collective

1401:    Input Arguments:
1402: +  ts - TS context obtained from TSCreate()
1403: -  reuse - PETSC_TRUE if the RHS Jacobian

1405:    Level: intermediate

1407: .seealso: TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
1408: @*/
1409: PetscErrorCode TSRHSJacobianSetReuse(TS ts,PetscBool reuse)
1410: {
1412:   ts->rhsjacobian.reuse = reuse;
1413:   return(0);
1414: }

1418: /*@C
1419:    TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.

1421:    Logically Collective on TS

1423:    Input Parameters:
1424: +  ts  - the TS context obtained from TSCreate()
1425: .  F   - vector to hold the residual (or NULL to have it created internally)
1426: .  fun - the function evaluation routine
1427: -  ctx - user-defined context for private data for the function evaluation routine (may be NULL)

1429:    Calling sequence of fun:
1430: $  fun(TS ts,PetscReal t,Vec U,Vec U_t,Vec U_tt,Vec F,ctx);

1432: +  t    - time at step/stage being solved
1433: .  U    - state vector
1434: .  U_t  - time derivative of state vector
1435: .  U_tt - second time derivative of state vector
1436: .  F    - function vector
1437: -  ctx  - [optional] user-defined context for matrix evaluation routine (may be NULL)

1439:    Level: beginner

1441: .keywords: TS, timestep, set, ODE, DAE, Function

1443: .seealso: TSSetI2Jacobian()
1444: @*/
1445: PetscErrorCode TSSetI2Function(TS ts,Vec F,TSI2Function fun,void *ctx)
1446: {
1447:   DM             dm;

1453:   TSSetIFunction(ts,F,NULL,NULL);
1454:   TSGetDM(ts,&dm);
1455:   DMTSSetI2Function(dm,fun,ctx);
1456:   return(0);
1457: }

1461: /*@C
1462:   TSGetI2Function - Returns the vector where the implicit residual is stored and the function/contex to compute it.

1464:   Not Collective

1466:   Input Parameter:
1467: . ts - the TS context

1469:   Output Parameter:
1470: + r - vector to hold residual (or NULL)
1471: . fun - the function to compute residual (or NULL)
1472: - ctx - the function context (or NULL)

1474:   Level: advanced

1476: .keywords: TS, nonlinear, get, function

1478: .seealso: TSSetI2Function(), SNESGetFunction()
1479: @*/
1480: PetscErrorCode TSGetI2Function(TS ts,Vec *r,TSI2Function *fun,void **ctx)
1481: {
1483:   SNES           snes;
1484:   DM             dm;

1488:   TSGetSNES(ts,&snes);
1489:   SNESGetFunction(snes,r,NULL,NULL);
1490:   TSGetDM(ts,&dm);
1491:   DMTSGetI2Function(dm,fun,ctx);
1492:   return(0);
1493: }

1497: /*@C
1498:    TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t  + a*dF/dU_tt
1499:         where F(t,U,U_t,U_tt) is the function you provided with TSSetI2Function().

1501:    Logically Collective on TS

1503:    Input Parameters:
1504: +  ts  - the TS context obtained from TSCreate()
1505: .  J   - Jacobian matrix
1506: .  P   - preconditioning matrix for J (may be same as J)
1507: .  jac - the Jacobian evaluation routine
1508: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)

1510:    Calling sequence of jac:
1511: $  jac(TS ts,PetscReal t,Vec U,Vec U_t,Vec U_tt,PetscReal v,PetscReal a,Mat J,Mat P,void *ctx);

1513: +  t    - time at step/stage being solved
1514: .  U    - state vector
1515: .  U_t  - time derivative of state vector
1516: .  U_tt - second time derivative of state vector
1517: .  v    - shift for U_t
1518: .  a    - shift for U_tt
1519: .  J    - Jacobian of G(U) = F(t,U,W+v*U,W'+a*U), equivalent to dF/dU + v*dF/dU_t  + a*dF/dU_tt
1520: .  P    - preconditioning matrix for J, may be same as J
1521: -  ctx  - [optional] user-defined context for matrix evaluation routine

1523:    Notes:
1524:    The matrices J and P are exactly the matrices that are used by SNES for the nonlinear solve.

1526:    The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1527:    the Jacobian of G(U) = F(t,U,W+v*U,W'+a*U) where F(t,U,U_t,U_tt) = 0 is the DAE to be solved.
1528:    The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U  where the positive "shift"
1529:    parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.

1531:    Level: beginner

1533: .keywords: TS, timestep, set, ODE, DAE, Jacobian

1535: .seealso: TSSetI2Function()
1536: @*/
1537: PetscErrorCode TSSetI2Jacobian(TS ts,Mat J,Mat P,TSI2Jacobian jac,void *ctx)
1538: {
1539:   DM             dm;

1546:   TSSetIJacobian(ts,J,P,NULL,NULL);
1547:   TSGetDM(ts,&dm);
1548:   DMTSSetI2Jacobian(dm,jac,ctx);
1549:   return(0);
1550: }

1554: /*@C
1555:   TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.

1557:   Not Collective, but parallel objects are returned if TS is parallel

1559:   Input Parameter:
1560: . ts  - The TS context obtained from TSCreate()

1562:   Output Parameters:
1563: + J  - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1564: . P - The matrix from which the preconditioner is constructed, often the same as J
1565: . jac - The function to compute the Jacobian matrices
1566: - ctx - User-defined context for Jacobian evaluation routine

1568:   Notes: You can pass in NULL for any return argument you do not need.

1570:   Level: advanced

1572: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber()

1574: .keywords: TS, timestep, get, matrix, Jacobian
1575: @*/
1576: PetscErrorCode  TSGetI2Jacobian(TS ts,Mat *J,Mat *P,TSI2Jacobian *jac,void **ctx)
1577: {
1579:   SNES           snes;
1580:   DM             dm;

1583:   TSGetSNES(ts,&snes);
1584:   SNESSetUpMatrices(snes);
1585:   SNESGetJacobian(snes,J,P,NULL,NULL);
1586:   TSGetDM(ts,&dm);
1587:   DMTSGetI2Jacobian(dm,jac,ctx);
1588:   return(0);
1589: }

1593: /*@
1594:   TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0

1596:   Collective on TS and Vec

1598:   Input Parameters:
1599: + ts - the TS context
1600: . t - current time
1601: . U - state vector
1602: . V - time derivative of state vector (U_t)
1603: - A - second time derivative of state vector (U_tt)

1605:   Output Parameter:
1606: . F - the residual vector

1608:   Note:
1609:   Most users should not need to explicitly call this routine, as it
1610:   is used internally within the nonlinear solvers.

1612:   Level: developer

1614: .keywords: TS, compute, function, vector

1616: .seealso: TSSetI2Function()
1617: @*/
1618: PetscErrorCode TSComputeI2Function(TS ts,PetscReal t,Vec U,Vec V,Vec A,Vec F)
1619: {
1620:   DM             dm;
1621:   TSI2Function   I2Function;
1622:   void           *ctx;
1623:   TSRHSFunction  rhsfunction;


1633:   TSGetDM(ts,&dm);
1634:   DMTSGetI2Function(dm,&I2Function,&ctx);
1635:   DMTSGetRHSFunction(dm,&rhsfunction,NULL);

1637:   if (!I2Function) {
1638:     TSComputeIFunction(ts,t,U,A,F,PETSC_FALSE);
1639:     return(0);
1640:   }

1642:   PetscLogEventBegin(TS_FunctionEval,ts,U,V,F);

1644:   PetscStackPush("TS user implicit function");
1645:   I2Function(ts,t,U,V,A,F,ctx);
1646:   PetscStackPop;

1648:   if (rhsfunction) {
1649:     Vec Frhs;
1650:     TSGetRHSVec_Private(ts,&Frhs);
1651:     TSComputeRHSFunction(ts,t,U,Frhs);
1652:     VecAXPY(F,-1,Frhs);
1653:   }

1655:   PetscLogEventEnd(TS_FunctionEval,ts,U,V,F);
1656:   return(0);
1657: }

1661: /*@
1662:   TSComputeI2Jacobian - Evaluates the Jacobian of the DAE

1664:   Collective on TS and Vec

1666:   Input Parameters:
1667: + ts - the TS context
1668: . t - current timestep
1669: . U - state vector
1670: . V - time derivative of state vector
1671: . A - second time derivative of state vector
1672: . shiftV - shift to apply, see note below
1673: - shiftA - shift to apply, see note below

1675:   Output Parameters:
1676: + J - Jacobian matrix
1677: - P - optional preconditioning matrix

1679:   Notes:
1680:   If F(t,U,V,A)=0 is the DAE, the required Jacobian is

1682:   dF/dU + shiftV*dF/dV + shiftA*dF/dA

1684:   Most users should not need to explicitly call this routine, as it
1685:   is used internally within the nonlinear solvers.

1687:   Level: developer

1689: .keywords: TS, compute, Jacobian, matrix

1691: .seealso:  TSSetI2Jacobian()
1692: @*/
1693: PetscErrorCode TSComputeI2Jacobian(TS ts,PetscReal t,Vec U,Vec V,Vec A,PetscReal shiftV,PetscReal shiftA,Mat J,Mat P)
1694: {
1695:   DM             dm;
1696:   TSI2Jacobian   I2Jacobian;
1697:   void           *ctx;
1698:   TSRHSJacobian  rhsjacobian;


1709:   TSGetDM(ts,&dm);
1710:   DMTSGetI2Jacobian(dm,&I2Jacobian,&ctx);
1711:   DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);

1713:   if (!I2Jacobian) {
1714:     TSComputeIJacobian(ts,t,U,A,shiftA,J,P,PETSC_FALSE);
1715:     return(0);
1716:   }

1718:   PetscLogEventBegin(TS_JacobianEval,ts,U,J,P);

1720:   PetscStackPush("TS user implicit Jacobian");
1721:   I2Jacobian(ts,t,U,V,A,shiftV,shiftA,J,P,ctx);
1722:   PetscStackPop;

1724:   if (rhsjacobian) {
1725:     Mat Jrhs,Prhs; MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
1726:     TSGetRHSMats_Private(ts,&Jrhs,&Prhs);
1727:     TSComputeRHSJacobian(ts,t,U,Jrhs,Prhs);
1728:     MatAXPY(J,-1,Jrhs,axpy);
1729:     if (P != J) {MatAXPY(P,-1,Prhs,axpy);}
1730:   }

1732:   PetscLogEventEnd(TS_JacobianEval,ts,U,J,P);
1733:   return(0);
1734: }

1738: /*@
1739:    TS2SetSolution - Sets the initial solution and time derivative vectors
1740:    for use by the TS routines handling second order equations.

1742:    Logically Collective on TS and Vec

1744:    Input Parameters:
1745: +  ts - the TS context obtained from TSCreate()
1746: .  u - the solution vector
1747: -  v - the time derivative vector

1749:    Level: beginner

1751: .keywords: TS, timestep, set, solution, initial conditions
1752: @*/
1753: PetscErrorCode  TS2SetSolution(TS ts,Vec u,Vec v)
1754: {

1761:   TSSetSolution(ts,u);
1762:   PetscObjectReference((PetscObject)v);
1763:   VecDestroy(&ts->vec_dot);
1764:   ts->vec_dot = v;
1765:   return(0);
1766: }

1770: /*@
1771:    TS2GetSolution - Returns the solution and time derivative at the present timestep
1772:    for second order equations. It is valid to call this routine inside the function
1773:    that you are evaluating in order to move to the new timestep. This vector not
1774:    changed until the solution at the next timestep has been calculated.

1776:    Not Collective, but Vec returned is parallel if TS is parallel

1778:    Input Parameter:
1779: .  ts - the TS context obtained from TSCreate()

1781:    Output Parameter:
1782: +  u - the vector containing the solution
1783: -  v - the vector containing the time derivative

1785:    Level: intermediate

1787: .seealso: TS2SetSolution(), TSGetTimeStep(), TSGetTime()

1789: .keywords: TS, timestep, get, solution
1790: @*/
1791: PetscErrorCode  TS2GetSolution(TS ts,Vec *u,Vec *v)
1792: {
1797:   if (u) *u = ts->vec_sol;
1798:   if (v) *v = ts->vec_dot;
1799:   return(0);
1800: }

1804: /*@C
1805:   TSLoad - Loads a KSP that has been stored in binary  with KSPView().

1807:   Collective on PetscViewer

1809:   Input Parameters:
1810: + newdm - the newly loaded TS, this needs to have been created with TSCreate() or
1811:            some related function before a call to TSLoad().
1812: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()

1814:    Level: intermediate

1816:   Notes:
1817:    The type is determined by the data in the file, any type set into the TS before this call is ignored.

1819:   Notes for advanced users:
1820:   Most users should not need to know the details of the binary storage
1821:   format, since TSLoad() and TSView() completely hide these details.
1822:   But for anyone who's interested, the standard binary matrix storage
1823:   format is
1824: .vb
1825:      has not yet been determined
1826: .ve

1828: .seealso: PetscViewerBinaryOpen(), TSView(), MatLoad(), VecLoad()
1829: @*/
1830: PetscErrorCode  TSLoad(TS ts, PetscViewer viewer)
1831: {
1833:   PetscBool      isbinary;
1834:   PetscInt       classid;
1835:   char           type[256];
1836:   DMTS           sdm;
1837:   DM             dm;

1842:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1843:   if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

1845:   PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1846:   if (classid != TS_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Not TS next in file");
1847:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1848:   TSSetType(ts, type);
1849:   if (ts->ops->load) {
1850:     (*ts->ops->load)(ts,viewer);
1851:   }
1852:   DMCreate(PetscObjectComm((PetscObject)ts),&dm);
1853:   DMLoad(dm,viewer);
1854:   TSSetDM(ts,dm);
1855:   DMCreateGlobalVector(ts->dm,&ts->vec_sol);
1856:   VecLoad(ts->vec_sol,viewer);
1857:   DMGetDMTS(ts->dm,&sdm);
1858:   DMTSLoad(sdm,viewer);
1859:   return(0);
1860: }

1862: #include <petscdraw.h>
1863: #if defined(PETSC_HAVE_SAWS)
1864: #include <petscviewersaws.h>
1865: #endif
1868: /*@C
1869:     TSView - Prints the TS data structure.

1871:     Collective on TS

1873:     Input Parameters:
1874: +   ts - the TS context obtained from TSCreate()
1875: -   viewer - visualization context

1877:     Options Database Key:
1878: .   -ts_view - calls TSView() at end of TSStep()

1880:     Notes:
1881:     The available visualization contexts include
1882: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1883: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1884:          output where only the first processor opens
1885:          the file.  All other processors send their
1886:          data to the first processor to print.

1888:     The user can open an alternative visualization context with
1889:     PetscViewerASCIIOpen() - output to a specified file.

1891:     Level: beginner

1893: .keywords: TS, timestep, view

1895: .seealso: PetscViewerASCIIOpen()
1896: @*/
1897: PetscErrorCode  TSView(TS ts,PetscViewer viewer)
1898: {
1900:   TSType         type;
1901:   PetscBool      iascii,isstring,isundials,isbinary,isdraw;
1902:   DMTS           sdm;
1903: #if defined(PETSC_HAVE_SAWS)
1904:   PetscBool      issaws;
1905: #endif

1909:   if (!viewer) {
1910:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts),&viewer);
1911:   }

1915:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1916:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1917:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1918:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1919: #if defined(PETSC_HAVE_SAWS)
1920:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1921: #endif
1922:   if (iascii) {
1923:     PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer);
1924:     PetscViewerASCIIPrintf(viewer,"  maximum steps=%D\n",ts->max_steps);
1925:     PetscViewerASCIIPrintf(viewer,"  maximum time=%g\n",(double)ts->max_time);
1926:     if (ts->problem_type == TS_NONLINEAR) {
1927:       PetscViewerASCIIPrintf(viewer,"  total number of nonlinear solver iterations=%D\n",ts->snes_its);
1928:       PetscViewerASCIIPrintf(viewer,"  total number of nonlinear solve failures=%D\n",ts->num_snes_failures);
1929:     }
1930:     PetscViewerASCIIPrintf(viewer,"  total number of linear solver iterations=%D\n",ts->ksp_its);
1931:     PetscViewerASCIIPrintf(viewer,"  total number of rejected steps=%D\n",ts->reject);
1932:     DMGetDMTS(ts->dm,&sdm);
1933:     DMTSView(sdm,viewer);
1934:     if (ts->ops->view) {
1935:       PetscViewerASCIIPushTab(viewer);
1936:       (*ts->ops->view)(ts,viewer);
1937:       PetscViewerASCIIPopTab(viewer);
1938:     }
1939:   } else if (isstring) {
1940:     TSGetType(ts,&type);
1941:     PetscViewerStringSPrintf(viewer," %-7.7s",type);
1942:   } else if (isbinary) {
1943:     PetscInt    classid = TS_FILE_CLASSID;
1944:     MPI_Comm    comm;
1945:     PetscMPIInt rank;
1946:     char        type[256];

1948:     PetscObjectGetComm((PetscObject)ts,&comm);
1949:     MPI_Comm_rank(comm,&rank);
1950:     if (!rank) {
1951:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1952:       PetscStrncpy(type,((PetscObject)ts)->type_name,256);
1953:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1954:     }
1955:     if (ts->ops->view) {
1956:       (*ts->ops->view)(ts,viewer);
1957:     }
1958:     DMView(ts->dm,viewer);
1959:     VecView(ts->vec_sol,viewer);
1960:     DMGetDMTS(ts->dm,&sdm);
1961:     DMTSView(sdm,viewer);
1962:   } else if (isdraw) {
1963:     PetscDraw draw;
1964:     char      str[36];
1965:     PetscReal x,y,bottom,h;

1967:     PetscViewerDrawGetDraw(viewer,0,&draw);
1968:     PetscDrawGetCurrentPoint(draw,&x,&y);
1969:     PetscStrcpy(str,"TS: ");
1970:     PetscStrcat(str,((PetscObject)ts)->type_name);
1971:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,str,NULL,&h);
1972:     bottom = y - h;
1973:     PetscDrawPushCurrentPoint(draw,x,bottom);
1974:     if (ts->ops->view) {
1975:       (*ts->ops->view)(ts,viewer);
1976:     }
1977:     PetscDrawPopCurrentPoint(draw);
1978: #if defined(PETSC_HAVE_SAWS)
1979:   } else if (issaws) {
1980:     PetscMPIInt rank;
1981:     const char  *name;

1983:     PetscObjectGetName((PetscObject)ts,&name);
1984:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1985:     if (!((PetscObject)ts)->amsmem && !rank) {
1986:       char       dir[1024];

1988:       PetscObjectViewSAWs((PetscObject)ts,viewer);
1989:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time_step",name);
1990:       PetscStackCallSAWs(SAWs_Register,(dir,&ts->steps,1,SAWs_READ,SAWs_INT));
1991:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time",name);
1992:       PetscStackCallSAWs(SAWs_Register,(dir,&ts->ptime,1,SAWs_READ,SAWs_DOUBLE));
1993:     }
1994:     if (ts->ops->view) {
1995:       (*ts->ops->view)(ts,viewer);
1996:     }
1997: #endif
1998:   }

2000:   PetscViewerASCIIPushTab(viewer);
2001:   PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);
2002:   PetscViewerASCIIPopTab(viewer);
2003:   return(0);
2004: }


2009: /*@
2010:    TSSetApplicationContext - Sets an optional user-defined context for
2011:    the timesteppers.

2013:    Logically Collective on TS

2015:    Input Parameters:
2016: +  ts - the TS context obtained from TSCreate()
2017: -  usrP - optional user context

2019:    Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
2020:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

2022:    Level: intermediate

2024: .keywords: TS, timestep, set, application, context

2026: .seealso: TSGetApplicationContext()
2027: @*/
2028: PetscErrorCode  TSSetApplicationContext(TS ts,void *usrP)
2029: {
2032:   ts->user = usrP;
2033:   return(0);
2034: }

2038: /*@
2039:     TSGetApplicationContext - Gets the user-defined context for the
2040:     timestepper.

2042:     Not Collective

2044:     Input Parameter:
2045: .   ts - the TS context obtained from TSCreate()

2047:     Output Parameter:
2048: .   usrP - user context

2050:    Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
2051:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

2053:     Level: intermediate

2055: .keywords: TS, timestep, get, application, context

2057: .seealso: TSSetApplicationContext()
2058: @*/
2059: PetscErrorCode  TSGetApplicationContext(TS ts,void *usrP)
2060: {
2063:   *(void**)usrP = ts->user;
2064:   return(0);
2065: }

2069: /*@
2070:    TSGetTimeStepNumber - Gets the number of time steps completed.

2072:    Not Collective

2074:    Input Parameter:
2075: .  ts - the TS context obtained from TSCreate()

2077:    Output Parameter:
2078: .  iter - number of steps completed so far

2080:    Level: intermediate

2082: .keywords: TS, timestep, get, iteration, number
2083: .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSSetPostStep()
2084: @*/
2085: PetscErrorCode  TSGetTimeStepNumber(TS ts,PetscInt *iter)
2086: {
2090:   *iter = ts->steps;
2091:   return(0);
2092: }

2096: /*@
2097:    TSSetInitialTimeStep - Sets the initial timestep to be used,
2098:    as well as the initial time.

2100:    Logically Collective on TS

2102:    Input Parameters:
2103: +  ts - the TS context obtained from TSCreate()
2104: .  initial_time - the initial time
2105: -  time_step - the size of the timestep

2107:    Level: intermediate

2109: .seealso: TSSetTimeStep(), TSGetTimeStep()

2111: .keywords: TS, set, initial, timestep
2112: @*/
2113: PetscErrorCode  TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step)
2114: {

2119:   TSSetTimeStep(ts,time_step);
2120:   TSSetTime(ts,initial_time);
2121:   return(0);
2122: }

2126: /*@
2127:    TSSetTimeStep - Allows one to reset the timestep at any time,
2128:    useful for simple pseudo-timestepping codes.

2130:    Logically Collective on TS

2132:    Input Parameters:
2133: +  ts - the TS context obtained from TSCreate()
2134: -  time_step - the size of the timestep

2136:    Level: intermediate

2138: .seealso: TSSetInitialTimeStep(), TSGetTimeStep()

2140: .keywords: TS, set, timestep
2141: @*/
2142: PetscErrorCode  TSSetTimeStep(TS ts,PetscReal time_step)
2143: {
2147:   ts->time_step = time_step;
2148:   return(0);
2149: }

2153: /*@
2154:    TSSetExactFinalTime - Determines whether to adapt the final time step to
2155:      match the exact final time, interpolate solution to the exact final time,
2156:      or just return at the final time TS computed.

2158:   Logically Collective on TS

2160:    Input Parameter:
2161: +   ts - the time-step context
2162: -   eftopt - exact final time option

2164: $  TS_EXACTFINALTIME_STEPOVER    - Don't do anything if final time is exceeded
2165: $  TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2166: $  TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time

2168:    Options Database:
2169: .   -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step at runtime

2171:    Warning: If you use the option TS_EXACTFINALTIME_STEPOVER the solution may be at a very different time
2172:     then the final time you selected.

2174:    Level: beginner

2176: .seealso: TSExactFinalTimeOption
2177: @*/
2178: PetscErrorCode  TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt)
2179: {
2183:   ts->exact_final_time = eftopt;
2184:   return(0);
2185: }

2189: /*@
2190:    TSGetTimeStep - Gets the current timestep size.

2192:    Not Collective

2194:    Input Parameter:
2195: .  ts - the TS context obtained from TSCreate()

2197:    Output Parameter:
2198: .  dt - the current timestep size

2200:    Level: intermediate

2202: .seealso: TSSetInitialTimeStep(), TSGetTimeStep()

2204: .keywords: TS, get, timestep
2205: @*/
2206: PetscErrorCode  TSGetTimeStep(TS ts,PetscReal *dt)
2207: {
2211:   *dt = ts->time_step;
2212:   return(0);
2213: }

2217: /*@
2218:    TSGetSolution - Returns the solution at the present timestep. It
2219:    is valid to call this routine inside the function that you are evaluating
2220:    in order to move to the new timestep. This vector not changed until
2221:    the solution at the next timestep has been calculated.

2223:    Not Collective, but Vec returned is parallel if TS is parallel

2225:    Input Parameter:
2226: .  ts - the TS context obtained from TSCreate()

2228:    Output Parameter:
2229: .  v - the vector containing the solution

2231:    Note: If you used TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); this does not return the solution at the requested
2232:    final time. It returns the solution at the next timestep.

2234:    Level: intermediate

2236: .seealso: TSGetTimeStep(), TSGetTime(), TSGetSolveTime()

2238: .keywords: TS, timestep, get, solution
2239: @*/
2240: PetscErrorCode  TSGetSolution(TS ts,Vec *v)
2241: {
2245:   *v = ts->vec_sol;
2246:   return(0);
2247: }

2251: /*@
2252:    TSGetCostGradients - Returns the gradients from the TSAdjointSolve()

2254:    Not Collective, but Vec returned is parallel if TS is parallel

2256:    Input Parameter:
2257: .  ts - the TS context obtained from TSCreate()

2259:    Output Parameter:
2260: +  lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
2261: -  mu - vectors containing the gradients of the cost functions with respect to the problem parameters

2263:    Level: intermediate

2265: .seealso: TSGetTimeStep()

2267: .keywords: TS, timestep, get, sensitivity
2268: @*/
2269: PetscErrorCode  TSGetCostGradients(TS ts,PetscInt *numcost,Vec **lambda,Vec **mu)
2270: {
2273:   if (numcost) *numcost = ts->numcost;
2274:   if (lambda)  *lambda  = ts->vecs_sensi;
2275:   if (mu)      *mu      = ts->vecs_sensip;
2276:   return(0);
2277: }

2279: /* ----- Routines to initialize and destroy a timestepper ---- */
2282: /*@
2283:   TSSetProblemType - Sets the type of problem to be solved.

2285:   Not collective

2287:   Input Parameters:
2288: + ts   - The TS
2289: - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2290: .vb
2291:          U_t - A U = 0      (linear)
2292:          U_t - A(t) U = 0   (linear)
2293:          F(t,U,U_t) = 0     (nonlinear)
2294: .ve

2296:    Level: beginner

2298: .keywords: TS, problem type
2299: .seealso: TSSetUp(), TSProblemType, TS
2300: @*/
2301: PetscErrorCode  TSSetProblemType(TS ts, TSProblemType type)
2302: {

2307:   ts->problem_type = type;
2308:   if (type == TS_LINEAR) {
2309:     SNES snes;
2310:     TSGetSNES(ts,&snes);
2311:     SNESSetType(snes,SNESKSPONLY);
2312:   }
2313:   return(0);
2314: }

2318: /*@C
2319:   TSGetProblemType - Gets the type of problem to be solved.

2321:   Not collective

2323:   Input Parameter:
2324: . ts   - The TS

2326:   Output Parameter:
2327: . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2328: .vb
2329:          M U_t = A U
2330:          M(t) U_t = A(t) U
2331:          F(t,U,U_t)
2332: .ve

2334:    Level: beginner

2336: .keywords: TS, problem type
2337: .seealso: TSSetUp(), TSProblemType, TS
2338: @*/
2339: PetscErrorCode  TSGetProblemType(TS ts, TSProblemType *type)
2340: {
2344:   *type = ts->problem_type;
2345:   return(0);
2346: }

2350: /*@
2351:    TSSetUp - Sets up the internal data structures for the later use
2352:    of a timestepper.

2354:    Collective on TS

2356:    Input Parameter:
2357: .  ts - the TS context obtained from TSCreate()

2359:    Notes:
2360:    For basic use of the TS solvers the user need not explicitly call
2361:    TSSetUp(), since these actions will automatically occur during
2362:    the call to TSStep().  However, if one wishes to control this
2363:    phase separately, TSSetUp() should be called after TSCreate()
2364:    and optional routines of the form TSSetXXX(), but before TSStep().

2366:    Level: advanced

2368: .keywords: TS, timestep, setup

2370: .seealso: TSCreate(), TSStep(), TSDestroy()
2371: @*/
2372: PetscErrorCode  TSSetUp(TS ts)
2373: {
2375:   DM             dm;
2376:   PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2377:   PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*);
2378:   TSIFunction    ifun;
2379:   TSIJacobian    ijac;
2380:   TSI2Jacobian   i2jac;
2381:   TSRHSJacobian  rhsjac;

2385:   if (ts->setupcalled) return(0);

2387:   ts->total_steps = 0;
2388:   if (!((PetscObject)ts)->type_name) {
2389:     TSGetIFunction(ts,NULL,&ifun,NULL);
2390:     TSSetType(ts,ifun ? TSBEULER : TSEULER);
2391:   }

2393:   if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first");

2395:   if (ts->rhsjacobian.reuse) {
2396:     Mat Amat,Pmat;
2397:     SNES snes;
2398:     TSGetSNES(ts,&snes);
2399:     SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);
2400:     /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2401:      * have displaced the RHS matrix */
2402:     if (Amat == ts->Arhs) {
2403:       MatDuplicate(ts->Arhs,MAT_DO_NOT_COPY_VALUES,&Amat);
2404:       SNESSetJacobian(snes,Amat,NULL,NULL,NULL);
2405:       MatDestroy(&Amat);
2406:     }
2407:     if (Pmat == ts->Brhs) {
2408:       MatDuplicate(ts->Brhs,MAT_DO_NOT_COPY_VALUES,&Pmat);
2409:       SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);
2410:       MatDestroy(&Pmat);
2411:     }
2412:   }
2413:   if (ts->ops->setup) {
2414:     (*ts->ops->setup)(ts);
2415:   }

2417:   /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2418:      to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2419:    */
2420:   TSGetDM(ts,&dm);
2421:   DMSNESGetFunction(dm,&func,NULL);
2422:   if (!func) {
2423:     DMSNESSetFunction(dm,SNESTSFormFunction,ts);
2424:   }
2425:   /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2426:      Otherwise, the SNES will use coloring internally to form the Jacobian.
2427:    */
2428:   DMSNESGetJacobian(dm,&jac,NULL);
2429:   DMTSGetIJacobian(dm,&ijac,NULL);
2430:   DMTSGetI2Jacobian(dm,&i2jac,NULL);
2431:   DMTSGetRHSJacobian(dm,&rhsjac,NULL);
2432:   if (!jac && (ijac || i2jac || rhsjac)) {
2433:     DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);
2434:   }
2435:   ts->setupcalled = PETSC_TRUE;
2436:   return(0);
2437: }

2441: /*@
2442:    TSAdjointSetUp - Sets up the internal data structures for the later use
2443:    of an adjoint solver

2445:    Collective on TS

2447:    Input Parameter:
2448: .  ts - the TS context obtained from TSCreate()

2450:    Level: advanced

2452: .keywords: TS, timestep, setup

2454: .seealso: TSCreate(), TSAdjointStep(), TSSetCostGradients()
2455: @*/
2456: PetscErrorCode  TSAdjointSetUp(TS ts)
2457: {

2462:   if (ts->adjointsetupcalled) return(0);
2463:   if (!ts->vecs_sensi) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetCostGradients() first");

2465:   if (ts->vec_costintegral) { /* if there is integral in the cost function*/
2466:     VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&ts->vecs_drdy);
2467:     if (ts->vecs_sensip){
2468:       VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&ts->vecs_drdp);
2469:     }
2470:   }

2472:   if (ts->ops->adjointsetup) {
2473:     (*ts->ops->adjointsetup)(ts);
2474:   }
2475:   ts->adjointsetupcalled = PETSC_TRUE;
2476:   return(0);
2477: }

2481: /*@
2482:    TSReset - Resets a TS context and removes any allocated Vecs and Mats.

2484:    Collective on TS

2486:    Input Parameter:
2487: .  ts - the TS context obtained from TSCreate()

2489:    Level: beginner

2491: .keywords: TS, timestep, reset

2493: .seealso: TSCreate(), TSSetup(), TSDestroy()
2494: @*/
2495: PetscErrorCode  TSReset(TS ts)
2496: {


2502:   if (ts->ops->reset) {
2503:     (*ts->ops->reset)(ts);
2504:   }
2505:   if (ts->snes) {SNESReset(ts->snes);}
2506:   if (ts->adapt) {TSAdaptReset(ts->adapt);}

2508:   MatDestroy(&ts->Arhs);
2509:   MatDestroy(&ts->Brhs);
2510:   VecDestroy(&ts->Frhs);
2511:   VecDestroy(&ts->vec_sol);
2512:   VecDestroy(&ts->vec_dot);
2513:   VecDestroy(&ts->vatol);
2514:   VecDestroy(&ts->vrtol);
2515:   VecDestroyVecs(ts->nwork,&ts->work);

2517:  if (ts->vec_costintegral) {
2518:     VecDestroyVecs(ts->numcost,&ts->vecs_drdy);
2519:     if (ts->vecs_drdp){
2520:       VecDestroyVecs(ts->numcost,&ts->vecs_drdp);
2521:     }
2522:   }
2523:   ts->vecs_sensi  = NULL;
2524:   ts->vecs_sensip = NULL;
2525:   MatDestroy(&ts->Jacp);
2526:   VecDestroy(&ts->vec_costintegral);
2527:   VecDestroy(&ts->vec_costintegrand);
2528:   ts->setupcalled = PETSC_FALSE;
2529:   return(0);
2530: }

2534: /*@
2535:    TSDestroy - Destroys the timestepper context that was created
2536:    with TSCreate().

2538:    Collective on TS

2540:    Input Parameter:
2541: .  ts - the TS context obtained from TSCreate()

2543:    Level: beginner

2545: .keywords: TS, timestepper, destroy

2547: .seealso: TSCreate(), TSSetUp(), TSSolve()
2548: @*/
2549: PetscErrorCode  TSDestroy(TS *ts)
2550: {

2554:   if (!*ts) return(0);
2556:   if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; return(0);}

2558:   TSReset((*ts));

2560:   /* if memory was published with SAWs then destroy it */
2561:   PetscObjectSAWsViewOff((PetscObject)*ts);
2562:   if ((*ts)->ops->destroy) {(*(*ts)->ops->destroy)((*ts));}

2564:   TSTrajectoryDestroy(&(*ts)->trajectory);

2566:   TSAdaptDestroy(&(*ts)->adapt);
2567:   TSEventDestroy(&(*ts)->event);

2569:   SNESDestroy(&(*ts)->snes);
2570:   DMDestroy(&(*ts)->dm);
2571:   TSMonitorCancel((*ts));
2572:   TSAdjointMonitorCancel((*ts));

2574:   PetscHeaderDestroy(ts);
2575:   return(0);
2576: }

2580: /*@
2581:    TSGetSNES - Returns the SNES (nonlinear solver) associated with
2582:    a TS (timestepper) context. Valid only for nonlinear problems.

2584:    Not Collective, but SNES is parallel if TS is parallel

2586:    Input Parameter:
2587: .  ts - the TS context obtained from TSCreate()

2589:    Output Parameter:
2590: .  snes - the nonlinear solver context

2592:    Notes:
2593:    The user can then directly manipulate the SNES context to set various
2594:    options, etc.  Likewise, the user can then extract and manipulate the
2595:    KSP, KSP, and PC contexts as well.

2597:    TSGetSNES() does not work for integrators that do not use SNES; in
2598:    this case TSGetSNES() returns NULL in snes.

2600:    Level: beginner

2602: .keywords: timestep, get, SNES
2603: @*/
2604: PetscErrorCode  TSGetSNES(TS ts,SNES *snes)
2605: {

2611:   if (!ts->snes) {
2612:     SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);
2613:     SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
2614:     PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->snes);
2615:     PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);
2616:     if (ts->dm) {SNESSetDM(ts->snes,ts->dm);}
2617:     if (ts->problem_type == TS_LINEAR) {
2618:       SNESSetType(ts->snes,SNESKSPONLY);
2619:     }
2620:   }
2621:   *snes = ts->snes;
2622:   return(0);
2623: }

2627: /*@
2628:    TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context

2630:    Collective

2632:    Input Parameter:
2633: +  ts - the TS context obtained from TSCreate()
2634: -  snes - the nonlinear solver context

2636:    Notes:
2637:    Most users should have the TS created by calling TSGetSNES()

2639:    Level: developer

2641: .keywords: timestep, set, SNES
2642: @*/
2643: PetscErrorCode TSSetSNES(TS ts,SNES snes)
2644: {
2646:   PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*);

2651:   PetscObjectReference((PetscObject)snes);
2652:   SNESDestroy(&ts->snes);

2654:   ts->snes = snes;

2656:   SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
2657:   SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);
2658:   if (func == SNESTSFormJacobian) {
2659:     SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);
2660:   }
2661:   return(0);
2662: }

2666: /*@
2667:    TSGetKSP - Returns the KSP (linear solver) associated with
2668:    a TS (timestepper) context.

2670:    Not Collective, but KSP is parallel if TS is parallel

2672:    Input Parameter:
2673: .  ts - the TS context obtained from TSCreate()

2675:    Output Parameter:
2676: .  ksp - the nonlinear solver context

2678:    Notes:
2679:    The user can then directly manipulate the KSP context to set various
2680:    options, etc.  Likewise, the user can then extract and manipulate the
2681:    KSP and PC contexts as well.

2683:    TSGetKSP() does not work for integrators that do not use KSP;
2684:    in this case TSGetKSP() returns NULL in ksp.

2686:    Level: beginner

2688: .keywords: timestep, get, KSP
2689: @*/
2690: PetscErrorCode  TSGetKSP(TS ts,KSP *ksp)
2691: {
2693:   SNES           snes;

2698:   if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first");
2699:   if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()");
2700:   TSGetSNES(ts,&snes);
2701:   SNESGetKSP(snes,ksp);
2702:   return(0);
2703: }

2705: /* ----------- Routines to set solver parameters ---------- */

2709: /*@
2710:    TSGetDuration - Gets the maximum number of timesteps to use and
2711:    maximum time for iteration.

2713:    Not Collective

2715:    Input Parameters:
2716: +  ts       - the TS context obtained from TSCreate()
2717: .  maxsteps - maximum number of iterations to use, or NULL
2718: -  maxtime  - final time to iterate to, or NULL

2720:    Level: intermediate

2722: .keywords: TS, timestep, get, maximum, iterations, time
2723: @*/
2724: PetscErrorCode  TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2725: {
2728:   if (maxsteps) {
2730:     *maxsteps = ts->max_steps;
2731:   }
2732:   if (maxtime) {
2734:     *maxtime = ts->max_time;
2735:   }
2736:   return(0);
2737: }

2741: /*@
2742:    TSSetDuration - Sets the maximum number of timesteps to use and
2743:    maximum time for iteration.

2745:    Logically Collective on TS

2747:    Input Parameters:
2748: +  ts - the TS context obtained from TSCreate()
2749: .  maxsteps - maximum number of iterations to use
2750: -  maxtime - final time to iterate to

2752:    Options Database Keys:
2753: .  -ts_max_steps <maxsteps> - Sets maxsteps
2754: .  -ts_final_time <maxtime> - Sets maxtime

2756:    Notes:
2757:    The default maximum number of iterations is 5000. Default time is 5.0

2759:    Level: intermediate

2761: .keywords: TS, timestep, set, maximum, iterations

2763: .seealso: TSSetExactFinalTime()
2764: @*/
2765: PetscErrorCode  TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime)
2766: {
2771:   if (maxsteps >= 0) ts->max_steps = maxsteps;
2772:   if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime;
2773:   return(0);
2774: }

2778: /*@
2779:    TSSetSolution - Sets the initial solution vector
2780:    for use by the TS routines.

2782:    Logically Collective on TS and Vec

2784:    Input Parameters:
2785: +  ts - the TS context obtained from TSCreate()
2786: -  u - the solution vector

2788:    Level: beginner

2790: .keywords: TS, timestep, set, solution, initial conditions
2791: @*/
2792: PetscErrorCode  TSSetSolution(TS ts,Vec u)
2793: {
2795:   DM             dm;

2800:   PetscObjectReference((PetscObject)u);
2801:   VecDestroy(&ts->vec_sol);
2802:   ts->vec_sol = u;

2804:   TSGetDM(ts,&dm);
2805:   DMShellSetGlobalVector(dm,u);
2806:   return(0);
2807: }

2811: /*@
2812:    TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time

2814:    Logically Collective on TS

2816:    Input Parameters:
2817: +  ts - the TS context obtained from TSCreate()
2818: .  steps - number of steps to use

2820:    Level: intermediate

2822:    Notes: Normally one does not call this and TSAdjointSolve() integrates back to the original timestep. One can call this
2823:           so as to integrate back to less than the original timestep

2825: .keywords: TS, timestep, set, maximum, iterations

2827: .seealso: TSSetExactFinalTime()
2828: @*/
2829: PetscErrorCode  TSAdjointSetSteps(TS ts,PetscInt steps)
2830: {
2834:   if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back a negative number of steps");
2835:   if (steps > ts->total_steps) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back more than the total number of forward steps");
2836:   ts->adjoint_max_steps = steps;
2837:   return(0);
2838: }

2842: /*@
2843:    TSSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial conditions and w.r.t. the problem parameters 
2844:       for use by the TSAdjoint routines.

2846:    Logically Collective on TS and Vec

2848:    Input Parameters:
2849: +  ts - the TS context obtained from TSCreate()
2850: .  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
2851: -  mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters

2853:    Level: beginner

2855:    Notes: the entries in these vectors must be correctly initialized with the values lamda_i = df/dy|finaltime  mu_i = df/dp|finaltime

2857: .keywords: TS, timestep, set, sensitivity, initial conditions
2858: @*/
2859: PetscErrorCode  TSSetCostGradients(TS ts,PetscInt numcost,Vec *lambda,Vec *mu)
2860: {
2864:   ts->vecs_sensi  = lambda;
2865:   ts->vecs_sensip = mu;
2866:   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");
2867:   ts->numcost  = numcost;
2868:   return(0);
2869: }

2873: /*@C
2874:   TSAdjointSetRHSJacobian - 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.

2876:   Logically Collective on TS

2878:   Input Parameters:
2879: + ts   - The TS context obtained from TSCreate()
2880: - func - The function

2882:   Calling sequence of func:
2883: $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx);
2884: +   t - current timestep
2885: .   y - input vector (current ODE solution)
2886: .   A - output matrix
2887: -   ctx - [optional] user-defined function context

2889:   Level: intermediate

2891:   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

2893: .keywords: TS, sensitivity
2894: .seealso:
2895: @*/
2896: PetscErrorCode  TSAdjointSetRHSJacobian(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx)
2897: {


2904:   ts->rhsjacobianp    = func;
2905:   ts->rhsjacobianpctx = ctx;
2906:   if(Amat) {
2907:     PetscObjectReference((PetscObject)Amat);
2908:     MatDestroy(&ts->Jacp);
2909:     ts->Jacp = Amat;
2910:   }
2911:   return(0);
2912: }

2916: /*@C
2917:   TSAdjointComputeRHSJacobian - Runs the user-defined Jacobian function.

2919:   Collective on TS

2921:   Input Parameters:
2922: . ts   - The TS context obtained from TSCreate()

2924:   Level: developer

2926: .keywords: TS, sensitivity
2927: .seealso: TSAdjointSetRHSJacobian()
2928: @*/
2929: PetscErrorCode  TSAdjointComputeRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat)
2930: {


2938:   PetscStackPush("TS user JacobianP function for sensitivity analysis");
2939:   (*ts->rhsjacobianp)(ts,t,X,Amat,ts->rhsjacobianpctx);
2940:   PetscStackPop;
2941:   return(0);
2942: }

2946: /*@C
2947:     TSSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions

2949:     Logically Collective on TS

2951:     Input Parameters:
2952: +   ts - the TS context obtained from TSCreate()
2953: .   numcost - number of gradients to be computed, this is the number of cost functions
2954: .   rf - routine for evaluating the integrand function
2955: .   drdyf - function that computes the gradients of the r's with respect to y,NULL if not a function y
2956: .   drdpf - function that computes the gradients of the r's with respect to p, NULL if not a function of p
2957: .   fwd ï flag indicating whether to evaluate cost integral in the forward run or the adjoint run
2958: -   ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

2960:     Calling sequence of rf:
2961: $     rf(TS ts,PetscReal t,Vec y,Vec f[],void *ctx);

2963: +   t - current timestep
2964: .   y - input vector
2965: .   f - function result; one vector entry for each cost function
2966: -   ctx - [optional] user-defined function context

2968:    Calling sequence of drdyf:
2969: $    PetscErroCode drdyf(TS ts,PetscReal t,Vec y,Vec *drdy,void *ctx);

2971:    Calling sequence of drdpf:
2972: $    PetscErroCode drdpf(TS ts,PetscReal t,Vec y,Vec *drdp,void *ctx);

2974:     Level: intermediate

2976:     Notes: For optimization there is generally a single cost function, numcost = 1. For sensitivities there may be multiple cost functions

2978: .keywords: TS, sensitivity analysis, timestep, set, quadrature, function

2980: .seealso: TSAdjointSetRHSJacobian(),TSGetCostGradients(), TSSetCostGradients()
2981: @*/
2982: PetscErrorCode  TSSetCostIntegrand(TS ts,PetscInt numcost,PetscErrorCode (*rf)(TS,PetscReal,Vec,Vec,void*),
2983:                                                           PetscErrorCode (*drdyf)(TS,PetscReal,Vec,Vec*,void*),
2984:                                                           PetscErrorCode (*drdpf)(TS,PetscReal,Vec,Vec*,void*),
2985:                                                           PetscBool fwd,void *ctx)
2986: {

2991:   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()");
2992:   if (!ts->numcost) ts->numcost=numcost;

2994:   ts->costintegralfwd  = fwd; /* Evaluate the cost integral in forward run if fwd is true */
2995:   VecCreateSeq(PETSC_COMM_SELF,numcost,&ts->vec_costintegral);
2996:   VecDuplicate(ts->vec_costintegral,&ts->vec_costintegrand);
2997:   ts->costintegrand    = rf;
2998:   ts->costintegrandctx = ctx;
2999:   ts->drdyfunction     = drdyf;
3000:   ts->drdpfunction     = drdpf;
3001:   return(0);
3002: }

3006: /*@
3007:    TSGetCostIntegral - Returns the values of the integral term in the cost functions.
3008:    It is valid to call the routine after a backward run.

3010:    Not Collective

3012:    Input Parameter:
3013: .  ts - the TS context obtained from TSCreate()

3015:    Output Parameter:
3016: .  v - the vector containing the integrals for each cost function

3018:    Level: intermediate

3020: .seealso: TSSetCostIntegrand()

3022: .keywords: TS, sensitivity analysis
3023: @*/
3024: PetscErrorCode  TSGetCostIntegral(TS ts,Vec *v)
3025: {
3029:   *v = ts->vec_costintegral;
3030:   return(0);
3031: }

3035: /*@
3036:    TSAdjointComputeCostIntegrand - Evaluates the integral function in the cost functions.

3038:    Input Parameters:
3039: +  ts - the TS context
3040: .  t - current time
3041: -  y - state vector, i.e. current solution

3043:    Output Parameter:
3044: .  q - vector of size numcost to hold the outputs

3046:    Note:
3047:    Most users should not need to explicitly call this routine, as it
3048:    is used internally within the sensitivity analysis context.

3050:    Level: developer

3052: .keywords: TS, compute

3054: .seealso: TSSetCostIntegrand()
3055: @*/
3056: PetscErrorCode TSAdjointComputeCostIntegrand(TS ts,PetscReal t,Vec y,Vec q)
3057: {


3065:   PetscLogEventBegin(TS_FunctionEval,ts,y,q,0);
3066:   if (ts->costintegrand) {
3067:     PetscStackPush("TS user integrand in the cost function");
3068:     (*ts->costintegrand)(ts,t,y,q,ts->costintegrandctx);
3069:     PetscStackPop;
3070:   } else {
3071:     VecZeroEntries(q);
3072:   }

3074:   PetscLogEventEnd(TS_FunctionEval,ts,y,q,0);
3075:   return(0);
3076: }

3080: /*@
3081:   TSAdjointComputeDRDYFunction - Runs the user-defined DRDY function.

3083:   Collective on TS

3085:   Input Parameters:
3086: . ts   - The TS context obtained from TSCreate()

3088:   Notes:
3089:   TSAdjointComputeDRDYFunction() is typically used for sensitivity implementation,
3090:   so most users would not generally call this routine themselves.

3092:   Level: developer

3094: .keywords: TS, sensitivity
3095: .seealso: TSAdjointComputeDRDYFunction()
3096: @*/
3097: PetscErrorCode  TSAdjointComputeDRDYFunction(TS ts,PetscReal t,Vec y,Vec *drdy)
3098: {


3105:   PetscStackPush("TS user DRDY function for sensitivity analysis");
3106:   (*ts->drdyfunction)(ts,t,y,drdy,ts->costintegrandctx);
3107:   PetscStackPop;
3108:   return(0);
3109: }

3113: /*@
3114:   TSAdjointComputeDRDPFunction - Runs the user-defined DRDP function.

3116:   Collective on TS

3118:   Input Parameters:
3119: . ts   - The TS context obtained from TSCreate()

3121:   Notes:
3122:   TSDRDPFunction() is typically used for sensitivity implementation,
3123:   so most users would not generally call this routine themselves.

3125:   Level: developer

3127: .keywords: TS, sensitivity
3128: .seealso: TSAdjointSetDRDPFunction()
3129: @*/
3130: PetscErrorCode  TSAdjointComputeDRDPFunction(TS ts,PetscReal t,Vec y,Vec *drdp)
3131: {


3138:   PetscStackPush("TS user DRDP function for sensitivity analysis");
3139:   (*ts->drdpfunction)(ts,t,y,drdp,ts->costintegrandctx);
3140:   PetscStackPop;
3141:   return(0);
3142: }

3146: /*@C
3147:   TSSetPreStep - Sets the general-purpose function
3148:   called once at the beginning of each time step.

3150:   Logically Collective on TS

3152:   Input Parameters:
3153: + ts   - The TS context obtained from TSCreate()
3154: - func - The function

3156:   Calling sequence of func:
3157: . func (TS ts);

3159:   Level: intermediate

3161: .keywords: TS, timestep
3162: .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep()
3163: @*/
3164: PetscErrorCode  TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
3165: {
3168:   ts->prestep = func;
3169:   return(0);
3170: }

3174: /*@
3175:   TSPreStep - Runs the user-defined pre-step function.

3177:   Collective on TS

3179:   Input Parameters:
3180: . ts   - The TS context obtained from TSCreate()

3182:   Notes:
3183:   TSPreStep() is typically used within time stepping implementations,
3184:   so most users would not generally call this routine themselves.

3186:   Level: developer

3188: .keywords: TS, timestep
3189: .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep()
3190: @*/
3191: PetscErrorCode  TSPreStep(TS ts)
3192: {

3197:   if (ts->prestep) {
3198:     PetscStackCallStandard((*ts->prestep),(ts));
3199:   }
3200:   return(0);
3201: }

3205: /*@C
3206:   TSSetPreStage - Sets the general-purpose function
3207:   called once at the beginning of each stage.

3209:   Logically Collective on TS

3211:   Input Parameters:
3212: + ts   - The TS context obtained from TSCreate()
3213: - func - The function

3215:   Calling sequence of func:
3216: . PetscErrorCode func(TS ts, PetscReal stagetime);

3218:   Level: intermediate

3220:   Note:
3221:   There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3222:   The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being
3223:   attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime().

3225: .keywords: TS, timestep
3226: .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3227: @*/
3228: PetscErrorCode  TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal))
3229: {
3232:   ts->prestage = func;
3233:   return(0);
3234: }

3238: /*@C
3239:   TSSetPostStage - Sets the general-purpose function
3240:   called once at the end of each stage.

3242:   Logically Collective on TS

3244:   Input Parameters:
3245: + ts   - The TS context obtained from TSCreate()
3246: - func - The function

3248:   Calling sequence of func:
3249: . PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y);

3251:   Level: intermediate

3253:   Note:
3254:   There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3255:   The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being
3256:   attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime().

3258: .keywords: TS, timestep
3259: .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3260: @*/
3261: PetscErrorCode  TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*))
3262: {
3265:   ts->poststage = func;
3266:   return(0);
3267: }

3271: /*@
3272:   TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage()

3274:   Collective on TS

3276:   Input Parameters:
3277: . ts          - The TS context obtained from TSCreate()
3278:   stagetime   - The absolute time of the current stage

3280:   Notes:
3281:   TSPreStage() is typically used within time stepping implementations,
3282:   most users would not generally call this routine themselves.

3284:   Level: developer

3286: .keywords: TS, timestep
3287: .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep()
3288: @*/
3289: PetscErrorCode  TSPreStage(TS ts, PetscReal stagetime)
3290: {

3295:   if (ts->prestage) {
3296:     PetscStackCallStandard((*ts->prestage),(ts,stagetime));
3297:   }
3298:   return(0);
3299: }

3303: /*@
3304:   TSPostStage - Runs the user-defined post-stage function set using TSSetPostStage()

3306:   Collective on TS

3308:   Input Parameters:
3309: . ts          - The TS context obtained from TSCreate()
3310:   stagetime   - The absolute time of the current stage
3311:   stageindex  - Stage number
3312:   Y           - Array of vectors (of size = total number
3313:                 of stages) with the stage solutions

3315:   Notes:
3316:   TSPostStage() is typically used within time stepping implementations,
3317:   most users would not generally call this routine themselves.

3319:   Level: developer

3321: .keywords: TS, timestep
3322: .seealso: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep()
3323: @*/
3324: PetscErrorCode  TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3325: {

3330:   if (ts->poststage) {
3331:     PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y));
3332:   }
3333:   return(0);
3334: }

3338: /*@C
3339:   TSSetPostStep - Sets the general-purpose function
3340:   called once at the end of each time step.

3342:   Logically Collective on TS

3344:   Input Parameters:
3345: + ts   - The TS context obtained from TSCreate()
3346: - func - The function

3348:   Calling sequence of func:
3349: $ func (TS ts);

3351:   Level: intermediate

3353: .keywords: TS, timestep
3354: .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime()
3355: @*/
3356: PetscErrorCode  TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
3357: {
3360:   ts->poststep = func;
3361:   return(0);
3362: }

3366: /*@
3367:   TSPostStep - Runs the user-defined post-step function.

3369:   Collective on TS

3371:   Input Parameters:
3372: . ts   - The TS context obtained from TSCreate()

3374:   Notes:
3375:   TSPostStep() is typically used within time stepping implementations,
3376:   so most users would not generally call this routine themselves.

3378:   Level: developer

3380: .keywords: TS, timestep
3381: @*/
3382: PetscErrorCode  TSPostStep(TS ts)
3383: {

3388:   if (ts->poststep) {
3389:     PetscStackCallStandard((*ts->poststep),(ts));
3390:   }
3391:   return(0);
3392: }

3394: /* ------------ Routines to set performance monitoring options ----------- */

3398: /*@C
3399:    TSMonitorSet - Sets an ADDITIONAL function that is to be used at every
3400:    timestep to display the iteration's  progress.

3402:    Logically Collective on TS

3404:    Input Parameters:
3405: +  ts - the TS context obtained from TSCreate()
3406: .  monitor - monitoring routine
3407: .  mctx - [optional] user-defined context for private data for the
3408:              monitor routine (use NULL if no context is desired)
3409: -  monitordestroy - [optional] routine that frees monitor context
3410:           (may be NULL)

3412:    Calling sequence of monitor:
3413: $    int monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)

3415: +    ts - the TS context
3416: .    steps - iteration number (after the final time step the monitor routine may be called with a step of -1, this indicates the solution has been interpolated to this time)
3417: .    time - current time
3418: .    u - current iterate
3419: -    mctx - [optional] monitoring context

3421:    Notes:
3422:    This routine adds an additional monitor to the list of monitors that
3423:    already has been loaded.

3425:    Fortran notes: Only a single monitor function can be set for each TS object

3427:    Level: intermediate

3429: .keywords: TS, timestep, set, monitor

3431: .seealso: TSMonitorDefault(), TSMonitorCancel()
3432: @*/
3433: PetscErrorCode  TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**))
3434: {
3436:   PetscInt       i;
3437:   PetscBool      identical;
3438: 
3441:   for (i=0; i<ts->numbermonitors;i++) {
3442:     PetscMonitorCompare((PetscErrorCode (*)(void))monitor,mctx,mdestroy,(PetscErrorCode (*)(void))ts->monitor[i],ts->monitorcontext[i],ts->monitordestroy[i],&identical);
3443:     if (identical) return(0);
3444:   }
3445:   if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3446:   ts->monitor[ts->numbermonitors]          = monitor;
3447:   ts->monitordestroy[ts->numbermonitors]   = mdestroy;
3448:   ts->monitorcontext[ts->numbermonitors++] = (void*)mctx;
3449:   return(0);
3450: }

3454: /*@C
3455:    TSMonitorCancel - Clears all the monitors that have been set on a time-step object.

3457:    Logically Collective on TS

3459:    Input Parameters:
3460: .  ts - the TS context obtained from TSCreate()

3462:    Notes:
3463:    There is no way to remove a single, specific monitor.

3465:    Level: intermediate

3467: .keywords: TS, timestep, set, monitor

3469: .seealso: TSMonitorDefault(), TSMonitorSet()
3470: @*/
3471: PetscErrorCode  TSMonitorCancel(TS ts)
3472: {
3474:   PetscInt       i;

3478:   for (i=0; i<ts->numbermonitors; i++) {
3479:     if (ts->monitordestroy[i]) {
3480:       (*ts->monitordestroy[i])(&ts->monitorcontext[i]);
3481:     }
3482:   }
3483:   ts->numbermonitors = 0;
3484:   return(0);
3485: }

3489: /*@C
3490:    TSMonitorDefault - The Default monitor, prints the timestep and time for each step

3492:    Level: intermediate

3494: .keywords: TS, set, monitor

3496: .seealso:  TSMonitorSet()
3497: @*/
3498: PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscViewerAndFormat *vf)
3499: {
3501:   PetscViewer    viewer =  vf->viewer;
3502:   PetscBool      iascii,ibinary;

3506:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
3507:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
3508:   PetscViewerPushFormat(viewer,vf->format);
3509:   if (iascii) {
3510:     PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
3511:     if (step == -1){ /* this indicates it is an interpolated solution */
3512:       PetscViewerASCIIPrintf(viewer,"Interpolated solution at time %g between steps %D and %D\n",(double)ptime,ts->steps-1,ts->steps);
3513:     } else {
3514:       PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
3515:     }
3516:     PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
3517:   } else if (ibinary) {
3518:     PetscMPIInt rank;
3519:     MPI_Comm_rank(PetscObjectComm((PetscObject)viewer),&rank);
3520:     if (!rank) {
3521:       PetscRealView(1,&ptime,viewer);
3522:     } else {
3523:       PetscRealView(0,&ptime,viewer);
3524:     }
3525:   }
3526:   PetscViewerPopFormat(viewer);
3527:   return(0);
3528: }

3532: /*@C
3533:    TSAdjointMonitorSet - Sets an ADDITIONAL function that is to be used at every
3534:    timestep to display the iteration's  progress.

3536:    Logically Collective on TS

3538:    Input Parameters:
3539: +  ts - the TS context obtained from TSCreate()
3540: .  adjointmonitor - monitoring routine
3541: .  adjointmctx - [optional] user-defined context for private data for the
3542:              monitor routine (use NULL if no context is desired)
3543: -  adjointmonitordestroy - [optional] routine that frees monitor context
3544:           (may be NULL)

3546:    Calling sequence of monitor:
3547: $    int adjointmonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *adjointmctx)

3549: +    ts - the TS context
3550: .    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
3551:                                been interpolated to)
3552: .    time - current time
3553: .    u - current iterate
3554: .    numcost - number of cost functionos
3555: .    lambda - sensitivities to initial conditions
3556: .    mu - sensitivities to parameters
3557: -    adjointmctx - [optional] adjoint monitoring context

3559:    Notes:
3560:    This routine adds an additional monitor to the list of monitors that
3561:    already has been loaded.

3563:    Fortran notes: Only a single monitor function can be set for each TS object

3565:    Level: intermediate

3567: .keywords: TS, timestep, set, adjoint, monitor

3569: .seealso: TSAdjointMonitorCancel()
3570: @*/
3571: PetscErrorCode  TSAdjointMonitorSet(TS ts,PetscErrorCode (*adjointmonitor)(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*),void *adjointmctx,PetscErrorCode (*adjointmdestroy)(void**))
3572: {
3574:   PetscInt       i;
3575:   PetscBool      identical;

3579:   for (i=0; i<ts->numbermonitors;i++) {
3580:     PetscMonitorCompare((PetscErrorCode (*)(void))adjointmonitor,adjointmctx,adjointmdestroy,(PetscErrorCode (*)(void))ts->adjointmonitor[i],ts->adjointmonitorcontext[i],ts->adjointmonitordestroy[i],&identical);
3581:     if (identical) return(0);
3582:   }
3583:   if (ts->numberadjointmonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many adjoint monitors set");
3584:   ts->adjointmonitor[ts->numberadjointmonitors]          = adjointmonitor;
3585:   ts->adjointmonitordestroy[ts->numberadjointmonitors]   = adjointmdestroy;
3586:   ts->adjointmonitorcontext[ts->numberadjointmonitors++] = (void*)adjointmctx;
3587:   return(0);
3588: }

3592: /*@C
3593:    TSAdjointMonitorCancel - Clears all the adjoint monitors that have been set on a time-step object.

3595:    Logically Collective on TS

3597:    Input Parameters:
3598: .  ts - the TS context obtained from TSCreate()

3600:    Notes:
3601:    There is no way to remove a single, specific monitor.

3603:    Level: intermediate

3605: .keywords: TS, timestep, set, adjoint, monitor

3607: .seealso: TSAdjointMonitorSet()
3608: @*/
3609: PetscErrorCode  TSAdjointMonitorCancel(TS ts)
3610: {
3612:   PetscInt       i;

3616:   for (i=0; i<ts->numberadjointmonitors; i++) {
3617:     if (ts->adjointmonitordestroy[i]) {
3618:       (*ts->adjointmonitordestroy[i])(&ts->adjointmonitorcontext[i]);
3619:     }
3620:   }
3621:   ts->numberadjointmonitors = 0;
3622:   return(0);
3623: }

3627: /*@C
3628:    TSAdjointMonitorDefault - the default monitor of adjoint computations

3630:    Level: intermediate

3632: .keywords: TS, set, monitor

3634: .seealso: TSAdjointMonitorSet()
3635: @*/
3636: PetscErrorCode TSAdjointMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscInt numcost,Vec *lambda,Vec *mu,PetscViewerAndFormat *vf)
3637: {
3639:   PetscViewer    viewer = vf->viewer;

3643:   PetscViewerPushFormat(viewer,vf->format);
3644:   PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
3645:   PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
3646:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
3647:   PetscViewerPopFormat(viewer);
3648:   return(0);
3649: }

3653: /*@
3654:    TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval

3656:    Collective on TS

3658:    Input Argument:
3659: +  ts - time stepping context
3660: -  t - time to interpolate to

3662:    Output Argument:
3663: .  U - state at given time

3665:    Level: intermediate

3667:    Developer Notes:
3668:    TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.

3670: .keywords: TS, set

3672: .seealso: TSSetExactFinalTime(), TSSolve()
3673: @*/
3674: PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U)
3675: {

3681:   if (t < ts->ptime_prev || t > ts->ptime) SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Requested time %g not in last time steps [%g,%g]",t,(double)ts->ptime_prev,(double)ts->ptime);
3682:   if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name);
3683:   (*ts->ops->interpolate)(ts,t,U);
3684:   return(0);
3685: }

3689: /*@
3690:    TSStep - Steps one time step

3692:    Collective on TS

3694:    Input Parameter:
3695: .  ts - the TS context obtained from TSCreate()

3697:    Level: developer

3699:    Notes:
3700:    The public interface for the ODE/DAE solvers is TSSolve(), you should almost for sure be using that routine and not this routine.

3702:    The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may
3703:    be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.

3705:    This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the
3706:    time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep.

3708: .keywords: TS, timestep, solve

3710: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate()
3711: @*/
3712: PetscErrorCode  TSStep(TS ts)
3713: {
3714:   PetscErrorCode   ierr;
3715:   static PetscBool cite = PETSC_FALSE;
3716:   PetscReal        ptime;

3720:   PetscCitationsRegister("@techreport{tspaper,\n"
3721:                                 "  title       = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3722:                                 "  author      = {Shrirang Abhyankar and Jed Brown and Emil Constantinescu and Debojyoti Ghosh and Barry F. Smith},\n"
3723:                                 "  type        = {Preprint},\n"
3724:                                 "  number      = {ANL/MCS-P5061-0114},\n"
3725:                                 "  institution = {Argonne National Laboratory},\n"
3726:                                 "  year        = {2014}\n}\n",&cite);

3728:   TSSetUp(ts);
3729:   TSTrajectorySetUp(ts->trajectory,ts);

3731:   if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"You must call TSSetExactFinalTime() or use -ts_exact_final_time <stepover,interpolate,matchstep> before calling TSStep()");
3732:   if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP && !ts->adapt) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Since TS is not adaptive you cannot use TS_EXACTFINALTIME_MATCHSTEP, suggest TS_EXACTFINALTIME_INTERPOLATE");

3734:   if (!ts->steps) ts->ptime_prev = ts->ptime;
3735:   ts->reason = TS_CONVERGED_ITERATING;
3736:   ptime = ts->ptime; ts->ptime_prev_rollback = ts->ptime_prev;
3737:   if (!ts->ops->step) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSStep not implemented for type '%s'",((PetscObject)ts)->type_name);
3738:   PetscLogEventBegin(TS_Step,ts,0,0,0);
3739:   (*ts->ops->step)(ts);
3740:   PetscLogEventEnd(TS_Step,ts,0,0,0);
3741:   ts->ptime_prev = ptime;
3742:   ts->steps++; ts->total_steps++;
3743:   ts->steprollback = PETSC_FALSE;
3744:   ts->steprestart  = PETSC_FALSE;

3746:   if (ts->reason < 0) {
3747:     if (ts->errorifstepfailed) {
3748:       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]);
3749:       else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
3750:     }
3751:   } else if (!ts->reason) {
3752:     if (ts->steps >= ts->max_steps)     ts->reason = TS_CONVERGED_ITS;
3753:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3754:   }
3755:   return(0);
3756: }

3760: /*@
3761:    TSAdjointStep - Steps one time step backward in the adjoint run

3763:    Collective on TS

3765:    Input Parameter:
3766: .  ts - the TS context obtained from TSCreate()

3768:    Level: intermediate

3770: .keywords: TS, adjoint, step

3772: .seealso: TSAdjointSetUp(), TSAdjointSolve()
3773: @*/
3774: PetscErrorCode  TSAdjointStep(TS ts)
3775: {
3776:   DM               dm;
3777:   PetscErrorCode   ierr;

3781:   TSGetDM(ts,&dm);
3782:   TSAdjointSetUp(ts);

3784:   VecViewFromOptions(ts->vec_sol,(PetscObject)ts,"-ts_view_solution");

3786:   ts->reason = TS_CONVERGED_ITERATING;
3787:   ts->ptime_prev = ts->ptime;
3788:   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);
3789:   PetscLogEventBegin(TS_AdjointStep,ts,0,0,0);
3790:   (*ts->ops->adjointstep)(ts);
3791:   PetscLogEventEnd(TS_AdjointStep,ts,0,0,0);
3792:   ts->steps++; ts->total_steps--;

3794:   if (ts->reason < 0) {
3795:     if (ts->errorifstepfailed) {
3796:       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]);
3797:       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]);
3798:       else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
3799:     }
3800:   } else if (!ts->reason) {
3801:     if (ts->steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS;
3802:   }
3803:   return(0);
3804: }

3808: /*@
3809:    TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3810:    at the end of a time step with a given order of accuracy.

3812:    Collective on TS

3814:    Input Arguments:
3815: +  ts - time stepping context
3816: .  wnormtype - norm type, either NORM_2 or NORM_INFINITY
3817: -  order - optional, desired order for the error evaluation or PETSC_DECIDE

3819:    Output Arguments:
3820: +  order - optional, the actual order of the error evaluation
3821: -  wlte - the weighted local truncation error norm

3823:    Level: advanced

3825:    Notes:
3826:    If the timestepper cannot evaluate the error in a particular step
3827:    (eg. in the first step or restart steps after event handling),
3828:    this routine returns wlte=-1.0 .

3830: .seealso: TSStep(), TSAdapt, TSErrorWeightedNorm()
3831: @*/
3832: PetscErrorCode TSEvaluateWLTE(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
3833: {

3843:   if (wnormtype != NORM_2 && wnormtype != NORM_INFINITY) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
3844:   if (!ts->ops->evaluatewlte) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateWLTE not implemented for type '%s'",((PetscObject)ts)->type_name);
3845:   (*ts->ops->evaluatewlte)(ts,wnormtype,order,wlte);
3846:   return(0);
3847: }

3851: /*@
3852:    TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.

3854:    Collective on TS

3856:    Input Arguments:
3857: +  ts - time stepping context
3858: .  order - desired order of accuracy
3859: -  done - whether the step was evaluated at this order (pass NULL to generate an error if not available)

3861:    Output Arguments:
3862: .  U - state at the end of the current step

3864:    Level: advanced

3866:    Notes:
3867:    This function cannot be called until all stages have been evaluated.
3868:    It is normally called by adaptive controllers before a step has been accepted and may also be called by the user after TSStep() has returned.

3870: .seealso: TSStep(), TSAdapt
3871: @*/
3872: PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done)
3873: {

3880:   if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name);
3881:   (*ts->ops->evaluatestep)(ts,order,U,done);
3882:   return(0);
3883: }

3887: /*@
3888:  TSForwardCostIntegral - Evaluate the cost integral in the forward run.
3889:  
3890:  Collective on TS
3891:  
3892:  Input Arguments:
3893:  .  ts - time stepping context
3894:  
3895:  Level: advanced
3896:  
3897:  Notes:
3898:  This function cannot be called until TSStep() has been completed.
3899:  
3900:  .seealso: TSSolve(), TSAdjointCostIntegral()
3901:  @*/
3902: PetscErrorCode TSForwardCostIntegral(TS ts)
3903: {
3906:     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);
3907:     (*ts->ops->forwardintegral)(ts);
3908:     return(0);
3909: }

3913: /*@
3914:    TSSolve - Steps the requested number of timesteps.

3916:    Collective on TS

3918:    Input Parameter:
3919: +  ts - the TS context obtained from TSCreate()
3920: -  u - the solution vector  (can be null if TSSetSolution() was used and TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP) was not used,
3921:                              otherwise must contain the initial conditions and will contain the solution at the final requested time

3923:    Level: beginner

3925:    Notes:
3926:    The final time returned by this function may be different from the time of the internally
3927:    held state accessible by TSGetSolution() and TSGetTime() because the method may have
3928:    stepped over the final time.

3930: .keywords: TS, timestep, solve

3932: .seealso: TSCreate(), TSSetSolution(), TSStep(), TSGetTime(), TSGetSolveTime()
3933: @*/
3934: PetscErrorCode TSSolve(TS ts,Vec u)
3935: {
3936:   Vec               solution;
3937:   PetscErrorCode    ierr;


3943:   if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE) {   /* Need ts->vec_sol to be distinct so it is not overwritten when we interpolate at the end */
3945:     if (!ts->vec_sol || u == ts->vec_sol) {
3946:       VecDuplicate(u,&solution);
3947:       TSSetSolution(ts,solution);
3948:       VecDestroy(&solution); /* grant ownership */
3949:     }
3950:     VecCopy(u,ts->vec_sol);
3951:   } else if (u) {
3952:     TSSetSolution(ts,u);
3953:   }
3954:   TSSetUp(ts);
3955:   TSTrajectorySetUp(ts->trajectory,ts);

3957:   if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"You must call TSSetExactFinalTime() or use -ts_exact_final_time <stepover,interpolate,matchstep> before calling TSSolve()");
3958:   if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP && !ts->adapt) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Since TS is not adaptive you cannot use TS_EXACTFINALTIME_MATCHSTEP, suggest TS_EXACTFINALTIME_INTERPOLATE");

3960:   /* reset time step and iteration counters */
3961:   ts->steps             = 0;
3962:   ts->ksp_its           = 0;
3963:   ts->snes_its          = 0;
3964:   ts->num_snes_failures = 0;
3965:   ts->reject            = 0;
3966:   ts->reason            = TS_CONVERGED_ITERATING;

3968:   TSViewFromOptions(ts,NULL,"-ts_view_pre");

3970:   if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
3971:     (*ts->ops->solve)(ts);
3972:     if (u) {VecCopy(ts->vec_sol,u);}
3973:     ts->solvetime = ts->ptime;
3974:     solution = ts->vec_sol;
3975:   } else { /* Step the requested number of timesteps. */
3976:     if (ts->steps >= ts->max_steps)     ts->reason = TS_CONVERGED_ITS;
3977:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3978:     TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);
3979:     TSEventInitialize(ts->event,ts,ts->ptime,ts->vec_sol);
3980:     ts->steprollback = PETSC_FALSE;
3981:     ts->steprestart  = PETSC_TRUE;

3983:     while (!ts->reason) {
3984:       TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
3985:       if (!ts->steprollback) {
3986:         TSPreStep(ts);
3987:       }
3988:       TSStep(ts);
3989:       if (ts->vec_costintegral && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
3990:         TSForwardCostIntegral(ts);
3991:       }
3992:       TSEventHandler(ts); /* The right-hand side may be changed due to event. Be careful with Any computation using the RHS information after this point. */
3993:       if (!ts->steprollback) {
3994:         TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);
3995:         TSPostStep(ts);
3996:       }
3997:     }
3998:     TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);

4000:     if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4001:       TSInterpolate(ts,ts->max_time,u);
4002:       ts->solvetime = ts->max_time;
4003:       solution = u;
4004:       TSMonitor(ts,-1,ts->solvetime,solution);
4005:     } else {
4006:       if (u) {VecCopy(ts->vec_sol,u);}
4007:       ts->solvetime = ts->ptime;
4008:       solution = ts->vec_sol;
4009:     }
4010:   }

4012:   TSViewFromOptions(ts,NULL,"-ts_view");
4013:   VecViewFromOptions(solution,NULL,"-ts_view_solution");
4014:   PetscObjectSAWsBlock((PetscObject)ts);
4015:   if (ts->adjoint_solve) {
4016:     TSAdjointSolve(ts);
4017:   }
4018:   return(0);
4019: }

4023: /*@
4024:  TSAdjointCostIntegral - Evaluate the cost integral in the adjoint run.
4025:  
4026:  Collective on TS
4027:  
4028:  Input Arguments:
4029:  .  ts - time stepping context
4030:  
4031:  Level: advanced
4032:  
4033:  Notes:
4034:  This function cannot be called until TSAdjointStep() has been completed.
4035:  
4036:  .seealso: TSAdjointSolve(), TSAdjointStep
4037:  @*/
4038: PetscErrorCode TSAdjointCostIntegral(TS ts)
4039: {
4042:     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);
4043:     (*ts->ops->adjointintegral)(ts);
4044:     return(0);
4045: }

4049: /*@
4050:    TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE

4052:    Collective on TS

4054:    Input Parameter:
4055: .  ts - the TS context obtained from TSCreate()

4057:    Options Database:
4058: . -ts_adjoint_view_solution <viewerinfo> - views the first gradient with respect to the initial conditions

4060:    Level: intermediate

4062:    Notes:
4063:    This must be called after a call to TSSolve() that solves the forward problem

4065:    By default this will integrate back to the initial time, one can use TSAdjointSetSteps() to step back to a later time

4067: .keywords: TS, timestep, solve

4069: .seealso: TSCreate(), TSSetCostGradients(), TSSetSolution(), TSAdjointStep()
4070: @*/
4071: PetscErrorCode TSAdjointSolve(TS ts)
4072: {
4073:   PetscErrorCode    ierr;

4077:   TSAdjointSetUp(ts);

4079:   /* reset time step and iteration counters */
4080:   ts->steps             = 0;
4081:   ts->ksp_its           = 0;
4082:   ts->snes_its          = 0;
4083:   ts->num_snes_failures = 0;
4084:   ts->reject            = 0;
4085:   ts->reason            = TS_CONVERGED_ITERATING;

4087:   if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->total_steps;

4089:   if (ts->steps >= ts->adjoint_max_steps)     ts->reason = TS_CONVERGED_ITS;
4090:   while (!ts->reason) {
4091:     TSTrajectoryGet(ts->trajectory,ts,ts->total_steps,&ts->ptime);
4092:     TSAdjointMonitor(ts,ts->total_steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
4093:     TSAdjointEventHandler(ts);
4094:     TSAdjointStep(ts);
4095:     if (ts->vec_costintegral && !ts->costintegralfwd) {
4096:       TSAdjointCostIntegral(ts);
4097:     }
4098:   }
4099:   TSTrajectoryGet(ts->trajectory,ts,ts->total_steps,&ts->ptime);
4100:   TSAdjointMonitor(ts,ts->total_steps,ts->ptime,ts->vec_sol,ts->numcost,ts->vecs_sensi,ts->vecs_sensip);
4101:   ts->solvetime = ts->ptime;
4102:   TSTrajectoryViewFromOptions(ts->trajectory,NULL,"-ts_trajectory_view");
4103:   VecViewFromOptions(ts->vecs_sensi[0],(PetscObject) ts, "-ts_adjoint_view_solution");
4104:   return(0);
4105: }

4109: /*@C
4110:    TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet()

4112:    Collective on TS

4114:    Input Parameters:
4115: +  ts - time stepping context obtained from TSCreate()
4116: .  step - step number that has just completed
4117: .  ptime - model time of the state
4118: -  u - state at the current model time

4120:    Notes:
4121:    TSMonitor() is typically used automatically within the time stepping implementations.
4122:    Users would almost never call this routine directly.

4124:    A step of -1 indicates that the monitor is being called on a solution obtained by interpolating from computed solutions

4126:    Level: developer

4128: .keywords: TS, timestep
4129: @*/
4130: PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u)
4131: {
4132:   DM             dm;
4133:   PetscInt       i,n = ts->numbermonitors;


4140:   TSGetDM(ts,&dm);
4141:   DMSetOutputSequenceNumber(dm,step,ptime);

4143:   VecLockPush(u);
4144:   for (i=0; i<n; i++) {
4145:     (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);
4146:   }
4147:   VecLockPop(u);
4148:   return(0);
4149: }

4153: /*@C
4154:    TSAdjointMonitor - Runs all user-provided adjoint monitor routines set using TSAdjointMonitorSet()

4156:    Collective on TS

4158:    Input Parameters:
4159: +  ts - time stepping context obtained from TSCreate()
4160: .  step - step number that has just completed
4161: .  ptime - model time of the state
4162: .  u - state at the current model time
4163: .  numcost - number of cost functions (dimension of lambda  or mu)
4164: .  lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables
4165: -  mu - vectors containing the gradients of the cost functions with respect to the problem parameters

4167:    Notes:
4168:    TSAdjointMonitor() is typically used automatically within the time stepping implementations.
4169:    Users would almost never call this routine directly.

4171:    Level: developer

4173: .keywords: TS, timestep
4174: @*/
4175: PetscErrorCode TSAdjointMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda, Vec *mu)
4176: {
4178:   PetscInt       i,n = ts->numberadjointmonitors;

4183:   VecLockPush(u);
4184:   for (i=0; i<n; i++) {
4185:     (*ts->adjointmonitor[i])(ts,step,ptime,u,numcost,lambda,mu,ts->adjointmonitorcontext[i]);
4186:   }
4187:   VecLockPop(u);
4188:   return(0);
4189: }

4191: /* ------------------------------------------------------------------------*/
4194: /*@C
4195:    TSMonitorLGCtxCreate - Creates a TSMonitorLGCtx context for use with
4196:    TS to monitor the solution process graphically in various ways

4198:    Collective on TS

4200:    Input Parameters:
4201: +  host - the X display to open, or null for the local machine
4202: .  label - the title to put in the title bar
4203: .  x, y - the screen coordinates of the upper left coordinate of the window
4204: .  m, n - the screen width and height in pixels
4205: -  howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time

4207:    Output Parameter:
4208: .  ctx - the context

4210:    Options Database Key:
4211: +  -ts_monitor_lg_timestep - automatically sets line graph monitor
4212: .  -ts_monitor_lg_solution - monitor the solution (or certain values of the solution by calling TSMonitorLGSetDisplayVariables() or TSMonitorLGCtxSetDisplayVariables())
4213: .  -ts_monitor_lg_error -  monitor the error
4214: .  -ts_monitor_lg_ksp_iterations - monitor the number of KSP iterations needed for each timestep
4215: .  -ts_monitor_lg_snes_iterations - monitor the number of SNES iterations needed for each timestep
4216: -  -lg_use_markers <true,false> - mark the data points (at each time step) on the plot; default is true

4218:    Notes:
4219:    Use TSMonitorLGCtxDestroy() to destroy.

4221:    One can provide a function that transforms the solution before plotting it with TSMonitorLGCtxSetTransform() or TSMonitorLGSetTransform()

4223:    Many of the functions that control the monitoring have two forms: TSMonitorLGSet/GetXXXX() and TSMonitorLGCtxSet/GetXXXX() the first take a TS object as the
4224:    first argument (if that TS object does not have a TSMonitorLGCtx associated with it the function call is ignored) and the second takes a TSMonitorLGCtx object
4225:    as the first argument.

4227:    One can control the names displayed for each solution or error variable with TSMonitorLGCtxSetVariableNames() or TSMonitorLGSetVariableNames()


4230:    Level: intermediate

4232: .keywords: TS, monitor, line graph, residual

4234: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError(), TSMonitorDefault(), VecView(), 
4235:            TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
4236:            TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
4237:            TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
4238:            TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()

4240: @*/
4241: PetscErrorCode  TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx)
4242: {
4243:   PetscDraw      draw;

4247:   PetscNew(ctx);
4248:   PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
4249:   PetscDrawSetFromOptions(draw);
4250:   PetscDrawLGCreate(draw,1,&(*ctx)->lg);
4251:   PetscDrawLGSetFromOptions((*ctx)->lg);
4252:   PetscDrawDestroy(&draw);
4253:   (*ctx)->howoften = howoften;
4254:   return(0);
4255: }

4259: PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
4260: {
4261:   TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
4262:   PetscReal      x   = ptime,y;

4266:   if (step < 0) return(0); /* -1 indicates an interpolated solution */
4267:   if (!step) {
4268:     PetscDrawAxis axis;
4269:     PetscDrawLGGetAxis(ctx->lg,&axis);
4270:     PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time Step");
4271:     PetscDrawLGReset(ctx->lg);
4272:   }
4273:   TSGetTimeStep(ts,&y);
4274:   PetscDrawLGAddPoint(ctx->lg,&x,&y);
4275:   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
4276:     PetscDrawLGDraw(ctx->lg);
4277:     PetscDrawLGSave(ctx->lg);
4278:   }
4279:   return(0);
4280: }

4284: /*@C
4285:    TSMonitorLGCtxDestroy - Destroys a line graph context that was created
4286:    with TSMonitorLGCtxCreate().

4288:    Collective on TSMonitorLGCtx

4290:    Input Parameter:
4291: .  ctx - the monitor context

4293:    Level: intermediate

4295: .keywords: TS, monitor, line graph, destroy

4297: .seealso: TSMonitorLGCtxCreate(),  TSMonitorSet(), TSMonitorLGTimeStep();
4298: @*/
4299: PetscErrorCode  TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx)
4300: {

4304:   if ((*ctx)->transformdestroy) {
4305:     ((*ctx)->transformdestroy)((*ctx)->transformctx);
4306:   }
4307:   PetscDrawLGDestroy(&(*ctx)->lg);
4308:   PetscStrArrayDestroy(&(*ctx)->names);
4309:   PetscStrArrayDestroy(&(*ctx)->displaynames);
4310:   PetscFree((*ctx)->displayvariables);
4311:   PetscFree((*ctx)->displayvalues);
4312:   PetscFree(*ctx);
4313:   return(0);
4314: }

4318: /*@
4319:    TSGetTime - Gets the time of the most recently completed step.

4321:    Not Collective

4323:    Input Parameter:
4324: .  ts - the TS context obtained from TSCreate()

4326:    Output Parameter:
4327: .  t  - the current time. This time may not corresponds to the final time set with TSSetDuration(), use TSGetSolveTime().

4329:    Level: beginner

4331:    Note:
4332:    When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(),
4333:    TSSetPreStage(), TSSetPostStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated.

4335: .seealso: TSSetInitialTimeStep(), TSGetTimeStep(), TSGetSolveTime()

4337: .keywords: TS, get, time
4338: @*/
4339: PetscErrorCode  TSGetTime(TS ts,PetscReal *t)
4340: {
4344:   *t = ts->ptime;
4345:   return(0);
4346: }

4350: /*@
4351:    TSGetPrevTime - Gets the starting time of the previously completed step.

4353:    Not Collective

4355:    Input Parameter:
4356: .  ts - the TS context obtained from TSCreate()

4358:    Output Parameter:
4359: .  t  - the previous time

4361:    Level: beginner

4363: .seealso: TSSetInitialTimeStep(), TSGetTimeStep()

4365: .keywords: TS, get, time
4366: @*/
4367: PetscErrorCode  TSGetPrevTime(TS ts,PetscReal *t)
4368: {
4372:   *t = ts->ptime_prev;
4373:   return(0);
4374: }

4378: /*@
4379:    TSSetTime - Allows one to reset the time.

4381:    Logically Collective on TS

4383:    Input Parameters:
4384: +  ts - the TS context obtained from TSCreate()
4385: -  time - the time

4387:    Level: intermediate

4389: .seealso: TSGetTime(), TSSetDuration()

4391: .keywords: TS, set, time
4392: @*/
4393: PetscErrorCode  TSSetTime(TS ts, PetscReal t)
4394: {
4398:   ts->ptime = t;
4399:   return(0);
4400: }

4404: /*@C
4405:    TSSetOptionsPrefix - Sets the prefix used for searching for all
4406:    TS options in the database.

4408:    Logically Collective on TS

4410:    Input Parameter:
4411: +  ts     - The TS context
4412: -  prefix - The prefix to prepend to all option names

4414:    Notes:
4415:    A hyphen (-) must NOT be given at the beginning of the prefix name.
4416:    The first character of all runtime options is AUTOMATICALLY the
4417:    hyphen.

4419:    Level: advanced

4421: .keywords: TS, set, options, prefix, database

4423: .seealso: TSSetFromOptions()

4425: @*/
4426: PetscErrorCode  TSSetOptionsPrefix(TS ts,const char prefix[])
4427: {
4429:   SNES           snes;

4433:   PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);
4434:   TSGetSNES(ts,&snes);
4435:   SNESSetOptionsPrefix(snes,prefix);
4436:   return(0);
4437: }


4442: /*@C
4443:    TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4444:    TS options in the database.

4446:    Logically Collective on TS

4448:    Input Parameter:
4449: +  ts     - The TS context
4450: -  prefix - The prefix to prepend to all option names

4452:    Notes:
4453:    A hyphen (-) must NOT be given at the beginning of the prefix name.
4454:    The first character of all runtime options is AUTOMATICALLY the
4455:    hyphen.

4457:    Level: advanced

4459: .keywords: TS, append, options, prefix, database

4461: .seealso: TSGetOptionsPrefix()

4463: @*/
4464: PetscErrorCode  TSAppendOptionsPrefix(TS ts,const char prefix[])
4465: {
4467:   SNES           snes;

4471:   PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);
4472:   TSGetSNES(ts,&snes);
4473:   SNESAppendOptionsPrefix(snes,prefix);
4474:   return(0);
4475: }

4479: /*@C
4480:    TSGetOptionsPrefix - Sets the prefix used for searching for all
4481:    TS options in the database.

4483:    Not Collective

4485:    Input Parameter:
4486: .  ts - The TS context

4488:    Output Parameter:
4489: .  prefix - A pointer to the prefix string used

4491:    Notes: On the fortran side, the user should pass in a string 'prifix' of
4492:    sufficient length to hold the prefix.

4494:    Level: intermediate

4496: .keywords: TS, get, options, prefix, database

4498: .seealso: TSAppendOptionsPrefix()
4499: @*/
4500: PetscErrorCode  TSGetOptionsPrefix(TS ts,const char *prefix[])
4501: {

4507:   PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);
4508:   return(0);
4509: }

4513: /*@C
4514:    TSGetRHSJacobian - Returns the Jacobian J at the present timestep.

4516:    Not Collective, but parallel objects are returned if TS is parallel

4518:    Input Parameter:
4519: .  ts  - The TS context obtained from TSCreate()

4521:    Output Parameters:
4522: +  Amat - The (approximate) Jacobian J of G, where U_t = G(U,t)  (or NULL)
4523: .  Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat  (or NULL)
4524: .  func - Function to compute the Jacobian of the RHS  (or NULL)
4525: -  ctx - User-defined context for Jacobian evaluation routine  (or NULL)

4527:    Notes: You can pass in NULL for any return argument you do not need.

4529:    Level: intermediate

4531: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber()

4533: .keywords: TS, timestep, get, matrix, Jacobian
4534: @*/
4535: PetscErrorCode  TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx)
4536: {
4538:   SNES           snes;
4539:   DM             dm;

4542:   TSGetSNES(ts,&snes);
4543:   SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
4544:   TSGetDM(ts,&dm);
4545:   DMTSGetRHSJacobian(dm,func,ctx);
4546:   return(0);
4547: }

4551: /*@C
4552:    TSGetIJacobian - Returns the implicit Jacobian at the present timestep.

4554:    Not Collective, but parallel objects are returned if TS is parallel

4556:    Input Parameter:
4557: .  ts  - The TS context obtained from TSCreate()

4559:    Output Parameters:
4560: +  Amat  - The (approximate) Jacobian of F(t,U,U_t)
4561: .  Pmat - The matrix from which the preconditioner is constructed, often the same as Amat
4562: .  f   - The function to compute the matrices
4563: - ctx - User-defined context for Jacobian evaluation routine

4565:    Notes: You can pass in NULL for any return argument you do not need.

4567:    Level: advanced

4569: .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber()

4571: .keywords: TS, timestep, get, matrix, Jacobian
4572: @*/
4573: PetscErrorCode  TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx)
4574: {
4576:   DM             dm;

4579:   if (Amat || Pmat) {
4580:     SNES snes;
4581:     TSGetSNES(ts,&snes);
4582:     SNESSetUpMatrices(snes);
4583:     SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
4584:   }
4585:   TSGetDM(ts,&dm);
4586:   DMTSGetIJacobian(dm,f,ctx);
4587:   return(0);
4588: }


4593: /*@C
4594:    TSMonitorDrawSolution - Monitors progress of the TS solvers by calling
4595:    VecView() for the solution at each timestep

4597:    Collective on TS

4599:    Input Parameters:
4600: +  ts - the TS context
4601: .  step - current time-step
4602: .  ptime - current time
4603: -  dummy - either a viewer or NULL

4605:    Options Database:
4606: .   -ts_monitor_draw_solution_initial - show initial solution as well as current solution

4608:    Notes: the initial solution and current solution are not display with a common axis scaling so generally the option -ts_monitor_draw_solution_initial
4609:        will look bad

4611:    Level: intermediate

4613: .keywords: TS,  vector, monitor, view

4615: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4616: @*/
4617: PetscErrorCode  TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4618: {
4619:   PetscErrorCode   ierr;
4620:   TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
4621:   PetscDraw        draw;

4624:   if (!step && ictx->showinitial) {
4625:     if (!ictx->initialsolution) {
4626:       VecDuplicate(u,&ictx->initialsolution);
4627:     }
4628:     VecCopy(u,ictx->initialsolution);
4629:   }
4630:   if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);

4632:   if (ictx->showinitial) {
4633:     PetscReal pause;
4634:     PetscViewerDrawGetPause(ictx->viewer,&pause);
4635:     PetscViewerDrawSetPause(ictx->viewer,0.0);
4636:     VecView(ictx->initialsolution,ictx->viewer);
4637:     PetscViewerDrawSetPause(ictx->viewer,pause);
4638:     PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);
4639:   }
4640:   VecView(u,ictx->viewer);
4641:   if (ictx->showtimestepandtime) {
4642:     PetscReal xl,yl,xr,yr,h;
4643:     char      time[32];

4645:     PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
4646:     PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
4647:     PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
4648:     h    = yl + .95*(yr - yl);
4649:     PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
4650:     PetscDrawFlush(draw);
4651:   }

4653:   if (ictx->showinitial) {
4654:     PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);
4655:   }
4656:   return(0);
4657: }

4661: /*@C
4662:    TSAdjointMonitorDrawSensi - Monitors progress of the adjoint TS solvers by calling
4663:    VecView() for the sensitivities to initial states at each timestep

4665:    Collective on TS

4667:    Input Parameters:
4668: +  ts - the TS context
4669: .  step - current time-step
4670: .  ptime - current time
4671: .  u - current state
4672: .  numcost - number of cost functions
4673: .  lambda - sensitivities to initial conditions
4674: .  mu - sensitivities to parameters
4675: -  dummy - either a viewer or NULL

4677:    Level: intermediate

4679: .keywords: TS,  vector, adjoint, monitor, view

4681: .seealso: TSAdjointMonitorSet(), TSAdjointMonitorDefault(), VecView()
4682: @*/
4683: PetscErrorCode  TSAdjointMonitorDrawSensi(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscInt numcost,Vec *lambda,Vec *mu,void *dummy)
4684: {
4685:   PetscErrorCode   ierr;
4686:   TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
4687:   PetscDraw        draw;
4688:   PetscReal        xl,yl,xr,yr,h;
4689:   char             time[32];

4692:   if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);

4694:   VecView(lambda[0],ictx->viewer);
4695:   PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
4696:   PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
4697:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
4698:   h    = yl + .95*(yr - yl);
4699:   PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
4700:   PetscDrawFlush(draw);
4701:   return(0);
4702: }

4706: /*@C
4707:    TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram

4709:    Collective on TS

4711:    Input Parameters:
4712: +  ts - the TS context
4713: .  step - current time-step
4714: .  ptime - current time
4715: -  dummy - either a viewer or NULL

4717:    Level: intermediate

4719: .keywords: TS,  vector, monitor, view

4721: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4722: @*/
4723: PetscErrorCode  TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4724: {
4725:   PetscErrorCode    ierr;
4726:   TSMonitorDrawCtx  ictx = (TSMonitorDrawCtx)dummy;
4727:   PetscDraw         draw;
4728:   PetscDrawAxis     axis;
4729:   PetscInt          n;
4730:   PetscMPIInt       size;
4731:   PetscReal         U0,U1,xl,yl,xr,yr,h;
4732:   char              time[32];
4733:   const PetscScalar *U;

4736:   MPI_Comm_size(PetscObjectComm((PetscObject)ts),&size);
4737:   if (size != 1) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Only allowed for sequential runs");
4738:   VecGetSize(u,&n);
4739:   if (n != 2) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Only for ODEs with two unknowns");

4741:   PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
4742:   PetscViewerDrawGetDrawAxis(ictx->viewer,0,&axis);
4743:   PetscDrawAxisGetLimits(axis,&xl,&xr,&yl,&yr);
4744:   if (!step) {
4745:     PetscDrawClear(draw);
4746:     PetscDrawAxisDraw(axis);
4747:   }

4749:   VecGetArrayRead(u,&U);
4750:   U0 = PetscRealPart(U[0]);
4751:   U1 = PetscRealPart(U[1]);
4752:   VecRestoreArrayRead(u,&U);
4753:   if ((U0 < xl) || (U1 < yl) || (U0 > xr) || (U1 > yr)) return(0);

4755:   PetscDrawCollectiveBegin(draw);
4756:   PetscDrawPoint(draw,U0,U1,PETSC_DRAW_BLACK);
4757:   if (ictx->showtimestepandtime) {
4758:     PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
4759:     PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
4760:     h    = yl + .95*(yr - yl);
4761:     PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
4762:   }
4763:   PetscDrawCollectiveEnd(draw);
4764:   PetscDrawFlush(draw);
4765:   PetscDrawSave(draw);
4766:   return(0);
4767: }


4772: /*@C
4773:    TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution()

4775:    Collective on TS

4777:    Input Parameters:
4778: .    ctx - the monitor context

4780:    Level: intermediate

4782: .keywords: TS,  vector, monitor, view

4784: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError()
4785: @*/
4786: PetscErrorCode  TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx)
4787: {

4791:   PetscViewerDestroy(&(*ictx)->viewer);
4792:   VecDestroy(&(*ictx)->initialsolution);
4793:   PetscFree(*ictx);
4794:   return(0);
4795: }

4799: /*@C
4800:    TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx

4802:    Collective on TS

4804:    Input Parameter:
4805: .    ts - time-step context

4807:    Output Patameter:
4808: .    ctx - the monitor context

4810:    Options Database:
4811: .   -ts_monitor_draw_solution_initial - show initial solution as well as current solution

4813:    Level: intermediate

4815: .keywords: TS,  vector, monitor, view

4817: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx()
4818: @*/
4819: PetscErrorCode  TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx)
4820: {
4821:   PetscErrorCode   ierr;

4824:   PetscNew(ctx);
4825:   PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);
4826:   PetscViewerSetFromOptions((*ctx)->viewer);

4828:   (*ctx)->howoften    = howoften;
4829:   (*ctx)->showinitial = PETSC_FALSE;
4830:   PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);

4832:   (*ctx)->showtimestepandtime = PETSC_FALSE;
4833:   PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);
4834:   return(0);
4835: }

4839: /*@C
4840:    TSMonitorDrawError - Monitors progress of the TS solvers by calling
4841:    VecView() for the error at each timestep

4843:    Collective on TS

4845:    Input Parameters:
4846: +  ts - the TS context
4847: .  step - current time-step
4848: .  ptime - current time
4849: -  dummy - either a viewer or NULL

4851:    Level: intermediate

4853: .keywords: TS,  vector, monitor, view

4855: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4856: @*/
4857: PetscErrorCode  TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4858: {
4859:   PetscErrorCode   ierr;
4860:   TSMonitorDrawCtx ctx    = (TSMonitorDrawCtx)dummy;
4861:   PetscViewer      viewer = ctx->viewer;
4862:   Vec              work;

4865:   if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) return(0);
4866:   VecDuplicate(u,&work);
4867:   TSComputeSolutionFunction(ts,ptime,work);
4868:   VecAXPY(work,-1.0,u);
4869:   VecView(work,viewer);
4870:   VecDestroy(&work);
4871:   return(0);
4872: }

4874: #include <petsc/private/dmimpl.h>
4877: /*@
4878:    TSSetDM - Sets the DM that may be used by some preconditioners

4880:    Logically Collective on TS and DM

4882:    Input Parameters:
4883: +  ts - the preconditioner context
4884: -  dm - the dm

4886:    Level: intermediate


4889: .seealso: TSGetDM(), SNESSetDM(), SNESGetDM()
4890: @*/
4891: PetscErrorCode  TSSetDM(TS ts,DM dm)
4892: {
4894:   SNES           snes;
4895:   DMTS           tsdm;

4899:   PetscObjectReference((PetscObject)dm);
4900:   if (ts->dm) {               /* Move the DMTS context over to the new DM unless the new DM already has one */
4901:     if (ts->dm->dmts && !dm->dmts) {
4902:       DMCopyDMTS(ts->dm,dm);
4903:       DMGetDMTS(ts->dm,&tsdm);
4904:       if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */
4905:         tsdm->originaldm = dm;
4906:       }
4907:     }
4908:     DMDestroy(&ts->dm);
4909:   }
4910:   ts->dm = dm;

4912:   TSGetSNES(ts,&snes);
4913:   SNESSetDM(snes,dm);
4914:   return(0);
4915: }

4919: /*@
4920:    TSGetDM - Gets the DM that may be used by some preconditioners

4922:    Not Collective

4924:    Input Parameter:
4925: . ts - the preconditioner context

4927:    Output Parameter:
4928: .  dm - the dm

4930:    Level: intermediate


4933: .seealso: TSSetDM(), SNESSetDM(), SNESGetDM()
4934: @*/
4935: PetscErrorCode  TSGetDM(TS ts,DM *dm)
4936: {

4941:   if (!ts->dm) {
4942:     DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);
4943:     if (ts->snes) {SNESSetDM(ts->snes,ts->dm);}
4944:   }
4945:   *dm = ts->dm;
4946:   return(0);
4947: }

4951: /*@
4952:    SNESTSFormFunction - Function to evaluate nonlinear residual

4954:    Logically Collective on SNES

4956:    Input Parameter:
4957: + snes - nonlinear solver
4958: . U - the current state at which to evaluate the residual
4959: - ctx - user context, must be a TS

4961:    Output Parameter:
4962: . F - the nonlinear residual

4964:    Notes:
4965:    This function is not normally called by users and is automatically registered with the SNES used by TS.
4966:    It is most frequently passed to MatFDColoringSetFunction().

4968:    Level: advanced

4970: .seealso: SNESSetFunction(), MatFDColoringSetFunction()
4971: @*/
4972: PetscErrorCode  SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx)
4973: {
4974:   TS             ts = (TS)ctx;

4982:   (ts->ops->snesfunction)(snes,U,F,ts);
4983:   return(0);
4984: }

4988: /*@
4989:    SNESTSFormJacobian - Function to evaluate the Jacobian

4991:    Collective on SNES

4993:    Input Parameter:
4994: + snes - nonlinear solver
4995: . U - the current state at which to evaluate the residual
4996: - ctx - user context, must be a TS

4998:    Output Parameter:
4999: + A - the Jacobian
5000: . B - the preconditioning matrix (may be the same as A)
5001: - flag - indicates any structure change in the matrix

5003:    Notes:
5004:    This function is not normally called by users and is automatically registered with the SNES used by TS.

5006:    Level: developer

5008: .seealso: SNESSetJacobian()
5009: @*/
5010: PetscErrorCode  SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx)
5011: {
5012:   TS             ts = (TS)ctx;

5023:   (ts->ops->snesjacobian)(snes,U,A,B,ts);
5024:   return(0);
5025: }

5029: /*@C
5030:    TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems Udot = A U only

5032:    Collective on TS

5034:    Input Arguments:
5035: +  ts - time stepping context
5036: .  t - time at which to evaluate
5037: .  U - state at which to evaluate
5038: -  ctx - context

5040:    Output Arguments:
5041: .  F - right hand side

5043:    Level: intermediate

5045:    Notes:
5046:    This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems.
5047:    The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian().

5049: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
5050: @*/
5051: PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
5052: {
5054:   Mat            Arhs,Brhs;

5057:   TSGetRHSMats_Private(ts,&Arhs,&Brhs);
5058:   TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);
5059:   MatMult(Arhs,U,F);
5060:   return(0);
5061: }

5065: /*@C
5066:    TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.

5068:    Collective on TS

5070:    Input Arguments:
5071: +  ts - time stepping context
5072: .  t - time at which to evaluate
5073: .  U - state at which to evaluate
5074: -  ctx - context

5076:    Output Arguments:
5077: +  A - pointer to operator
5078: .  B - pointer to preconditioning matrix
5079: -  flg - matrix structure flag

5081:    Level: intermediate

5083:    Notes:
5084:    This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems.

5086: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear()
5087: @*/
5088: PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
5089: {
5091:   return(0);
5092: }

5096: /*@C
5097:    TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only

5099:    Collective on TS

5101:    Input Arguments:
5102: +  ts - time stepping context
5103: .  t - time at which to evaluate
5104: .  U - state at which to evaluate
5105: .  Udot - time derivative of state vector
5106: -  ctx - context

5108:    Output Arguments:
5109: .  F - left hand side

5111:    Level: intermediate

5113:    Notes:
5114:    The assumption here is that the left hand side is of the form A*Udot (and not A*Udot + B*U). For other cases, the
5115:    user is required to write their own TSComputeIFunction.
5116:    This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems.
5117:    The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian().

5119:    Note that using this function is NOT equivalent to using TSComputeRHSFunctionLinear() since that solves Udot = A U

5121: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant(), TSComputeRHSFunctionLinear()
5122: @*/
5123: PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
5124: {
5126:   Mat            A,B;

5129:   TSGetIJacobian(ts,&A,&B,NULL,NULL);
5130:   TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);
5131:   MatMult(A,Udot,F);
5132:   return(0);
5133: }

5137: /*@C
5138:    TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE

5140:    Collective on TS

5142:    Input Arguments:
5143: +  ts - time stepping context
5144: .  t - time at which to evaluate
5145: .  U - state at which to evaluate
5146: .  Udot - time derivative of state vector
5147: .  shift - shift to apply
5148: -  ctx - context

5150:    Output Arguments:
5151: +  A - pointer to operator
5152: .  B - pointer to preconditioning matrix
5153: -  flg - matrix structure flag

5155:    Level: advanced

5157:    Notes:
5158:    This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems.

5160:    It is only appropriate for problems of the form

5162: $     M Udot = F(U,t)

5164:   where M is constant and F is non-stiff.  The user must pass M to TSSetIJacobian().  The current implementation only
5165:   works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing
5166:   an implicit operator of the form

5168: $    shift*M + J

5170:   where J is the Jacobian of -F(U).  Support may be added in a future version of PETSc, but for now, the user must store
5171:   a copy of M or reassemble it when requested.

5173: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear()
5174: @*/
5175: PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx)
5176: {

5180:   MatScale(A, shift / ts->ijacobian.shift);
5181:   ts->ijacobian.shift = shift;
5182:   return(0);
5183: }

5187: /*@
5188:    TSGetEquationType - Gets the type of the equation that TS is solving.

5190:    Not Collective

5192:    Input Parameter:
5193: .  ts - the TS context

5195:    Output Parameter:
5196: .  equation_type - see TSEquationType

5198:    Level: beginner

5200: .keywords: TS, equation type

5202: .seealso: TSSetEquationType(), TSEquationType
5203: @*/
5204: PetscErrorCode  TSGetEquationType(TS ts,TSEquationType *equation_type)
5205: {
5209:   *equation_type = ts->equation_type;
5210:   return(0);
5211: }

5215: /*@
5216:    TSSetEquationType - Sets the type of the equation that TS is solving.

5218:    Not Collective

5220:    Input Parameter:
5221: +  ts - the TS context
5222: -  equation_type - see TSEquationType

5224:    Level: advanced

5226: .keywords: TS, equation type

5228: .seealso: TSGetEquationType(), TSEquationType
5229: @*/
5230: PetscErrorCode  TSSetEquationType(TS ts,TSEquationType equation_type)
5231: {
5234:   ts->equation_type = equation_type;
5235:   return(0);
5236: }

5240: /*@
5241:    TSGetConvergedReason - Gets the reason the TS iteration was stopped.

5243:    Not Collective

5245:    Input Parameter:
5246: .  ts - the TS context

5248:    Output Parameter:
5249: .  reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
5250:             manual pages for the individual convergence tests for complete lists

5252:    Level: beginner

5254:    Notes:
5255:    Can only be called after the call to TSSolve() is complete.

5257: .keywords: TS, nonlinear, set, convergence, test

5259: .seealso: TSSetConvergenceTest(), TSConvergedReason
5260: @*/
5261: PetscErrorCode  TSGetConvergedReason(TS ts,TSConvergedReason *reason)
5262: {
5266:   *reason = ts->reason;
5267:   return(0);
5268: }

5272: /*@
5273:    TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve.

5275:    Not Collective

5277:    Input Parameter:
5278: +  ts - the TS context
5279: .  reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
5280:             manual pages for the individual convergence tests for complete lists

5282:    Level: advanced

5284:    Notes:
5285:    Can only be called during TSSolve() is active.

5287: .keywords: TS, nonlinear, set, convergence, test

5289: .seealso: TSConvergedReason
5290: @*/
5291: PetscErrorCode  TSSetConvergedReason(TS ts,TSConvergedReason reason)
5292: {
5295:   ts->reason = reason;
5296:   return(0);
5297: }

5301: /*@
5302:    TSGetSolveTime - Gets the time after a call to TSSolve()

5304:    Not Collective

5306:    Input Parameter:
5307: .  ts - the TS context

5309:    Output Parameter:
5310: .  ftime - the final time. This time corresponds to the final time set with TSSetDuration()

5312:    Level: beginner

5314:    Notes:
5315:    Can only be called after the call to TSSolve() is complete.

5317: .keywords: TS, nonlinear, set, convergence, test

5319: .seealso: TSSetConvergenceTest(), TSConvergedReason
5320: @*/
5321: PetscErrorCode  TSGetSolveTime(TS ts,PetscReal *ftime)
5322: {
5326:   *ftime = ts->solvetime;
5327:   return(0);
5328: }

5332: /*@
5333:    TSGetTotalSteps - Gets the total number of steps done since the last call to TSSetUp() or TSCreate()

5335:    Not Collective

5337:    Input Parameter:
5338: .  ts - the TS context

5340:    Output Parameter:
5341: .  steps - the number of steps

5343:    Level: beginner

5345:    Notes:
5346:    Includes the number of steps for all calls to TSSolve() since TSSetUp() was called

5348: .keywords: TS, nonlinear, set, convergence, test

5350: .seealso: TSSetConvergenceTest(), TSConvergedReason
5351: @*/
5352: PetscErrorCode  TSGetTotalSteps(TS ts,PetscInt *steps)
5353: {
5357:   *steps = ts->total_steps;
5358:   return(0);
5359: }

5363: /*@
5364:    TSGetSNESIterations - Gets the total number of nonlinear iterations
5365:    used by the time integrator.

5367:    Not Collective

5369:    Input Parameter:
5370: .  ts - TS context

5372:    Output Parameter:
5373: .  nits - number of nonlinear iterations

5375:    Notes:
5376:    This counter is reset to zero for each successive call to TSSolve().

5378:    Level: intermediate

5380: .keywords: TS, get, number, nonlinear, iterations

5382: .seealso:  TSGetKSPIterations()
5383: @*/
5384: PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits)
5385: {
5389:   *nits = ts->snes_its;
5390:   return(0);
5391: }

5395: /*@
5396:    TSGetKSPIterations - Gets the total number of linear iterations
5397:    used by the time integrator.

5399:    Not Collective

5401:    Input Parameter:
5402: .  ts - TS context

5404:    Output Parameter:
5405: .  lits - number of linear iterations

5407:    Notes:
5408:    This counter is reset to zero for each successive call to TSSolve().

5410:    Level: intermediate

5412: .keywords: TS, get, number, linear, iterations

5414: .seealso:  TSGetSNESIterations(), SNESGetKSPIterations()
5415: @*/
5416: PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits)
5417: {
5421:   *lits = ts->ksp_its;
5422:   return(0);
5423: }

5427: /*@
5428:    TSGetStepRejections - Gets the total number of rejected steps.

5430:    Not Collective

5432:    Input Parameter:
5433: .  ts - TS context

5435:    Output Parameter:
5436: .  rejects - number of steps rejected

5438:    Notes:
5439:    This counter is reset to zero for each successive call to TSSolve().

5441:    Level: intermediate

5443: .keywords: TS, get, number

5445: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails()
5446: @*/
5447: PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects)
5448: {
5452:   *rejects = ts->reject;
5453:   return(0);
5454: }

5458: /*@
5459:    TSGetSNESFailures - Gets the total number of failed SNES solves

5461:    Not Collective

5463:    Input Parameter:
5464: .  ts - TS context

5466:    Output Parameter:
5467: .  fails - number of failed nonlinear solves

5469:    Notes:
5470:    This counter is reset to zero for each successive call to TSSolve().

5472:    Level: intermediate

5474: .keywords: TS, get, number

5476: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures()
5477: @*/
5478: PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails)
5479: {
5483:   *fails = ts->num_snes_failures;
5484:   return(0);
5485: }

5489: /*@
5490:    TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails

5492:    Not Collective

5494:    Input Parameter:
5495: +  ts - TS context
5496: -  rejects - maximum number of rejected steps, pass -1 for unlimited

5498:    Notes:
5499:    The counter is reset to zero for each step

5501:    Options Database Key:
5502:  .  -ts_max_reject - Maximum number of step rejections before a step fails

5504:    Level: intermediate

5506: .keywords: TS, set, maximum, number

5508: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
5509: @*/
5510: PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects)
5511: {
5514:   ts->max_reject = rejects;
5515:   return(0);
5516: }

5520: /*@
5521:    TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves

5523:    Not Collective

5525:    Input Parameter:
5526: +  ts - TS context
5527: -  fails - maximum number of failed nonlinear solves, pass -1 for unlimited

5529:    Notes:
5530:    The counter is reset to zero for each successive call to TSSolve().

5532:    Options Database Key:
5533:  .  -ts_max_snes_failures - Maximum number of nonlinear solve failures

5535:    Level: intermediate

5537: .keywords: TS, set, maximum, number

5539: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason()
5540: @*/
5541: PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails)
5542: {
5545:   ts->max_snes_failures = fails;
5546:   return(0);
5547: }

5551: /*@
5552:    TSSetErrorIfStepFails - Error if no step succeeds

5554:    Not Collective

5556:    Input Parameter:
5557: +  ts - TS context
5558: -  err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure

5560:    Options Database Key:
5561:  .  -ts_error_if_step_fails - Error if no step succeeds

5563:    Level: intermediate

5565: .keywords: TS, set, error

5567: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
5568: @*/
5569: PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err)
5570: {
5573:   ts->errorifstepfailed = err;
5574:   return(0);
5575: }

5579: /*@C
5580:    TSMonitorSolution - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file or a PetscDraw object

5582:    Collective on TS

5584:    Input Parameters:
5585: +  ts - the TS context
5586: .  step - current time-step
5587: .  ptime - current time
5588: .  u - current state
5589: -  vf - viewer and its format

5591:    Level: intermediate

5593: .keywords: TS,  vector, monitor, view

5595: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5596: @*/
5597: PetscErrorCode  TSMonitorSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscViewerAndFormat *vf)
5598: {

5602:   PetscViewerPushFormat(vf->viewer,vf->format);
5603:   VecView(u,vf->viewer);
5604:   PetscViewerPopFormat(vf->viewer);
5605:   return(0);
5606: }

5610: /*@C
5611:    TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep.

5613:    Collective on TS

5615:    Input Parameters:
5616: +  ts - the TS context
5617: .  step - current time-step
5618: .  ptime - current time
5619: .  u - current state
5620: -  filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)

5622:    Level: intermediate

5624:    Notes:
5625:    The VTK format does not allow writing multiple time steps in the same file, therefore a different file will be written for each time step.
5626:    These are named according to the file name template.

5628:    This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy().

5630: .keywords: TS,  vector, monitor, view

5632: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5633: @*/
5634: PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate)
5635: {
5637:   char           filename[PETSC_MAX_PATH_LEN];
5638:   PetscViewer    viewer;

5641:   if (step < 0) return(0); /* -1 indicates interpolated solution */
5642:   PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);
5643:   PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);
5644:   VecView(u,viewer);
5645:   PetscViewerDestroy(&viewer);
5646:   return(0);
5647: }

5651: /*@C
5652:    TSMonitorSolutionVTKDestroy - Destroy context for monitoring

5654:    Collective on TS

5656:    Input Parameters:
5657: .  filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)

5659:    Level: intermediate

5661:    Note:
5662:    This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK().

5664: .keywords: TS,  vector, monitor, view

5666: .seealso: TSMonitorSet(), TSMonitorSolutionVTK()
5667: @*/
5668: PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate)
5669: {

5673:   PetscFree(*(char**)filenametemplate);
5674:   return(0);
5675: }

5679: /*@
5680:    TSGetAdapt - Get the adaptive controller context for the current method

5682:    Collective on TS if controller has not been created yet

5684:    Input Arguments:
5685: .  ts - time stepping context

5687:    Output Arguments:
5688: .  adapt - adaptive controller

5690:    Level: intermediate

5692: .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose()
5693: @*/
5694: PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt)
5695: {

5701:   if (!ts->adapt) {
5702:     TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);
5703:     PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);
5704:     PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);
5705:   }
5706:   *adapt = ts->adapt;
5707:   return(0);
5708: }

5712: /*@
5713:    TSSetTolerances - Set tolerances for local truncation error when using adaptive controller

5715:    Logically Collective

5717:    Input Arguments:
5718: +  ts - time integration context
5719: .  atol - scalar absolute tolerances, PETSC_DECIDE to leave current value
5720: .  vatol - vector of absolute tolerances or NULL, used in preference to atol if present
5721: .  rtol - scalar relative tolerances, PETSC_DECIDE to leave current value
5722: -  vrtol - vector of relative tolerances or NULL, used in preference to atol if present

5724:    Options Database keys:
5725: +  -ts_rtol <rtol> - relative tolerance for local truncation error
5726: -  -ts_atol <atol> Absolute tolerance for local truncation error

5728:    Notes:
5729:    With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5730:    (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5731:    computed only for the differential or the algebraic part then this can be done using the vector of
5732:    tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the 
5733:    differential part and infinity for the algebraic part, the LTE calculation will include only the
5734:    differential variables.

5736:    Level: beginner

5738: .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances()
5739: @*/
5740: PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol)
5741: {

5745:   if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol;
5746:   if (vatol) {
5747:     PetscObjectReference((PetscObject)vatol);
5748:     VecDestroy(&ts->vatol);
5749:     ts->vatol = vatol;
5750:   }
5751:   if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol;
5752:   if (vrtol) {
5753:     PetscObjectReference((PetscObject)vrtol);
5754:     VecDestroy(&ts->vrtol);
5755:     ts->vrtol = vrtol;
5756:   }
5757:   return(0);
5758: }

5762: /*@
5763:    TSGetTolerances - Get tolerances for local truncation error when using adaptive controller

5765:    Logically Collective

5767:    Input Arguments:
5768: .  ts - time integration context

5770:    Output Arguments:
5771: +  atol - scalar absolute tolerances, NULL to ignore
5772: .  vatol - vector of absolute tolerances, NULL to ignore
5773: .  rtol - scalar relative tolerances, NULL to ignore
5774: -  vrtol - vector of relative tolerances, NULL to ignore

5776:    Level: beginner

5778: .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances()
5779: @*/
5780: PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol)
5781: {
5783:   if (atol)  *atol  = ts->atol;
5784:   if (vatol) *vatol = ts->vatol;
5785:   if (rtol)  *rtol  = ts->rtol;
5786:   if (vrtol) *vrtol = ts->vrtol;
5787:   return(0);
5788: }

5792: /*@
5793:    TSErrorWeightedNorm2 - compute a weighted 2-norm of the difference between two state vectors

5795:    Collective on TS

5797:    Input Arguments:
5798: +  ts - time stepping context
5799: .  U - state vector, usually ts->vec_sol
5800: -  Y - state vector to be compared to U

5802:    Output Arguments:
5803: .  norm - weighted norm, a value of 1.0 is considered small

5805:    Level: developer

5807: .seealso: TSErrorWeightedNorm(), TSErrorWeightedNormInfinity()
5808: @*/
5809: PetscErrorCode TSErrorWeightedNorm2(TS ts,Vec U,Vec Y,PetscReal *norm)
5810: {
5811:   PetscErrorCode    ierr;
5812:   PetscInt          i,n,N,rstart;
5813:   const PetscScalar *u,*y;
5814:   PetscReal         sum,gsum;
5815:   PetscReal         tol;

5825:   if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector");

5827:   VecGetSize(U,&N);
5828:   VecGetLocalSize(U,&n);
5829:   VecGetOwnershipRange(U,&rstart,NULL);
5830:   VecGetArrayRead(U,&u);
5831:   VecGetArrayRead(Y,&y);
5832:   sum  = 0.;
5833:   if (ts->vatol && ts->vrtol) {
5834:     const PetscScalar *atol,*rtol;
5835:     VecGetArrayRead(ts->vatol,&atol);
5836:     VecGetArrayRead(ts->vrtol,&rtol);
5837:     for (i=0; i<n; i++) {
5838:       tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5839:       sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol);
5840:     }
5841:     VecRestoreArrayRead(ts->vatol,&atol);
5842:     VecRestoreArrayRead(ts->vrtol,&rtol);
5843:   } else if (ts->vatol) {       /* vector atol, scalar rtol */
5844:     const PetscScalar *atol;
5845:     VecGetArrayRead(ts->vatol,&atol);
5846:     for (i=0; i<n; i++) {
5847:       tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5848:       sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol);
5849:     }
5850:     VecRestoreArrayRead(ts->vatol,&atol);
5851:   } else if (ts->vrtol) {       /* scalar atol, vector rtol */
5852:     const PetscScalar *rtol;
5853:     VecGetArrayRead(ts->vrtol,&rtol);
5854:     for (i=0; i<n; i++) {
5855:       tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5856:       sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol);
5857:     }
5858:     VecRestoreArrayRead(ts->vrtol,&rtol);
5859:   } else {                      /* scalar atol, scalar rtol */
5860:     for (i=0; i<n; i++) {
5861:       tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5862:       sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol);
5863:     }
5864:   }
5865:   VecRestoreArrayRead(U,&u);
5866:   VecRestoreArrayRead(Y,&y);

5868:   MPIU_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));
5869:   *norm = PetscSqrtReal(gsum / N);

5871:   if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
5872:   return(0);
5873: }

5877: /*@
5878:    TSErrorWeightedNormInfinity - compute a weighted infinity-norm of the difference between two state vectors

5880:    Collective on TS

5882:    Input Arguments:
5883: +  ts - time stepping context
5884: .  U - state vector, usually ts->vec_sol
5885: -  Y - state vector to be compared to U

5887:    Output Arguments:
5888: .  norm - weighted norm, a value of 1.0 is considered small

5890:    Level: developer

5892: .seealso: TSErrorWeightedNorm(), TSErrorWeightedNorm2()
5893: @*/
5894: PetscErrorCode TSErrorWeightedNormInfinity(TS ts,Vec U,Vec Y,PetscReal *norm)
5895: {
5896:   PetscErrorCode    ierr;
5897:   PetscInt          i,n,N,rstart,k;
5898:   const PetscScalar *u,*y;
5899:   PetscReal         max,gmax;
5900:   PetscReal         tol;

5910:   if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector");

5912:   VecGetSize(U,&N);
5913:   VecGetLocalSize(U,&n);
5914:   VecGetOwnershipRange(U,&rstart,NULL);
5915:   VecGetArrayRead(U,&u);
5916:   VecGetArrayRead(Y,&y);
5917:   if (ts->vatol && ts->vrtol) {
5918:     const PetscScalar *atol,*rtol;
5919:     VecGetArrayRead(ts->vatol,&atol);
5920:     VecGetArrayRead(ts->vrtol,&rtol);
5921:     k = 0;
5922:     tol = PetscRealPart(atol[k]) + PetscRealPart(rtol[k]) * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k]));
5923:     max = PetscAbsScalar(y[k] - u[k]) / tol;
5924:     for (i=1; i<n; i++) {
5925:       tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5926:       max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol);
5927:     }
5928:     VecRestoreArrayRead(ts->vatol,&atol);
5929:     VecRestoreArrayRead(ts->vrtol,&rtol);
5930:   } else if (ts->vatol) {       /* vector atol, scalar rtol */
5931:     const PetscScalar *atol;
5932:     VecGetArrayRead(ts->vatol,&atol);
5933:     k = 0;
5934:     tol = PetscRealPart(atol[k]) + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k]));
5935:     max = PetscAbsScalar(y[k] - u[k]) / tol;
5936:     for (i=1; i<n; i++) {
5937:       tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5938:       max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol);
5939:     }
5940:     VecRestoreArrayRead(ts->vatol,&atol);
5941:   } else if (ts->vrtol) {       /* scalar atol, vector rtol */
5942:     const PetscScalar *rtol;
5943:     VecGetArrayRead(ts->vrtol,&rtol);
5944:     k = 0;
5945:     tol = ts->atol + PetscRealPart(rtol[k]) * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k]));
5946:     max = PetscAbsScalar(y[k] - u[k]) / tol;
5947:     for (i=1; i<n; i++) {
5948:       tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5949:       max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol);
5950:     }
5951:     VecRestoreArrayRead(ts->vrtol,&rtol);
5952:   } else {                      /* scalar atol, scalar rtol */
5953:     k = 0;
5954:     tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[k]),PetscAbsScalar(y[k]));
5955:     max = PetscAbsScalar(y[k] - u[k]) / tol;
5956:     for (i=1; i<n; i++) {
5957:       tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5958:       max = PetscMax(max,PetscAbsScalar(y[i] - u[i]) / tol);
5959:     }
5960:   }
5961:   VecRestoreArrayRead(U,&u);
5962:   VecRestoreArrayRead(Y,&y);

5964:   MPIU_Allreduce(&max,&gmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)ts));
5965:   *norm = gmax;

5967:   if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
5968:   return(0);
5969: }

5973: /*@
5974:    TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors

5976:    Collective on TS

5978:    Input Arguments:
5979: +  ts - time stepping context
5980: .  U - state vector, usually ts->vec_sol
5981: .  Y - state vector to be compared to U
5982: -  wnormtype - norm type, either NORM_2 or NORM_INFINITY

5984:    Output Arguments:
5985: .  norm - weighted norm, a value of 1.0 is considered small


5988:    Options Database Keys:
5989: .  -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5991:    Level: developer

5993: .seealso: TSErrorWeightedNormInfinity(), TSErrorWeightedNorm2()
5994: @*/
5995: PetscErrorCode TSErrorWeightedNorm(TS ts,Vec U,Vec Y,NormType wnormtype,PetscReal *norm)
5996: {

6000:   if (wnormtype == NORM_2) {
6001:     TSErrorWeightedNorm2(ts,U,Y,norm);
6002:   } else if(wnormtype == NORM_INFINITY) {
6003:     TSErrorWeightedNormInfinity(ts,U,Y,norm);
6004:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
6005:   return(0);
6006: }

6010: /*@
6011:    TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler

6013:    Logically Collective on TS

6015:    Input Arguments:
6016: +  ts - time stepping context
6017: -  cfltime - maximum stable time step if using forward Euler (value can be different on each process)

6019:    Note:
6020:    After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()

6022:    Level: intermediate

6024: .seealso: TSGetCFLTime(), TSADAPTCFL
6025: @*/
6026: PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime)
6027: {
6030:   ts->cfltime_local = cfltime;
6031:   ts->cfltime       = -1.;
6032:   return(0);
6033: }

6037: /*@
6038:    TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler

6040:    Collective on TS

6042:    Input Arguments:
6043: .  ts - time stepping context

6045:    Output Arguments:
6046: .  cfltime - maximum stable time step for forward Euler

6048:    Level: advanced

6050: .seealso: TSSetCFLTimeLocal()
6051: @*/
6052: PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime)
6053: {

6057:   if (ts->cfltime < 0) {
6058:     MPIU_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));
6059:   }
6060:   *cfltime = ts->cfltime;
6061:   return(0);
6062: }

6066: /*@
6067:    TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu

6069:    Input Parameters:
6070: .  ts   - the TS context.
6071: .  xl   - lower bound.
6072: .  xu   - upper bound.

6074:    Notes:
6075:    If this routine is not called then the lower and upper bounds are set to
6076:    PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().

6078:    Level: advanced

6080: @*/
6081: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
6082: {
6084:   SNES           snes;

6087:   TSGetSNES(ts,&snes);
6088:   SNESVISetVariableBounds(snes,xl,xu);
6089:   return(0);
6090: }

6092: #if defined(PETSC_HAVE_MATLAB_ENGINE)
6093: #include <mex.h>

6095: typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext;

6099: /*
6100:    TSComputeFunction_Matlab - Calls the function that has been set with
6101:                          TSSetFunctionMatlab().

6103:    Collective on TS

6105:    Input Parameters:
6106: +  snes - the TS context
6107: -  u - input vector

6109:    Output Parameter:
6110: .  y - function vector, as set by TSSetFunction()

6112:    Notes:
6113:    TSComputeFunction() is typically used within nonlinear solvers
6114:    implementations, so most users would not generally call this routine
6115:    themselves.

6117:    Level: developer

6119: .keywords: TS, nonlinear, compute, function

6121: .seealso: TSSetFunction(), TSGetFunction()
6122: */
6123: PetscErrorCode  TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx)
6124: {
6125:   PetscErrorCode  ierr;
6126:   TSMatlabContext *sctx = (TSMatlabContext*)ctx;
6127:   int             nlhs  = 1,nrhs = 7;
6128:   mxArray         *plhs[1],*prhs[7];
6129:   long long int   lx = 0,lxdot = 0,ly = 0,ls = 0;


6139:   PetscMemcpy(&ls,&snes,sizeof(snes));
6140:   PetscMemcpy(&lx,&u,sizeof(u));
6141:   PetscMemcpy(&lxdot,&udot,sizeof(udot));
6142:   PetscMemcpy(&ly,&y,sizeof(u));

6144:   prhs[0] =  mxCreateDoubleScalar((double)ls);
6145:   prhs[1] =  mxCreateDoubleScalar(time);
6146:   prhs[2] =  mxCreateDoubleScalar((double)lx);
6147:   prhs[3] =  mxCreateDoubleScalar((double)lxdot);
6148:   prhs[4] =  mxCreateDoubleScalar((double)ly);
6149:   prhs[5] =  mxCreateString(sctx->funcname);
6150:   prhs[6] =  sctx->ctx;
6151:    mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");
6152:    mxGetScalar(plhs[0]);
6153:   mxDestroyArray(prhs[0]);
6154:   mxDestroyArray(prhs[1]);
6155:   mxDestroyArray(prhs[2]);
6156:   mxDestroyArray(prhs[3]);
6157:   mxDestroyArray(prhs[4]);
6158:   mxDestroyArray(prhs[5]);
6159:   mxDestroyArray(plhs[0]);
6160:   return(0);
6161: }


6166: /*
6167:    TSSetFunctionMatlab - Sets the function evaluation routine and function
6168:    vector for use by the TS routines in solving ODEs
6169:    equations from MATLAB. Here the function is a string containing the name of a MATLAB function

6171:    Logically Collective on TS

6173:    Input Parameters:
6174: +  ts - the TS context
6175: -  func - function evaluation routine

6177:    Calling sequence of func:
6178: $    func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx);

6180:    Level: beginner

6182: .keywords: TS, nonlinear, set, function

6184: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
6185: */
6186: PetscErrorCode  TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx)
6187: {
6188:   PetscErrorCode  ierr;
6189:   TSMatlabContext *sctx;

6192:   /* currently sctx is memory bleed */
6193:   PetscMalloc(sizeof(TSMatlabContext),&sctx);
6194:   PetscStrallocpy(func,&sctx->funcname);
6195:   /*
6196:      This should work, but it doesn't
6197:   sctx->ctx = ctx;
6198:   mexMakeArrayPersistent(sctx->ctx);
6199:   */
6200:   sctx->ctx = mxDuplicateArray(ctx);

6202:   TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);
6203:   return(0);
6204: }

6208: /*
6209:    TSComputeJacobian_Matlab - Calls the function that has been set with
6210:                          TSSetJacobianMatlab().

6212:    Collective on TS

6214:    Input Parameters:
6215: +  ts - the TS context
6216: .  u - input vector
6217: .  A, B - the matrices
6218: -  ctx - user context

6220:    Level: developer

6222: .keywords: TS, nonlinear, compute, function

6224: .seealso: TSSetFunction(), TSGetFunction()
6225: @*/
6226: PetscErrorCode  TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat A,Mat B,void *ctx)
6227: {
6228:   PetscErrorCode  ierr;
6229:   TSMatlabContext *sctx = (TSMatlabContext*)ctx;
6230:   int             nlhs  = 2,nrhs = 9;
6231:   mxArray         *plhs[2],*prhs[9];
6232:   long long int   lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0;


6238:   /* call Matlab function in ctx with arguments u and y */

6240:   PetscMemcpy(&ls,&ts,sizeof(ts));
6241:   PetscMemcpy(&lx,&u,sizeof(u));
6242:   PetscMemcpy(&lxdot,&udot,sizeof(u));
6243:   PetscMemcpy(&lA,A,sizeof(u));
6244:   PetscMemcpy(&lB,B,sizeof(u));

6246:   prhs[0] =  mxCreateDoubleScalar((double)ls);
6247:   prhs[1] =  mxCreateDoubleScalar((double)time);
6248:   prhs[2] =  mxCreateDoubleScalar((double)lx);
6249:   prhs[3] =  mxCreateDoubleScalar((double)lxdot);
6250:   prhs[4] =  mxCreateDoubleScalar((double)shift);
6251:   prhs[5] =  mxCreateDoubleScalar((double)lA);
6252:   prhs[6] =  mxCreateDoubleScalar((double)lB);
6253:   prhs[7] =  mxCreateString(sctx->funcname);
6254:   prhs[8] =  sctx->ctx;
6255:    mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");
6256:    mxGetScalar(plhs[0]);
6257:   mxDestroyArray(prhs[0]);
6258:   mxDestroyArray(prhs[1]);
6259:   mxDestroyArray(prhs[2]);
6260:   mxDestroyArray(prhs[3]);
6261:   mxDestroyArray(prhs[4]);
6262:   mxDestroyArray(prhs[5]);
6263:   mxDestroyArray(prhs[6]);
6264:   mxDestroyArray(prhs[7]);
6265:   mxDestroyArray(plhs[0]);
6266:   mxDestroyArray(plhs[1]);
6267:   return(0);
6268: }


6273: /*
6274:    TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
6275:    vector for use by the TS routines in solving ODEs from MATLAB. Here the function is a string containing the name of a MATLAB function

6277:    Logically Collective on TS

6279:    Input Parameters:
6280: +  ts - the TS context
6281: .  A,B - Jacobian matrices
6282: .  func - function evaluation routine
6283: -  ctx - user context

6285:    Calling sequence of func:
6286: $    flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx);


6289:    Level: developer

6291: .keywords: TS, nonlinear, set, function

6293: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
6294: */
6295: PetscErrorCode  TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx)
6296: {
6297:   PetscErrorCode  ierr;
6298:   TSMatlabContext *sctx;

6301:   /* currently sctx is memory bleed */
6302:   PetscMalloc(sizeof(TSMatlabContext),&sctx);
6303:   PetscStrallocpy(func,&sctx->funcname);
6304:   /*
6305:      This should work, but it doesn't
6306:   sctx->ctx = ctx;
6307:   mexMakeArrayPersistent(sctx->ctx);
6308:   */
6309:   sctx->ctx = mxDuplicateArray(ctx);

6311:   TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);
6312:   return(0);
6313: }

6317: /*
6318:    TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab().

6320:    Collective on TS

6322: .seealso: TSSetFunction(), TSGetFunction()
6323: @*/
6324: PetscErrorCode  TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx)
6325: {
6326:   PetscErrorCode  ierr;
6327:   TSMatlabContext *sctx = (TSMatlabContext*)ctx;
6328:   int             nlhs  = 1,nrhs = 6;
6329:   mxArray         *plhs[1],*prhs[6];
6330:   long long int   lx = 0,ls = 0;


6336:   PetscMemcpy(&ls,&ts,sizeof(ts));
6337:   PetscMemcpy(&lx,&u,sizeof(u));

6339:   prhs[0] =  mxCreateDoubleScalar((double)ls);
6340:   prhs[1] =  mxCreateDoubleScalar((double)it);
6341:   prhs[2] =  mxCreateDoubleScalar((double)time);
6342:   prhs[3] =  mxCreateDoubleScalar((double)lx);
6343:   prhs[4] =  mxCreateString(sctx->funcname);
6344:   prhs[5] =  sctx->ctx;
6345:    mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");
6346:    mxGetScalar(plhs[0]);
6347:   mxDestroyArray(prhs[0]);
6348:   mxDestroyArray(prhs[1]);
6349:   mxDestroyArray(prhs[2]);
6350:   mxDestroyArray(prhs[3]);
6351:   mxDestroyArray(prhs[4]);
6352:   mxDestroyArray(plhs[0]);
6353:   return(0);
6354: }


6359: /*
6360:    TSMonitorSetMatlab - Sets the monitor function from Matlab

6362:    Level: developer

6364: .keywords: TS, nonlinear, set, function

6366: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
6367: */
6368: PetscErrorCode  TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx)
6369: {
6370:   PetscErrorCode  ierr;
6371:   TSMatlabContext *sctx;

6374:   /* currently sctx is memory bleed */
6375:   PetscMalloc(sizeof(TSMatlabContext),&sctx);
6376:   PetscStrallocpy(func,&sctx->funcname);
6377:   /*
6378:      This should work, but it doesn't
6379:   sctx->ctx = ctx;
6380:   mexMakeArrayPersistent(sctx->ctx);
6381:   */
6382:   sctx->ctx = mxDuplicateArray(ctx);

6384:   TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);
6385:   return(0);
6386: }
6387: #endif

6391: /*@C
6392:    TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector
6393:        in a time based line graph

6395:    Collective on TS

6397:    Input Parameters:
6398: +  ts - the TS context
6399: .  step - current time-step
6400: .  ptime - current time
6401: .  u - current solution
6402: -  dctx - the TSMonitorLGCtx object that contains all the options for the monitoring, this is created with TSMonitorLGCtxCreate()

6404:    Options Database:
6405: .   -ts_monitor_lg_solution_variables

6407:    Level: intermediate

6409:    Notes: Each process in a parallel run displays its component solutions in a separate window

6411: .keywords: TS,  vector, monitor, view

6413: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
6414:            TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
6415:            TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
6416:            TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()
6417: @*/
6418: PetscErrorCode  TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
6419: {
6420:   PetscErrorCode    ierr;
6421:   TSMonitorLGCtx    ctx = (TSMonitorLGCtx)dctx;
6422:   const PetscScalar *yy;
6423:   Vec               v;

6426:   if (step < 0) return(0); /* -1 indicates interpolated solution */
6427:   if (!step) {
6428:     PetscDrawAxis axis;
6429:     PetscInt      dim;
6430:     PetscDrawLGGetAxis(ctx->lg,&axis);
6431:     PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");
6432:     if (ctx->names && !ctx->displaynames) {
6433:       char      **displaynames;
6434:       PetscBool flg;
6435:       VecGetLocalSize(u,&dim);
6436:       PetscMalloc((dim+1)*sizeof(char*),&displaynames);
6437:       PetscMemzero(displaynames,(dim+1)*sizeof(char*));
6438:       PetscOptionsGetStringArray(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",displaynames,&dim,&flg);
6439:       if (flg) {
6440:         TSMonitorLGCtxSetDisplayVariables(ctx,(const char *const *)displaynames);
6441:       }
6442:       PetscStrArrayDestroy(&displaynames);
6443:     }
6444:     if (ctx->displaynames) {
6445:       PetscDrawLGSetDimension(ctx->lg,ctx->ndisplayvariables);
6446:       PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->displaynames);
6447:     } else if (ctx->names) {
6448:       VecGetLocalSize(u,&dim);
6449:       PetscDrawLGSetDimension(ctx->lg,dim);
6450:       PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->names);
6451:     } else {
6452:       VecGetLocalSize(u,&dim);
6453:       PetscDrawLGSetDimension(ctx->lg,dim);
6454:     }
6455:     PetscDrawLGReset(ctx->lg);
6456:   }

6458:   if (!ctx->transform) v = u;
6459:   else {(*ctx->transform)(ctx->transformctx,u,&v);}
6460:   VecGetArrayRead(v,&yy);
6461:   if (ctx->displaynames) {
6462:     PetscInt i;
6463:     for (i=0; i<ctx->ndisplayvariables; i++)
6464:       ctx->displayvalues[i] = PetscRealPart(yy[ctx->displayvariables[i]]);
6465:     PetscDrawLGAddCommonPoint(ctx->lg,ptime,ctx->displayvalues);
6466:   } else {
6467: #if defined(PETSC_USE_COMPLEX)
6468:     PetscInt  i,n;
6469:     PetscReal *yreal;
6470:     VecGetLocalSize(v,&n);
6471:     PetscMalloc1(n,&yreal);
6472:     for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
6473:     PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
6474:     PetscFree(yreal);
6475: #else
6476:     PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
6477: #endif
6478:   }
6479:   VecRestoreArrayRead(v,&yy);
6480:   if (ctx->transform) {VecDestroy(&v);}

6482:   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
6483:     PetscDrawLGDraw(ctx->lg);
6484:     PetscDrawLGSave(ctx->lg);
6485:   }
6486:   return(0);
6487: }


6492: /*@C
6493:    TSMonitorLGSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot

6495:    Collective on TS

6497:    Input Parameters:
6498: +  ts - the TS context
6499: -  names - the names of the components, final string must be NULL

6501:    Level: intermediate

6503:    Notes: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored

6505: .keywords: TS,  vector, monitor, view

6507: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetVariableNames()
6508: @*/
6509: PetscErrorCode  TSMonitorLGSetVariableNames(TS ts,const char * const *names)
6510: {
6511:   PetscErrorCode    ierr;
6512:   PetscInt          i;

6515:   for (i=0; i<ts->numbermonitors; i++) {
6516:     if (ts->monitor[i] == TSMonitorLGSolution) {
6517:       TSMonitorLGCtxSetVariableNames((TSMonitorLGCtx)ts->monitorcontext[i],names);
6518:       break;
6519:     }
6520:   }
6521:   return(0);
6522: }

6526: /*@C
6527:    TSMonitorLGCtxSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot

6529:    Collective on TS

6531:    Input Parameters:
6532: +  ts - the TS context
6533: -  names - the names of the components, final string must be NULL

6535:    Level: intermediate

6537: .keywords: TS,  vector, monitor, view

6539: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGSetVariableNames()
6540: @*/
6541: PetscErrorCode  TSMonitorLGCtxSetVariableNames(TSMonitorLGCtx ctx,const char * const *names)
6542: {
6543:   PetscErrorCode    ierr;

6546:   PetscStrArrayDestroy(&ctx->names);
6547:   PetscStrArrayallocpy(names,&ctx->names);
6548:   return(0);
6549: }

6553: /*@C
6554:    TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot

6556:    Collective on TS

6558:    Input Parameter:
6559: .  ts - the TS context

6561:    Output Parameter:
6562: .  names - the names of the components, final string must be NULL

6564:    Level: intermediate

6566:    Notes: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored

6568: .keywords: TS,  vector, monitor, view

6570: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
6571: @*/
6572: PetscErrorCode  TSMonitorLGGetVariableNames(TS ts,const char *const **names)
6573: {
6574:   PetscInt       i;

6577:   *names = NULL;
6578:   for (i=0; i<ts->numbermonitors; i++) {
6579:     if (ts->monitor[i] == TSMonitorLGSolution) {
6580:       TSMonitorLGCtx  ctx = (TSMonitorLGCtx) ts->monitorcontext[i];
6581:       *names = (const char *const *)ctx->names;
6582:       break;
6583:     }
6584:   }
6585:   return(0);
6586: }

6590: /*@C
6591:    TSMonitorLGCtxSetDisplayVariables - Sets the variables that are to be display in the monitor

6593:    Collective on TS

6595:    Input Parameters:
6596: +  ctx - the TSMonitorLG context
6597: .  displaynames - the names of the components, final string must be NULL

6599:    Level: intermediate

6601: .keywords: TS,  vector, monitor, view

6603: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
6604: @*/
6605: PetscErrorCode  TSMonitorLGCtxSetDisplayVariables(TSMonitorLGCtx ctx,const char * const *displaynames)
6606: {
6607:   PetscInt          j = 0,k;
6608:   PetscErrorCode    ierr;

6611:   if (!ctx->names) return(0);
6612:   PetscStrArrayDestroy(&ctx->displaynames);
6613:   PetscStrArrayallocpy(displaynames,&ctx->displaynames);
6614:   while (displaynames[j]) j++;
6615:   ctx->ndisplayvariables = j;
6616:   PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvariables);
6617:   PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvalues);
6618:   j = 0;
6619:   while (displaynames[j]) {
6620:     k = 0;
6621:     while (ctx->names[k]) {
6622:       PetscBool flg;
6623:       PetscStrcmp(displaynames[j],ctx->names[k],&flg);
6624:       if (flg) {
6625:         ctx->displayvariables[j] = k;
6626:         break;
6627:       }
6628:       k++;
6629:     }
6630:     j++;
6631:   }
6632:   return(0);
6633: }


6638: /*@C
6639:    TSMonitorLGSetDisplayVariables - Sets the variables that are to be display in the monitor

6641:    Collective on TS

6643:    Input Parameters:
6644: +  ts - the TS context
6645: .  displaynames - the names of the components, final string must be NULL

6647:    Notes: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored

6649:    Level: intermediate

6651: .keywords: TS,  vector, monitor, view

6653: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
6654: @*/
6655: PetscErrorCode  TSMonitorLGSetDisplayVariables(TS ts,const char * const *displaynames)
6656: {
6657:   PetscInt          i;
6658:   PetscErrorCode    ierr;

6661:   for (i=0; i<ts->numbermonitors; i++) {
6662:     if (ts->monitor[i] == TSMonitorLGSolution) {
6663:       TSMonitorLGCtxSetDisplayVariables((TSMonitorLGCtx)ts->monitorcontext[i],displaynames);
6664:       break;
6665:     }
6666:   }
6667:   return(0);
6668: }

6672: /*@C
6673:    TSMonitorLGSetTransform - Solution vector will be transformed by provided function before being displayed

6675:    Collective on TS

6677:    Input Parameters:
6678: +  ts - the TS context
6679: .  transform - the transform function
6680: .  destroy - function to destroy the optional context
6681: -  ctx - optional context used by transform function

6683:    Notes: If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored

6685:    Level: intermediate

6687: .keywords: TS,  vector, monitor, view

6689: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGCtxSetTransform()
6690: @*/
6691: PetscErrorCode  TSMonitorLGSetTransform(TS ts,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
6692: {
6693:   PetscInt          i;
6694:   PetscErrorCode    ierr;

6697:   for (i=0; i<ts->numbermonitors; i++) {
6698:     if (ts->monitor[i] == TSMonitorLGSolution) {
6699:       TSMonitorLGCtxSetTransform((TSMonitorLGCtx)ts->monitorcontext[i],transform,destroy,tctx);
6700:     }
6701:   }
6702:   return(0);
6703: }

6707: /*@C
6708:    TSMonitorLGCtxSetTransform - Solution vector will be transformed by provided function before being displayed

6710:    Collective on TSLGCtx

6712:    Input Parameters:
6713: +  ts - the TS context
6714: .  transform - the transform function
6715: .  destroy - function to destroy the optional context
6716: -  ctx - optional context used by transform function

6718:    Level: intermediate

6720: .keywords: TS,  vector, monitor, view

6722: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGSetTransform()
6723: @*/
6724: PetscErrorCode  TSMonitorLGCtxSetTransform(TSMonitorLGCtx ctx,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
6725: {
6727:   ctx->transform    = transform;
6728:   ctx->transformdestroy = destroy;
6729:   ctx->transformctx = tctx;
6730:   return(0);
6731: }

6735: /*@C
6736:    TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector
6737:        in a time based line graph

6739:    Collective on TS

6741:    Input Parameters:
6742: +  ts - the TS context
6743: .  step - current time-step
6744: .  ptime - current time
6745: .  u - current solution
6746: -  dctx - TSMonitorLGCtx object created with TSMonitorLGCtxCreate()

6748:    Level: intermediate

6750:    Notes: Each process in a parallel run displays its component errors in a separate window

6752:    The user must provide the solution using TSSetSolutionFunction() to use this monitor.

6754:    Options Database Keys:
6755: .  -ts_monitor_lg_error - create a graphical monitor of error history

6757: .keywords: TS,  vector, monitor, view

6759: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
6760: @*/
6761: PetscErrorCode  TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
6762: {
6763:   PetscErrorCode    ierr;
6764:   TSMonitorLGCtx    ctx = (TSMonitorLGCtx)dummy;
6765:   const PetscScalar *yy;
6766:   Vec               y;

6769:   if (!step) {
6770:     PetscDrawAxis axis;
6771:     PetscInt      dim;
6772:     PetscDrawLGGetAxis(ctx->lg,&axis);
6773:     PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");
6774:     VecGetLocalSize(u,&dim);
6775:     PetscDrawLGSetDimension(ctx->lg,dim);
6776:     PetscDrawLGReset(ctx->lg);
6777:   }
6778:   VecDuplicate(u,&y);
6779:   TSComputeSolutionFunction(ts,ptime,y);
6780:   VecAXPY(y,-1.0,u);
6781:   VecGetArrayRead(y,&yy);
6782: #if defined(PETSC_USE_COMPLEX)
6783:   {
6784:     PetscReal *yreal;
6785:     PetscInt  i,n;
6786:     VecGetLocalSize(y,&n);
6787:     PetscMalloc1(n,&yreal);
6788:     for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
6789:     PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
6790:     PetscFree(yreal);
6791:   }
6792: #else
6793:   PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
6794: #endif
6795:   VecRestoreArrayRead(y,&yy);
6796:   VecDestroy(&y);
6797:   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
6798:     PetscDrawLGDraw(ctx->lg);
6799:     PetscDrawLGSave(ctx->lg);
6800:   }
6801:   return(0);
6802: }

6806: PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
6807: {
6808:   TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
6809:   PetscReal      x   = ptime,y;
6811:   PetscInt       its;

6814:   if (n < 0) return(0); /* -1 indicates interpolated solution */
6815:   if (!n) {
6816:     PetscDrawAxis axis;
6817:     PetscDrawLGGetAxis(ctx->lg,&axis);
6818:     PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");
6819:     PetscDrawLGReset(ctx->lg);
6820:     ctx->snes_its = 0;
6821:   }
6822:   TSGetSNESIterations(ts,&its);
6823:   y    = its - ctx->snes_its;
6824:   PetscDrawLGAddPoint(ctx->lg,&x,&y);
6825:   if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
6826:     PetscDrawLGDraw(ctx->lg);
6827:     PetscDrawLGSave(ctx->lg);
6828:   }
6829:   ctx->snes_its = its;
6830:   return(0);
6831: }

6835: PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
6836: {
6837:   TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
6838:   PetscReal      x   = ptime,y;
6840:   PetscInt       its;

6843:   if (n < 0) return(0); /* -1 indicates interpolated solution */
6844:   if (!n) {
6845:     PetscDrawAxis axis;
6846:     PetscDrawLGGetAxis(ctx->lg,&axis);
6847:     PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");
6848:     PetscDrawLGReset(ctx->lg);
6849:     ctx->ksp_its = 0;
6850:   }
6851:   TSGetKSPIterations(ts,&its);
6852:   y    = its - ctx->ksp_its;
6853:   PetscDrawLGAddPoint(ctx->lg,&x,&y);
6854:   if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
6855:     PetscDrawLGDraw(ctx->lg);
6856:     PetscDrawLGSave(ctx->lg);
6857:   }
6858:   ctx->ksp_its = its;
6859:   return(0);
6860: }

6864: /*@
6865:    TSComputeLinearStability - computes the linear stability function at a point

6867:    Collective on TS and Vec

6869:    Input Parameters:
6870: +  ts - the TS context
6871: -  xr,xi - real and imaginary part of input arguments

6873:    Output Parameters:
6874: .  yr,yi - real and imaginary part of function value

6876:    Level: developer

6878: .keywords: TS, compute

6880: .seealso: TSSetRHSFunction(), TSComputeIFunction()
6881: @*/
6882: PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
6883: {

6888:   if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method");
6889:   (*ts->ops->linearstability)(ts,xr,xi,yr,yi);
6890:   return(0);
6891: }

6893: /* ------------------------------------------------------------------------*/
6896: /*@C
6897:    TSMonitorEnvelopeCtxCreate - Creates a context for use with TSMonitorEnvelope()

6899:    Collective on TS

6901:    Input Parameters:
6902: .  ts  - the ODE solver object

6904:    Output Parameter:
6905: .  ctx - the context

6907:    Level: intermediate

6909: .keywords: TS, monitor, line graph, residual, seealso

6911: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()

6913: @*/
6914: PetscErrorCode  TSMonitorEnvelopeCtxCreate(TS ts,TSMonitorEnvelopeCtx *ctx)
6915: {

6919:   PetscNew(ctx);
6920:   return(0);
6921: }

6925: /*@C
6926:    TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution

6928:    Collective on TS

6930:    Input Parameters:
6931: +  ts - the TS context
6932: .  step - current time-step
6933: .  ptime - current time
6934: .  u  - current solution
6935: -  dctx - the envelope context

6937:    Options Database:
6938: .  -ts_monitor_envelope

6940:    Level: intermediate

6942:    Notes: after a solve you can use TSMonitorEnvelopeGetBounds() to access the envelope

6944: .keywords: TS,  vector, monitor, view

6946: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxCreate()
6947: @*/
6948: PetscErrorCode  TSMonitorEnvelope(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
6949: {
6950:   PetscErrorCode       ierr;
6951:   TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dctx;

6954:   if (!ctx->max) {
6955:     VecDuplicate(u,&ctx->max);
6956:     VecDuplicate(u,&ctx->min);
6957:     VecCopy(u,ctx->max);
6958:     VecCopy(u,ctx->min);
6959:   } else {
6960:     VecPointwiseMax(ctx->max,u,ctx->max);
6961:     VecPointwiseMin(ctx->min,u,ctx->min);
6962:   }
6963:   return(0);
6964: }


6969: /*@C
6970:    TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution

6972:    Collective on TS

6974:    Input Parameter:
6975: .  ts - the TS context

6977:    Output Parameter:
6978: +  max - the maximum values
6979: -  min - the minimum values

6981:    Notes: If the TS does not have a TSMonitorEnvelopeCtx associated with it then this function is ignored

6983:    Level: intermediate

6985: .keywords: TS,  vector, monitor, view

6987: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
6988: @*/
6989: PetscErrorCode  TSMonitorEnvelopeGetBounds(TS ts,Vec *max,Vec *min)
6990: {
6991:   PetscInt i;

6994:   if (max) *max = NULL;
6995:   if (min) *min = NULL;
6996:   for (i=0; i<ts->numbermonitors; i++) {
6997:     if (ts->monitor[i] == TSMonitorEnvelope) {
6998:       TSMonitorEnvelopeCtx  ctx = (TSMonitorEnvelopeCtx) ts->monitorcontext[i];
6999:       if (max) *max = ctx->max;
7000:       if (min) *min = ctx->min;
7001:       break;
7002:     }
7003:   }
7004:   return(0);
7005: }

7009: /*@C
7010:    TSMonitorEnvelopeCtxDestroy - Destroys a context that was created  with TSMonitorEnvelopeCtxCreate().

7012:    Collective on TSMonitorEnvelopeCtx

7014:    Input Parameter:
7015: .  ctx - the monitor context

7017:    Level: intermediate

7019: .keywords: TS, monitor, line graph, destroy

7021: .seealso: TSMonitorLGCtxCreate(),  TSMonitorSet(), TSMonitorLGTimeStep()
7022: @*/
7023: PetscErrorCode  TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx)
7024: {

7028:   VecDestroy(&(*ctx)->min);
7029:   VecDestroy(&(*ctx)->max);
7030:   PetscFree(*ctx);
7031:   return(0);
7032: }

7036: /*@
7037:    TSRollBack - Rolls back one time step

7039:    Collective on TS

7041:    Input Parameter:
7042: .  ts - the TS context obtained from TSCreate()

7044:    Level: advanced

7046: .keywords: TS, timestep, rollback

7048: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate()
7049: @*/
7050: PetscErrorCode  TSRollBack(TS ts)
7051: {

7056:   if (ts->steprollback) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TSRollBack already called");
7057:   if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name);
7058:   (*ts->ops->rollback)(ts);
7059:   ts->time_step = ts->ptime - ts->ptime_prev;
7060:   ts->ptime = ts->ptime_prev;
7061:   ts->ptime_prev = ts->ptime_prev_rollback;
7062:   ts->steps--; ts->total_steps--;
7063:   ts->steprollback = PETSC_TRUE;
7064:   return(0);
7065: }

7069: /*@
7070:    TSGetStages - Get the number of stages and stage values

7072:    Input Parameter:
7073: .  ts - the TS context obtained from TSCreate()

7075:    Level: advanced

7077: .keywords: TS, getstages

7079: .seealso: TSCreate()
7080: @*/
7081: PetscErrorCode  TSGetStages(TS ts,PetscInt *ns,Vec **Y)
7082: {


7089:   if (!ts->ops->getstages) *ns=0;
7090:   else {
7091:     (*ts->ops->getstages)(ts,ns,Y);
7092:   }
7093:   return(0);
7094: }

7098: /*@C
7099:   TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.

7101:   Collective on SNES

7103:   Input Parameters:
7104: + ts - the TS context
7105: . t - current timestep
7106: . U - state vector
7107: . Udot - time derivative of state vector
7108: . shift - shift to apply, see note below
7109: - ctx - an optional user context

7111:   Output Parameters:
7112: + J - Jacobian matrix (not altered in this routine)
7113: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)

7115:   Level: intermediate

7117:   Notes:
7118:   If F(t,U,Udot)=0 is the DAE, the required Jacobian is

7120:   dF/dU + shift*dF/dUdot

7122:   Most users should not need to explicitly call this routine, as it
7123:   is used internally within the nonlinear solvers.

7125:   This will first try to get the coloring from the DM.  If the DM type has no coloring
7126:   routine, then it will try to get the coloring from the matrix.  This requires that the
7127:   matrix have nonzero entries precomputed.

7129: .keywords: TS, finite differences, Jacobian, coloring, sparse
7130: .seealso: TSSetIJacobian(), MatFDColoringCreate(), MatFDColoringSetFunction()
7131: @*/
7132: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat J,Mat B,void *ctx)
7133: {
7134:   SNES           snes;
7135:   MatFDColoring  color;
7136:   PetscBool      hascolor, matcolor = PETSC_FALSE;

7140:   PetscOptionsGetBool(((PetscObject)ts)->options,((PetscObject) ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL);
7141:   PetscObjectQuery((PetscObject) B, "TSMatFDColoring", (PetscObject *) &color);
7142:   if (!color) {
7143:     DM         dm;
7144:     ISColoring iscoloring;

7146:     TSGetDM(ts, &dm);
7147:     DMHasColoring(dm, &hascolor);
7148:     if (hascolor && !matcolor) {
7149:       DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring);
7150:       MatFDColoringCreate(B, iscoloring, &color);
7151:       MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);
7152:       MatFDColoringSetFromOptions(color);
7153:       MatFDColoringSetUp(B, iscoloring, color);
7154:       ISColoringDestroy(&iscoloring);
7155:     } else {
7156:       MatColoring mc;

7158:       MatColoringCreate(B, &mc);
7159:       MatColoringSetDistance(mc, 2);
7160:       MatColoringSetType(mc, MATCOLORINGSL);
7161:       MatColoringSetFromOptions(mc);
7162:       MatColoringApply(mc, &iscoloring);
7163:       MatColoringDestroy(&mc);
7164:       MatFDColoringCreate(B, iscoloring, &color);
7165:       MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);
7166:       MatFDColoringSetFromOptions(color);
7167:       MatFDColoringSetUp(B, iscoloring, color);
7168:       ISColoringDestroy(&iscoloring);
7169:     }
7170:     PetscObjectCompose((PetscObject) B, "TSMatFDColoring", (PetscObject) color);
7171:     PetscObjectDereference((PetscObject) color);
7172:   }
7173:   TSGetSNES(ts, &snes);
7174:   MatFDColoringApply(B, color, U, snes);
7175:   if (J != B) {
7176:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
7177:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
7178:   }
7179:   return(0);
7180: }

7184: /*@
7185:     TSSetFunctionDomainError - Set the function testing if the current state vector is valid

7187:     Input Parameters:
7188:     ts - the TS context
7189:     func - function called within TSFunctionDomainError

7191:     Level: intermediate

7193: .keywords: TS, state, domain
7194: .seealso: TSAdaptCheckStage(), TSFunctionDomainError()
7195: @*/

7197: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS,PetscReal,Vec,PetscBool*))
7198: {
7201:   ts->functiondomainerror = func;
7202:   return(0);
7203: }

7207: /*@
7208:     TSFunctionDomainError - Check if the current state is valid

7210:     Input Parameters:
7211:     ts - the TS context
7212:     stagetime - time of the simulation
7213:     Y - state vector to check.

7215:     Output Parameter:
7216:     accept - Set to PETSC_FALSE if the current state vector is valid.

7218:     Note:
7219:     This function should be used to ensure the state is in a valid part of the space.
7220:     For example, one can ensure here all values are positive.

7222:     Level: advanced
7223: @*/
7224: PetscErrorCode TSFunctionDomainError(TS ts,PetscReal stagetime,Vec Y,PetscBool* accept)
7225: {


7231:   *accept = PETSC_TRUE;
7232:   if (ts->functiondomainerror) {
7233:     PetscStackCallStandard((*ts->functiondomainerror),(ts,stagetime,Y,accept));
7234:   }
7235:   return(0);
7236: }

7238: #undef  __FUNCT__
7240: /*@C
7241:   TSClone - This function clones a time step object. 

7243:   Collective on MPI_Comm

7245:   Input Parameter:
7246: . tsin    - The input TS

7248:   Output Parameter:
7249: . tsout   - The output TS (cloned)

7251:   Notes:
7252:   This function is used to create a clone of a TS object. It is used in ARKIMEX for initializing the slope for first stage explicit methods. It will likely be replaced in the future with a mechanism of switching methods on the fly. 

7254:   When using TSDestroy() on a clone the user has to first reset the correct TS reference in the embedded SNES object: e.g.: by running SNES snes_dup=NULL; TSGetSNES(ts,&snes_dup); TSSetSNES(ts,snes_dup);

7256:   Level: developer

7258: .keywords: TS, clone
7259: .seealso: TSCreate(), TSSetType(), TSSetUp(), TSDestroy(), TSSetProblemType()
7260: @*/
7261: PetscErrorCode  TSClone(TS tsin, TS *tsout)
7262: {
7263:   TS             t;
7265:   SNES           snes_start;
7266:   DM             dm;
7267:   TSType         type;

7271:   *tsout = NULL;

7273:   PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView);

7275:   /* General TS description */
7276:   t->numbermonitors    = 0;
7277:   t->setupcalled       = 0;
7278:   t->ksp_its           = 0;
7279:   t->snes_its          = 0;
7280:   t->nwork             = 0;
7281:   t->rhsjacobian.time  = -1e20;
7282:   t->rhsjacobian.scale = 1.;
7283:   t->ijacobian.shift   = 1.;

7285:   TSGetSNES(tsin,&snes_start);
7286:   TSSetSNES(t,snes_start);

7288:   TSGetDM(tsin,&dm);
7289:   TSSetDM(t,dm);

7291:   t->adapt = tsin->adapt;
7292:   PetscObjectReference((PetscObject)t->adapt);

7294:   t->problem_type      = tsin->problem_type;
7295:   t->ptime             = tsin->ptime;
7296:   t->time_step         = tsin->time_step;
7297:   t->max_time          = tsin->max_time;
7298:   t->steps             = tsin->steps;
7299:   t->max_steps         = tsin->max_steps;
7300:   t->equation_type     = tsin->equation_type;
7301:   t->atol              = tsin->atol;
7302:   t->rtol              = tsin->rtol;
7303:   t->max_snes_failures = tsin->max_snes_failures;
7304:   t->max_reject        = tsin->max_reject;
7305:   t->errorifstepfailed = tsin->errorifstepfailed;

7307:   TSGetType(tsin,&type);
7308:   TSSetType(t,type);

7310:   t->vec_sol           = NULL;

7312:   t->cfltime          = tsin->cfltime;
7313:   t->cfltime_local    = tsin->cfltime_local;
7314:   t->exact_final_time = tsin->exact_final_time;

7316:   PetscMemcpy(t->ops,tsin->ops,sizeof(struct _TSOps));

7318:   if (((PetscObject)tsin)->fortran_func_pointers) {
7319:     PetscInt i;
7320:     PetscMalloc((10)*sizeof(void(*)(void)),&((PetscObject)t)->fortran_func_pointers);
7321:     for (i=0; i<10; i++) {
7322:       ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
7323:     }
7324:   }
7325:   *tsout = t;
7326:   return(0);
7327: }