Actual source code: ts.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1:  #include <petsc/private/tsimpl.h>
  2:  #include <petscdmshell.h>
  3:  #include <petscdmda.h>
  4:  #include <petscviewer.h>
  5:  #include <petscdraw.h>

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

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

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

 16:    Collective on TS

 18:    Input Parameters:
 19: +  ts - TS object you wish to monitor
 20: .  name - the monitor type one is seeking
 21: .  help - message indicating what monitoring is done
 22: .  manual - manual page for the monitor
 23: .  monitor - the monitor function
 24: -  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

 26:    Level: developer

 28: .seealso: PetscOptionsGetViewer(), PetscOptionsGetReal(), PetscOptionsHasName(), PetscOptionsGetString(),
 29:           PetscOptionsGetIntArray(), PetscOptionsGetRealArray(), PetscOptionsBool()
 30:           PetscOptionsInt(), PetscOptionsString(), PetscOptionsReal(), PetscOptionsBool(),
 31:           PetscOptionsName(), PetscOptionsBegin(), PetscOptionsEnd(), PetscOptionsHead(),
 32:           PetscOptionsStringArray(),PetscOptionsRealArray(), PetscOptionsScalar(),
 33:           PetscOptionsBoolGroupBegin(), PetscOptionsBoolGroup(), PetscOptionsBoolGroupEnd(),
 34:           PetscOptionsFList(), PetscOptionsEList()
 35: @*/
 36: PetscErrorCode  TSMonitorSetFromOptions(TS ts,const char name[],const char help[], const char manual[],PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,PetscViewerAndFormat*),PetscErrorCode (*monitorsetup)(TS,PetscViewerAndFormat*))
 37: {
 38:   PetscErrorCode    ierr;
 39:   PetscViewer       viewer;
 40:   PetscViewerFormat format;
 41:   PetscBool         flg;

 44:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject) ts)->options,((PetscObject)ts)->prefix,name,&viewer,&format,&flg);
 45:   if (flg) {
 46:     PetscViewerAndFormat *vf;
 47:     PetscViewerAndFormatCreate(viewer,format,&vf);
 48:     PetscObjectDereference((PetscObject)viewer);
 49:     if (monitorsetup) {
 50:       (*monitorsetup)(ts,vf);
 51:     }
 52:     TSMonitorSet(ts,(PetscErrorCode (*)(TS,PetscInt,PetscReal,Vec,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
 53:   }
 54:   return(0);
 55: }

 57: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt,TSAdaptType default_type)
 58: {

 64:   if (!((PetscObject)adapt)->type_name) {
 65:     TSAdaptSetType(adapt,default_type);
 66:   }
 67:   return(0);
 68: }

 70: /*@
 71:    TSSetFromOptions - Sets various TS parameters from user options.

 73:    Collective on TS

 75:    Input Parameter:
 76: .  ts - the TS context obtained from TSCreate()

 78:    Options Database Keys:
 79: +  -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSALPHA, TSGLLE, TSSSP, TSGLEE, TSBSYMP
 80: .  -ts_save_trajectory - checkpoint the solution at each time-step
 81: .  -ts_max_time <time> - maximum time to compute to
 82: .  -ts_max_steps <steps> - maximum number of time-steps to take
 83: .  -ts_init_time <time> - initial time to start computation
 84: .  -ts_final_time <time> - final time to compute to (deprecated: use -ts_max_time)
 85: .  -ts_dt <dt> - initial time step
 86: .  -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
 87: .  -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
 88: .  -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
 89: .  -ts_error_if_step_fails <true,false> - Error if no step succeeds
 90: .  -ts_rtol <rtol> - relative tolerance for local truncation error
 91: .  -ts_atol <atol> Absolute tolerance for local truncation error
 92: .  -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function
 93: .  -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - test the Jacobian at each iteration against finite difference with RHS function
 94: .  -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
 95: .  -ts_fd_color - Use finite differences with coloring to compute IJacobian
 96: .  -ts_monitor - print information at each timestep
 97: .  -ts_monitor_lg_solution - Monitor solution graphically
 98: .  -ts_monitor_lg_error - Monitor error graphically
 99: .  -ts_monitor_error - Monitors norm of error
100: .  -ts_monitor_lg_timestep - Monitor timestep size graphically
101: .  -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
102: .  -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
103: .  -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
104: .  -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
105: .  -ts_monitor_draw_solution - Monitor solution graphically
106: .  -ts_monitor_draw_solution_phase  <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
107: .  -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
108: .  -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
109: .  -ts_monitor_solution_vtk <filename.vts,filename.vtu> - Save each time step to a binary file, use filename-%%03D.vts (filename-%%03D.vtu)
110: .  -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time

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

114:    Level: beginner

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

118: .seealso: TSGetType()
119: @*/
120: PetscErrorCode  TSSetFromOptions(TS ts)
121: {
122:   PetscBool              opt,flg,tflg;
123:   PetscErrorCode         ierr;
124:   char                   monfilename[PETSC_MAX_PATH_LEN];
125:   PetscReal              time_step;
126:   TSExactFinalTimeOption eftopt;
127:   char                   dir[16];
128:   TSIFunction            ifun;
129:   const char             *defaultType;
130:   char                   typeName[256];


135:   TSRegisterAll();
136:   TSGetIFunction(ts,NULL,&ifun,NULL);

138:   PetscObjectOptionsBegin((PetscObject)ts);
139:   if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
140:   else defaultType = ifun ? TSBEULER : TSEULER;
141:   PetscOptionsFList("-ts_type","TS method","TSSetType",TSList,defaultType,typeName,256,&opt);
142:   if (opt) {
143:     TSSetType(ts,typeName);
144:   } else {
145:     TSSetType(ts,defaultType);
146:   }

148:   /* Handle generic TS options */
149:   PetscOptionsDeprecated("-ts_final_time","-ts_max_time","3.10",NULL);
150:   PetscOptionsReal("-ts_max_time","Maximum time to run to","TSSetMaxTime",ts->max_time,&ts->max_time,NULL);
151:   PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetMaxSteps",ts->max_steps,&ts->max_steps,NULL);
152:   PetscOptionsReal("-ts_init_time","Initial time","TSSetTime",ts->ptime,&ts->ptime,NULL);
153:   PetscOptionsReal("-ts_dt","Initial time step","TSSetTimeStep",ts->time_step,&time_step,&flg);
154:   if (flg) {TSSetTimeStep(ts,time_step);}
155:   PetscOptionsEnum("-ts_exact_final_time","Option for handling of final time step","TSSetExactFinalTime",TSExactFinalTimeOptions,(PetscEnum)ts->exact_final_time,(PetscEnum*)&eftopt,&flg);
156:   if (flg) {TSSetExactFinalTime(ts,eftopt);}
157:   PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","TSSetMaxSNESFailures",ts->max_snes_failures,&ts->max_snes_failures,NULL);
158:   PetscOptionsInt("-ts_max_reject","Maximum number of step rejections before step fails","TSSetMaxStepRejections",ts->max_reject,&ts->max_reject,NULL);
159:   PetscOptionsBool("-ts_error_if_step_fails","Error if no step succeeds","TSSetErrorIfStepFails",ts->errorifstepfailed,&ts->errorifstepfailed,NULL);
160:   PetscOptionsReal("-ts_rtol","Relative tolerance for local truncation error","TSSetTolerances",ts->rtol,&ts->rtol,NULL);
161:   PetscOptionsReal("-ts_atol","Absolute tolerance for local truncation error","TSSetTolerances",ts->atol,&ts->atol,NULL);

163:   PetscOptionsBool("-ts_rhs_jacobian_test_mult","Test the RHS Jacobian for consistency with RHS at each solve ","None",ts->testjacobian,&ts->testjacobian,NULL);
164:   PetscOptionsBool("-ts_rhs_jacobian_test_mult_transpose","Test the RHS Jacobian transpose for consistency with RHS at each solve ","None",ts->testjacobiantranspose,&ts->testjacobiantranspose,NULL);
165: #if defined(PETSC_HAVE_SAWS)
166:   {
167:   PetscBool set;
168:   flg  = PETSC_FALSE;
169:   PetscOptionsBool("-ts_saws_block","Block for SAWs memory snooper at end of TSSolve","PetscObjectSAWsBlock",((PetscObject)ts)->amspublishblock,&flg,&set);
170:   if (set) {
171:     PetscObjectSAWsSetBlock((PetscObject)ts,flg);
172:   }
173:   }
174: #endif

176:   /* Monitor options */
177:   TSMonitorSetFromOptions(ts,"-ts_monitor","Monitor time and timestep size","TSMonitorDefault",TSMonitorDefault,NULL);
178:   TSMonitorSetFromOptions(ts,"-ts_monitor_extreme","Monitor extreme values of the solution","TSMonitorExtreme",TSMonitorExtreme,NULL);
179:   TSMonitorSetFromOptions(ts,"-ts_monitor_solution","View the solution at each timestep","TSMonitorSolution",TSMonitorSolution,NULL);

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

184:   PetscOptionsName("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",&opt);
185:   if (opt) {
186:     TSMonitorLGCtx ctx;
187:     PetscInt       howoften = 1;

189:     PetscOptionsInt("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",howoften,&howoften,NULL);
190:     TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
191:     TSMonitorSet(ts,TSMonitorLGSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
192:   }

194:   PetscOptionsName("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",&opt);
195:   if (opt) {
196:     TSMonitorLGCtx ctx;
197:     PetscInt       howoften = 1;

199:     PetscOptionsInt("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",howoften,&howoften,NULL);
200:     TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
201:     TSMonitorSet(ts,TSMonitorLGError,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
202:   }
203:   TSMonitorSetFromOptions(ts,"-ts_monitor_error","View the error at each timestep","TSMonitorError",TSMonitorError,NULL);

205:   PetscOptionsName("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",&opt);
206:   if (opt) {
207:     TSMonitorLGCtx ctx;
208:     PetscInt       howoften = 1;

210:     PetscOptionsInt("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);
211:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
212:     TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
213:   }
214:   PetscOptionsName("-ts_monitor_lg_timestep_log","Monitor log timestep size graphically","TSMonitorLGTimeStep",&opt);
215:   if (opt) {
216:     TSMonitorLGCtx ctx;
217:     PetscInt       howoften = 1;

219:     PetscOptionsInt("-ts_monitor_lg_timestep_log","Monitor log timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);
220:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
221:     TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
222:     ctx->semilogy = PETSC_TRUE;
223:   }

225:   PetscOptionsName("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",&opt);
226:   if (opt) {
227:     TSMonitorLGCtx ctx;
228:     PetscInt       howoften = 1;

230:     PetscOptionsInt("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",howoften,&howoften,NULL);
231:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
232:     TSMonitorSet(ts,TSMonitorLGSNESIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
233:   }
234:   PetscOptionsName("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",&opt);
235:   if (opt) {
236:     TSMonitorLGCtx ctx;
237:     PetscInt       howoften = 1;

239:     PetscOptionsInt("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",howoften,&howoften,NULL);
240:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),NULL,NULL,PETSC_DECIDE,PETSC_DECIDE,400,300,howoften,&ctx);
241:     TSMonitorSet(ts,TSMonitorLGKSPIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);
242:   }
243:   PetscOptionsName("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",&opt);
244:   if (opt) {
245:     TSMonitorSPEigCtx ctx;
246:     PetscInt          howoften = 1;

248:     PetscOptionsInt("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",howoften,&howoften,NULL);
249:     TSMonitorSPEigCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
250:     TSMonitorSet(ts,TSMonitorSPEig,ctx,(PetscErrorCode (*)(void**))TSMonitorSPEigCtxDestroy);
251:   }
252:   PetscOptionsName("-ts_monitor_sp_swarm","Display particle phase from the DMSwarm","TSMonitorSPSwarm",&opt);
253:   if (opt) {
254:     TSMonitorSPCtx  ctx;
255:     PetscInt        howoften = 1;
256:     PetscOptionsInt("-ts_monitor_sp_swarm","Display particles phase from the DMSwarm","TSMonitorSPSwarm",howoften,&howoften,NULL);
257:     TSMonitorSPCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
258:     TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode (*)(void**))TSMonitorSPCtxDestroy);
259:   }
260:   opt  = PETSC_FALSE;
261:   PetscOptionsName("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",&opt);
262:   if (opt) {
263:     TSMonitorDrawCtx ctx;
264:     PetscInt         howoften = 1;

266:     PetscOptionsInt("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",howoften,&howoften,NULL);
267:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,"Computed Solution",PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
268:     TSMonitorSet(ts,TSMonitorDrawSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
269:   }
270:   opt  = PETSC_FALSE;
271:   PetscOptionsName("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",&opt);
272:   if (opt) {
273:     TSMonitorDrawCtx ctx;
274:     PetscReal        bounds[4];
275:     PetscInt         n = 4;
276:     PetscDraw        draw;
277:     PetscDrawAxis    axis;

279:     PetscOptionsRealArray("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",bounds,&n,NULL);
280:     if (n != 4) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Must provide bounding box of phase field");
281:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,1,&ctx);
282:     PetscViewerDrawGetDraw(ctx->viewer,0,&draw);
283:     PetscViewerDrawGetDrawAxis(ctx->viewer,0,&axis);
284:     PetscDrawAxisSetLimits(axis,bounds[0],bounds[2],bounds[1],bounds[3]);
285:     PetscDrawAxisSetLabels(axis,"Phase Diagram","Variable 1","Variable 2");
286:     TSMonitorSet(ts,TSMonitorDrawSolutionPhase,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
287:   }
288:   opt  = PETSC_FALSE;
289:   PetscOptionsName("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",&opt);
290:   if (opt) {
291:     TSMonitorDrawCtx ctx;
292:     PetscInt         howoften = 1;

294:     PetscOptionsInt("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",howoften,&howoften,NULL);
295:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,"Error",PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
296:     TSMonitorSet(ts,TSMonitorDrawError,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
297:   }
298:   opt  = PETSC_FALSE;
299:   PetscOptionsName("-ts_monitor_draw_solution_function","Monitor solution provided by TSMonitorSetSolutionFunction() graphically","TSMonitorDrawSolutionFunction",&opt);
300:   if (opt) {
301:     TSMonitorDrawCtx ctx;
302:     PetscInt         howoften = 1;

304:     PetscOptionsInt("-ts_monitor_draw_solution_function","Monitor solution provided by TSMonitorSetSolutionFunction() graphically","TSMonitorDrawSolutionFunction",howoften,&howoften,NULL);
305:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,"Solution provided by user function",PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);
306:     TSMonitorSet(ts,TSMonitorDrawSolutionFunction,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);
307:   }

309:   opt  = PETSC_FALSE;
310:   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);
311:   if (flg) {
312:     const char *ptr,*ptr2;
313:     char       *filetemplate;
314:     if (!monfilename[0]) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
315:     /* Do some cursory validation of the input. */
316:     PetscStrstr(monfilename,"%",(char**)&ptr);
317:     if (!ptr) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts");
318:     for (ptr++; ptr && *ptr; ptr++) {
319:       PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);
320:       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");
321:       if (ptr2) break;
322:     }
323:     PetscStrallocpy(monfilename,&filetemplate);
324:     TSMonitorSet(ts,TSMonitorSolutionVTK,filetemplate,(PetscErrorCode (*)(void**))TSMonitorSolutionVTKDestroy);
325:   }

327:   PetscOptionsString("-ts_monitor_dmda_ray","Display a ray of the solution","None","y=0",dir,16,&flg);
328:   if (flg) {
329:     TSMonitorDMDARayCtx *rayctx;
330:     int                  ray = 0;
331:     DMDADirection        ddir;
332:     DM                   da;
333:     PetscMPIInt          rank;

335:     if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
336:     if (dir[0] == 'x') ddir = DMDA_X;
337:     else if (dir[0] == 'y') ddir = DMDA_Y;
338:     else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir);
339:     sscanf(dir+2,"%d",&ray);

341:     PetscInfo2(((PetscObject)ts),"Displaying DMDA ray %c = %D\n",dir[0],ray);
342:     PetscNew(&rayctx);
343:     TSGetDM(ts,&da);
344:     DMDAGetRay(da,ddir,ray,&rayctx->ray,&rayctx->scatter);
345:     MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);
346:     if (!rank) {
347:       PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,600,300,&rayctx->viewer);
348:     }
349:     rayctx->lgctx = NULL;
350:     TSMonitorSet(ts,TSMonitorDMDARay,rayctx,TSMonitorDMDARayDestroy);
351:   }
352:   PetscOptionsString("-ts_monitor_lg_dmda_ray","Display a ray of the solution","None","x=0",dir,16,&flg);
353:   if (flg) {
354:     TSMonitorDMDARayCtx *rayctx;
355:     int                 ray = 0;
356:     DMDADirection       ddir;
357:     DM                  da;
358:     PetscInt            howoften = 1;

360:     if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
361:     if      (dir[0] == 'x') ddir = DMDA_X;
362:     else if (dir[0] == 'y') ddir = DMDA_Y;
363:     else SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
364:     sscanf(dir+2, "%d", &ray);

366:     PetscInfo2(((PetscObject) ts),"Displaying LG DMDA ray %c = %D\n", dir[0], ray);
367:     PetscNew(&rayctx);
368:     TSGetDM(ts, &da);
369:     DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);
370:     TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&rayctx->lgctx);
371:     TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy);
372:   }

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

378:     TSMonitorEnvelopeCtxCreate(ts,&ctx);
379:     TSMonitorSet(ts,TSMonitorEnvelope,ctx,(PetscErrorCode (*)(void**))TSMonitorEnvelopeCtxDestroy);
380:   }

382:   flg  = PETSC_FALSE;
383:   PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeJacobianDefaultColor", flg, &flg, NULL);
384:   if (flg) {
385:     DM   dm;
386:     DMTS tdm;

388:     TSGetDM(ts, &dm);
389:     DMGetDMTS(dm, &tdm);
390:     tdm->ijacobianctx = NULL;
391:     TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, 0);
392:     PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n");
393:   }

395:   /* Handle specific TS options */
396:   if (ts->ops->setfromoptions) {
397:     (*ts->ops->setfromoptions)(PetscOptionsObject,ts);
398:   }

400:   /* Handle TSAdapt options */
401:   TSGetAdapt(ts,&ts->adapt);
402:   TSAdaptSetDefaultType(ts->adapt,ts->default_adapt_type);
403:   TSAdaptSetFromOptions(PetscOptionsObject,ts->adapt);

405:   /* TS trajectory must be set after TS, since it may use some TS options above */
406:   tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
407:   PetscOptionsBool("-ts_save_trajectory","Save the solution at each timestep","TSSetSaveTrajectory",tflg,&tflg,NULL);
408:   if (tflg) {
409:     TSSetSaveTrajectory(ts);
410:   }

412:   TSAdjointSetFromOptions(PetscOptionsObject,ts);

414:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
415:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)ts);
416:   PetscOptionsEnd();

418:   if (ts->trajectory) {
419:     TSTrajectorySetFromOptions(ts->trajectory,ts);
420:   }

422:   /* why do we have to do this here and not during TSSetUp? */
423:   TSGetSNES(ts,&ts->snes);
424:   if (ts->problem_type == TS_LINEAR) {
425:     PetscObjectTypeCompareAny((PetscObject)ts->snes,&flg,SNESKSPONLY,SNESKSPTRANSPOSEONLY,"");
426:     if (!flg) { SNESSetType(ts->snes,SNESKSPONLY); }
427:   }
428:   SNESSetFromOptions(ts->snes);
429:   return(0);
430: }

432: /*@
433:    TSGetTrajectory - Gets the trajectory from a TS if it exists

435:    Collective on TS

437:    Input Parameters:
438: .  ts - the TS context obtained from TSCreate()

440:    Output Parameters;
441: .  tr - the TSTrajectory object, if it exists

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

445:    Level: advanced

447: .seealso: TSGetTrajectory(), TSAdjointSolve(), TSTrajectory, TSTrajectoryCreate()

449: .keywords: TS, set, checkpoint,
450: @*/
451: PetscErrorCode  TSGetTrajectory(TS ts,TSTrajectory *tr)
452: {
455:   *tr = ts->trajectory;
456:   return(0);
457: }

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

462:    Collective on TS

464:    Input Parameters:
465: .  ts - the TS context obtained from TSCreate()

467:    Options Database:
468: +  -ts_save_trajectory - saves the trajectory to a file
469: -  -ts_trajectory_type type

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

473:     The TSTRAJECTORYVISUALIZATION files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
474:    MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m

476:    Level: intermediate

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

480: .keywords: TS, set, checkpoint,
481: @*/
482: PetscErrorCode  TSSetSaveTrajectory(TS ts)
483: {

488:   if (!ts->trajectory) {
489:     TSTrajectoryCreate(PetscObjectComm((PetscObject)ts),&ts->trajectory);
490:   }
491:   return(0);
492: }

494: /*@
495:    TSResetTrajectory - Destroys and recreates the internal TSTrajectory object

497:    Collective on TS

499:    Input Parameters:
500: .  ts - the TS context obtained from TSCreate()

502:    Level: intermediate

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

506: .keywords: TS, set, checkpoint,
507: @*/
508: PetscErrorCode  TSResetTrajectory(TS ts)
509: {

514:   if (ts->trajectory) {
515:     TSTrajectoryDestroy(&ts->trajectory);
516:     TSTrajectoryCreate(PetscObjectComm((PetscObject)ts),&ts->trajectory);
517:   }
518:   return(0);
519: }

521: /*@
522:    TSComputeRHSJacobian - Computes the Jacobian matrix that has been
523:       set with TSSetRHSJacobian().

525:    Collective on TS and Vec

527:    Input Parameters:
528: +  ts - the TS context
529: .  t - current timestep
530: -  U - input vector

532:    Output Parameters:
533: +  A - Jacobian matrix
534: .  B - optional preconditioning matrix
535: -  flag - flag indicating matrix structure

537:    Notes:
538:    Most users should not need to explicitly call this routine, as it
539:    is used internally within the nonlinear solvers.

541:    See KSPSetOperators() for important information about setting the
542:    flag parameter.

544:    Level: developer

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

548: .seealso:  TSSetRHSJacobian(), KSPSetOperators()
549: @*/
550: PetscErrorCode  TSComputeRHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B)
551: {
552:   PetscErrorCode   ierr;
553:   PetscObjectState Ustate;
554:   PetscObjectId    Uid;
555:   DM               dm;
556:   DMTS             tsdm;
557:   TSRHSJacobian    rhsjacobianfunc;
558:   void             *ctx;
559:   TSIJacobian      ijacobianfunc;
560:   TSRHSFunction    rhsfunction;

566:   TSGetDM(ts,&dm);
567:   DMGetDMTS(dm,&tsdm);
568:   DMTSGetRHSJacobian(dm,&rhsjacobianfunc,&ctx);
569:   DMTSGetIJacobian(dm,&ijacobianfunc,NULL);
570:   DMTSGetRHSFunction(dm,&rhsfunction,&ctx);
571:   PetscObjectStateGet((PetscObject)U,&Ustate);
572:   PetscObjectGetId((PetscObject)U,&Uid);

574:   if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) {
575:     /* restore back RHS Jacobian matrices if they have been shifted/scaled */
576:     if (A == ts->Arhs) {
577:       if (ts->rhsjacobian.shift != 0) {
578:         MatShift(A,-ts->rhsjacobian.shift);
579:       }
580:       if (ts->rhsjacobian.scale != 1.) {
581:         MatScale(A,1./ts->rhsjacobian.scale);
582:       }
583:     }
584:     if (B && B == ts->Brhs && A != B) {
585:       if (ts->rhsjacobian.shift != 0) {
586:         MatShift(B,-ts->rhsjacobian.shift);
587:       }
588:       if (ts->rhsjacobian.scale != 1.) {
589:         MatScale(B,1./ts->rhsjacobian.scale);
590:       }
591:     }
592:     ts->rhsjacobian.shift = 0;
593:     ts->rhsjacobian.scale = 1.;
594:     return(0);
595:   }

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

599:   if (ts->rhsjacobian.reuse) {
600:     if (A == ts->Arhs) {
601:       /* MatScale has a short path for this case.
602:          However, this code path is taken the first time TSComputeRHSJacobian is called
603:          and the matrices have not assembled yet */
604:       if (ts->rhsjacobian.shift != 0) {
605:         MatShift(A,-ts->rhsjacobian.shift);
606:       }
607:       if (ts->rhsjacobian.scale != 1.) {
608:         MatScale(A,1./ts->rhsjacobian.scale);
609:       }
610:     }
611:     if (B && B == ts->Brhs && A != B) {
612:       if (ts->rhsjacobian.shift != 0) {
613:         MatShift(B,-ts->rhsjacobian.shift);
614:       }
615:       if (ts->rhsjacobian.scale != 1.) {
616:         MatScale(B,1./ts->rhsjacobian.scale);
617:       }
618:     }
619:   }

621:   if (rhsjacobianfunc) {
622:     PetscBool missing;
623:     PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);
624:     PetscStackPush("TS user Jacobian function");
625:     (*rhsjacobianfunc)(ts,t,U,A,B,ctx);
626:     PetscStackPop;
627:     PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);
628:     if (A) {
629:       MatMissingDiagonal(A,&missing,NULL);
630:       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");
631:     }
632:     if (B && B != A) {
633:       MatMissingDiagonal(B,&missing,NULL);
634:       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");
635:     }
636:   } else {
637:     MatZeroEntries(A);
638:     if (B && A != B) {MatZeroEntries(B);}
639:   }
640:   ts->rhsjacobian.time  = t;
641:   ts->rhsjacobian.shift = 0;
642:   ts->rhsjacobian.scale = 1.;
643:   PetscObjectGetId((PetscObject)U,&ts->rhsjacobian.Xid);
644:   PetscObjectStateGet((PetscObject)U,&ts->rhsjacobian.Xstate);
645:   return(0);
646: }

648: /*@
649:    TSComputeRHSFunction - Evaluates the right-hand-side function.

651:    Collective on TS and Vec

653:    Input Parameters:
654: +  ts - the TS context
655: .  t - current time
656: -  U - state vector

658:    Output Parameter:
659: .  y - right hand side

661:    Note:
662:    Most users should not need to explicitly call this routine, as it
663:    is used internally within the nonlinear solvers.

665:    Level: developer

667: .keywords: TS, compute

669: .seealso: TSSetRHSFunction(), TSComputeIFunction()
670: @*/
671: PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec U,Vec y)
672: {
674:   TSRHSFunction  rhsfunction;
675:   TSIFunction    ifunction;
676:   void           *ctx;
677:   DM             dm;

683:   TSGetDM(ts,&dm);
684:   DMTSGetRHSFunction(dm,&rhsfunction,&ctx);
685:   DMTSGetIFunction(dm,&ifunction,NULL);

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

689:   PetscLogEventBegin(TS_FunctionEval,ts,U,y,0);
690:   if (rhsfunction) {
691:     PetscStackPush("TS user right-hand-side function");
692:     (*rhsfunction)(ts,t,U,y,ctx);
693:     PetscStackPop;
694:   } else {
695:     VecZeroEntries(y);
696:   }

698:   PetscLogEventEnd(TS_FunctionEval,ts,U,y,0);
699:   return(0);
700: }

702: /*@
703:    TSComputeSolutionFunction - Evaluates the solution function.

705:    Collective on TS and Vec

707:    Input Parameters:
708: +  ts - the TS context
709: -  t - current time

711:    Output Parameter:
712: .  U - the solution

714:    Note:
715:    Most users should not need to explicitly call this routine, as it
716:    is used internally within the nonlinear solvers.

718:    Level: developer

720: .keywords: TS, compute

722: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
723: @*/
724: PetscErrorCode TSComputeSolutionFunction(TS ts,PetscReal t,Vec U)
725: {
726:   PetscErrorCode     ierr;
727:   TSSolutionFunction solutionfunction;
728:   void               *ctx;
729:   DM                 dm;

734:   TSGetDM(ts,&dm);
735:   DMTSGetSolutionFunction(dm,&solutionfunction,&ctx);

737:   if (solutionfunction) {
738:     PetscStackPush("TS user solution function");
739:     (*solutionfunction)(ts,t,U,ctx);
740:     PetscStackPop;
741:   }
742:   return(0);
743: }
744: /*@
745:    TSComputeForcingFunction - Evaluates the forcing function.

747:    Collective on TS and Vec

749:    Input Parameters:
750: +  ts - the TS context
751: -  t - current time

753:    Output Parameter:
754: .  U - the function value

756:    Note:
757:    Most users should not need to explicitly call this routine, as it
758:    is used internally within the nonlinear solvers.

760:    Level: developer

762: .keywords: TS, compute

764: .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction()
765: @*/
766: PetscErrorCode TSComputeForcingFunction(TS ts,PetscReal t,Vec U)
767: {
768:   PetscErrorCode     ierr, (*forcing)(TS,PetscReal,Vec,void*);
769:   void               *ctx;
770:   DM                 dm;

775:   TSGetDM(ts,&dm);
776:   DMTSGetForcingFunction(dm,&forcing,&ctx);

778:   if (forcing) {
779:     PetscStackPush("TS user forcing function");
780:     (*forcing)(ts,t,U,ctx);
781:     PetscStackPop;
782:   }
783:   return(0);
784: }

786: static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs)
787: {
788:   Vec            F;

792:   *Frhs = NULL;
793:   TSGetIFunction(ts,&F,NULL,NULL);
794:   if (!ts->Frhs) {
795:     VecDuplicate(F,&ts->Frhs);
796:   }
797:   *Frhs = ts->Frhs;
798:   return(0);
799: }

801: PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs)
802: {
803:   Mat            A,B;
805:   TSIJacobian    ijacobian;

808:   if (Arhs) *Arhs = NULL;
809:   if (Brhs) *Brhs = NULL;
810:   TSGetIJacobian(ts,&A,&B,&ijacobian,NULL);
811:   if (Arhs) {
812:     if (!ts->Arhs) {
813:       if (ijacobian) {
814:         MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);
815:       } else {
816:         ts->Arhs = A;
817:         PetscObjectReference((PetscObject)A);
818:       }
819:     } else {
820:       PetscBool flg;
821:       SNESGetUseMatrixFree(ts->snes,NULL,&flg);
822:       /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
823:       if (flg && !ijacobian && ts->Arhs == ts->Brhs){
824:         PetscObjectDereference((PetscObject)ts->Arhs);
825:         ts->Arhs = A;
826:         PetscObjectReference((PetscObject)A);
827:       }
828:     }
829:     *Arhs = ts->Arhs;
830:   }
831:   if (Brhs) {
832:     if (!ts->Brhs) {
833:       if (A != B) {
834:         if (ijacobian) {
835:           MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);
836:         } else {
837:           ts->Brhs = B;
838:           PetscObjectReference((PetscObject)B);
839:         }
840:       } else {
841:         PetscObjectReference((PetscObject)ts->Arhs);
842:         ts->Brhs = ts->Arhs;
843:       }
844:     }
845:     *Brhs = ts->Brhs;
846:   }
847:   return(0);
848: }

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

853:    Collective on TS and Vec

855:    Input Parameters:
856: +  ts - the TS context
857: .  t - current time
858: .  U - state vector
859: .  Udot - time derivative of state vector
860: -  imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate

862:    Output Parameter:
863: .  Y - right hand side

865:    Note:
866:    Most users should not need to explicitly call this routine, as it
867:    is used internally within the nonlinear solvers.

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

872:    Level: developer

874: .keywords: TS, compute

876: .seealso: TSSetIFunction(), TSComputeRHSFunction()
877: @*/
878: PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec Y,PetscBool imex)
879: {
881:   TSIFunction    ifunction;
882:   TSRHSFunction  rhsfunction;
883:   void           *ctx;
884:   DM             dm;


892:   TSGetDM(ts,&dm);
893:   DMTSGetIFunction(dm,&ifunction,&ctx);
894:   DMTSGetRHSFunction(dm,&rhsfunction,NULL);

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

898:   PetscLogEventBegin(TS_FunctionEval,ts,U,Udot,Y);
899:   if (ifunction) {
900:     PetscStackPush("TS user implicit function");
901:     (*ifunction)(ts,t,U,Udot,Y,ctx);
902:     PetscStackPop;
903:   }
904:   if (imex) {
905:     if (!ifunction) {
906:       VecCopy(Udot,Y);
907:     }
908:   } else if (rhsfunction) {
909:     if (ifunction) {
910:       Vec Frhs;
911:       TSGetRHSVec_Private(ts,&Frhs);
912:       TSComputeRHSFunction(ts,t,U,Frhs);
913:       VecAXPY(Y,-1,Frhs);
914:     } else {
915:       TSComputeRHSFunction(ts,t,U,Y);
916:       VecAYPX(Y,-1,Udot);
917:     }
918:   }
919:   PetscLogEventEnd(TS_FunctionEval,ts,U,Udot,Y);
920:   return(0);
921: }

923: /*@
924:    TSComputeIJacobian - Evaluates the Jacobian of the DAE

926:    Collective on TS and Vec

928:    Input
929:       Input Parameters:
930: +  ts - the TS context
931: .  t - current timestep
932: .  U - state vector
933: .  Udot - time derivative of state vector
934: .  shift - shift to apply, see note below
935: -  imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate

937:    Output Parameters:
938: +  A - Jacobian matrix
939: -  B - matrix from which the preconditioner is constructed; often the same as A

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

944:    dF/dU + shift*dF/dUdot

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

949:    Level: developer

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

953: .seealso:  TSSetIJacobian()
954: @*/
955: PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,PetscBool imex)
956: {
958:   TSIJacobian    ijacobian;
959:   TSRHSJacobian  rhsjacobian;
960:   DM             dm;
961:   void           *ctx;


972:   TSGetDM(ts,&dm);
973:   DMTSGetIJacobian(dm,&ijacobian,&ctx);
974:   DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);

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

978:   PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);
979:   if (ijacobian) {
980:     PetscBool missing;
981:     PetscStackPush("TS user implicit Jacobian");
982:     (*ijacobian)(ts,t,U,Udot,shift,A,B,ctx);
983:     PetscStackPop;
984:     MatMissingDiagonal(A,&missing,NULL);
985:     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");
986:     if (B != A) {
987:       MatMissingDiagonal(B,&missing,NULL);
988:       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");
989:     }
990:   }
991:   if (imex) {
992:     if (!ijacobian) {  /* system was written as Udot = G(t,U) */
993:       PetscBool assembled;
994:       if (rhsjacobian) {
995:         Mat Arhs = NULL;
996:         TSGetRHSMats_Private(ts,&Arhs,NULL);
997:         if (A == Arhs) {
998:           if (rhsjacobian == TSComputeRHSJacobianConstant) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Unsupported operation! cannot use TSComputeRHSJacobianConstant");
999:           ts->rhsjacobian.time = PETSC_MIN_REAL;
1000:         }
1001:       }
1002:       MatZeroEntries(A);
1003:       MatAssembled(A,&assembled);
1004:       if (!assembled) {
1005:         MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1006:         MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1007:       }
1008:       MatShift(A,shift);
1009:       if (A != B) {
1010:         MatZeroEntries(B);
1011:         MatAssembled(B,&assembled);
1012:         if (!assembled) {
1013:           MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1014:           MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1015:         }
1016:         MatShift(B,shift);
1017:       }
1018:     }
1019:   } else {
1020:     Mat Arhs = NULL,Brhs = NULL;
1021:     if (rhsjacobian) {
1022:       TSGetRHSMats_Private(ts,&Arhs,&Brhs);
1023:       TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);
1024:     }
1025:     if (Arhs == A) {           /* No IJacobian, so we only have the RHS matrix */
1026:       PetscBool flg;
1027:       ts->rhsjacobian.scale = -1;
1028:       ts->rhsjacobian.shift = shift;
1029:       SNESGetUseMatrixFree(ts->snes,NULL,&flg);
1030:       /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
1031:       if (!flg) {
1032:         MatScale(A,-1);
1033:         MatShift(A,shift);
1034:       }
1035:       if (A != B) {
1036:         MatScale(B,-1);
1037:         MatShift(B,shift);
1038:       }
1039:     } else if (Arhs) {          /* Both IJacobian and RHSJacobian */
1040:       MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
1041:       if (!ijacobian) {         /* No IJacobian provided, but we have a separate RHS matrix */
1042:         MatZeroEntries(A);
1043:         MatShift(A,shift);
1044:         if (A != B) {
1045:           MatZeroEntries(B);
1046:           MatShift(B,shift);
1047:         }
1048:       }
1049:       MatAXPY(A,-1,Arhs,axpy);
1050:       if (A != B) {
1051:         MatAXPY(B,-1,Brhs,axpy);
1052:       }
1053:     }
1054:   }
1055:   PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);
1056:   return(0);
1057: }

1059: /*@C
1060:     TSSetRHSFunction - Sets the routine for evaluating the function,
1061:     where U_t = G(t,u).

1063:     Logically Collective on TS

1065:     Input Parameters:
1066: +   ts - the TS context obtained from TSCreate()
1067: .   r - vector to put the computed right hand side (or NULL to have it created)
1068: .   f - routine for evaluating the right-hand-side function
1069: -   ctx - [optional] user-defined context for private data for the
1070:           function evaluation routine (may be NULL)

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

1075: +   t - current timestep
1076: .   u - input vector
1077: .   F - function vector
1078: -   ctx - [optional] user-defined function context

1080:     Level: beginner

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

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

1087: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSSetIFunction()
1088: @*/
1089: PetscErrorCode  TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx)
1090: {
1092:   SNES           snes;
1093:   Vec            ralloc = NULL;
1094:   DM             dm;


1100:   TSGetDM(ts,&dm);
1101:   DMTSSetRHSFunction(dm,f,ctx);
1102:   TSGetSNES(ts,&snes);
1103:   if (!r && !ts->dm && ts->vec_sol) {
1104:     VecDuplicate(ts->vec_sol,&ralloc);
1105:     r = ralloc;
1106:   }
1107:   SNESSetFunction(snes,r,SNESTSFormFunction,ts);
1108:   VecDestroy(&ralloc);
1109:   return(0);
1110: }

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

1115:     Logically Collective on TS

1117:     Input Parameters:
1118: +   ts - the TS context obtained from TSCreate()
1119: .   f - routine for evaluating the solution
1120: -   ctx - [optional] user-defined context for private data for the
1121:           function evaluation routine (may be NULL)

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

1126: +   t - current timestep
1127: .   u - output vector
1128: -   ctx - [optional] user-defined function context

1130:     Options Database:
1131: +  -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided TSSetSolutionFunction()
1132: -  -ts_monitor_draw_error - Monitor error graphically, requires user to have provided TSSetSolutionFunction()

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

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

1141:     Level: beginner

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

1145: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetForcingFunction(), TSSetSolution(), TSGetSolution(), TSMonitorLGError(), TSMonitorDrawError()
1146: @*/
1147: PetscErrorCode  TSSetSolutionFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx)
1148: {
1150:   DM             dm;

1154:   TSGetDM(ts,&dm);
1155:   DMTSSetSolutionFunction(dm,f,ctx);
1156:   return(0);
1157: }

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

1162:     Logically Collective on TS

1164:     Input Parameters:
1165: +   ts - the TS context obtained from TSCreate()
1166: .   func - routine for evaluating the forcing function
1167: -   ctx - [optional] user-defined context for private data for the
1168:           function evaluation routine (may be NULL)

1170:     Calling sequence of func:
1171: $     func (TS ts,PetscReal t,Vec f,void *ctx);

1173: +   t - current timestep
1174: .   f - output vector
1175: -   ctx - [optional] user-defined function context

1177:     Notes:
1178:     This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1179:     create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1180:     definition of the problem you are solving and hence possibly introducing bugs.

1182:     This replaces the ODE F(u,u_t,t) = 0 the TS is solving with F(u,u_t,t) - func(t) = 0

1184:     This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1185:     parameters can be passed in the ctx variable.

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

1189:     Level: beginner

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

1193: .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetSolutionFunction()
1194: @*/
1195: PetscErrorCode  TSSetForcingFunction(TS ts,TSForcingFunction func,void *ctx)
1196: {
1198:   DM             dm;

1202:   TSGetDM(ts,&dm);
1203:   DMTSSetForcingFunction(dm,func,ctx);
1204:   return(0);
1205: }

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

1211:    Logically Collective on TS

1213:    Input Parameters:
1214: +  ts  - the TS context obtained from TSCreate()
1215: .  Amat - (approximate) Jacobian matrix
1216: .  Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1217: .  f   - the Jacobian evaluation routine
1218: -  ctx - [optional] user-defined context for private data for the
1219:          Jacobian evaluation routine (may be NULL)

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

1224: +  t - current timestep
1225: .  u - input vector
1226: .  Amat - (approximate) Jacobian matrix
1227: .  Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1228: -  ctx - [optional] user-defined context for matrix evaluation routine

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

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

1236:    Level: beginner

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

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

1242: @*/
1243: PetscErrorCode  TSSetRHSJacobian(TS ts,Mat Amat,Mat Pmat,TSRHSJacobian f,void *ctx)
1244: {
1246:   SNES           snes;
1247:   DM             dm;
1248:   TSIJacobian    ijacobian;


1257:   TSGetDM(ts,&dm);
1258:   DMTSSetRHSJacobian(dm,f,ctx);
1259:   if (f == TSComputeRHSJacobianConstant) {
1260:     /* Handle this case automatically for the user; otherwise user should call themselves. */
1261:     TSRHSJacobianSetReuse(ts,PETSC_TRUE);
1262:   }
1263:   DMTSGetIJacobian(dm,&ijacobian,NULL);
1264:   TSGetSNES(ts,&snes);
1265:   if (!ijacobian) {
1266:     SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
1267:   }
1268:   if (Amat) {
1269:     PetscObjectReference((PetscObject)Amat);
1270:     MatDestroy(&ts->Arhs);
1271:     ts->Arhs = Amat;
1272:   }
1273:   if (Pmat) {
1274:     PetscObjectReference((PetscObject)Pmat);
1275:     MatDestroy(&ts->Brhs);
1276:     ts->Brhs = Pmat;
1277:   }
1278:   return(0);
1279: }

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

1284:    Logically Collective on TS

1286:    Input Parameters:
1287: +  ts  - the TS context obtained from TSCreate()
1288: .  r   - vector to hold the residual (or NULL to have it created internally)
1289: .  f   - the function evaluation routine
1290: -  ctx - user-defined context for private data for the function evaluation routine (may be NULL)

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

1295: +  t   - time at step/stage being solved
1296: .  u   - state vector
1297: .  u_t - time derivative of state vector
1298: .  F   - function vector
1299: -  ctx - [optional] user-defined context for matrix evaluation routine

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

1304:    Level: beginner

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

1308: .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian()
1309: @*/
1310: PetscErrorCode  TSSetIFunction(TS ts,Vec r,TSIFunction f,void *ctx)
1311: {
1313:   SNES           snes;
1314:   Vec            ralloc = NULL;
1315:   DM             dm;


1321:   TSGetDM(ts,&dm);
1322:   DMTSSetIFunction(dm,f,ctx);

1324:   TSGetSNES(ts,&snes);
1325:   if (!r && !ts->dm && ts->vec_sol) {
1326:     VecDuplicate(ts->vec_sol,&ralloc);
1327:     r  = ralloc;
1328:   }
1329:   SNESSetFunction(snes,r,SNESTSFormFunction,ts);
1330:   VecDestroy(&ralloc);
1331:   return(0);
1332: }

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

1337:    Not Collective

1339:    Input Parameter:
1340: .  ts - the TS context

1342:    Output Parameter:
1343: +  r - vector to hold residual (or NULL)
1344: .  func - the function to compute residual (or NULL)
1345: -  ctx - the function context (or NULL)

1347:    Level: advanced

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

1351: .seealso: TSSetIFunction(), SNESGetFunction()
1352: @*/
1353: PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx)
1354: {
1356:   SNES           snes;
1357:   DM             dm;

1361:   TSGetSNES(ts,&snes);
1362:   SNESGetFunction(snes,r,NULL,NULL);
1363:   TSGetDM(ts,&dm);
1364:   DMTSGetIFunction(dm,func,ctx);
1365:   return(0);
1366: }

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

1371:    Not Collective

1373:    Input Parameter:
1374: .  ts - the TS context

1376:    Output Parameter:
1377: +  r - vector to hold computed right hand side (or NULL)
1378: .  func - the function to compute right hand side (or NULL)
1379: -  ctx - the function context (or NULL)

1381:    Level: advanced

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

1385: .seealso: TSSetRHSFunction(), SNESGetFunction()
1386: @*/
1387: PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx)
1388: {
1390:   SNES           snes;
1391:   DM             dm;

1395:   TSGetSNES(ts,&snes);
1396:   SNESGetFunction(snes,r,NULL,NULL);
1397:   TSGetDM(ts,&dm);
1398:   DMTSGetRHSFunction(dm,func,ctx);
1399:   return(0);
1400: }

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

1406:    Logically Collective on TS

1408:    Input Parameters:
1409: +  ts  - the TS context obtained from TSCreate()
1410: .  Amat - (approximate) Jacobian matrix
1411: .  Pmat - matrix used to compute preconditioner (usually the same as Amat)
1412: .  f   - the Jacobian evaluation routine
1413: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)

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

1418: +  t    - time at step/stage being solved
1419: .  U    - state vector
1420: .  U_t  - time derivative of state vector
1421: .  a    - shift
1422: .  Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1423: .  Pmat - matrix used for constructing preconditioner, usually the same as Amat
1424: -  ctx  - [optional] user-defined context for matrix evaluation routine

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

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

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

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

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

1444:    Level: beginner

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

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

1450: @*/
1451: PetscErrorCode  TSSetIJacobian(TS ts,Mat Amat,Mat Pmat,TSIJacobian f,void *ctx)
1452: {
1454:   SNES           snes;
1455:   DM             dm;


1464:   TSGetDM(ts,&dm);
1465:   DMTSSetIJacobian(dm,f,ctx);

1467:   TSGetSNES(ts,&snes);
1468:   SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);
1469:   return(0);
1470: }

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

1478:    Logically Collective

1480:    Input Arguments:
1481: +  ts - TS context obtained from TSCreate()
1482: -  reuse - PETSC_TRUE if the RHS Jacobian

1484:    Level: intermediate

1486: .seealso: TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
1487: @*/
1488: PetscErrorCode TSRHSJacobianSetReuse(TS ts,PetscBool reuse)
1489: {
1491:   ts->rhsjacobian.reuse = reuse;
1492:   return(0);
1493: }

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

1498:    Logically Collective on TS

1500:    Input Parameters:
1501: +  ts  - the TS context obtained from TSCreate()
1502: .  F   - vector to hold the residual (or NULL to have it created internally)
1503: .  fun - the function evaluation routine
1504: -  ctx - user-defined context for private data for the function evaluation routine (may be NULL)

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

1509: +  t    - time at step/stage being solved
1510: .  U    - state vector
1511: .  U_t  - time derivative of state vector
1512: .  U_tt - second time derivative of state vector
1513: .  F    - function vector
1514: -  ctx  - [optional] user-defined context for matrix evaluation routine (may be NULL)

1516:    Level: beginner

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

1520: .seealso: TSSetI2Jacobian()
1521: @*/
1522: PetscErrorCode TSSetI2Function(TS ts,Vec F,TSI2Function fun,void *ctx)
1523: {
1524:   DM             dm;

1530:   TSSetIFunction(ts,F,NULL,NULL);
1531:   TSGetDM(ts,&dm);
1532:   DMTSSetI2Function(dm,fun,ctx);
1533:   return(0);
1534: }

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

1539:   Not Collective

1541:   Input Parameter:
1542: . ts - the TS context

1544:   Output Parameter:
1545: + r - vector to hold residual (or NULL)
1546: . fun - the function to compute residual (or NULL)
1547: - ctx - the function context (or NULL)

1549:   Level: advanced

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

1553: .seealso: TSSetI2Function(), SNESGetFunction()
1554: @*/
1555: PetscErrorCode TSGetI2Function(TS ts,Vec *r,TSI2Function *fun,void **ctx)
1556: {
1558:   SNES           snes;
1559:   DM             dm;

1563:   TSGetSNES(ts,&snes);
1564:   SNESGetFunction(snes,r,NULL,NULL);
1565:   TSGetDM(ts,&dm);
1566:   DMTSGetI2Function(dm,fun,ctx);
1567:   return(0);
1568: }

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

1574:    Logically Collective on TS

1576:    Input Parameters:
1577: +  ts  - the TS context obtained from TSCreate()
1578: .  J   - Jacobian matrix
1579: .  P   - preconditioning matrix for J (may be same as J)
1580: .  jac - the Jacobian evaluation routine
1581: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)

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

1586: +  t    - time at step/stage being solved
1587: .  U    - state vector
1588: .  U_t  - time derivative of state vector
1589: .  U_tt - second time derivative of state vector
1590: .  v    - shift for U_t
1591: .  a    - shift for U_tt
1592: .  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
1593: .  P    - preconditioning matrix for J, may be same as J
1594: -  ctx  - [optional] user-defined context for matrix evaluation routine

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

1599:    The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1600:    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.
1601:    The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U  where the positive "shift"
1602:    parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.

1604:    Level: beginner

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

1608: .seealso: TSSetI2Function()
1609: @*/
1610: PetscErrorCode TSSetI2Jacobian(TS ts,Mat J,Mat P,TSI2Jacobian jac,void *ctx)
1611: {
1612:   DM             dm;

1619:   TSSetIJacobian(ts,J,P,NULL,NULL);
1620:   TSGetDM(ts,&dm);
1621:   DMTSSetI2Jacobian(dm,jac,ctx);
1622:   return(0);
1623: }

1625: /*@C
1626:   TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.

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

1630:   Input Parameter:
1631: . ts  - The TS context obtained from TSCreate()

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

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

1642:   Level: advanced

1644: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetStepNumber()

1646: .keywords: TS, timestep, get, matrix, Jacobian
1647: @*/
1648: PetscErrorCode  TSGetI2Jacobian(TS ts,Mat *J,Mat *P,TSI2Jacobian *jac,void **ctx)
1649: {
1651:   SNES           snes;
1652:   DM             dm;

1655:   TSGetSNES(ts,&snes);
1656:   SNESSetUpMatrices(snes);
1657:   SNESGetJacobian(snes,J,P,NULL,NULL);
1658:   TSGetDM(ts,&dm);
1659:   DMTSGetI2Jacobian(dm,jac,ctx);
1660:   return(0);
1661: }

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

1666:   Collective on TS and Vec

1668:   Input Parameters:
1669: + ts - the TS context
1670: . t - current time
1671: . U - state vector
1672: . V - time derivative of state vector (U_t)
1673: - A - second time derivative of state vector (U_tt)

1675:   Output Parameter:
1676: . F - the residual vector

1678:   Note:
1679:   Most users should not need to explicitly call this routine, as it
1680:   is used internally within the nonlinear solvers.

1682:   Level: developer

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

1686: .seealso: TSSetI2Function()
1687: @*/
1688: PetscErrorCode TSComputeI2Function(TS ts,PetscReal t,Vec U,Vec V,Vec A,Vec F)
1689: {
1690:   DM             dm;
1691:   TSI2Function   I2Function;
1692:   void           *ctx;
1693:   TSRHSFunction  rhsfunction;


1703:   TSGetDM(ts,&dm);
1704:   DMTSGetI2Function(dm,&I2Function,&ctx);
1705:   DMTSGetRHSFunction(dm,&rhsfunction,NULL);

1707:   if (!I2Function) {
1708:     TSComputeIFunction(ts,t,U,A,F,PETSC_FALSE);
1709:     return(0);
1710:   }

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

1714:   PetscStackPush("TS user implicit function");
1715:   I2Function(ts,t,U,V,A,F,ctx);
1716:   PetscStackPop;

1718:   if (rhsfunction) {
1719:     Vec Frhs;
1720:     TSGetRHSVec_Private(ts,&Frhs);
1721:     TSComputeRHSFunction(ts,t,U,Frhs);
1722:     VecAXPY(F,-1,Frhs);
1723:   }

1725:   PetscLogEventEnd(TS_FunctionEval,ts,U,V,F);
1726:   return(0);
1727: }

1729: /*@
1730:   TSComputeI2Jacobian - Evaluates the Jacobian of the DAE

1732:   Collective on TS and Vec

1734:   Input Parameters:
1735: + ts - the TS context
1736: . t - current timestep
1737: . U - state vector
1738: . V - time derivative of state vector
1739: . A - second time derivative of state vector
1740: . shiftV - shift to apply, see note below
1741: - shiftA - shift to apply, see note below

1743:   Output Parameters:
1744: + J - Jacobian matrix
1745: - P - optional preconditioning matrix

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

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

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

1755:   Level: developer

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

1759: .seealso:  TSSetI2Jacobian()
1760: @*/
1761: PetscErrorCode TSComputeI2Jacobian(TS ts,PetscReal t,Vec U,Vec V,Vec A,PetscReal shiftV,PetscReal shiftA,Mat J,Mat P)
1762: {
1763:   DM             dm;
1764:   TSI2Jacobian   I2Jacobian;
1765:   void           *ctx;
1766:   TSRHSJacobian  rhsjacobian;


1777:   TSGetDM(ts,&dm);
1778:   DMTSGetI2Jacobian(dm,&I2Jacobian,&ctx);
1779:   DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);

1781:   if (!I2Jacobian) {
1782:     TSComputeIJacobian(ts,t,U,A,shiftA,J,P,PETSC_FALSE);
1783:     return(0);
1784:   }

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

1788:   PetscStackPush("TS user implicit Jacobian");
1789:   I2Jacobian(ts,t,U,V,A,shiftV,shiftA,J,P,ctx);
1790:   PetscStackPop;

1792:   if (rhsjacobian) {
1793:     Mat Jrhs,Prhs; MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
1794:     TSGetRHSMats_Private(ts,&Jrhs,&Prhs);
1795:     TSComputeRHSJacobian(ts,t,U,Jrhs,Prhs);
1796:     MatAXPY(J,-1,Jrhs,axpy);
1797:     if (P != J) {MatAXPY(P,-1,Prhs,axpy);}
1798:   }

1800:   PetscLogEventEnd(TS_JacobianEval,ts,U,J,P);
1801:   return(0);
1802: }

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

1808:    Logically Collective on TS and Vec

1810:    Input Parameters:
1811: +  ts - the TS context obtained from TSCreate()
1812: .  u - the solution vector
1813: -  v - the time derivative vector

1815:    Level: beginner

1817: .keywords: TS, timestep, set, solution, initial conditions
1818: @*/
1819: PetscErrorCode  TS2SetSolution(TS ts,Vec u,Vec v)
1820: {

1827:   TSSetSolution(ts,u);
1828:   PetscObjectReference((PetscObject)v);
1829:   VecDestroy(&ts->vec_dot);
1830:   ts->vec_dot = v;
1831:   return(0);
1832: }

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

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

1842:    Input Parameter:
1843: .  ts - the TS context obtained from TSCreate()

1845:    Output Parameter:
1846: +  u - the vector containing the solution
1847: -  v - the vector containing the time derivative

1849:    Level: intermediate

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

1853: .keywords: TS, timestep, get, solution
1854: @*/
1855: PetscErrorCode  TS2GetSolution(TS ts,Vec *u,Vec *v)
1856: {
1861:   if (u) *u = ts->vec_sol;
1862:   if (v) *v = ts->vec_dot;
1863:   return(0);
1864: }

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

1869:   Collective on PetscViewer

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

1876:    Level: intermediate

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

1881:   Notes for advanced users:
1882:   Most users should not need to know the details of the binary storage
1883:   format, since TSLoad() and TSView() completely hide these details.
1884:   But for anyone who's interested, the standard binary matrix storage
1885:   format is
1886: .vb
1887:      has not yet been determined
1888: .ve

1890: .seealso: PetscViewerBinaryOpen(), TSView(), MatLoad(), VecLoad()
1891: @*/
1892: PetscErrorCode  TSLoad(TS ts, PetscViewer viewer)
1893: {
1895:   PetscBool      isbinary;
1896:   PetscInt       classid;
1897:   char           type[256];
1898:   DMTS           sdm;
1899:   DM             dm;

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

1907:   PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1908:   if (classid != TS_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Not TS next in file");
1909:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1910:   TSSetType(ts, type);
1911:   if (ts->ops->load) {
1912:     (*ts->ops->load)(ts,viewer);
1913:   }
1914:   DMCreate(PetscObjectComm((PetscObject)ts),&dm);
1915:   DMLoad(dm,viewer);
1916:   TSSetDM(ts,dm);
1917:   DMCreateGlobalVector(ts->dm,&ts->vec_sol);
1918:   VecLoad(ts->vec_sol,viewer);
1919:   DMGetDMTS(ts->dm,&sdm);
1920:   DMTSLoad(sdm,viewer);
1921:   return(0);
1922: }

1924:  #include <petscdraw.h>
1925: #if defined(PETSC_HAVE_SAWS)
1926:  #include <petscviewersaws.h>
1927: #endif
1928: /*@C
1929:     TSView - Prints the TS data structure.

1931:     Collective on TS

1933:     Input Parameters:
1934: +   ts - the TS context obtained from TSCreate()
1935: -   viewer - visualization context

1937:     Options Database Key:
1938: .   -ts_view - calls TSView() at end of TSStep()

1940:     Notes:
1941:     The available visualization contexts include
1942: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1943: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1944:          output where only the first processor opens
1945:          the file.  All other processors send their
1946:          data to the first processor to print.

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

1951:     Level: beginner

1953: .keywords: TS, timestep, view

1955: .seealso: PetscViewerASCIIOpen()
1956: @*/
1957: PetscErrorCode  TSView(TS ts,PetscViewer viewer)
1958: {
1960:   TSType         type;
1961:   PetscBool      iascii,isstring,isundials,isbinary,isdraw;
1962:   DMTS           sdm;
1963: #if defined(PETSC_HAVE_SAWS)
1964:   PetscBool      issaws;
1965: #endif

1969:   if (!viewer) {
1970:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts),&viewer);
1971:   }

1975:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1976:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1977:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1978:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1979: #if defined(PETSC_HAVE_SAWS)
1980:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1981: #endif
1982:   if (iascii) {
1983:     PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer);
1984:     if (ts->ops->view) {
1985:       PetscViewerASCIIPushTab(viewer);
1986:       (*ts->ops->view)(ts,viewer);
1987:       PetscViewerASCIIPopTab(viewer);
1988:     }
1989:     if (ts->max_steps < PETSC_MAX_INT) {
1990:       PetscViewerASCIIPrintf(viewer,"  maximum steps=%D\n",ts->max_steps);
1991:     }
1992:     if (ts->max_time < PETSC_MAX_REAL) {
1993:       PetscViewerASCIIPrintf(viewer,"  maximum time=%g\n",(double)ts->max_time);
1994:     }
1995:     if (ts->usessnes) {
1996:       PetscBool lin;
1997:       if (ts->problem_type == TS_NONLINEAR) {
1998:         PetscViewerASCIIPrintf(viewer,"  total number of nonlinear solver iterations=%D\n",ts->snes_its);
1999:       }
2000:       PetscViewerASCIIPrintf(viewer,"  total number of linear solver iterations=%D\n",ts->ksp_its);
2001:       PetscObjectTypeCompareAny((PetscObject)ts->snes,&lin,SNESKSPONLY,SNESKSPTRANSPOSEONLY,"");
2002:       PetscViewerASCIIPrintf(viewer,"  total number of %slinear solve failures=%D\n",lin ? "" : "non",ts->num_snes_failures);
2003:     }
2004:     PetscViewerASCIIPrintf(viewer,"  total number of rejected steps=%D\n",ts->reject);
2005:     if (ts->vrtol) {
2006:       PetscViewerASCIIPrintf(viewer,"  using vector of relative error tolerances, ");
2007:     } else {
2008:       PetscViewerASCIIPrintf(viewer,"  using relative error tolerance of %g, ",(double)ts->rtol);
2009:     }
2010:     if (ts->vatol) {
2011:       PetscViewerASCIIPrintf(viewer,"  using vector of absolute error tolerances\n");
2012:     } else {
2013:       PetscViewerASCIIPrintf(viewer,"  using absolute error tolerance of %g\n",(double)ts->atol);
2014:     }
2015:     PetscViewerASCIIPushTab(viewer);
2016:     TSAdaptView(ts->adapt,viewer);
2017:     PetscViewerASCIIPopTab(viewer);
2018:     if (ts->snes && ts->usessnes)  {
2019:       PetscViewerASCIIPushTab(viewer);
2020:       SNESView(ts->snes,viewer);
2021:       PetscViewerASCIIPopTab(viewer);
2022:     }
2023:     DMGetDMTS(ts->dm,&sdm);
2024:     DMTSView(sdm,viewer);
2025:   } else if (isstring) {
2026:     TSGetType(ts,&type);
2027:     PetscViewerStringSPrintf(viewer," %-7.7s",type);
2028:   } else if (isbinary) {
2029:     PetscInt    classid = TS_FILE_CLASSID;
2030:     MPI_Comm    comm;
2031:     PetscMPIInt rank;
2032:     char        type[256];

2034:     PetscObjectGetComm((PetscObject)ts,&comm);
2035:     MPI_Comm_rank(comm,&rank);
2036:     if (!rank) {
2037:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
2038:       PetscStrncpy(type,((PetscObject)ts)->type_name,256);
2039:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
2040:     }
2041:     if (ts->ops->view) {
2042:       (*ts->ops->view)(ts,viewer);
2043:     }
2044:     if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
2045:     DMView(ts->dm,viewer);
2046:     VecView(ts->vec_sol,viewer);
2047:     DMGetDMTS(ts->dm,&sdm);
2048:     DMTSView(sdm,viewer);
2049:   } else if (isdraw) {
2050:     PetscDraw draw;
2051:     char      str[36];
2052:     PetscReal x,y,bottom,h;

2054:     PetscViewerDrawGetDraw(viewer,0,&draw);
2055:     PetscDrawGetCurrentPoint(draw,&x,&y);
2056:     PetscStrcpy(str,"TS: ");
2057:     PetscStrcat(str,((PetscObject)ts)->type_name);
2058:     PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,str,NULL,&h);
2059:     bottom = y - h;
2060:     PetscDrawPushCurrentPoint(draw,x,bottom);
2061:     if (ts->ops->view) {
2062:       (*ts->ops->view)(ts,viewer);
2063:     }
2064:     if (ts->adapt) {TSAdaptView(ts->adapt,viewer);}
2065:     if (ts->snes)  {SNESView(ts->snes,viewer);}
2066:     PetscDrawPopCurrentPoint(draw);
2067: #if defined(PETSC_HAVE_SAWS)
2068:   } else if (issaws) {
2069:     PetscMPIInt rank;
2070:     const char  *name;

2072:     PetscObjectGetName((PetscObject)ts,&name);
2073:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
2074:     if (!((PetscObject)ts)->amsmem && !rank) {
2075:       char       dir[1024];

2077:       PetscObjectViewSAWs((PetscObject)ts,viewer);
2078:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time_step",name);
2079:       PetscStackCallSAWs(SAWs_Register,(dir,&ts->steps,1,SAWs_READ,SAWs_INT));
2080:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time",name);
2081:       PetscStackCallSAWs(SAWs_Register,(dir,&ts->ptime,1,SAWs_READ,SAWs_DOUBLE));
2082:     }
2083:     if (ts->ops->view) {
2084:       (*ts->ops->view)(ts,viewer);
2085:     }
2086: #endif
2087:   }

2089:   PetscViewerASCIIPushTab(viewer);
2090:   PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);
2091:   PetscViewerASCIIPopTab(viewer);
2092:   return(0);
2093: }

2095: /*@
2096:    TSSetApplicationContext - Sets an optional user-defined context for
2097:    the timesteppers.

2099:    Logically Collective on TS

2101:    Input Parameters:
2102: +  ts - the TS context obtained from TSCreate()
2103: -  usrP - optional user context

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

2109:    Level: intermediate

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

2113: .seealso: TSGetApplicationContext()
2114: @*/
2115: PetscErrorCode  TSSetApplicationContext(TS ts,void *usrP)
2116: {
2119:   ts->user = usrP;
2120:   return(0);
2121: }

2123: /*@
2124:     TSGetApplicationContext - Gets the user-defined context for the
2125:     timestepper.

2127:     Not Collective

2129:     Input Parameter:
2130: .   ts - the TS context obtained from TSCreate()

2132:     Output Parameter:
2133: .   usrP - user context

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

2139:     Level: intermediate

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

2143: .seealso: TSSetApplicationContext()
2144: @*/
2145: PetscErrorCode  TSGetApplicationContext(TS ts,void *usrP)
2146: {
2149:   *(void**)usrP = ts->user;
2150:   return(0);
2151: }

2153: /*@
2154:    TSGetStepNumber - Gets the number of steps completed.

2156:    Not Collective

2158:    Input Parameter:
2159: .  ts - the TS context obtained from TSCreate()

2161:    Output Parameter:
2162: .  steps - number of steps completed so far

2164:    Level: intermediate

2166: .keywords: TS, timestep, get, iteration, number
2167: .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSSetPostStep()
2168: @*/
2169: PetscErrorCode TSGetStepNumber(TS ts,PetscInt *steps)
2170: {
2174:   *steps = ts->steps;
2175:   return(0);
2176: }

2178: /*@
2179:    TSSetStepNumber - Sets the number of steps completed.

2181:    Logically Collective on TS

2183:    Input Parameters:
2184: +  ts - the TS context
2185: -  steps - number of steps completed so far

2187:    Notes:
2188:    For most uses of the TS solvers the user need not explicitly call
2189:    TSSetStepNumber(), as the step counter is appropriately updated in
2190:    TSSolve()/TSStep()/TSRollBack(). Power users may call this routine to
2191:    reinitialize timestepping by setting the step counter to zero (and time
2192:    to the initial time) to solve a similar problem with different initial
2193:    conditions or parameters. Other possible use case is to continue
2194:    timestepping from a previously interrupted run in such a way that TS
2195:    monitors will be called with a initial nonzero step counter.

2197:    Level: advanced

2199: .keywords: TS, timestep, set, iteration, number
2200: .seealso: TSGetStepNumber(), TSSetTime(), TSSetTimeStep(), TSSetSolution()
2201: @*/
2202: PetscErrorCode TSSetStepNumber(TS ts,PetscInt steps)
2203: {
2207:   if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Step number must be non-negative");
2208:   ts->steps = steps;
2209:   return(0);
2210: }

2212: /*@
2213:    TSSetTimeStep - Allows one to reset the timestep at any time,
2214:    useful for simple pseudo-timestepping codes.

2216:    Logically Collective on TS

2218:    Input Parameters:
2219: +  ts - the TS context obtained from TSCreate()
2220: -  time_step - the size of the timestep

2222:    Level: intermediate

2224: .seealso: TSGetTimeStep(), TSSetTime()

2226: .keywords: TS, set, timestep
2227: @*/
2228: PetscErrorCode  TSSetTimeStep(TS ts,PetscReal time_step)
2229: {
2233:   ts->time_step = time_step;
2234:   return(0);
2235: }

2237: /*@
2238:    TSSetExactFinalTime - Determines whether to adapt the final time step to
2239:      match the exact final time, interpolate solution to the exact final time,
2240:      or just return at the final time TS computed.

2242:   Logically Collective on TS

2244:    Input Parameter:
2245: +   ts - the time-step context
2246: -   eftopt - exact final time option

2248: $  TS_EXACTFINALTIME_STEPOVER    - Don't do anything if final time is exceeded
2249: $  TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2250: $  TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time

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

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

2258:    Level: beginner

2260: .seealso: TSExactFinalTimeOption, TSGetExactFinalTime()
2261: @*/
2262: PetscErrorCode TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt)
2263: {
2267:   ts->exact_final_time = eftopt;
2268:   return(0);
2269: }

2271: /*@
2272:    TSGetExactFinalTime - Gets the exact final time option.

2274:    Not Collective

2276:    Input Parameter:
2277: .  ts - the TS context

2279:    Output Parameter:
2280: .  eftopt - exact final time option

2282:    Level: beginner

2284: .seealso: TSExactFinalTimeOption, TSSetExactFinalTime()
2285: @*/
2286: PetscErrorCode TSGetExactFinalTime(TS ts,TSExactFinalTimeOption *eftopt)
2287: {
2291:   *eftopt = ts->exact_final_time;
2292:   return(0);
2293: }

2295: /*@
2296:    TSGetTimeStep - Gets the current timestep size.

2298:    Not Collective

2300:    Input Parameter:
2301: .  ts - the TS context obtained from TSCreate()

2303:    Output Parameter:
2304: .  dt - the current timestep size

2306:    Level: intermediate

2308: .seealso: TSSetTimeStep(), TSGetTime()

2310: .keywords: TS, get, timestep
2311: @*/
2312: PetscErrorCode  TSGetTimeStep(TS ts,PetscReal *dt)
2313: {
2317:   *dt = ts->time_step;
2318:   return(0);
2319: }

2321: /*@
2322:    TSGetSolution - Returns the solution at the present timestep. It
2323:    is valid to call this routine inside the function that you are evaluating
2324:    in order to move to the new timestep. This vector not changed until
2325:    the solution at the next timestep has been calculated.

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

2329:    Input Parameter:
2330: .  ts - the TS context obtained from TSCreate()

2332:    Output Parameter:
2333: .  v - the vector containing the solution

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

2338:    Level: intermediate

2340: .seealso: TSGetTimeStep(), TSGetTime(), TSGetSolveTime(), TSGetSolutionComponents(), TSSetSolutionFunction()

2342: .keywords: TS, timestep, get, solution
2343: @*/
2344: PetscErrorCode  TSGetSolution(TS ts,Vec *v)
2345: {
2349:   *v = ts->vec_sol;
2350:   return(0);
2351: }

2353: /*@
2354:    TSGetSolutionComponents - Returns any solution components at the present
2355:    timestep, if available for the time integration method being used.
2356:    Solution components are quantities that share the same size and
2357:    structure as the solution vector.

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

2361:    Parameters :
2362: .  ts - the TS context obtained from TSCreate() (input parameter).
2363: .  n - If v is PETSC_NULL, then the number of solution components is
2364:        returned through n, else the n-th solution component is
2365:        returned in v.
2366: .  v - the vector containing the n-th solution component
2367:        (may be PETSC_NULL to use this function to find out
2368:         the number of solutions components).

2370:    Level: advanced

2372: .seealso: TSGetSolution()

2374: .keywords: TS, timestep, get, solution
2375: @*/
2376: PetscErrorCode  TSGetSolutionComponents(TS ts,PetscInt *n,Vec *v)
2377: {

2382:   if (!ts->ops->getsolutioncomponents) *n = 0;
2383:   else {
2384:     (*ts->ops->getsolutioncomponents)(ts,n,v);
2385:   }
2386:   return(0);
2387: }

2389: /*@
2390:    TSGetAuxSolution - Returns an auxiliary solution at the present
2391:    timestep, if available for the time integration method being used.

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

2395:    Parameters :
2396: .  ts - the TS context obtained from TSCreate() (input parameter).
2397: .  v - the vector containing the auxiliary solution

2399:    Level: intermediate

2401: .seealso: TSGetSolution()

2403: .keywords: TS, timestep, get, solution
2404: @*/
2405: PetscErrorCode  TSGetAuxSolution(TS ts,Vec *v)
2406: {

2411:   if (ts->ops->getauxsolution) {
2412:     (*ts->ops->getauxsolution)(ts,v);
2413:   } else {
2414:     VecZeroEntries(*v);
2415:   }
2416:   return(0);
2417: }

2419: /*@
2420:    TSGetTimeError - Returns the estimated error vector, if the chosen
2421:    TSType has an error estimation functionality.

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

2425:    Note: MUST call after TSSetUp()

2427:    Parameters :
2428: .  ts - the TS context obtained from TSCreate() (input parameter).
2429: .  n - current estimate (n=0) or previous one (n=-1)
2430: .  v - the vector containing the error (same size as the solution).

2432:    Level: intermediate

2434: .seealso: TSGetSolution(), TSSetTimeError()

2436: .keywords: TS, timestep, get, error
2437: @*/
2438: PetscErrorCode  TSGetTimeError(TS ts,PetscInt n,Vec *v)
2439: {

2444:   if (ts->ops->gettimeerror) {
2445:     (*ts->ops->gettimeerror)(ts,n,v);
2446:   } else {
2447:     VecZeroEntries(*v);
2448:   }
2449:   return(0);
2450: }

2452: /*@
2453:    TSSetTimeError - Sets the estimated error vector, if the chosen
2454:    TSType has an error estimation functionality. This can be used
2455:    to restart such a time integrator with a given error vector.

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

2459:    Parameters :
2460: .  ts - the TS context obtained from TSCreate() (input parameter).
2461: .  v - the vector containing the error (same size as the solution).

2463:    Level: intermediate

2465: .seealso: TSSetSolution(), TSGetTimeError)

2467: .keywords: TS, timestep, get, error
2468: @*/
2469: PetscErrorCode  TSSetTimeError(TS ts,Vec v)
2470: {

2475:   if (!ts->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetUp() first");
2476:   if (ts->ops->settimeerror) {
2477:     (*ts->ops->settimeerror)(ts,v);
2478:   }
2479:   return(0);
2480: }

2482: /* ----- Routines to initialize and destroy a timestepper ---- */
2483: /*@
2484:   TSSetProblemType - Sets the type of problem to be solved.

2486:   Not collective

2488:   Input Parameters:
2489: + ts   - The TS
2490: - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2491: .vb
2492:          U_t - A U = 0      (linear)
2493:          U_t - A(t) U = 0   (linear)
2494:          F(t,U,U_t) = 0     (nonlinear)
2495: .ve

2497:    Level: beginner

2499: .keywords: TS, problem type
2500: .seealso: TSSetUp(), TSProblemType, TS
2501: @*/
2502: PetscErrorCode  TSSetProblemType(TS ts, TSProblemType type)
2503: {

2508:   ts->problem_type = type;
2509:   if (type == TS_LINEAR) {
2510:     SNES snes;
2511:     TSGetSNES(ts,&snes);
2512:     SNESSetType(snes,SNESKSPONLY);
2513:   }
2514:   return(0);
2515: }

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

2520:   Not collective

2522:   Input Parameter:
2523: . ts   - The TS

2525:   Output Parameter:
2526: . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2527: .vb
2528:          M U_t = A U
2529:          M(t) U_t = A(t) U
2530:          F(t,U,U_t)
2531: .ve

2533:    Level: beginner

2535: .keywords: TS, problem type
2536: .seealso: TSSetUp(), TSProblemType, TS
2537: @*/
2538: PetscErrorCode  TSGetProblemType(TS ts, TSProblemType *type)
2539: {
2543:   *type = ts->problem_type;
2544:   return(0);
2545: }

2547: /*@
2548:    TSSetUp - Sets up the internal data structures for the later use
2549:    of a timestepper.

2551:    Collective on TS

2553:    Input Parameter:
2554: .  ts - the TS context obtained from TSCreate()

2556:    Notes:
2557:    For basic use of the TS solvers the user need not explicitly call
2558:    TSSetUp(), since these actions will automatically occur during
2559:    the call to TSStep() or TSSolve().  However, if one wishes to control this
2560:    phase separately, TSSetUp() should be called after TSCreate()
2561:    and optional routines of the form TSSetXXX(), but before TSStep() and TSSolve().

2563:    Level: advanced

2565: .keywords: TS, timestep, setup

2567: .seealso: TSCreate(), TSStep(), TSDestroy(), TSSolve()
2568: @*/
2569: PetscErrorCode  TSSetUp(TS ts)
2570: {
2572:   DM             dm;
2573:   PetscErrorCode (*func)(SNES,Vec,Vec,void*);
2574:   PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*);
2575:   TSIFunction    ifun;
2576:   TSIJacobian    ijac;
2577:   TSI2Jacobian   i2jac;
2578:   TSRHSJacobian  rhsjac;
2579:   PetscBool      isnone;

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

2585:   if (!((PetscObject)ts)->type_name) {
2586:     TSGetIFunction(ts,NULL,&ifun,NULL);
2587:     TSSetType(ts,ifun ? TSBEULER : TSEULER);
2588:   }

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

2592:   TSGetRHSJacobian(ts,NULL,NULL,&rhsjac,NULL);
2593:   if (ts->rhsjacobian.reuse && rhsjac == TSComputeRHSJacobianConstant) {
2594:     Mat Amat,Pmat;
2595:     SNES snes;
2596:     TSGetSNES(ts,&snes);
2597:     SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);
2598:     /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2599:      * have displaced the RHS matrix */
2600:     if (Amat && Amat == ts->Arhs) {
2601:       /* we need to copy the values of the matrix because for the constant Jacobian case the user will never set the numerical values in this new location */
2602:       MatDuplicate(ts->Arhs,MAT_COPY_VALUES,&Amat);
2603:       SNESSetJacobian(snes,Amat,NULL,NULL,NULL);
2604:       MatDestroy(&Amat);
2605:     }
2606:     if (Pmat && Pmat == ts->Brhs) {
2607:       MatDuplicate(ts->Brhs,MAT_COPY_VALUES,&Pmat);
2608:       SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);
2609:       MatDestroy(&Pmat);
2610:     }
2611:   }

2613:   TSGetAdapt(ts,&ts->adapt);
2614:   TSAdaptSetDefaultType(ts->adapt,ts->default_adapt_type);

2616:   if (ts->ops->setup) {
2617:     (*ts->ops->setup)(ts);
2618:   }

2620:   /* Attempt to check/preset a default value for the exact final time option */
2621:   PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&isnone);
2622:   if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED)
2623:     ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;

2625:   /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2626:      to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2627:    */
2628:   TSGetDM(ts,&dm);
2629:   DMSNESGetFunction(dm,&func,NULL);
2630:   if (!func) {
2631:     DMSNESSetFunction(dm,SNESTSFormFunction,ts);
2632:   }
2633:   /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2634:      Otherwise, the SNES will use coloring internally to form the Jacobian.
2635:    */
2636:   DMSNESGetJacobian(dm,&jac,NULL);
2637:   DMTSGetIJacobian(dm,&ijac,NULL);
2638:   DMTSGetI2Jacobian(dm,&i2jac,NULL);
2639:   DMTSGetRHSJacobian(dm,&rhsjac,NULL);
2640:   if (!jac && (ijac || i2jac || rhsjac)) {
2641:     DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);
2642:   }

2644:   /* if time integration scheme has a starting method, call it */
2645:   if (ts->ops->startingmethod) {
2646:     (*ts->ops->startingmethod)(ts);
2647:   }

2649:   ts->setupcalled = PETSC_TRUE;
2650:   return(0);
2651: }

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

2656:    Collective on TS

2658:    Input Parameter:
2659: .  ts - the TS context obtained from TSCreate()

2661:    Level: beginner

2663: .keywords: TS, timestep, reset

2665: .seealso: TSCreate(), TSSetup(), TSDestroy()
2666: @*/
2667: PetscErrorCode  TSReset(TS ts)
2668: {
2669:   TS_RHSSplitLink ilink = ts->tsrhssplit,next;
2670:   PetscErrorCode  ierr;


2675:   if (ts->ops->reset) {
2676:     (*ts->ops->reset)(ts);
2677:   }
2678:   if (ts->snes) {SNESReset(ts->snes);}
2679:   if (ts->adapt) {TSAdaptReset(ts->adapt);}

2681:   MatDestroy(&ts->Arhs);
2682:   MatDestroy(&ts->Brhs);
2683:   VecDestroy(&ts->Frhs);
2684:   VecDestroy(&ts->vec_sol);
2685:   VecDestroy(&ts->vec_dot);
2686:   VecDestroy(&ts->vatol);
2687:   VecDestroy(&ts->vrtol);
2688:   VecDestroyVecs(ts->nwork,&ts->work);

2690:   VecDestroyVecs(ts->numcost,&ts->vecs_drdy);
2691:   VecDestroyVecs(ts->numcost,&ts->vecs_drdp);

2693:   MatDestroy(&ts->Jacp);
2694:   VecDestroy(&ts->vec_costintegral);
2695:   VecDestroy(&ts->vec_costintegrand);
2696:   MatDestroy(&ts->mat_sensip);

2698:   while (ilink) {
2699:     next = ilink->next;
2700:     TSDestroy(&ilink->ts);
2701:     PetscFree(ilink->splitname);
2702:     ISDestroy(&ilink->is);
2703:     PetscFree(ilink);
2704:     ilink = next;
2705:   }
2706:   ts->num_rhs_splits = 0;
2707:   ts->setupcalled = PETSC_FALSE;
2708:   return(0);
2709: }

2711: /*@
2712:    TSDestroy - Destroys the timestepper context that was created
2713:    with TSCreate().

2715:    Collective on TS

2717:    Input Parameter:
2718: .  ts - the TS context obtained from TSCreate()

2720:    Level: beginner

2722: .keywords: TS, timestepper, destroy

2724: .seealso: TSCreate(), TSSetUp(), TSSolve()
2725: @*/
2726: PetscErrorCode  TSDestroy(TS *ts)
2727: {

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

2735:   TSReset((*ts));

2737:   /* if memory was published with SAWs then destroy it */
2738:   PetscObjectSAWsViewOff((PetscObject)*ts);
2739:   if ((*ts)->ops->destroy) {(*(*ts)->ops->destroy)((*ts));}

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

2743:   TSAdaptDestroy(&(*ts)->adapt);
2744:   TSEventDestroy(&(*ts)->event);

2746:   SNESDestroy(&(*ts)->snes);
2747:   DMDestroy(&(*ts)->dm);
2748:   TSMonitorCancel((*ts));
2749:   TSAdjointMonitorCancel((*ts));

2751:   PetscHeaderDestroy(ts);
2752:   return(0);
2753: }

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

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

2761:    Input Parameter:
2762: .  ts - the TS context obtained from TSCreate()

2764:    Output Parameter:
2765: .  snes - the nonlinear solver context

2767:    Notes:
2768:    The user can then directly manipulate the SNES context to set various
2769:    options, etc.  Likewise, the user can then extract and manipulate the
2770:    KSP, KSP, and PC contexts as well.

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

2775:    Level: beginner

2777: .keywords: timestep, get, SNES
2778: @*/
2779: PetscErrorCode  TSGetSNES(TS ts,SNES *snes)
2780: {

2786:   if (!ts->snes) {
2787:     SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);
2788:     PetscObjectSetOptions((PetscObject)ts->snes,((PetscObject)ts)->options);
2789:     SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
2790:     PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->snes);
2791:     PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);
2792:     if (ts->dm) {SNESSetDM(ts->snes,ts->dm);}
2793:     if (ts->problem_type == TS_LINEAR) {
2794:       SNESSetType(ts->snes,SNESKSPONLY);
2795:     }
2796:   }
2797:   *snes = ts->snes;
2798:   return(0);
2799: }

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

2804:    Collective

2806:    Input Parameter:
2807: +  ts - the TS context obtained from TSCreate()
2808: -  snes - the nonlinear solver context

2810:    Notes:
2811:    Most users should have the TS created by calling TSGetSNES()

2813:    Level: developer

2815: .keywords: timestep, set, SNES
2816: @*/
2817: PetscErrorCode TSSetSNES(TS ts,SNES snes)
2818: {
2820:   PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*);

2825:   PetscObjectReference((PetscObject)snes);
2826:   SNESDestroy(&ts->snes);

2828:   ts->snes = snes;

2830:   SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);
2831:   SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);
2832:   if (func == SNESTSFormJacobian) {
2833:     SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);
2834:   }
2835:   return(0);
2836: }

2838: /*@
2839:    TSGetKSP - Returns the KSP (linear solver) associated with
2840:    a TS (timestepper) context.

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

2844:    Input Parameter:
2845: .  ts - the TS context obtained from TSCreate()

2847:    Output Parameter:
2848: .  ksp - the nonlinear solver context

2850:    Notes:
2851:    The user can then directly manipulate the KSP context to set various
2852:    options, etc.  Likewise, the user can then extract and manipulate the
2853:    KSP and PC contexts as well.

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

2858:    Level: beginner

2860: .keywords: timestep, get, KSP
2861: @*/
2862: PetscErrorCode  TSGetKSP(TS ts,KSP *ksp)
2863: {
2865:   SNES           snes;

2870:   if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first");
2871:   if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()");
2872:   TSGetSNES(ts,&snes);
2873:   SNESGetKSP(snes,ksp);
2874:   return(0);
2875: }

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

2879: /*@
2880:    TSSetMaxSteps - Sets the maximum number of steps to use.

2882:    Logically Collective on TS

2884:    Input Parameters:
2885: +  ts - the TS context obtained from TSCreate()
2886: -  maxsteps - maximum number of steps to use

2888:    Options Database Keys:
2889: .  -ts_max_steps <maxsteps> - Sets maxsteps

2891:    Notes:
2892:    The default maximum number of steps is 5000

2894:    Level: intermediate

2896: .keywords: TS, timestep, set, maximum, steps

2898: .seealso: TSGetMaxSteps(), TSSetMaxTime(), TSSetExactFinalTime()
2899: @*/
2900: PetscErrorCode TSSetMaxSteps(TS ts,PetscInt maxsteps)
2901: {
2905:   if (maxsteps < 0 ) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of steps must be non-negative");
2906:   ts->max_steps = maxsteps;
2907:   return(0);
2908: }

2910: /*@
2911:    TSGetMaxSteps - Gets the maximum number of steps to use.

2913:    Not Collective

2915:    Input Parameters:
2916: .  ts - the TS context obtained from TSCreate()

2918:    Output Parameter:
2919: .  maxsteps - maximum number of steps to use

2921:    Level: advanced

2923: .keywords: TS, timestep, get, maximum, steps

2925: .seealso: TSSetMaxSteps(), TSGetMaxTime(), TSSetMaxTime()
2926: @*/
2927: PetscErrorCode TSGetMaxSteps(TS ts,PetscInt *maxsteps)
2928: {
2932:   *maxsteps = ts->max_steps;
2933:   return(0);
2934: }

2936: /*@
2937:    TSSetMaxTime - Sets the maximum (or final) time for timestepping.

2939:    Logically Collective on TS

2941:    Input Parameters:
2942: +  ts - the TS context obtained from TSCreate()
2943: -  maxtime - final time to step to

2945:    Options Database Keys:
2946: .  -ts_max_time <maxtime> - Sets maxtime

2948:    Notes:
2949:    The default maximum time is 5.0

2951:    Level: intermediate

2953: .keywords: TS, timestep, set, maximum, time

2955: .seealso: TSGetMaxTime(), TSSetMaxSteps(), TSSetExactFinalTime()
2956: @*/
2957: PetscErrorCode TSSetMaxTime(TS ts,PetscReal maxtime)
2958: {
2962:   ts->max_time = maxtime;
2963:   return(0);
2964: }

2966: /*@
2967:    TSGetMaxTime - Gets the maximum (or final) time for timestepping.

2969:    Not Collective

2971:    Input Parameters:
2972: .  ts - the TS context obtained from TSCreate()

2974:    Output Parameter:
2975: .  maxtime - final time to step to

2977:    Level: advanced

2979: .keywords: TS, timestep, get, maximum, time

2981: .seealso: TSSetMaxTime(), TSGetMaxSteps(), TSSetMaxSteps()
2982: @*/
2983: PetscErrorCode TSGetMaxTime(TS ts,PetscReal *maxtime)
2984: {
2988:   *maxtime = ts->max_time;
2989:   return(0);
2990: }

2992: /*@
2993:    TSSetInitialTimeStep - Deprecated, use TSSetTime() and TSSetTimeStep().

2995:    Level: deprecated

2997: @*/
2998: PetscErrorCode  TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step)
2999: {
3003:   TSSetTime(ts,initial_time);
3004:   TSSetTimeStep(ts,time_step);
3005:   return(0);
3006: }

3008: /*@
3009:    TSGetDuration - Deprecated, use TSGetMaxSteps() and TSGetMaxTime().

3011:    Level: deprecated

3013: @*/
3014: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
3015: {
3018:   if (maxsteps) {
3020:     *maxsteps = ts->max_steps;
3021:   }
3022:   if (maxtime) {
3024:     *maxtime = ts->max_time;
3025:   }
3026:   return(0);
3027: }

3029: /*@
3030:    TSSetDuration - Deprecated, use TSSetMaxSteps() and TSSetMaxTime().

3032:    Level: deprecated

3034: @*/
3035: PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime)
3036: {
3041:   if (maxsteps >= 0) ts->max_steps = maxsteps;
3042:   if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime;
3043:   return(0);
3044: }

3046: /*@
3047:    TSGetTimeStepNumber - Deprecated, use TSGetStepNumber().

3049:    Level: deprecated

3051: @*/
3052: PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt *steps) { return TSGetStepNumber(ts,steps); }

3054: /*@
3055:    TSGetTotalSteps - Deprecated, use TSGetStepNumber().

3057:    Level: deprecated

3059: @*/
3060: PetscErrorCode TSGetTotalSteps(TS ts,PetscInt *steps) { return TSGetStepNumber(ts,steps); }

3062: /*@
3063:    TSSetSolution - Sets the initial solution vector
3064:    for use by the TS routines.

3066:    Logically Collective on TS and Vec

3068:    Input Parameters:
3069: +  ts - the TS context obtained from TSCreate()
3070: -  u - the solution vector

3072:    Level: beginner

3074: .keywords: TS, timestep, set, solution, initial values

3076: .seealso: TSSetSolutionFunction(), TSGetSolution(), TSCreate()
3077: @*/
3078: PetscErrorCode  TSSetSolution(TS ts,Vec u)
3079: {
3081:   DM             dm;

3086:   PetscObjectReference((PetscObject)u);
3087:   VecDestroy(&ts->vec_sol);
3088:   ts->vec_sol = u;

3090:   TSGetDM(ts,&dm);
3091:   DMShellSetGlobalVector(dm,u);
3092:   return(0);
3093: }

3095: /*@C
3096:   TSSetPreStep - Sets the general-purpose function
3097:   called once at the beginning of each time step.

3099:   Logically Collective on TS

3101:   Input Parameters:
3102: + ts   - The TS context obtained from TSCreate()
3103: - func - The function

3105:   Calling sequence of func:
3106: . func (TS ts);

3108:   Level: intermediate

3110: .keywords: TS, timestep
3111: .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep(), TSRestartStep()
3112: @*/
3113: PetscErrorCode  TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
3114: {
3117:   ts->prestep = func;
3118:   return(0);
3119: }

3121: /*@
3122:   TSPreStep - Runs the user-defined pre-step function.

3124:   Collective on TS

3126:   Input Parameters:
3127: . ts   - The TS context obtained from TSCreate()

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

3133:   Level: developer

3135: .keywords: TS, timestep
3136: .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep()
3137: @*/
3138: PetscErrorCode  TSPreStep(TS ts)
3139: {

3144:   if (ts->prestep) {
3145:     Vec              U;
3146:     PetscObjectState sprev,spost;

3148:     TSGetSolution(ts,&U);
3149:     PetscObjectStateGet((PetscObject)U,&sprev);
3150:     PetscStackCallStandard((*ts->prestep),(ts));
3151:     PetscObjectStateGet((PetscObject)U,&spost);
3152:     if (sprev != spost) {TSRestartStep(ts);}
3153:   }
3154:   return(0);
3155: }

3157: /*@C
3158:   TSSetPreStage - Sets the general-purpose function
3159:   called once at the beginning of each stage.

3161:   Logically Collective on TS

3163:   Input Parameters:
3164: + ts   - The TS context obtained from TSCreate()
3165: - func - The function

3167:   Calling sequence of func:
3168: . PetscErrorCode func(TS ts, PetscReal stagetime);

3170:   Level: intermediate

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

3177: .keywords: TS, timestep
3178: .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3179: @*/
3180: PetscErrorCode  TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal))
3181: {
3184:   ts->prestage = func;
3185:   return(0);
3186: }

3188: /*@C
3189:   TSSetPostStage - Sets the general-purpose function
3190:   called once at the end of each stage.

3192:   Logically Collective on TS

3194:   Input Parameters:
3195: + ts   - The TS context obtained from TSCreate()
3196: - func - The function

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

3201:   Level: intermediate

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

3208: .keywords: TS, timestep
3209: .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3210: @*/
3211: PetscErrorCode  TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*))
3212: {
3215:   ts->poststage = func;
3216:   return(0);
3217: }

3219: /*@C
3220:   TSSetPostEvaluate - Sets the general-purpose function
3221:   called once at the end of each step evaluation.

3223:   Logically Collective on TS

3225:   Input Parameters:
3226: + ts   - The TS context obtained from TSCreate()
3227: - func - The function

3229:   Calling sequence of func:
3230: . PetscErrorCode func(TS ts);

3232:   Level: intermediate

3234:   Note:
3235:   Semantically, TSSetPostEvaluate() differs from TSSetPostStep() since the function it sets is called before event-handling
3236:   thus guaranteeing the same solution (computed by the time-stepper) will be passed to it. On the other hand, TSPostStep()
3237:   may be passed a different solution, possibly changed by the event handler. TSPostEvaluate() is called after the next step
3238:   solution is evaluated allowing to modify it, if need be. The solution can be obtained with TSGetSolution(), the time step
3239:   with TSGetTimeStep(), and the time at the start of the step is available via TSGetTime()

3241: .keywords: TS, timestep
3242: .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext()
3243: @*/
3244: PetscErrorCode  TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS))
3245: {
3248:   ts->postevaluate = func;
3249:   return(0);
3250: }

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

3255:   Collective on TS

3257:   Input Parameters:
3258: . ts          - The TS context obtained from TSCreate()
3259:   stagetime   - The absolute time of the current stage

3261:   Notes:
3262:   TSPreStage() is typically used within time stepping implementations,
3263:   most users would not generally call this routine themselves.

3265:   Level: developer

3267: .keywords: TS, timestep
3268: .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep()
3269: @*/
3270: PetscErrorCode  TSPreStage(TS ts, PetscReal stagetime)
3271: {
3274:   if (ts->prestage) {
3275:     PetscStackCallStandard((*ts->prestage),(ts,stagetime));
3276:   }
3277:   return(0);
3278: }

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

3283:   Collective on TS

3285:   Input Parameters:
3286: . ts          - The TS context obtained from TSCreate()
3287:   stagetime   - The absolute time of the current stage
3288:   stageindex  - Stage number
3289:   Y           - Array of vectors (of size = total number
3290:                 of stages) with the stage solutions

3292:   Notes:
3293:   TSPostStage() is typically used within time stepping implementations,
3294:   most users would not generally call this routine themselves.

3296:   Level: developer

3298: .keywords: TS, timestep
3299: .seealso: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep()
3300: @*/
3301: PetscErrorCode  TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3302: {
3305:   if (ts->poststage) {
3306:     PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y));
3307:   }
3308:   return(0);
3309: }

3311: /*@
3312:   TSPostEvaluate - Runs the user-defined post-evaluate function set using TSSetPostEvaluate()

3314:   Collective on TS

3316:   Input Parameters:
3317: . ts          - The TS context obtained from TSCreate()

3319:   Notes:
3320:   TSPostEvaluate() is typically used within time stepping implementations,
3321:   most users would not generally call this routine themselves.

3323:   Level: developer

3325: .keywords: TS, timestep
3326: .seealso: TSSetPostEvaluate(), TSSetPreStep(), TSPreStep(), TSPostStep()
3327: @*/
3328: PetscErrorCode  TSPostEvaluate(TS ts)
3329: {

3334:   if (ts->postevaluate) {
3335:     Vec              U;
3336:     PetscObjectState sprev,spost;

3338:     TSGetSolution(ts,&U);
3339:     PetscObjectStateGet((PetscObject)U,&sprev);
3340:     PetscStackCallStandard((*ts->postevaluate),(ts));
3341:     PetscObjectStateGet((PetscObject)U,&spost);
3342:     if (sprev != spost) {TSRestartStep(ts);}
3343:   }
3344:   return(0);
3345: }

3347: /*@C
3348:   TSSetPostStep - Sets the general-purpose function
3349:   called once at the end of each time step.

3351:   Logically Collective on TS

3353:   Input Parameters:
3354: + ts   - The TS context obtained from TSCreate()
3355: - func - The function

3357:   Calling sequence of func:
3358: $ func (TS ts);

3360:   Notes:
3361:   The function set by TSSetPostStep() is called after each successful step. The solution vector X
3362:   obtained by TSGetSolution() may be different than that computed at the step end if the event handler
3363:   locates an event and TSPostEvent() modifies it. Use TSSetPostEvaluate() if an unmodified solution is needed instead.

3365:   Level: intermediate

3367: .keywords: TS, timestep
3368: .seealso: TSSetPreStep(), TSSetPreStage(), TSSetPostEvaluate(), TSGetTimeStep(), TSGetStepNumber(), TSGetTime(), TSRestartStep()
3369: @*/
3370: PetscErrorCode  TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
3371: {
3374:   ts->poststep = func;
3375:   return(0);
3376: }

3378: /*@
3379:   TSPostStep - Runs the user-defined post-step function.

3381:   Collective on TS

3383:   Input Parameters:
3384: . ts   - The TS context obtained from TSCreate()

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

3390:   Level: developer

3392: .keywords: TS, timestep
3393: @*/
3394: PetscErrorCode  TSPostStep(TS ts)
3395: {

3400:   if (ts->poststep) {
3401:     Vec              U;
3402:     PetscObjectState sprev,spost;

3404:     TSGetSolution(ts,&U);
3405:     PetscObjectStateGet((PetscObject)U,&sprev);
3406:     PetscStackCallStandard((*ts->poststep),(ts));
3407:     PetscObjectStateGet((PetscObject)U,&spost);
3408:     if (sprev != spost) {TSRestartStep(ts);}
3409:   }
3410:   return(0);
3411: }

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

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

3419:    Logically Collective on TS

3421:    Input Parameters:
3422: +  ts - the TS context obtained from TSCreate()
3423: .  monitor - monitoring routine
3424: .  mctx - [optional] user-defined context for private data for the
3425:              monitor routine (use NULL if no context is desired)
3426: -  monitordestroy - [optional] routine that frees monitor context
3427:           (may be NULL)

3429:    Calling sequence of monitor:
3430: $    PetscErrorCode monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)

3432: +    ts - the TS context
3433: .    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)
3434: .    time - current time
3435: .    u - current iterate
3436: -    mctx - [optional] monitoring context

3438:    Notes:
3439:    This routine adds an additional monitor to the list of monitors that
3440:    already has been loaded.

3442:    Fortran Notes:
3443:     Only a single monitor function can be set for each TS object

3445:    Level: intermediate

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

3449: .seealso: TSMonitorDefault(), TSMonitorCancel()
3450: @*/
3451: PetscErrorCode  TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**))
3452: {
3454:   PetscInt       i;
3455:   PetscBool      identical;

3459:   for (i=0; i<ts->numbermonitors;i++) {
3460:     PetscMonitorCompare((PetscErrorCode (*)(void))monitor,mctx,mdestroy,(PetscErrorCode (*)(void))ts->monitor[i],ts->monitorcontext[i],ts->monitordestroy[i],&identical);
3461:     if (identical) return(0);
3462:   }
3463:   if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
3464:   ts->monitor[ts->numbermonitors]          = monitor;
3465:   ts->monitordestroy[ts->numbermonitors]   = mdestroy;
3466:   ts->monitorcontext[ts->numbermonitors++] = (void*)mctx;
3467:   return(0);
3468: }

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

3473:    Logically Collective on TS

3475:    Input Parameters:
3476: .  ts - the TS context obtained from TSCreate()

3478:    Notes:
3479:    There is no way to remove a single, specific monitor.

3481:    Level: intermediate

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

3485: .seealso: TSMonitorDefault(), TSMonitorSet()
3486: @*/
3487: PetscErrorCode  TSMonitorCancel(TS ts)
3488: {
3490:   PetscInt       i;

3494:   for (i=0; i<ts->numbermonitors; i++) {
3495:     if (ts->monitordestroy[i]) {
3496:       (*ts->monitordestroy[i])(&ts->monitorcontext[i]);
3497:     }
3498:   }
3499:   ts->numbermonitors = 0;
3500:   return(0);
3501: }

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

3506:    Level: intermediate

3508: .keywords: TS, set, monitor

3510: .seealso:  TSMonitorSet()
3511: @*/
3512: PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscViewerAndFormat *vf)
3513: {
3515:   PetscViewer    viewer =  vf->viewer;
3516:   PetscBool      iascii,ibinary;

3520:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
3521:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
3522:   PetscViewerPushFormat(viewer,vf->format);
3523:   if (iascii) {
3524:     PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
3525:     if (step == -1){ /* this indicates it is an interpolated solution */
3526:       PetscViewerASCIIPrintf(viewer,"Interpolated solution at time %g between steps %D and %D\n",(double)ptime,ts->steps-1,ts->steps);
3527:     } else {
3528:       PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)\n" : "\n");
3529:     }
3530:     PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
3531:   } else if (ibinary) {
3532:     PetscMPIInt rank;
3533:     MPI_Comm_rank(PetscObjectComm((PetscObject)viewer),&rank);
3534:     if (!rank) {
3535:       PetscBool skipHeader;
3536:       PetscInt  classid = REAL_FILE_CLASSID;

3538:       PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
3539:       if (!skipHeader) {
3540:          PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
3541:        }
3542:       PetscRealView(1,&ptime,viewer);
3543:     } else {
3544:       PetscRealView(0,&ptime,viewer);
3545:     }
3546:   }
3547:   PetscViewerPopFormat(viewer);
3548:   return(0);
3549: }

3551: /*@C
3552:    TSMonitorExtreme - Prints the extreme values of the solution at each timestep

3554:    Level: intermediate

3556: .keywords: TS, set, monitor

3558: .seealso:  TSMonitorSet()
3559: @*/
3560: PetscErrorCode TSMonitorExtreme(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscViewerAndFormat *vf)
3561: {
3563:   PetscViewer    viewer =  vf->viewer;
3564:   PetscBool      iascii;
3565:   PetscReal      max,min;

3567: 
3570:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
3571:   PetscViewerPushFormat(viewer,vf->format);
3572:   if (iascii) {
3573:     VecMax(v,NULL,&max);
3574:     VecMin(v,NULL,&min);
3575:     PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);
3576:     PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g%s max %g min %g\n",step,(double)ts->time_step,(double)ptime,ts->steprollback ? " (r)" : "",(double)max,(double)min);
3577:     PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);
3578:   }
3579:   PetscViewerPopFormat(viewer);
3580:   return(0);
3581: }

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

3586:    Collective on TS

3588:    Input Argument:
3589: +  ts - time stepping context
3590: -  t - time to interpolate to

3592:    Output Argument:
3593: .  U - state at given time

3595:    Level: intermediate

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

3600: .keywords: TS, set

3602: .seealso: TSSetExactFinalTime(), TSSolve()
3603: @*/
3604: PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U)
3605: {

3611:   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);
3612:   if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name);
3613:   (*ts->ops->interpolate)(ts,t,U);
3614:   return(0);
3615: }

3617: /*@
3618:    TSStep - Steps one time step

3620:    Collective on TS

3622:    Input Parameter:
3623: .  ts - the TS context obtained from TSCreate()

3625:    Level: developer

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

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

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

3636: .keywords: TS, timestep, solve

3638: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate()
3639: @*/
3640: PetscErrorCode  TSStep(TS ts)
3641: {
3642:   PetscErrorCode   ierr;
3643:   static PetscBool cite = PETSC_FALSE;
3644:   PetscReal        ptime;

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

3656:   TSSetUp(ts);
3657:   TSTrajectorySetUp(ts->trajectory,ts);

3659:   if (ts->max_time >= PETSC_MAX_REAL && ts->max_steps == PETSC_MAX_INT) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
3660:   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()");
3661:   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");

3663:   if (!ts->steps) ts->ptime_prev = ts->ptime;
3664:   ptime = ts->ptime; ts->ptime_prev_rollback = ts->ptime_prev;
3665:   ts->reason = TS_CONVERGED_ITERATING;
3666:   if (!ts->ops->step) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSStep not implemented for type '%s'",((PetscObject)ts)->type_name);
3667:   PetscLogEventBegin(TS_Step,ts,0,0,0);
3668:   (*ts->ops->step)(ts);
3669:   PetscLogEventEnd(TS_Step,ts,0,0,0);
3670:   ts->ptime_prev = ptime;
3671:   ts->steps++;
3672:   ts->steprollback = PETSC_FALSE;
3673:   ts->steprestart  = PETSC_FALSE;

3675:   if (ts->reason < 0) {
3676:     if (ts->errorifstepfailed) {
3677:       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]);
3678:       else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]);
3679:     }
3680:   } else if (!ts->reason) {
3681:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3682:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3683:   }
3684:   return(0);
3685: }

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

3691:    Collective on TS

3693:    Input Arguments:
3694: +  ts - time stepping context
3695: .  wnormtype - norm type, either NORM_2 or NORM_INFINITY
3696: -  order - optional, desired order for the error evaluation or PETSC_DECIDE

3698:    Output Arguments:
3699: +  order - optional, the actual order of the error evaluation
3700: -  wlte - the weighted local truncation error norm

3702:    Level: advanced

3704:    Notes:
3705:    If the timestepper cannot evaluate the error in a particular step
3706:    (eg. in the first step or restart steps after event handling),
3707:    this routine returns wlte=-1.0 .

3709: .seealso: TSStep(), TSAdapt, TSErrorWeightedNorm()
3710: @*/
3711: PetscErrorCode TSEvaluateWLTE(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
3712: {

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

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

3731:    Collective on TS

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

3738:    Output Arguments:
3739: .  U - state at the end of the current step

3741:    Level: advanced

3743:    Notes:
3744:    This function cannot be called until all stages have been evaluated.
3745:    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.

3747: .seealso: TSStep(), TSAdapt
3748: @*/
3749: PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done)
3750: {

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

3762: /*@
3763:    TSSolve - Steps the requested number of timesteps.

3765:    Collective on TS

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

3772:    Level: beginner

3774:    Notes:
3775:    The final time returned by this function may be different from the time of the internally
3776:    held state accessible by TSGetSolution() and TSGetTime() because the method may have
3777:    stepped over the final time.

3779: .keywords: TS, timestep, solve

3781: .seealso: TSCreate(), TSSetSolution(), TSStep(), TSGetTime(), TSGetSolveTime()
3782: @*/
3783: PetscErrorCode TSSolve(TS ts,Vec u)
3784: {
3785:   Vec               solution;
3786:   PetscErrorCode    ierr;


3792:   if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && u) {   /* Need ts->vec_sol to be distinct so it is not overwritten when we interpolate at the end */
3793:     if (!ts->vec_sol || u == ts->vec_sol) {
3794:       VecDuplicate(u,&solution);
3795:       TSSetSolution(ts,solution);
3796:       VecDestroy(&solution); /* grant ownership */
3797:     }
3798:     VecCopy(u,ts->vec_sol);
3799:     if (ts->forward_solve) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
3800:   } else if (u) {
3801:     TSSetSolution(ts,u);
3802:   }
3803:   TSSetUp(ts);
3804:   TSTrajectorySetUp(ts->trajectory,ts);

3806:   if (ts->max_time >= PETSC_MAX_REAL && ts->max_steps == PETSC_MAX_INT) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
3807:   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()");
3808:   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");

3810:   if (ts->forward_solve) {
3811:     TSForwardSetUp(ts);
3812:   }

3814:   /* reset number of steps only when the step is not restarted. ARKIMEX
3815:      restarts the step after an event. Resetting these counters in such case causes
3816:      TSTrajectory to incorrectly save the output files
3817:   */
3818:   /* reset time step and iteration counters */
3819:   if (!ts->steps) {
3820:     ts->ksp_its           = 0;
3821:     ts->snes_its          = 0;
3822:     ts->num_snes_failures = 0;
3823:     ts->reject            = 0;
3824:     ts->steprestart       = PETSC_TRUE;
3825:     ts->steprollback      = PETSC_FALSE;
3826:   }
3827:   if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP && ts->ptime + ts->time_step > ts->max_time) ts->time_step = ts->max_time - ts->ptime;
3828:   ts->reason = TS_CONVERGED_ITERATING;

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

3832:   if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
3833:     (*ts->ops->solve)(ts);
3834:     if (u) {VecCopy(ts->vec_sol,u);}
3835:     ts->solvetime = ts->ptime;
3836:     solution = ts->vec_sol;
3837:   } else { /* Step the requested number of timesteps. */
3838:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3839:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;

3841:     if (!ts->steps) {
3842:       TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);
3843:       TSEventInitialize(ts->event,ts,ts->ptime,ts->vec_sol);
3844:     }

3846:     while (!ts->reason) {
3847:       TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);
3848:       if (!ts->steprollback) {
3849:         TSPreStep(ts);
3850:       }
3851:       TSStep(ts);
3852:       if (ts->testjacobian) {
3853:         TSRHSJacobianTest(ts,NULL);
3854:       }
3855:       if (ts->testjacobiantranspose) {
3856:         TSRHSJacobianTestTranspose(ts,NULL);
3857:       }
3858:       if (ts->vec_costintegral && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
3859:         TSForwardCostIntegral(ts);
3860:       }
3861:       if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
3862:         TSForwardStep(ts);
3863:       }
3864:       TSPostEvaluate(ts);
3865:       TSEventHandler(ts); /* The right-hand side may be changed due to event. Be careful with Any computation using the RHS information after this point. */
3866:       if (ts->steprollback) {
3867:         TSPostEvaluate(ts);
3868:       }
3869:       if (!ts->steprollback) {
3870:         TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);
3871:         TSPostStep(ts);
3872:       }
3873:     }
3874:     TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);

3876:     if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
3877:       TSInterpolate(ts,ts->max_time,u);
3878:       ts->solvetime = ts->max_time;
3879:       solution = u;
3880:       TSMonitor(ts,-1,ts->solvetime,solution);
3881:     } else {
3882:       if (u) {VecCopy(ts->vec_sol,u);}
3883:       ts->solvetime = ts->ptime;
3884:       solution = ts->vec_sol;
3885:     }
3886:   }

3888:   TSViewFromOptions(ts,NULL,"-ts_view");
3889:   VecViewFromOptions(solution,NULL,"-ts_view_solution");
3890:   PetscObjectSAWsBlock((PetscObject)ts);
3891:   if (ts->adjoint_solve) {
3892:     TSAdjointSolve(ts);
3893:   }
3894:   return(0);
3895: }

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

3900:    Collective on TS

3902:    Input Parameters:
3903: +  ts - time stepping context obtained from TSCreate()
3904: .  step - step number that has just completed
3905: .  ptime - model time of the state
3906: -  u - state at the current model time

3908:    Notes:
3909:    TSMonitor() is typically used automatically within the time stepping implementations.
3910:    Users would almost never call this routine directly.

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

3914:    Level: developer

3916: .keywords: TS, timestep
3917: @*/
3918: PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u)
3919: {
3920:   DM             dm;
3921:   PetscInt       i,n = ts->numbermonitors;


3928:   TSGetDM(ts,&dm);
3929:   DMSetOutputSequenceNumber(dm,step,ptime);

3931:   VecLockReadPush(u);
3932:   for (i=0; i<n; i++) {
3933:     (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);
3934:   }
3935:   VecLockReadPop(u);
3936:   return(0);
3937: }

3939: /* ------------------------------------------------------------------------*/
3940: /*@C
3941:    TSMonitorLGCtxCreate - Creates a TSMonitorLGCtx context for use with
3942:    TS to monitor the solution process graphically in various ways

3944:    Collective on TS

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

3953:    Output Parameter:
3954: .  ctx - the context

3956:    Options Database Key:
3957: +  -ts_monitor_lg_timestep - automatically sets line graph monitor
3958: +  -ts_monitor_lg_timestep_log - automatically sets line graph monitor
3959: .  -ts_monitor_lg_solution - monitor the solution (or certain values of the solution by calling TSMonitorLGSetDisplayVariables() or TSMonitorLGCtxSetDisplayVariables())
3960: .  -ts_monitor_lg_error -  monitor the error
3961: .  -ts_monitor_lg_ksp_iterations - monitor the number of KSP iterations needed for each timestep
3962: .  -ts_monitor_lg_snes_iterations - monitor the number of SNES iterations needed for each timestep
3963: -  -lg_use_markers <true,false> - mark the data points (at each time step) on the plot; default is true

3965:    Notes:
3966:    Use TSMonitorLGCtxDestroy() to destroy.

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

3970:    Many of the functions that control the monitoring have two forms: TSMonitorLGSet/GetXXXX() and TSMonitorLGCtxSet/GetXXXX() the first take a TS object as the
3971:    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
3972:    as the first argument.

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

3976:    Level: intermediate

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

3980: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError(), TSMonitorDefault(), VecView(),
3981:            TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
3982:            TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
3983:            TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
3984:            TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()

3986: @*/
3987: PetscErrorCode  TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx)
3988: {
3989:   PetscDraw      draw;

3993:   PetscNew(ctx);
3994:   PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
3995:   PetscDrawSetFromOptions(draw);
3996:   PetscDrawLGCreate(draw,1,&(*ctx)->lg);
3997:   PetscDrawLGSetFromOptions((*ctx)->lg);
3998:   PetscDrawDestroy(&draw);
3999:   (*ctx)->howoften = howoften;
4000:   return(0);
4001: }

4003: PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx)
4004: {
4005:   TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
4006:   PetscReal      x   = ptime,y;

4010:   if (step < 0) return(0); /* -1 indicates an interpolated solution */
4011:   if (!step) {
4012:     PetscDrawAxis axis;
4013:     const char *ylabel = ctx->semilogy ? "Log Time Step" : "Time Step";
4014:     PetscDrawLGGetAxis(ctx->lg,&axis);
4015:     PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time",ylabel);
4016:     PetscDrawLGReset(ctx->lg);
4017:   }
4018:   TSGetTimeStep(ts,&y);
4019:   if (ctx->semilogy) y = PetscLog10Real(y);
4020:   PetscDrawLGAddPoint(ctx->lg,&x,&y);
4021:   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
4022:     PetscDrawLGDraw(ctx->lg);
4023:     PetscDrawLGSave(ctx->lg);
4024:   }
4025:   return(0);
4026: }

4028: /*@C
4029:    TSMonitorLGCtxDestroy - Destroys a line graph context that was created
4030:    with TSMonitorLGCtxCreate().

4032:    Collective on TSMonitorLGCtx

4034:    Input Parameter:
4035: .  ctx - the monitor context

4037:    Level: intermediate

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

4041: .seealso: TSMonitorLGCtxCreate(),  TSMonitorSet(), TSMonitorLGTimeStep();
4042: @*/
4043: PetscErrorCode  TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx)
4044: {

4048:   if ((*ctx)->transformdestroy) {
4049:     ((*ctx)->transformdestroy)((*ctx)->transformctx);
4050:   }
4051:   PetscDrawLGDestroy(&(*ctx)->lg);
4052:   PetscStrArrayDestroy(&(*ctx)->names);
4053:   PetscStrArrayDestroy(&(*ctx)->displaynames);
4054:   PetscFree((*ctx)->displayvariables);
4055:   PetscFree((*ctx)->displayvalues);
4056:   PetscFree(*ctx);
4057:   return(0);
4058: }

4060: /*
4061:   
4062:   Creates a TS Monitor SPCtx for use with DM Swarm particle visualizations

4064: */
4065: PetscErrorCode TSMonitorSPCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorSPCtx *ctx)
4066: {
4067:   PetscDraw      draw;

4071:   PetscNew(ctx);
4072:   PetscDrawCreate(comm,host,label,x,y,m,n,&draw);
4073:   PetscDrawSetFromOptions(draw);
4074:   PetscDrawSPCreate(draw,1,&(*ctx)->sp);
4075:   PetscDrawDestroy(&draw);
4076:   (*ctx)->howoften = howoften;
4077:   return(0);

4079: }

4081: /* 
4082:   Destroys a TSMonitorSPCtx that was created with TSMonitorSPCtxCreate 
4083: */
4084: PetscErrorCode TSMonitorSPCtxDestroy(TSMonitorSPCtx *ctx)
4085: {

4089: 
4090:   PetscDrawSPDestroy(&(*ctx)->sp);
4091:   PetscFree(*ctx);
4092: 
4093:   return(0);

4095: }

4097: /*@
4098:    TSGetTime - Gets the time of the most recently completed step.

4100:    Not Collective

4102:    Input Parameter:
4103: .  ts - the TS context obtained from TSCreate()

4105:    Output Parameter:
4106: .  t  - the current time. This time may not corresponds to the final time set with TSSetMaxTime(), use TSGetSolveTime().

4108:    Level: beginner

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

4114: .seealso:  TSGetSolveTime(), TSSetTime(), TSGetTimeStep()

4116: .keywords: TS, get, time
4117: @*/
4118: PetscErrorCode  TSGetTime(TS ts,PetscReal *t)
4119: {
4123:   *t = ts->ptime;
4124:   return(0);
4125: }

4127: /*@
4128:    TSGetPrevTime - Gets the starting time of the previously completed step.

4130:    Not Collective

4132:    Input Parameter:
4133: .  ts - the TS context obtained from TSCreate()

4135:    Output Parameter:
4136: .  t  - the previous time

4138:    Level: beginner

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

4142: .keywords: TS, get, time
4143: @*/
4144: PetscErrorCode  TSGetPrevTime(TS ts,PetscReal *t)
4145: {
4149:   *t = ts->ptime_prev;
4150:   return(0);
4151: }

4153: /*@
4154:    TSSetTime - Allows one to reset the time.

4156:    Logically Collective on TS

4158:    Input Parameters:
4159: +  ts - the TS context obtained from TSCreate()
4160: -  time - the time

4162:    Level: intermediate

4164: .seealso: TSGetTime(), TSSetMaxSteps()

4166: .keywords: TS, set, time
4167: @*/
4168: PetscErrorCode  TSSetTime(TS ts, PetscReal t)
4169: {
4173:   ts->ptime = t;
4174:   return(0);
4175: }

4177: /*@C
4178:    TSSetOptionsPrefix - Sets the prefix used for searching for all
4179:    TS options in the database.

4181:    Logically Collective on TS

4183:    Input Parameter:
4184: +  ts     - The TS context
4185: -  prefix - The prefix to prepend to all option names

4187:    Notes:
4188:    A hyphen (-) must NOT be given at the beginning of the prefix name.
4189:    The first character of all runtime options is AUTOMATICALLY the
4190:    hyphen.

4192:    Level: advanced

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

4196: .seealso: TSSetFromOptions()

4198: @*/
4199: PetscErrorCode  TSSetOptionsPrefix(TS ts,const char prefix[])
4200: {
4202:   SNES           snes;

4206:   PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);
4207:   TSGetSNES(ts,&snes);
4208:   SNESSetOptionsPrefix(snes,prefix);
4209:   return(0);
4210: }

4212: /*@C
4213:    TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4214:    TS options in the database.

4216:    Logically Collective on TS

4218:    Input Parameter:
4219: +  ts     - The TS context
4220: -  prefix - The prefix to prepend to all option names

4222:    Notes:
4223:    A hyphen (-) must NOT be given at the beginning of the prefix name.
4224:    The first character of all runtime options is AUTOMATICALLY the
4225:    hyphen.

4227:    Level: advanced

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

4231: .seealso: TSGetOptionsPrefix()

4233: @*/
4234: PetscErrorCode  TSAppendOptionsPrefix(TS ts,const char prefix[])
4235: {
4237:   SNES           snes;

4241:   PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);
4242:   TSGetSNES(ts,&snes);
4243:   SNESAppendOptionsPrefix(snes,prefix);
4244:   return(0);
4245: }

4247: /*@C
4248:    TSGetOptionsPrefix - Sets the prefix used for searching for all
4249:    TS options in the database.

4251:    Not Collective

4253:    Input Parameter:
4254: .  ts - The TS context

4256:    Output Parameter:
4257: .  prefix - A pointer to the prefix string used

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

4263:    Level: intermediate

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

4267: .seealso: TSAppendOptionsPrefix()
4268: @*/
4269: PetscErrorCode  TSGetOptionsPrefix(TS ts,const char *prefix[])
4270: {

4276:   PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);
4277:   return(0);
4278: }

4280: /*@C
4281:    TSGetRHSJacobian - Returns the Jacobian J at the present timestep.

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

4285:    Input Parameter:
4286: .  ts  - The TS context obtained from TSCreate()

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

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

4297:    Level: intermediate

4299: .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetStepNumber()

4301: .keywords: TS, timestep, get, matrix, Jacobian
4302: @*/
4303: PetscErrorCode  TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx)
4304: {
4306:   DM             dm;

4309:   if (Amat || Pmat) {
4310:     SNES snes;
4311:     TSGetSNES(ts,&snes);
4312:     SNESSetUpMatrices(snes);
4313:     SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
4314:   }
4315:   TSGetDM(ts,&dm);
4316:   DMTSGetRHSJacobian(dm,func,ctx);
4317:   return(0);
4318: }

4320: /*@C
4321:    TSGetIJacobian - Returns the implicit Jacobian at the present timestep.

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

4325:    Input Parameter:
4326: .  ts  - The TS context obtained from TSCreate()

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

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

4337:    Level: advanced

4339: .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetStepNumber()

4341: .keywords: TS, timestep, get, matrix, Jacobian
4342: @*/
4343: PetscErrorCode  TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx)
4344: {
4346:   DM             dm;

4349:   if (Amat || Pmat) {
4350:     SNES snes;
4351:     TSGetSNES(ts,&snes);
4352:     SNESSetUpMatrices(snes);
4353:     SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);
4354:   }
4355:   TSGetDM(ts,&dm);
4356:   DMTSGetIJacobian(dm,f,ctx);
4357:   return(0);
4358: }

4360: /*@C
4361:    TSMonitorDrawSolution - Monitors progress of the TS solvers by calling
4362:    VecView() for the solution at each timestep

4364:    Collective on TS

4366:    Input Parameters:
4367: +  ts - the TS context
4368: .  step - current time-step
4369: .  ptime - current time
4370: -  dummy - either a viewer or NULL

4372:    Options Database:
4373: .   -ts_monitor_draw_solution_initial - show initial solution as well as current solution

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

4379:    Level: intermediate

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

4383: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4384: @*/
4385: PetscErrorCode  TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4386: {
4387:   PetscErrorCode   ierr;
4388:   TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy;
4389:   PetscDraw        draw;

4392:   if (!step && ictx->showinitial) {
4393:     if (!ictx->initialsolution) {
4394:       VecDuplicate(u,&ictx->initialsolution);
4395:     }
4396:     VecCopy(u,ictx->initialsolution);
4397:   }
4398:   if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) return(0);

4400:   if (ictx->showinitial) {
4401:     PetscReal pause;
4402:     PetscViewerDrawGetPause(ictx->viewer,&pause);
4403:     PetscViewerDrawSetPause(ictx->viewer,0.0);
4404:     VecView(ictx->initialsolution,ictx->viewer);
4405:     PetscViewerDrawSetPause(ictx->viewer,pause);
4406:     PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);
4407:   }
4408:   VecView(u,ictx->viewer);
4409:   if (ictx->showtimestepandtime) {
4410:     PetscReal xl,yl,xr,yr,h;
4411:     char      time[32];

4413:     PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
4414:     PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
4415:     PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
4416:     h    = yl + .95*(yr - yl);
4417:     PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
4418:     PetscDrawFlush(draw);
4419:   }

4421:   if (ictx->showinitial) {
4422:     PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);
4423:   }
4424:   return(0);
4425: }

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

4430:    Collective on TS

4432:    Input Parameters:
4433: +  ts - the TS context
4434: .  step - current time-step
4435: .  ptime - current time
4436: -  dummy - either a viewer or NULL

4438:    Level: intermediate

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

4442: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
4443: @*/
4444: PetscErrorCode  TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4445: {
4446:   PetscErrorCode    ierr;
4447:   TSMonitorDrawCtx  ictx = (TSMonitorDrawCtx)dummy;
4448:   PetscDraw         draw;
4449:   PetscDrawAxis     axis;
4450:   PetscInt          n;
4451:   PetscMPIInt       size;
4452:   PetscReal         U0,U1,xl,yl,xr,yr,h;
4453:   char              time[32];
4454:   const PetscScalar *U;

4457:   MPI_Comm_size(PetscObjectComm((PetscObject)ts),&size);
4458:   if (size != 1) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Only allowed for sequential runs");
4459:   VecGetSize(u,&n);
4460:   if (n != 2) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Only for ODEs with two unknowns");

4462:   PetscViewerDrawGetDraw(ictx->viewer,0,&draw);
4463:   PetscViewerDrawGetDrawAxis(ictx->viewer,0,&axis);
4464:   PetscDrawAxisGetLimits(axis,&xl,&xr,&yl,&yr);
4465:   if (!step) {
4466:     PetscDrawClear(draw);
4467:     PetscDrawAxisDraw(axis);
4468:   }

4470:   VecGetArrayRead(u,&U);
4471:   U0 = PetscRealPart(U[0]);
4472:   U1 = PetscRealPart(U[1]);
4473:   VecRestoreArrayRead(u,&U);
4474:   if ((U0 < xl) || (U1 < yl) || (U0 > xr) || (U1 > yr)) return(0);

4476:   PetscDrawCollectiveBegin(draw);
4477:   PetscDrawPoint(draw,U0,U1,PETSC_DRAW_BLACK);
4478:   if (ictx->showtimestepandtime) {
4479:     PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
4480:     PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);
4481:     h    = yl + .95*(yr - yl);
4482:     PetscDrawStringCentered(draw,.5*(xl+xr),h,PETSC_DRAW_BLACK,time);
4483:   }
4484:   PetscDrawCollectiveEnd(draw);
4485:   PetscDrawFlush(draw);
4486:   PetscDrawPause(draw);
4487:   PetscDrawSave(draw);
4488:   return(0);
4489: }

4491: /*@C
4492:    TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution()

4494:    Collective on TS

4496:    Input Parameters:
4497: .    ctx - the monitor context

4499:    Level: intermediate

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

4503: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError()
4504: @*/
4505: PetscErrorCode  TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx)
4506: {

4510:   PetscViewerDestroy(&(*ictx)->viewer);
4511:   VecDestroy(&(*ictx)->initialsolution);
4512:   PetscFree(*ictx);
4513:   return(0);
4514: }

4516: /*@C
4517:    TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx

4519:    Collective on TS

4521:    Input Parameter:
4522: .    ts - time-step context

4524:    Output Patameter:
4525: .    ctx - the monitor context

4527:    Options Database:
4528: .   -ts_monitor_draw_solution_initial - show initial solution as well as current solution

4530:    Level: intermediate

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

4534: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx()
4535: @*/
4536: PetscErrorCode  TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx)
4537: {
4538:   PetscErrorCode   ierr;

4541:   PetscNew(ctx);
4542:   PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);
4543:   PetscViewerSetFromOptions((*ctx)->viewer);

4545:   (*ctx)->howoften    = howoften;
4546:   (*ctx)->showinitial = PETSC_FALSE;
4547:   PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);

4549:   (*ctx)->showtimestepandtime = PETSC_FALSE;
4550:   PetscOptionsGetBool(NULL,NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);
4551:   return(0);
4552: }

4554: /*@C
4555:    TSMonitorDrawSolutionFunction - Monitors progress of the TS solvers by calling
4556:    VecView() for the solution provided by TSSetSolutionFunction() at each timestep

4558:    Collective on TS

4560:    Input Parameters:
4561: +  ts - the TS context
4562: .  step - current time-step
4563: .  ptime - current time
4564: -  dummy - either a viewer or NULL

4566:    Options Database:
4567: .  -ts_monitor_draw_solution_function - Monitor error graphically, requires user to have provided TSSetSolutionFunction()

4569:    Level: intermediate

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

4573: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
4574: @*/
4575: PetscErrorCode  TSMonitorDrawSolutionFunction(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4576: {
4577:   PetscErrorCode   ierr;
4578:   TSMonitorDrawCtx ctx    = (TSMonitorDrawCtx)dummy;
4579:   PetscViewer      viewer = ctx->viewer;
4580:   Vec              work;

4583:   if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) return(0);
4584:   VecDuplicate(u,&work);
4585:   TSComputeSolutionFunction(ts,ptime,work);
4586:   VecView(work,viewer);
4587:   VecDestroy(&work);
4588:   return(0);
4589: }

4591: /*@C
4592:    TSMonitorDrawError - Monitors progress of the TS solvers by calling
4593:    VecView() for the error at each timestep

4595:    Collective on TS

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

4603:    Options Database:
4604: .  -ts_monitor_draw_error - Monitor error graphically, requires user to have provided TSSetSolutionFunction()

4606:    Level: intermediate

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

4610: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
4611: @*/
4612: PetscErrorCode  TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
4613: {
4614:   PetscErrorCode   ierr;
4615:   TSMonitorDrawCtx ctx    = (TSMonitorDrawCtx)dummy;
4616:   PetscViewer      viewer = ctx->viewer;
4617:   Vec              work;

4620:   if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) return(0);
4621:   VecDuplicate(u,&work);
4622:   TSComputeSolutionFunction(ts,ptime,work);
4623:   VecAXPY(work,-1.0,u);
4624:   VecView(work,viewer);
4625:   VecDestroy(&work);
4626:   return(0);
4627: }

4629:  #include <petsc/private/dmimpl.h>
4630: /*@
4631:    TSSetDM - Sets the DM that may be used by some nonlinear solvers or preconditioners under the TS

4633:    Logically Collective on TS and DM

4635:    Input Parameters:
4636: +  ts - the ODE integrator object
4637: -  dm - the dm, cannot be NULL

4639:    Notes:
4640:    A DM can only be used for solving one problem at a time because information about the problem is stored on the DM,
4641:    even when not using interfaces like DMTSSetIFunction().  Use DMClone() to get a distinct DM when solving
4642:    different problems using the same function space.

4644:    Level: intermediate

4646: .seealso: TSGetDM(), SNESSetDM(), SNESGetDM()
4647: @*/
4648: PetscErrorCode  TSSetDM(TS ts,DM dm)
4649: {
4651:   SNES           snes;
4652:   DMTS           tsdm;

4657:   PetscObjectReference((PetscObject)dm);
4658:   if (ts->dm) {               /* Move the DMTS context over to the new DM unless the new DM already has one */
4659:     if (ts->dm->dmts && !dm->dmts) {
4660:       DMCopyDMTS(ts->dm,dm);
4661:       DMGetDMTS(ts->dm,&tsdm);
4662:       if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */
4663:         tsdm->originaldm = dm;
4664:       }
4665:     }
4666:     DMDestroy(&ts->dm);
4667:   }
4668:   ts->dm = dm;

4670:   TSGetSNES(ts,&snes);
4671:   SNESSetDM(snes,dm);
4672:   return(0);
4673: }

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

4678:    Not Collective

4680:    Input Parameter:
4681: . ts - the preconditioner context

4683:    Output Parameter:
4684: .  dm - the dm

4686:    Level: intermediate

4688: .seealso: TSSetDM(), SNESSetDM(), SNESGetDM()
4689: @*/
4690: PetscErrorCode  TSGetDM(TS ts,DM *dm)
4691: {

4696:   if (!ts->dm) {
4697:     DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);
4698:     if (ts->snes) {SNESSetDM(ts->snes,ts->dm);}
4699:   }
4700:   *dm = ts->dm;
4701:   return(0);
4702: }

4704: /*@
4705:    SNESTSFormFunction - Function to evaluate nonlinear residual

4707:    Logically Collective on SNES

4709:    Input Parameter:
4710: + snes - nonlinear solver
4711: . U - the current state at which to evaluate the residual
4712: - ctx - user context, must be a TS

4714:    Output Parameter:
4715: . F - the nonlinear residual

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

4721:    Level: advanced

4723: .seealso: SNESSetFunction(), MatFDColoringSetFunction()
4724: @*/
4725: PetscErrorCode  SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx)
4726: {
4727:   TS             ts = (TS)ctx;

4735:   (ts->ops->snesfunction)(snes,U,F,ts);
4736:   return(0);
4737: }

4739: /*@
4740:    SNESTSFormJacobian - Function to evaluate the Jacobian

4742:    Collective on SNES

4744:    Input Parameter:
4745: + snes - nonlinear solver
4746: . U - the current state at which to evaluate the residual
4747: - ctx - user context, must be a TS

4749:    Output Parameter:
4750: + A - the Jacobian
4751: . B - the preconditioning matrix (may be the same as A)
4752: - flag - indicates any structure change in the matrix

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

4757:    Level: developer

4759: .seealso: SNESSetJacobian()
4760: @*/
4761: PetscErrorCode  SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx)
4762: {
4763:   TS             ts = (TS)ctx;

4774:   (ts->ops->snesjacobian)(snes,U,A,B,ts);
4775:   return(0);
4776: }

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

4781:    Collective on TS

4783:    Input Arguments:
4784: +  ts - time stepping context
4785: .  t - time at which to evaluate
4786: .  U - state at which to evaluate
4787: -  ctx - context

4789:    Output Arguments:
4790: .  F - right hand side

4792:    Level: intermediate

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

4798: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant()
4799: @*/
4800: PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
4801: {
4803:   Mat            Arhs,Brhs;

4806:   TSGetRHSMats_Private(ts,&Arhs,&Brhs);
4807:   TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);
4808:   MatMult(Arhs,U,F);
4809:   return(0);
4810: }

4812: /*@C
4813:    TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.

4815:    Collective on TS

4817:    Input Arguments:
4818: +  ts - time stepping context
4819: .  t - time at which to evaluate
4820: .  U - state at which to evaluate
4821: -  ctx - context

4823:    Output Arguments:
4824: +  A - pointer to operator
4825: .  B - pointer to preconditioning matrix
4826: -  flg - matrix structure flag

4828:    Level: intermediate

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

4833: .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear()
4834: @*/
4835: PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
4836: {
4838:   return(0);
4839: }

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

4844:    Collective on TS

4846:    Input Arguments:
4847: +  ts - time stepping context
4848: .  t - time at which to evaluate
4849: .  U - state at which to evaluate
4850: .  Udot - time derivative of state vector
4851: -  ctx - context

4853:    Output Arguments:
4854: .  F - left hand side

4856:    Level: intermediate

4858:    Notes:
4859:    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
4860:    user is required to write their own TSComputeIFunction.
4861:    This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems.
4862:    The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian().

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

4866: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant(), TSComputeRHSFunctionLinear()
4867: @*/
4868: PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
4869: {
4871:   Mat            A,B;

4874:   TSGetIJacobian(ts,&A,&B,NULL,NULL);
4875:   TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);
4876:   MatMult(A,Udot,F);
4877:   return(0);
4878: }

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

4883:    Collective on TS

4885:    Input Arguments:
4886: +  ts - time stepping context
4887: .  t - time at which to evaluate
4888: .  U - state at which to evaluate
4889: .  Udot - time derivative of state vector
4890: .  shift - shift to apply
4891: -  ctx - context

4893:    Output Arguments:
4894: +  A - pointer to operator
4895: .  B - pointer to preconditioning matrix
4896: -  flg - matrix structure flag

4898:    Level: advanced

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

4903:    It is only appropriate for problems of the form

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

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

4911: $    shift*M + J

4913:   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
4914:   a copy of M or reassemble it when requested.

4916: .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear()
4917: @*/
4918: PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx)
4919: {

4923:   MatScale(A, shift / ts->ijacobian.shift);
4924:   ts->ijacobian.shift = shift;
4925:   return(0);
4926: }

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

4931:    Not Collective

4933:    Input Parameter:
4934: .  ts - the TS context

4936:    Output Parameter:
4937: .  equation_type - see TSEquationType

4939:    Level: beginner

4941: .keywords: TS, equation type

4943: .seealso: TSSetEquationType(), TSEquationType
4944: @*/
4945: PetscErrorCode  TSGetEquationType(TS ts,TSEquationType *equation_type)
4946: {
4950:   *equation_type = ts->equation_type;
4951:   return(0);
4952: }

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

4957:    Not Collective

4959:    Input Parameter:
4960: +  ts - the TS context
4961: -  equation_type - see TSEquationType

4963:    Level: advanced

4965: .keywords: TS, equation type

4967: .seealso: TSGetEquationType(), TSEquationType
4968: @*/
4969: PetscErrorCode  TSSetEquationType(TS ts,TSEquationType equation_type)
4970: {
4973:   ts->equation_type = equation_type;
4974:   return(0);
4975: }

4977: /*@
4978:    TSGetConvergedReason - Gets the reason the TS iteration was stopped.

4980:    Not Collective

4982:    Input Parameter:
4983: .  ts - the TS context

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

4989:    Level: beginner

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

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

4996: .seealso: TSSetConvergenceTest(), TSConvergedReason
4997: @*/
4998: PetscErrorCode  TSGetConvergedReason(TS ts,TSConvergedReason *reason)
4999: {
5003:   *reason = ts->reason;
5004:   return(0);
5005: }

5007: /*@
5008:    TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve.

5010:    Not Collective

5012:    Input Parameter:
5013: +  ts - the TS context
5014: .  reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
5015:             manual pages for the individual convergence tests for complete lists

5017:    Level: advanced

5019:    Notes:
5020:    Can only be called during TSSolve() is active.

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

5024: .seealso: TSConvergedReason
5025: @*/
5026: PetscErrorCode  TSSetConvergedReason(TS ts,TSConvergedReason reason)
5027: {
5030:   ts->reason = reason;
5031:   return(0);
5032: }

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

5037:    Not Collective

5039:    Input Parameter:
5040: .  ts - the TS context

5042:    Output Parameter:
5043: .  ftime - the final time. This time corresponds to the final time set with TSSetMaxTime()

5045:    Level: beginner

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

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

5052: .seealso: TSSetConvergenceTest(), TSConvergedReason
5053: @*/
5054: PetscErrorCode  TSGetSolveTime(TS ts,PetscReal *ftime)
5055: {
5059:   *ftime = ts->solvetime;
5060:   return(0);
5061: }

5063: /*@
5064:    TSGetSNESIterations - Gets the total number of nonlinear iterations
5065:    used by the time integrator.

5067:    Not Collective

5069:    Input Parameter:
5070: .  ts - TS context

5072:    Output Parameter:
5073: .  nits - number of nonlinear iterations

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

5078:    Level: intermediate

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

5082: .seealso:  TSGetKSPIterations()
5083: @*/
5084: PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits)
5085: {
5089:   *nits = ts->snes_its;
5090:   return(0);
5091: }

5093: /*@
5094:    TSGetKSPIterations - Gets the total number of linear iterations
5095:    used by the time integrator.

5097:    Not Collective

5099:    Input Parameter:
5100: .  ts - TS context

5102:    Output Parameter:
5103: .  lits - number of linear iterations

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

5108:    Level: intermediate

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

5112: .seealso:  TSGetSNESIterations(), SNESGetKSPIterations()
5113: @*/
5114: PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits)
5115: {
5119:   *lits = ts->ksp_its;
5120:   return(0);
5121: }

5123: /*@
5124:    TSGetStepRejections - Gets the total number of rejected steps.

5126:    Not Collective

5128:    Input Parameter:
5129: .  ts - TS context

5131:    Output Parameter:
5132: .  rejects - number of steps rejected

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

5137:    Level: intermediate

5139: .keywords: TS, get, number

5141: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails()
5142: @*/
5143: PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects)
5144: {
5148:   *rejects = ts->reject;
5149:   return(0);
5150: }

5152: /*@
5153:    TSGetSNESFailures - Gets the total number of failed SNES solves

5155:    Not Collective

5157:    Input Parameter:
5158: .  ts - TS context

5160:    Output Parameter:
5161: .  fails - number of failed nonlinear solves

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

5166:    Level: intermediate

5168: .keywords: TS, get, number

5170: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures()
5171: @*/
5172: PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails)
5173: {
5177:   *fails = ts->num_snes_failures;
5178:   return(0);
5179: }

5181: /*@
5182:    TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails

5184:    Not Collective

5186:    Input Parameter:
5187: +  ts - TS context
5188: -  rejects - maximum number of rejected steps, pass -1 for unlimited

5190:    Notes:
5191:    The counter is reset to zero for each step

5193:    Options Database Key:
5194:  .  -ts_max_reject - Maximum number of step rejections before a step fails

5196:    Level: intermediate

5198: .keywords: TS, set, maximum, number

5200: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
5201: @*/
5202: PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects)
5203: {
5206:   ts->max_reject = rejects;
5207:   return(0);
5208: }

5210: /*@
5211:    TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves

5213:    Not Collective

5215:    Input Parameter:
5216: +  ts - TS context
5217: -  fails - maximum number of failed nonlinear solves, pass -1 for unlimited

5219:    Notes:
5220:    The counter is reset to zero for each successive call to TSSolve().

5222:    Options Database Key:
5223:  .  -ts_max_snes_failures - Maximum number of nonlinear solve failures

5225:    Level: intermediate

5227: .keywords: TS, set, maximum, number

5229: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason()
5230: @*/
5231: PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails)
5232: {
5235:   ts->max_snes_failures = fails;
5236:   return(0);
5237: }

5239: /*@
5240:    TSSetErrorIfStepFails - Error if no step succeeds

5242:    Not Collective

5244:    Input Parameter:
5245: +  ts - TS context
5246: -  err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure

5248:    Options Database Key:
5249:  .  -ts_error_if_step_fails - Error if no step succeeds

5251:    Level: intermediate

5253: .keywords: TS, set, error

5255: .seealso:  TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason()
5256: @*/
5257: PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err)
5258: {
5261:   ts->errorifstepfailed = err;
5262:   return(0);
5263: }

5265: /*@C
5266:    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

5268:    Collective on TS

5270:    Input Parameters:
5271: +  ts - the TS context
5272: .  step - current time-step
5273: .  ptime - current time
5274: .  u - current state
5275: -  vf - viewer and its format

5277:    Level: intermediate

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

5281: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5282: @*/
5283: PetscErrorCode  TSMonitorSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscViewerAndFormat *vf)
5284: {

5288:   PetscViewerPushFormat(vf->viewer,vf->format);
5289:   VecView(u,vf->viewer);
5290:   PetscViewerPopFormat(vf->viewer);
5291:   return(0);
5292: }

5294: /*@C
5295:    TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep.

5297:    Collective on TS

5299:    Input Parameters:
5300: +  ts - the TS context
5301: .  step - current time-step
5302: .  ptime - current time
5303: .  u - current state
5304: -  filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)

5306:    Level: intermediate

5308:    Notes:
5309:    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.
5310:    These are named according to the file name template.

5312:    This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy().

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

5316: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView()
5317: @*/
5318: PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate)
5319: {
5321:   char           filename[PETSC_MAX_PATH_LEN];
5322:   PetscViewer    viewer;

5325:   if (step < 0) return(0); /* -1 indicates interpolated solution */
5326:   PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);
5327:   PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);
5328:   VecView(u,viewer);
5329:   PetscViewerDestroy(&viewer);
5330:   return(0);
5331: }

5333: /*@C
5334:    TSMonitorSolutionVTKDestroy - Destroy context for monitoring

5336:    Collective on TS

5338:    Input Parameters:
5339: .  filenametemplate - string containing a format specifier for the integer time step (e.g. %03D)

5341:    Level: intermediate

5343:    Note:
5344:    This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK().

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

5348: .seealso: TSMonitorSet(), TSMonitorSolutionVTK()
5349: @*/
5350: PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate)
5351: {

5355:   PetscFree(*(char**)filenametemplate);
5356:   return(0);
5357: }

5359: /*@
5360:    TSGetAdapt - Get the adaptive controller context for the current method

5362:    Collective on TS if controller has not been created yet

5364:    Input Arguments:
5365: .  ts - time stepping context

5367:    Output Arguments:
5368: .  adapt - adaptive controller

5370:    Level: intermediate

5372: .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose()
5373: @*/
5374: PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt)
5375: {

5381:   if (!ts->adapt) {
5382:     TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);
5383:     PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);
5384:     PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);
5385:   }
5386:   *adapt = ts->adapt;
5387:   return(0);
5388: }

5390: /*@
5391:    TSSetTolerances - Set tolerances for local truncation error when using adaptive controller

5393:    Logically Collective

5395:    Input Arguments:
5396: +  ts - time integration context
5397: .  atol - scalar absolute tolerances, PETSC_DECIDE to leave current value
5398: .  vatol - vector of absolute tolerances or NULL, used in preference to atol if present
5399: .  rtol - scalar relative tolerances, PETSC_DECIDE to leave current value
5400: -  vrtol - vector of relative tolerances or NULL, used in preference to atol if present

5402:    Options Database keys:
5403: +  -ts_rtol <rtol> - relative tolerance for local truncation error
5404: -  -ts_atol <atol> Absolute tolerance for local truncation error

5406:    Notes:
5407:    With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5408:    (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5409:    computed only for the differential or the algebraic part then this can be done using the vector of
5410:    tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
5411:    differential part and infinity for the algebraic part, the LTE calculation will include only the
5412:    differential variables.

5414:    Level: beginner

5416: .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances()
5417: @*/
5418: PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol)
5419: {

5423:   if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol;
5424:   if (vatol) {
5425:     PetscObjectReference((PetscObject)vatol);
5426:     VecDestroy(&ts->vatol);
5427:     ts->vatol = vatol;
5428:   }
5429:   if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol;
5430:   if (vrtol) {
5431:     PetscObjectReference((PetscObject)vrtol);
5432:     VecDestroy(&ts->vrtol);
5433:     ts->vrtol = vrtol;
5434:   }
5435:   return(0);
5436: }

5438: /*@
5439:    TSGetTolerances - Get tolerances for local truncation error when using adaptive controller

5441:    Logically Collective

5443:    Input Arguments:
5444: .  ts - time integration context

5446:    Output Arguments:
5447: +  atol - scalar absolute tolerances, NULL to ignore
5448: .  vatol - vector of absolute tolerances, NULL to ignore
5449: .  rtol - scalar relative tolerances, NULL to ignore
5450: -  vrtol - vector of relative tolerances, NULL to ignore

5452:    Level: beginner

5454: .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances()
5455: @*/
5456: PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol)
5457: {
5459:   if (atol)  *atol  = ts->atol;
5460:   if (vatol) *vatol = ts->vatol;
5461:   if (rtol)  *rtol  = ts->rtol;
5462:   if (vrtol) *vrtol = ts->vrtol;
5463:   return(0);
5464: }

5466: /*@
5467:    TSErrorWeightedNorm2 - compute a weighted 2-norm of the difference between two state vectors

5469:    Collective on TS

5471:    Input Arguments:
5472: +  ts - time stepping context
5473: .  U - state vector, usually ts->vec_sol
5474: -  Y - state vector to be compared to U

5476:    Output Arguments:
5477: .  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5478: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5479: .  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5481:    Level: developer

5483: .seealso: TSErrorWeightedNorm(), TSErrorWeightedNormInfinity()
5484: @*/
5485: PetscErrorCode TSErrorWeightedNorm2(TS ts,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
5486: {
5487:   PetscErrorCode    ierr;
5488:   PetscInt          i,n,N,rstart;
5489:   PetscInt          n_loc,na_loc,nr_loc;
5490:   PetscReal         n_glb,na_glb,nr_glb;
5491:   const PetscScalar *u,*y;
5492:   PetscReal         sum,suma,sumr,gsum,gsuma,gsumr,diff;
5493:   PetscReal         tol,tola,tolr;
5494:   PetscReal         err_loc[6],err_glb[6];

5506:   if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector");

5508:   VecGetSize(U,&N);
5509:   VecGetLocalSize(U,&n);
5510:   VecGetOwnershipRange(U,&rstart,NULL);
5511:   VecGetArrayRead(U,&u);
5512:   VecGetArrayRead(Y,&y);
5513:   sum  = 0.; n_loc  = 0;
5514:   suma = 0.; na_loc = 0;
5515:   sumr = 0.; nr_loc = 0;
5516:   if (ts->vatol && ts->vrtol) {
5517:     const PetscScalar *atol,*rtol;
5518:     VecGetArrayRead(ts->vatol,&atol);
5519:     VecGetArrayRead(ts->vrtol,&rtol);
5520:     for (i=0; i<n; i++) {
5521:       diff = PetscAbsScalar(y[i] - u[i]);
5522:       tola = PetscRealPart(atol[i]);
5523:       if(tola>0.){
5524:         suma  += PetscSqr(diff/tola);
5525:         na_loc++;
5526:       }
5527:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5528:       if(tolr>0.){
5529:         sumr  += PetscSqr(diff/tolr);
5530:         nr_loc++;
5531:       }
5532:       tol=tola+tolr;
5533:       if(tol>0.){
5534:         sum  += PetscSqr(diff/tol);
5535:         n_loc++;
5536:       }
5537:     }
5538:     VecRestoreArrayRead(ts->vatol,&atol);
5539:     VecRestoreArrayRead(ts->vrtol,&rtol);
5540:   } else if (ts->vatol) {       /* vector atol, scalar rtol */
5541:     const PetscScalar *atol;
5542:     VecGetArrayRead(ts->vatol,&atol);
5543:     for (i=0; i<n; i++) {
5544:       diff = PetscAbsScalar(y[i] - u[i]);
5545:       tola = PetscRealPart(atol[i]);
5546:       if(tola>0.){
5547:         suma  += PetscSqr(diff/tola);
5548:         na_loc++;
5549:       }
5550:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5551:       if(tolr>0.){
5552:         sumr  += PetscSqr(diff/tolr);
5553:         nr_loc++;
5554:       }
5555:       tol=tola+tolr;
5556:       if(tol>0.){
5557:         sum  += PetscSqr(diff/tol);
5558:         n_loc++;
5559:       }
5560:     }
5561:     VecRestoreArrayRead(ts->vatol,&atol);
5562:   } else if (ts->vrtol) {       /* scalar atol, vector rtol */
5563:     const PetscScalar *rtol;
5564:     VecGetArrayRead(ts->vrtol,&rtol);
5565:     for (i=0; i<n; i++) {
5566:       diff = PetscAbsScalar(y[i] - u[i]);
5567:       tola = ts->atol;
5568:       if(tola>0.){
5569:         suma  += PetscSqr(diff/tola);
5570:         na_loc++;
5571:       }
5572:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5573:       if(tolr>0.){
5574:         sumr  += PetscSqr(diff/tolr);
5575:         nr_loc++;
5576:       }
5577:       tol=tola+tolr;
5578:       if(tol>0.){
5579:         sum  += PetscSqr(diff/tol);
5580:         n_loc++;
5581:       }
5582:     }
5583:     VecRestoreArrayRead(ts->vrtol,&rtol);
5584:   } else {                      /* scalar atol, scalar rtol */
5585:     for (i=0; i<n; i++) {
5586:       diff = PetscAbsScalar(y[i] - u[i]);
5587:      tola = ts->atol;
5588:       if(tola>0.){
5589:         suma  += PetscSqr(diff/tola);
5590:         na_loc++;
5591:       }
5592:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5593:       if(tolr>0.){
5594:         sumr  += PetscSqr(diff/tolr);
5595:         nr_loc++;
5596:       }
5597:       tol=tola+tolr;
5598:       if(tol>0.){
5599:         sum  += PetscSqr(diff/tol);
5600:         n_loc++;
5601:       }
5602:     }
5603:   }
5604:   VecRestoreArrayRead(U,&u);
5605:   VecRestoreArrayRead(Y,&y);

5607:   err_loc[0] = sum;
5608:   err_loc[1] = suma;
5609:   err_loc[2] = sumr;
5610:   err_loc[3] = (PetscReal)n_loc;
5611:   err_loc[4] = (PetscReal)na_loc;
5612:   err_loc[5] = (PetscReal)nr_loc;

5614:   MPIU_Allreduce(err_loc,err_glb,6,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));

5616:   gsum   = err_glb[0];
5617:   gsuma  = err_glb[1];
5618:   gsumr  = err_glb[2];
5619:   n_glb  = err_glb[3];
5620:   na_glb = err_glb[4];
5621:   nr_glb = err_glb[5];

5623:   *norm  = 0.;
5624:   if(n_glb>0. ){*norm  = PetscSqrtReal(gsum  / n_glb );}
5625:   *norma = 0.;
5626:   if(na_glb>0.){*norma = PetscSqrtReal(gsuma / na_glb);}
5627:   *normr = 0.;
5628:   if(nr_glb>0.){*normr = PetscSqrtReal(gsumr / nr_glb);}

5630:   if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
5631:   if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
5632:   if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
5633:   return(0);
5634: }

5636: /*@
5637:    TSErrorWeightedNormInfinity - compute a weighted infinity-norm of the difference between two state vectors

5639:    Collective on TS

5641:    Input Arguments:
5642: +  ts - time stepping context
5643: .  U - state vector, usually ts->vec_sol
5644: -  Y - state vector to be compared to U

5646:    Output Arguments:
5647: .  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5648: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5649: .  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5651:    Level: developer

5653: .seealso: TSErrorWeightedNorm(), TSErrorWeightedNorm2()
5654: @*/
5655: PetscErrorCode TSErrorWeightedNormInfinity(TS ts,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
5656: {
5657:   PetscErrorCode    ierr;
5658:   PetscInt          i,n,N,rstart;
5659:   const PetscScalar *u,*y;
5660:   PetscReal         max,gmax,maxa,gmaxa,maxr,gmaxr;
5661:   PetscReal         tol,tola,tolr,diff;
5662:   PetscReal         err_loc[3],err_glb[3];

5674:   if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"U and Y cannot be the same vector");

5676:   VecGetSize(U,&N);
5677:   VecGetLocalSize(U,&n);
5678:   VecGetOwnershipRange(U,&rstart,NULL);
5679:   VecGetArrayRead(U,&u);
5680:   VecGetArrayRead(Y,&y);

5682:   max=0.;
5683:   maxa=0.;
5684:   maxr=0.;

5686:   if (ts->vatol && ts->vrtol) {     /* vector atol, vector rtol */
5687:     const PetscScalar *atol,*rtol;
5688:     VecGetArrayRead(ts->vatol,&atol);
5689:     VecGetArrayRead(ts->vrtol,&rtol);

5691:     for (i=0; i<n; i++) {
5692:       diff = PetscAbsScalar(y[i] - u[i]);
5693:       tola = PetscRealPart(atol[i]);
5694:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5695:       tol  = tola+tolr;
5696:       if(tola>0.){
5697:         maxa = PetscMax(maxa,diff / tola);
5698:       }
5699:       if(tolr>0.){
5700:         maxr = PetscMax(maxr,diff / tolr);
5701:       }
5702:       if(tol>0.){
5703:         max = PetscMax(max,diff / tol);
5704:       }
5705:     }
5706:     VecRestoreArrayRead(ts->vatol,&atol);
5707:     VecRestoreArrayRead(ts->vrtol,&rtol);
5708:   } else if (ts->vatol) {       /* vector atol, scalar rtol */
5709:     const PetscScalar *atol;
5710:     VecGetArrayRead(ts->vatol,&atol);
5711:     for (i=0; i<n; i++) {
5712:       diff = PetscAbsScalar(y[i] - u[i]);
5713:       tola = PetscRealPart(atol[i]);
5714:       tolr = ts->rtol  * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5715:       tol  = tola+tolr;
5716:       if(tola>0.){
5717:         maxa = PetscMax(maxa,diff / tola);
5718:       }
5719:       if(tolr>0.){
5720:         maxr = PetscMax(maxr,diff / tolr);
5721:       }
5722:       if(tol>0.){
5723:         max = PetscMax(max,diff / tol);
5724:       }
5725:     }
5726:     VecRestoreArrayRead(ts->vatol,&atol);
5727:   } else if (ts->vrtol) {       /* scalar atol, vector rtol */
5728:     const PetscScalar *rtol;
5729:     VecGetArrayRead(ts->vrtol,&rtol);

5731:     for (i=0; i<n; i++) {
5732:       diff = PetscAbsScalar(y[i] - u[i]);
5733:       tola = ts->atol;
5734:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5735:       tol  = tola+tolr;
5736:       if(tola>0.){
5737:         maxa = PetscMax(maxa,diff / tola);
5738:       }
5739:       if(tolr>0.){
5740:         maxr = PetscMax(maxr,diff / tolr);
5741:       }
5742:       if(tol>0.){
5743:         max = PetscMax(max,diff / tol);
5744:       }
5745:     }
5746:     VecRestoreArrayRead(ts->vrtol,&rtol);
5747:   } else {                      /* scalar atol, scalar rtol */

5749:     for (i=0; i<n; i++) {
5750:       diff = PetscAbsScalar(y[i] - u[i]);
5751:       tola = ts->atol;
5752:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5753:       tol  = tola+tolr;
5754:       if(tola>0.){
5755:         maxa = PetscMax(maxa,diff / tola);
5756:       }
5757:       if(tolr>0.){
5758:         maxr = PetscMax(maxr,diff / tolr);
5759:       }
5760:       if(tol>0.){
5761:         max = PetscMax(max,diff / tol);
5762:       }
5763:     }
5764:   }
5765:   VecRestoreArrayRead(U,&u);
5766:   VecRestoreArrayRead(Y,&y);
5767:   err_loc[0] = max;
5768:   err_loc[1] = maxa;
5769:   err_loc[2] = maxr;
5770:   MPIU_Allreduce(err_loc,err_glb,3,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)ts));
5771:   gmax   = err_glb[0];
5772:   gmaxa  = err_glb[1];
5773:   gmaxr  = err_glb[2];

5775:   *norm = gmax;
5776:   *norma = gmaxa;
5777:   *normr = gmaxr;
5778:   if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
5779:     if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
5780:     if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
5781:   return(0);
5782: }

5784: /*@
5785:    TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances

5787:    Collective on TS

5789:    Input Arguments:
5790: +  ts - time stepping context
5791: .  U - state vector, usually ts->vec_sol
5792: .  Y - state vector to be compared to U
5793: -  wnormtype - norm type, either NORM_2 or NORM_INFINITY

5795:    Output Arguments:
5796: .  norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5797: .  norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5798: .  normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

5800:    Options Database Keys:
5801: .  -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5803:    Level: developer

5805: .seealso: TSErrorWeightedNormInfinity(), TSErrorWeightedNorm2(), TSErrorWeightedENorm
5806: @*/
5807: PetscErrorCode TSErrorWeightedNorm(TS ts,Vec U,Vec Y,NormType wnormtype,PetscReal *norm,PetscReal *norma,PetscReal *normr)
5808: {

5812:   if (wnormtype == NORM_2) {
5813:     TSErrorWeightedNorm2(ts,U,Y,norm,norma,normr);
5814:   } else if(wnormtype == NORM_INFINITY) {
5815:     TSErrorWeightedNormInfinity(ts,U,Y,norm,norma,normr);
5816:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
5817:   return(0);
5818: }


5821: /*@
5822:    TSErrorWeightedENorm2 - compute a weighted 2 error norm based on supplied absolute and relative tolerances

5824:    Collective on TS

5826:    Input Arguments:
5827: +  ts - time stepping context
5828: .  E - error vector
5829: .  U - state vector, usually ts->vec_sol
5830: -  Y - state vector, previous time step

5832:    Output Arguments:
5833: .  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5834: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5835: .  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5837:    Level: developer

5839: .seealso: TSErrorWeightedENorm(), TSErrorWeightedENormInfinity()
5840: @*/
5841: PetscErrorCode TSErrorWeightedENorm2(TS ts,Vec E,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
5842: {
5843:   PetscErrorCode    ierr;
5844:   PetscInt          i,n,N,rstart;
5845:   PetscInt          n_loc,na_loc,nr_loc;
5846:   PetscReal         n_glb,na_glb,nr_glb;
5847:   const PetscScalar *e,*u,*y;
5848:   PetscReal         err,sum,suma,sumr,gsum,gsuma,gsumr;
5849:   PetscReal         tol,tola,tolr;
5850:   PetscReal         err_loc[6],err_glb[6];


5866:   VecGetSize(E,&N);
5867:   VecGetLocalSize(E,&n);
5868:   VecGetOwnershipRange(E,&rstart,NULL);
5869:   VecGetArrayRead(E,&e);
5870:   VecGetArrayRead(U,&u);
5871:   VecGetArrayRead(Y,&y);
5872:   sum  = 0.; n_loc  = 0;
5873:   suma = 0.; na_loc = 0;
5874:   sumr = 0.; nr_loc = 0;
5875:   if (ts->vatol && ts->vrtol) {
5876:     const PetscScalar *atol,*rtol;
5877:     VecGetArrayRead(ts->vatol,&atol);
5878:     VecGetArrayRead(ts->vrtol,&rtol);
5879:     for (i=0; i<n; i++) {
5880:       err = PetscAbsScalar(e[i]);
5881:       tola = PetscRealPart(atol[i]);
5882:       if(tola>0.){
5883:         suma  += PetscSqr(err/tola);
5884:         na_loc++;
5885:       }
5886:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5887:       if(tolr>0.){
5888:         sumr  += PetscSqr(err/tolr);
5889:         nr_loc++;
5890:       }
5891:       tol=tola+tolr;
5892:       if(tol>0.){
5893:         sum  += PetscSqr(err/tol);
5894:         n_loc++;
5895:       }
5896:     }
5897:     VecRestoreArrayRead(ts->vatol,&atol);
5898:     VecRestoreArrayRead(ts->vrtol,&rtol);
5899:   } else if (ts->vatol) {       /* vector atol, scalar rtol */
5900:     const PetscScalar *atol;
5901:     VecGetArrayRead(ts->vatol,&atol);
5902:     for (i=0; i<n; i++) {
5903:       err = PetscAbsScalar(e[i]);
5904:       tola = PetscRealPart(atol[i]);
5905:       if(tola>0.){
5906:         suma  += PetscSqr(err/tola);
5907:         na_loc++;
5908:       }
5909:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5910:       if(tolr>0.){
5911:         sumr  += PetscSqr(err/tolr);
5912:         nr_loc++;
5913:       }
5914:       tol=tola+tolr;
5915:       if(tol>0.){
5916:         sum  += PetscSqr(err/tol);
5917:         n_loc++;
5918:       }
5919:     }
5920:     VecRestoreArrayRead(ts->vatol,&atol);
5921:   } else if (ts->vrtol) {       /* scalar atol, vector rtol */
5922:     const PetscScalar *rtol;
5923:     VecGetArrayRead(ts->vrtol,&rtol);
5924:     for (i=0; i<n; i++) {
5925:       err = PetscAbsScalar(e[i]);
5926:       tola = ts->atol;
5927:       if(tola>0.){
5928:         suma  += PetscSqr(err/tola);
5929:         na_loc++;
5930:       }
5931:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5932:       if(tolr>0.){
5933:         sumr  += PetscSqr(err/tolr);
5934:         nr_loc++;
5935:       }
5936:       tol=tola+tolr;
5937:       if(tol>0.){
5938:         sum  += PetscSqr(err/tol);
5939:         n_loc++;
5940:       }
5941:     }
5942:     VecRestoreArrayRead(ts->vrtol,&rtol);
5943:   } else {                      /* scalar atol, scalar rtol */
5944:     for (i=0; i<n; i++) {
5945:       err = PetscAbsScalar(e[i]);
5946:      tola = ts->atol;
5947:       if(tola>0.){
5948:         suma  += PetscSqr(err/tola);
5949:         na_loc++;
5950:       }
5951:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
5952:       if(tolr>0.){
5953:         sumr  += PetscSqr(err/tolr);
5954:         nr_loc++;
5955:       }
5956:       tol=tola+tolr;
5957:       if(tol>0.){
5958:         sum  += PetscSqr(err/tol);
5959:         n_loc++;
5960:       }
5961:     }
5962:   }
5963:   VecRestoreArrayRead(E,&e);
5964:   VecRestoreArrayRead(U,&u);
5965:   VecRestoreArrayRead(Y,&y);

5967:   err_loc[0] = sum;
5968:   err_loc[1] = suma;
5969:   err_loc[2] = sumr;
5970:   err_loc[3] = (PetscReal)n_loc;
5971:   err_loc[4] = (PetscReal)na_loc;
5972:   err_loc[5] = (PetscReal)nr_loc;

5974:   MPIU_Allreduce(err_loc,err_glb,6,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));

5976:   gsum   = err_glb[0];
5977:   gsuma  = err_glb[1];
5978:   gsumr  = err_glb[2];
5979:   n_glb  = err_glb[3];
5980:   na_glb = err_glb[4];
5981:   nr_glb = err_glb[5];

5983:   *norm  = 0.;
5984:   if(n_glb>0. ){*norm  = PetscSqrtReal(gsum  / n_glb );}
5985:   *norma = 0.;
5986:   if(na_glb>0.){*norma = PetscSqrtReal(gsuma / na_glb);}
5987:   *normr = 0.;
5988:   if(nr_glb>0.){*normr = PetscSqrtReal(gsumr / nr_glb);}

5990:   if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
5991:   if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
5992:   if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
5993:   return(0);
5994: }

5996: /*@
5997:    TSErrorWeightedENormInfinity - compute a weighted infinity error norm based on supplied absolute and relative tolerances
5998:    Collective on TS

6000:    Input Arguments:
6001: +  ts - time stepping context
6002: .  E - error vector
6003: .  U - state vector, usually ts->vec_sol
6004: -  Y - state vector, previous time step

6006:    Output Arguments:
6007: .  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
6008: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
6009: .  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

6011:    Level: developer

6013: .seealso: TSErrorWeightedENorm(), TSErrorWeightedENorm2()
6014: @*/
6015: PetscErrorCode TSErrorWeightedENormInfinity(TS ts,Vec E,Vec U,Vec Y,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6016: {
6017:   PetscErrorCode    ierr;
6018:   PetscInt          i,n,N,rstart;
6019:   const PetscScalar *e,*u,*y;
6020:   PetscReal         err,max,gmax,maxa,gmaxa,maxr,gmaxr;
6021:   PetscReal         tol,tola,tolr;
6022:   PetscReal         err_loc[3],err_glb[3];


6038:   VecGetSize(E,&N);
6039:   VecGetLocalSize(E,&n);
6040:   VecGetOwnershipRange(E,&rstart,NULL);
6041:   VecGetArrayRead(E,&e);
6042:   VecGetArrayRead(U,&u);
6043:   VecGetArrayRead(Y,&y);

6045:   max=0.;
6046:   maxa=0.;
6047:   maxr=0.;

6049:   if (ts->vatol && ts->vrtol) {     /* vector atol, vector rtol */
6050:     const PetscScalar *atol,*rtol;
6051:     VecGetArrayRead(ts->vatol,&atol);
6052:     VecGetArrayRead(ts->vrtol,&rtol);

6054:     for (i=0; i<n; i++) {
6055:       err = PetscAbsScalar(e[i]);
6056:       tola = PetscRealPart(atol[i]);
6057:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6058:       tol  = tola+tolr;
6059:       if(tola>0.){
6060:         maxa = PetscMax(maxa,err / tola);
6061:       }
6062:       if(tolr>0.){
6063:         maxr = PetscMax(maxr,err / tolr);
6064:       }
6065:       if(tol>0.){
6066:         max = PetscMax(max,err / tol);
6067:       }
6068:     }
6069:     VecRestoreArrayRead(ts->vatol,&atol);
6070:     VecRestoreArrayRead(ts->vrtol,&rtol);
6071:   } else if (ts->vatol) {       /* vector atol, scalar rtol */
6072:     const PetscScalar *atol;
6073:     VecGetArrayRead(ts->vatol,&atol);
6074:     for (i=0; i<n; i++) {
6075:       err = PetscAbsScalar(e[i]);
6076:       tola = PetscRealPart(atol[i]);
6077:       tolr = ts->rtol  * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6078:       tol  = tola+tolr;
6079:       if(tola>0.){
6080:         maxa = PetscMax(maxa,err / tola);
6081:       }
6082:       if(tolr>0.){
6083:         maxr = PetscMax(maxr,err / tolr);
6084:       }
6085:       if(tol>0.){
6086:         max = PetscMax(max,err / tol);
6087:       }
6088:     }
6089:     VecRestoreArrayRead(ts->vatol,&atol);
6090:   } else if (ts->vrtol) {       /* scalar atol, vector rtol */
6091:     const PetscScalar *rtol;
6092:     VecGetArrayRead(ts->vrtol,&rtol);

6094:     for (i=0; i<n; i++) {
6095:       err = PetscAbsScalar(e[i]);
6096:       tola = ts->atol;
6097:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6098:       tol  = tola+tolr;
6099:       if(tola>0.){
6100:         maxa = PetscMax(maxa,err / tola);
6101:       }
6102:       if(tolr>0.){
6103:         maxr = PetscMax(maxr,err / tolr);
6104:       }
6105:       if(tol>0.){
6106:         max = PetscMax(max,err / tol);
6107:       }
6108:     }
6109:     VecRestoreArrayRead(ts->vrtol,&rtol);
6110:   } else {                      /* scalar atol, scalar rtol */

6112:     for (i=0; i<n; i++) {
6113:       err = PetscAbsScalar(e[i]);
6114:       tola = ts->atol;
6115:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i]));
6116:       tol  = tola+tolr;
6117:       if(tola>0.){
6118:         maxa = PetscMax(maxa,err / tola);
6119:       }
6120:       if(tolr>0.){
6121:         maxr = PetscMax(maxr,err / tolr);
6122:       }
6123:       if(tol>0.){
6124:         max = PetscMax(max,err / tol);
6125:       }
6126:     }
6127:   }
6128:   VecRestoreArrayRead(E,&e);
6129:   VecRestoreArrayRead(U,&u);
6130:   VecRestoreArrayRead(Y,&y);
6131:   err_loc[0] = max;
6132:   err_loc[1] = maxa;
6133:   err_loc[2] = maxr;
6134:   MPIU_Allreduce(err_loc,err_glb,3,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)ts));
6135:   gmax   = err_glb[0];
6136:   gmaxa  = err_glb[1];
6137:   gmaxr  = err_glb[2];

6139:   *norm = gmax;
6140:   *norma = gmaxa;
6141:   *normr = gmaxr;
6142:   if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
6143:     if (PetscIsInfOrNanScalar(*norma)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norma");
6144:     if (PetscIsInfOrNanScalar(*normr)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in normr");
6145:   return(0);
6146: }

6148: /*@
6149:    TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances

6151:    Collective on TS

6153:    Input Arguments:
6154: +  ts - time stepping context
6155: .  E - error vector
6156: .  U - state vector, usually ts->vec_sol
6157: .  Y - state vector, previous time step
6158: -  wnormtype - norm type, either NORM_2 or NORM_INFINITY

6160:    Output Arguments:
6161: .  norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
6162: .  norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
6163: .  normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

6165:    Options Database Keys:
6166: .  -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

6168:    Level: developer

6170: .seealso: TSErrorWeightedENormInfinity(), TSErrorWeightedENorm2(), TSErrorWeightedNormInfinity(), TSErrorWeightedNorm2()
6171: @*/
6172: PetscErrorCode TSErrorWeightedENorm(TS ts,Vec E,Vec U,Vec Y,NormType wnormtype,PetscReal *norm,PetscReal *norma,PetscReal *normr)
6173: {

6177:   if (wnormtype == NORM_2) {
6178:     TSErrorWeightedENorm2(ts,E,U,Y,norm,norma,normr);
6179:   } else if(wnormtype == NORM_INFINITY) {
6180:     TSErrorWeightedENormInfinity(ts,E,U,Y,norm,norma,normr);
6181:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for norm type %s",NormTypes[wnormtype]);
6182:   return(0);
6183: }


6186: /*@
6187:    TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler

6189:    Logically Collective on TS

6191:    Input Arguments:
6192: +  ts - time stepping context
6193: -  cfltime - maximum stable time step if using forward Euler (value can be different on each process)

6195:    Note:
6196:    After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()

6198:    Level: intermediate

6200: .seealso: TSGetCFLTime(), TSADAPTCFL
6201: @*/
6202: PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime)
6203: {
6206:   ts->cfltime_local = cfltime;
6207:   ts->cfltime       = -1.;
6208:   return(0);
6209: }

6211: /*@
6212:    TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler

6214:    Collective on TS

6216:    Input Arguments:
6217: .  ts - time stepping context

6219:    Output Arguments:
6220: .  cfltime - maximum stable time step for forward Euler

6222:    Level: advanced

6224: .seealso: TSSetCFLTimeLocal()
6225: @*/
6226: PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime)
6227: {

6231:   if (ts->cfltime < 0) {
6232:     MPIU_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));
6233:   }
6234:   *cfltime = ts->cfltime;
6235:   return(0);
6236: }

6238: /*@
6239:    TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu

6241:    Input Parameters:
6242: .  ts   - the TS context.
6243: .  xl   - lower bound.
6244: .  xu   - upper bound.

6246:    Notes:
6247:    If this routine is not called then the lower and upper bounds are set to
6248:    PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().

6250:    Level: advanced

6252: @*/
6253: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
6254: {
6256:   SNES           snes;

6259:   TSGetSNES(ts,&snes);
6260:   SNESVISetVariableBounds(snes,xl,xu);
6261:   return(0);
6262: }

6264: #if defined(PETSC_HAVE_MATLAB_ENGINE)
6265: #include <mex.h>

6267: typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext;

6269: /*
6270:    TSComputeFunction_Matlab - Calls the function that has been set with
6271:                          TSSetFunctionMatlab().

6273:    Collective on TS

6275:    Input Parameters:
6276: +  snes - the TS context
6277: -  u - input vector

6279:    Output Parameter:
6280: .  y - function vector, as set by TSSetFunction()

6282:    Notes:
6283:    TSComputeFunction() is typically used within nonlinear solvers
6284:    implementations, so most users would not generally call this routine
6285:    themselves.

6287:    Level: developer

6289: .keywords: TS, nonlinear, compute, function

6291: .seealso: TSSetFunction(), TSGetFunction()
6292: */
6293: PetscErrorCode  TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx)
6294: {
6295:   PetscErrorCode  ierr;
6296:   TSMatlabContext *sctx = (TSMatlabContext*)ctx;
6297:   int             nlhs  = 1,nrhs = 7;
6298:   mxArray         *plhs[1],*prhs[7];
6299:   long long int   lx = 0,lxdot = 0,ly = 0,ls = 0;


6309:   PetscMemcpy(&ls,&snes,sizeof(snes));
6310:   PetscMemcpy(&lx,&u,sizeof(u));
6311:   PetscMemcpy(&lxdot,&udot,sizeof(udot));
6312:   PetscMemcpy(&ly,&y,sizeof(u));

6314:   prhs[0] =  mxCreateDoubleScalar((double)ls);
6315:   prhs[1] =  mxCreateDoubleScalar(time);
6316:   prhs[2] =  mxCreateDoubleScalar((double)lx);
6317:   prhs[3] =  mxCreateDoubleScalar((double)lxdot);
6318:   prhs[4] =  mxCreateDoubleScalar((double)ly);
6319:   prhs[5] =  mxCreateString(sctx->funcname);
6320:   prhs[6] =  sctx->ctx;
6321:    mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");
6322:    mxGetScalar(plhs[0]);
6323:   mxDestroyArray(prhs[0]);
6324:   mxDestroyArray(prhs[1]);
6325:   mxDestroyArray(prhs[2]);
6326:   mxDestroyArray(prhs[3]);
6327:   mxDestroyArray(prhs[4]);
6328:   mxDestroyArray(prhs[5]);
6329:   mxDestroyArray(plhs[0]);
6330:   return(0);
6331: }

6333: /*
6334:    TSSetFunctionMatlab - Sets the function evaluation routine and function
6335:    vector for use by the TS routines in solving ODEs
6336:    equations from MATLAB. Here the function is a string containing the name of a MATLAB function

6338:    Logically Collective on TS

6340:    Input Parameters:
6341: +  ts - the TS context
6342: -  func - function evaluation routine

6344:    Calling sequence of func:
6345: $    func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx);

6347:    Level: beginner

6349: .keywords: TS, nonlinear, set, function

6351: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
6352: */
6353: PetscErrorCode  TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx)
6354: {
6355:   PetscErrorCode  ierr;
6356:   TSMatlabContext *sctx;

6359:   /* currently sctx is memory bleed */
6360:   PetscNew(&sctx);
6361:   PetscStrallocpy(func,&sctx->funcname);
6362:   /*
6363:      This should work, but it doesn't
6364:   sctx->ctx = ctx;
6365:   mexMakeArrayPersistent(sctx->ctx);
6366:   */
6367:   sctx->ctx = mxDuplicateArray(ctx);

6369:   TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);
6370:   return(0);
6371: }

6373: /*
6374:    TSComputeJacobian_Matlab - Calls the function that has been set with
6375:                          TSSetJacobianMatlab().

6377:    Collective on TS

6379:    Input Parameters:
6380: +  ts - the TS context
6381: .  u - input vector
6382: .  A, B - the matrices
6383: -  ctx - user context

6385:    Level: developer

6387: .keywords: TS, nonlinear, compute, function

6389: .seealso: TSSetFunction(), TSGetFunction()
6390: @*/
6391: PetscErrorCode  TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat A,Mat B,void *ctx)
6392: {
6393:   PetscErrorCode  ierr;
6394:   TSMatlabContext *sctx = (TSMatlabContext*)ctx;
6395:   int             nlhs  = 2,nrhs = 9;
6396:   mxArray         *plhs[2],*prhs[9];
6397:   long long int   lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0;


6403:   /* call Matlab function in ctx with arguments u and y */

6405:   PetscMemcpy(&ls,&ts,sizeof(ts));
6406:   PetscMemcpy(&lx,&u,sizeof(u));
6407:   PetscMemcpy(&lxdot,&udot,sizeof(u));
6408:   PetscMemcpy(&lA,A,sizeof(u));
6409:   PetscMemcpy(&lB,B,sizeof(u));

6411:   prhs[0] =  mxCreateDoubleScalar((double)ls);
6412:   prhs[1] =  mxCreateDoubleScalar((double)time);
6413:   prhs[2] =  mxCreateDoubleScalar((double)lx);
6414:   prhs[3] =  mxCreateDoubleScalar((double)lxdot);
6415:   prhs[4] =  mxCreateDoubleScalar((double)shift);
6416:   prhs[5] =  mxCreateDoubleScalar((double)lA);
6417:   prhs[6] =  mxCreateDoubleScalar((double)lB);
6418:   prhs[7] =  mxCreateString(sctx->funcname);
6419:   prhs[8] =  sctx->ctx;
6420:    mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");
6421:    mxGetScalar(plhs[0]);
6422:   mxDestroyArray(prhs[0]);
6423:   mxDestroyArray(prhs[1]);
6424:   mxDestroyArray(prhs[2]);
6425:   mxDestroyArray(prhs[3]);
6426:   mxDestroyArray(prhs[4]);
6427:   mxDestroyArray(prhs[5]);
6428:   mxDestroyArray(prhs[6]);
6429:   mxDestroyArray(prhs[7]);
6430:   mxDestroyArray(plhs[0]);
6431:   mxDestroyArray(plhs[1]);
6432:   return(0);
6433: }

6435: /*
6436:    TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices
6437:    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

6439:    Logically Collective on TS

6441:    Input Parameters:
6442: +  ts - the TS context
6443: .  A,B - Jacobian matrices
6444: .  func - function evaluation routine
6445: -  ctx - user context

6447:    Calling sequence of func:
6448: $    flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx);

6450:    Level: developer

6452: .keywords: TS, nonlinear, set, function

6454: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
6455: */
6456: PetscErrorCode  TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx)
6457: {
6458:   PetscErrorCode  ierr;
6459:   TSMatlabContext *sctx;

6462:   /* currently sctx is memory bleed */
6463:   PetscNew(&sctx);
6464:   PetscStrallocpy(func,&sctx->funcname);
6465:   /*
6466:      This should work, but it doesn't
6467:   sctx->ctx = ctx;
6468:   mexMakeArrayPersistent(sctx->ctx);
6469:   */
6470:   sctx->ctx = mxDuplicateArray(ctx);

6472:   TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);
6473:   return(0);
6474: }

6476: /*
6477:    TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab().

6479:    Collective on TS

6481: .seealso: TSSetFunction(), TSGetFunction()
6482: @*/
6483: PetscErrorCode  TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx)
6484: {
6485:   PetscErrorCode  ierr;
6486:   TSMatlabContext *sctx = (TSMatlabContext*)ctx;
6487:   int             nlhs  = 1,nrhs = 6;
6488:   mxArray         *plhs[1],*prhs[6];
6489:   long long int   lx = 0,ls = 0;


6495:   PetscMemcpy(&ls,&ts,sizeof(ts));
6496:   PetscMemcpy(&lx,&u,sizeof(u));

6498:   prhs[0] =  mxCreateDoubleScalar((double)ls);
6499:   prhs[1] =  mxCreateDoubleScalar((double)it);
6500:   prhs[2] =  mxCreateDoubleScalar((double)time);
6501:   prhs[3] =  mxCreateDoubleScalar((double)lx);
6502:   prhs[4] =  mxCreateString(sctx->funcname);
6503:   prhs[5] =  sctx->ctx;
6504:    mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");
6505:    mxGetScalar(plhs[0]);
6506:   mxDestroyArray(prhs[0]);
6507:   mxDestroyArray(prhs[1]);
6508:   mxDestroyArray(prhs[2]);
6509:   mxDestroyArray(prhs[3]);
6510:   mxDestroyArray(prhs[4]);
6511:   mxDestroyArray(plhs[0]);
6512:   return(0);
6513: }

6515: /*
6516:    TSMonitorSetMatlab - Sets the monitor function from Matlab

6518:    Level: developer

6520: .keywords: TS, nonlinear, set, function

6522: .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction()
6523: */
6524: PetscErrorCode  TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx)
6525: {
6526:   PetscErrorCode  ierr;
6527:   TSMatlabContext *sctx;

6530:   /* currently sctx is memory bleed */
6531:   PetscNew(&sctx);
6532:   PetscStrallocpy(func,&sctx->funcname);
6533:   /*
6534:      This should work, but it doesn't
6535:   sctx->ctx = ctx;
6536:   mexMakeArrayPersistent(sctx->ctx);
6537:   */
6538:   sctx->ctx = mxDuplicateArray(ctx);

6540:   TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);
6541:   return(0);
6542: }
6543: #endif

6545: /*@C
6546:    TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector
6547:        in a time based line graph

6549:    Collective on TS

6551:    Input Parameters:
6552: +  ts - the TS context
6553: .  step - current time-step
6554: .  ptime - current time
6555: .  u - current solution
6556: -  dctx - the TSMonitorLGCtx object that contains all the options for the monitoring, this is created with TSMonitorLGCtxCreate()

6558:    Options Database:
6559: .   -ts_monitor_lg_solution_variables

6561:    Level: intermediate

6563:    Notes:
6564:     Each process in a parallel run displays its component solutions in a separate window

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

6568: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGCtxCreate(), TSMonitorLGCtxSetVariableNames(), TSMonitorLGCtxGetVariableNames(),
6569:            TSMonitorLGSetVariableNames(), TSMonitorLGGetVariableNames(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetDisplayVariables(),
6570:            TSMonitorLGCtxSetTransform(), TSMonitorLGSetTransform(), TSMonitorLGError(), TSMonitorLGSNESIterations(), TSMonitorLGKSPIterations(),
6571:            TSMonitorEnvelopeCtxCreate(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxDestroy(), TSMonitorEnvelop()
6572: @*/
6573: PetscErrorCode  TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
6574: {
6575:   PetscErrorCode    ierr;
6576:   TSMonitorLGCtx    ctx = (TSMonitorLGCtx)dctx;
6577:   const PetscScalar *yy;
6578:   Vec               v;

6581:   if (step < 0) return(0); /* -1 indicates interpolated solution */
6582:   if (!step) {
6583:     PetscDrawAxis axis;
6584:     PetscInt      dim;
6585:     PetscDrawLGGetAxis(ctx->lg,&axis);
6586:     PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");
6587:     if (!ctx->names) {
6588:       PetscBool flg;
6589:       /* user provides names of variables to plot but no names has been set so assume names are integer values */
6590:       PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",&flg);
6591:       if (flg) {
6592:         PetscInt i,n;
6593:         char     **names;
6594:         VecGetSize(u,&n);
6595:         PetscMalloc1(n+1,&names);
6596:         for (i=0; i<n; i++) {
6597:           PetscMalloc1(5,&names[i]);
6598:           PetscSNPrintf(names[i],5,"%D",i);
6599:         }
6600:         names[n] = NULL;
6601:         ctx->names = names;
6602:       }
6603:     }
6604:     if (ctx->names && !ctx->displaynames) {
6605:       char      **displaynames;
6606:       PetscBool flg;
6607:       VecGetLocalSize(u,&dim);
6608:       PetscMalloc1(dim+1,&displaynames);
6609:       PetscMemzero(displaynames,(dim+1)*sizeof(char*));
6610:       PetscOptionsGetStringArray(((PetscObject)ts)->options,((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",displaynames,&dim,&flg);
6611:       if (flg) {
6612:         TSMonitorLGCtxSetDisplayVariables(ctx,(const char *const *)displaynames);
6613:       }
6614:       PetscStrArrayDestroy(&displaynames);
6615:     }
6616:     if (ctx->displaynames) {
6617:       PetscDrawLGSetDimension(ctx->lg,ctx->ndisplayvariables);
6618:       PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->displaynames);
6619:     } else if (ctx->names) {
6620:       VecGetLocalSize(u,&dim);
6621:       PetscDrawLGSetDimension(ctx->lg,dim);
6622:       PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->names);
6623:     } else {
6624:       VecGetLocalSize(u,&dim);
6625:       PetscDrawLGSetDimension(ctx->lg,dim);
6626:     }
6627:     PetscDrawLGReset(ctx->lg);
6628:   }

6630:   if (!ctx->transform) v = u;
6631:   else {(*ctx->transform)(ctx->transformctx,u,&v);}
6632:   VecGetArrayRead(v,&yy);
6633:   if (ctx->displaynames) {
6634:     PetscInt i;
6635:     for (i=0; i<ctx->ndisplayvariables; i++)
6636:       ctx->displayvalues[i] = PetscRealPart(yy[ctx->displayvariables[i]]);
6637:     PetscDrawLGAddCommonPoint(ctx->lg,ptime,ctx->displayvalues);
6638:   } else {
6639: #if defined(PETSC_USE_COMPLEX)
6640:     PetscInt  i,n;
6641:     PetscReal *yreal;
6642:     VecGetLocalSize(v,&n);
6643:     PetscMalloc1(n,&yreal);
6644:     for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
6645:     PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
6646:     PetscFree(yreal);
6647: #else
6648:     PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
6649: #endif
6650:   }
6651:   VecRestoreArrayRead(v,&yy);
6652:   if (ctx->transform) {VecDestroy(&v);}

6654:   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
6655:     PetscDrawLGDraw(ctx->lg);
6656:     PetscDrawLGSave(ctx->lg);
6657:   }
6658:   return(0);
6659: }

6661: /*@C
6662:    TSMonitorLGSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot

6664:    Collective on TS

6666:    Input Parameters:
6667: +  ts - the TS context
6668: -  names - the names of the components, final string must be NULL

6670:    Level: intermediate

6672:    Notes:
6673:     If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored

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

6677: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetVariableNames()
6678: @*/
6679: PetscErrorCode  TSMonitorLGSetVariableNames(TS ts,const char * const *names)
6680: {
6681:   PetscErrorCode    ierr;
6682:   PetscInt          i;

6685:   for (i=0; i<ts->numbermonitors; i++) {
6686:     if (ts->monitor[i] == TSMonitorLGSolution) {
6687:       TSMonitorLGCtxSetVariableNames((TSMonitorLGCtx)ts->monitorcontext[i],names);
6688:       break;
6689:     }
6690:   }
6691:   return(0);
6692: }

6694: /*@C
6695:    TSMonitorLGCtxSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot

6697:    Collective on TS

6699:    Input Parameters:
6700: +  ts - the TS context
6701: -  names - the names of the components, final string must be NULL

6703:    Level: intermediate

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

6707: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGSetVariableNames()
6708: @*/
6709: PetscErrorCode  TSMonitorLGCtxSetVariableNames(TSMonitorLGCtx ctx,const char * const *names)
6710: {
6711:   PetscErrorCode    ierr;

6714:   PetscStrArrayDestroy(&ctx->names);
6715:   PetscStrArrayallocpy(names,&ctx->names);
6716:   return(0);
6717: }

6719: /*@C
6720:    TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot

6722:    Collective on TS

6724:    Input Parameter:
6725: .  ts - the TS context

6727:    Output Parameter:
6728: .  names - the names of the components, final string must be NULL

6730:    Level: intermediate

6732:    Notes:
6733:     If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored

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

6737: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
6738: @*/
6739: PetscErrorCode  TSMonitorLGGetVariableNames(TS ts,const char *const **names)
6740: {
6741:   PetscInt       i;

6744:   *names = NULL;
6745:   for (i=0; i<ts->numbermonitors; i++) {
6746:     if (ts->monitor[i] == TSMonitorLGSolution) {
6747:       TSMonitorLGCtx  ctx = (TSMonitorLGCtx) ts->monitorcontext[i];
6748:       *names = (const char *const *)ctx->names;
6749:       break;
6750:     }
6751:   }
6752:   return(0);
6753: }

6755: /*@C
6756:    TSMonitorLGCtxSetDisplayVariables - Sets the variables that are to be display in the monitor

6758:    Collective on TS

6760:    Input Parameters:
6761: +  ctx - the TSMonitorLG context
6762: .  displaynames - the names of the components, final string must be NULL

6764:    Level: intermediate

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

6768: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
6769: @*/
6770: PetscErrorCode  TSMonitorLGCtxSetDisplayVariables(TSMonitorLGCtx ctx,const char * const *displaynames)
6771: {
6772:   PetscInt          j = 0,k;
6773:   PetscErrorCode    ierr;

6776:   if (!ctx->names) return(0);
6777:   PetscStrArrayDestroy(&ctx->displaynames);
6778:   PetscStrArrayallocpy(displaynames,&ctx->displaynames);
6779:   while (displaynames[j]) j++;
6780:   ctx->ndisplayvariables = j;
6781:   PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvariables);
6782:   PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvalues);
6783:   j = 0;
6784:   while (displaynames[j]) {
6785:     k = 0;
6786:     while (ctx->names[k]) {
6787:       PetscBool flg;
6788:       PetscStrcmp(displaynames[j],ctx->names[k],&flg);
6789:       if (flg) {
6790:         ctx->displayvariables[j] = k;
6791:         break;
6792:       }
6793:       k++;
6794:     }
6795:     j++;
6796:   }
6797:   return(0);
6798: }

6800: /*@C
6801:    TSMonitorLGSetDisplayVariables - Sets the variables that are to be display in the monitor

6803:    Collective on TS

6805:    Input Parameters:
6806: +  ts - the TS context
6807: .  displaynames - the names of the components, final string must be NULL

6809:    Notes:
6810:     If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored

6812:    Level: intermediate

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

6816: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames()
6817: @*/
6818: PetscErrorCode  TSMonitorLGSetDisplayVariables(TS ts,const char * const *displaynames)
6819: {
6820:   PetscInt          i;
6821:   PetscErrorCode    ierr;

6824:   for (i=0; i<ts->numbermonitors; i++) {
6825:     if (ts->monitor[i] == TSMonitorLGSolution) {
6826:       TSMonitorLGCtxSetDisplayVariables((TSMonitorLGCtx)ts->monitorcontext[i],displaynames);
6827:       break;
6828:     }
6829:   }
6830:   return(0);
6831: }

6833: /*@C
6834:    TSMonitorLGSetTransform - Solution vector will be transformed by provided function before being displayed

6836:    Collective on TS

6838:    Input Parameters:
6839: +  ts - the TS context
6840: .  transform - the transform function
6841: .  destroy - function to destroy the optional context
6842: -  ctx - optional context used by transform function

6844:    Notes:
6845:     If the TS object does not have a TSMonitorLGCtx associated with it then this function is ignored

6847:    Level: intermediate

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

6851: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGCtxSetTransform()
6852: @*/
6853: PetscErrorCode  TSMonitorLGSetTransform(TS ts,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
6854: {
6855:   PetscInt          i;
6856:   PetscErrorCode    ierr;

6859:   for (i=0; i<ts->numbermonitors; i++) {
6860:     if (ts->monitor[i] == TSMonitorLGSolution) {
6861:       TSMonitorLGCtxSetTransform((TSMonitorLGCtx)ts->monitorcontext[i],transform,destroy,tctx);
6862:     }
6863:   }
6864:   return(0);
6865: }

6867: /*@C
6868:    TSMonitorLGCtxSetTransform - Solution vector will be transformed by provided function before being displayed

6870:    Collective on TSLGCtx

6872:    Input Parameters:
6873: +  ts - the TS context
6874: .  transform - the transform function
6875: .  destroy - function to destroy the optional context
6876: -  ctx - optional context used by transform function

6878:    Level: intermediate

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

6882: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGSetTransform()
6883: @*/
6884: PetscErrorCode  TSMonitorLGCtxSetTransform(TSMonitorLGCtx ctx,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx)
6885: {
6887:   ctx->transform    = transform;
6888:   ctx->transformdestroy = destroy;
6889:   ctx->transformctx = tctx;
6890:   return(0);
6891: }

6893: /*@C
6894:    TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the error
6895:        in a time based line graph

6897:    Collective on TS

6899:    Input Parameters:
6900: +  ts - the TS context
6901: .  step - current time-step
6902: .  ptime - current time
6903: .  u - current solution
6904: -  dctx - TSMonitorLGCtx object created with TSMonitorLGCtxCreate()

6906:    Level: intermediate

6908:    Notes:
6909:     Each process in a parallel run displays its component errors in a separate window

6911:    The user must provide the solution using TSSetSolutionFunction() to use this monitor.

6913:    Options Database Keys:
6914: .  -ts_monitor_lg_error - create a graphical monitor of error history

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

6918: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
6919: @*/
6920: PetscErrorCode  TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy)
6921: {
6922:   PetscErrorCode    ierr;
6923:   TSMonitorLGCtx    ctx = (TSMonitorLGCtx)dummy;
6924:   const PetscScalar *yy;
6925:   Vec               y;

6928:   if (!step) {
6929:     PetscDrawAxis axis;
6930:     PetscInt      dim;
6931:     PetscDrawLGGetAxis(ctx->lg,&axis);
6932:     PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Error");
6933:     VecGetLocalSize(u,&dim);
6934:     PetscDrawLGSetDimension(ctx->lg,dim);
6935:     PetscDrawLGReset(ctx->lg);
6936:   }
6937:   VecDuplicate(u,&y);
6938:   TSComputeSolutionFunction(ts,ptime,y);
6939:   VecAXPY(y,-1.0,u);
6940:   VecGetArrayRead(y,&yy);
6941: #if defined(PETSC_USE_COMPLEX)
6942:   {
6943:     PetscReal *yreal;
6944:     PetscInt  i,n;
6945:     VecGetLocalSize(y,&n);
6946:     PetscMalloc1(n,&yreal);
6947:     for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]);
6948:     PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);
6949:     PetscFree(yreal);
6950:   }
6951: #else
6952:   PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);
6953: #endif
6954:   VecRestoreArrayRead(y,&yy);
6955:   VecDestroy(&y);
6956:   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
6957:     PetscDrawLGDraw(ctx->lg);
6958:     PetscDrawLGSave(ctx->lg);
6959:   }
6960:   return(0);
6961: }

6963: /*@C
6964:    TSMonitorSPSwarmSolution - Graphically displays phase plots of DMSwarm particles on a scatter plot

6966:    Input Parameters:
6967: +  ts - the TS context
6968: .  step - current time-step
6969: .  ptime - current time
6970: .  u - current solution
6971: -  dctx - the TSMonitorSPCtx object that contains all the options for the monitoring, this is created with TSMonitorSPCtxCreate()

6973:    Options Database:
6974: .   -ts_monitor_sp_swarm

6976:    Level: intermediate

6978: .keywords: TS,  vector, monitor, view, swarm
6979: @*/
6980: PetscErrorCode TSMonitorSPSwarmSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
6981: {
6982:   PetscErrorCode    ierr;
6983:   TSMonitorSPCtx    ctx = (TSMonitorSPCtx)dctx;
6984:   const PetscScalar *yy;
6985:   PetscReal       *y,*x;
6986:   PetscInt          Np, p, dim=2;
6987:   DM                dm;

6990: 
6991:   if (step < 0) return(0); /* -1 indicates interpolated solution */
6992:   if (!step) {
6993:     PetscDrawAxis axis;
6994:     PetscDrawSPGetAxis(ctx->sp,&axis);
6995:     PetscDrawAxisSetLabels(axis,"Particles","X","Y");
6996:     PetscDrawAxisSetLimits(axis, -5, 5, -5, 5);
6997:     PetscDrawAxisSetHoldLimits(axis, PETSC_TRUE);
6998:     TSGetDM(ts, &dm);
6999:     DMGetDimension(dm, &dim);
7000:     if(dim!=2) SETERRQ(PETSC_COMM_SELF, ierr, "Dimensions improper for monitor arguments! Current support: two dimensions.");
7001:     VecGetLocalSize(u, &Np);
7002:     Np /= 2*dim;
7003:     PetscDrawSPSetDimension(ctx->sp, Np);
7004:     PetscDrawSPReset(ctx->sp);
7005:   }
7006: 
7007:   VecGetLocalSize(u, &Np);
7008:   Np /= 2*dim;
7009:   VecGetArrayRead(u,&yy);
7010:   PetscMalloc2(Np, &x, Np, &y);
7011:   /* get points from solution vector */
7012:   for (p=0; p<Np; ++p){
7013:     x[p] = PetscRealPart(yy[2*dim*p]);
7014:     y[p] = PetscRealPart(yy[2*dim*p+1]);
7015:   }
7016:   VecRestoreArrayRead(u,&yy);
7017: 
7018:   if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) {
7019:     PetscDrawSPAddPoint(ctx->sp,x,y);
7020:     PetscDrawSPDraw(ctx->sp,PETSC_FALSE);
7021:     PetscDrawSPSave(ctx->sp);
7022:   }

7024:   PetscFree2(x, y);

7026:   return(0);
7027: }



7031: /*@C
7032:    TSMonitorError - Monitors progress of the TS solvers by printing the 2 norm of the error at each timestep

7034:    Collective on TS

7036:    Input Parameters:
7037: +  ts - the TS context
7038: .  step - current time-step
7039: .  ptime - current time
7040: .  u - current solution
7041: -  dctx - unused context

7043:    Level: intermediate

7045:    The user must provide the solution using TSSetSolutionFunction() to use this monitor.

7047:    Options Database Keys:
7048: .  -ts_monitor_error - create a graphical monitor of error history

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

7052: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction()
7053: @*/
7054: PetscErrorCode  TSMonitorError(TS ts,PetscInt step,PetscReal ptime,Vec u,PetscViewerAndFormat *vf)
7055: {
7056:   PetscErrorCode    ierr;
7057:   Vec               y;
7058:   PetscReal         nrm;
7059:   PetscBool         flg;

7062:   VecDuplicate(u,&y);
7063:   TSComputeSolutionFunction(ts,ptime,y);
7064:   VecAXPY(y,-1.0,u);
7065:   PetscObjectTypeCompare((PetscObject)vf->viewer,PETSCVIEWERASCII,&flg);
7066:   if (flg) {
7067:     VecNorm(y,NORM_2,&nrm);
7068:     PetscViewerASCIIPrintf(vf->viewer,"2-norm of error %g\n",(double)nrm);
7069:   }
7070:   PetscObjectTypeCompare((PetscObject)vf->viewer,PETSCVIEWERDRAW,&flg);
7071:   if (flg) {
7072:     VecView(y,vf->viewer);
7073:   }
7074:   VecDestroy(&y);
7075:   return(0);
7076: }

7078: PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
7079: {
7080:   TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
7081:   PetscReal      x   = ptime,y;
7083:   PetscInt       its;

7086:   if (n < 0) return(0); /* -1 indicates interpolated solution */
7087:   if (!n) {
7088:     PetscDrawAxis axis;
7089:     PetscDrawLGGetAxis(ctx->lg,&axis);
7090:     PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");
7091:     PetscDrawLGReset(ctx->lg);
7092:     ctx->snes_its = 0;
7093:   }
7094:   TSGetSNESIterations(ts,&its);
7095:   y    = its - ctx->snes_its;
7096:   PetscDrawLGAddPoint(ctx->lg,&x,&y);
7097:   if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
7098:     PetscDrawLGDraw(ctx->lg);
7099:     PetscDrawLGSave(ctx->lg);
7100:   }
7101:   ctx->snes_its = its;
7102:   return(0);
7103: }

7105: PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx)
7106: {
7107:   TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx;
7108:   PetscReal      x   = ptime,y;
7110:   PetscInt       its;

7113:   if (n < 0) return(0); /* -1 indicates interpolated solution */
7114:   if (!n) {
7115:     PetscDrawAxis axis;
7116:     PetscDrawLGGetAxis(ctx->lg,&axis);
7117:     PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");
7118:     PetscDrawLGReset(ctx->lg);
7119:     ctx->ksp_its = 0;
7120:   }
7121:   TSGetKSPIterations(ts,&its);
7122:   y    = its - ctx->ksp_its;
7123:   PetscDrawLGAddPoint(ctx->lg,&x,&y);
7124:   if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) {
7125:     PetscDrawLGDraw(ctx->lg);
7126:     PetscDrawLGSave(ctx->lg);
7127:   }
7128:   ctx->ksp_its = its;
7129:   return(0);
7130: }

7132: /*@
7133:    TSComputeLinearStability - computes the linear stability function at a point

7135:    Collective on TS and Vec

7137:    Input Parameters:
7138: +  ts - the TS context
7139: -  xr,xi - real and imaginary part of input arguments

7141:    Output Parameters:
7142: .  yr,yi - real and imaginary part of function value

7144:    Level: developer

7146: .keywords: TS, compute

7148: .seealso: TSSetRHSFunction(), TSComputeIFunction()
7149: @*/
7150: PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
7151: {

7156:   if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method");
7157:   (*ts->ops->linearstability)(ts,xr,xi,yr,yi);
7158:   return(0);
7159: }

7161: /* ------------------------------------------------------------------------*/
7162: /*@C
7163:    TSMonitorEnvelopeCtxCreate - Creates a context for use with TSMonitorEnvelope()

7165:    Collective on TS

7167:    Input Parameters:
7168: .  ts  - the ODE solver object

7170:    Output Parameter:
7171: .  ctx - the context

7173:    Level: intermediate

7175: .keywords: TS, monitor, line graph, residual, seealso

7177: .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()

7179: @*/
7180: PetscErrorCode  TSMonitorEnvelopeCtxCreate(TS ts,TSMonitorEnvelopeCtx *ctx)
7181: {

7185:   PetscNew(ctx);
7186:   return(0);
7187: }

7189: /*@C
7190:    TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution

7192:    Collective on TS

7194:    Input Parameters:
7195: +  ts - the TS context
7196: .  step - current time-step
7197: .  ptime - current time
7198: .  u  - current solution
7199: -  dctx - the envelope context

7201:    Options Database:
7202: .  -ts_monitor_envelope

7204:    Level: intermediate

7206:    Notes:
7207:     after a solve you can use TSMonitorEnvelopeGetBounds() to access the envelope

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

7211: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorEnvelopeGetBounds(), TSMonitorEnvelopeCtxCreate()
7212: @*/
7213: PetscErrorCode  TSMonitorEnvelope(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dctx)
7214: {
7215:   PetscErrorCode       ierr;
7216:   TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dctx;

7219:   if (!ctx->max) {
7220:     VecDuplicate(u,&ctx->max);
7221:     VecDuplicate(u,&ctx->min);
7222:     VecCopy(u,ctx->max);
7223:     VecCopy(u,ctx->min);
7224:   } else {
7225:     VecPointwiseMax(ctx->max,u,ctx->max);
7226:     VecPointwiseMin(ctx->min,u,ctx->min);
7227:   }
7228:   return(0);
7229: }

7231: /*@C
7232:    TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution

7234:    Collective on TS

7236:    Input Parameter:
7237: .  ts - the TS context

7239:    Output Parameter:
7240: +  max - the maximum values
7241: -  min - the minimum values

7243:    Notes:
7244:     If the TS does not have a TSMonitorEnvelopeCtx associated with it then this function is ignored

7246:    Level: intermediate

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

7250: .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables()
7251: @*/
7252: PetscErrorCode  TSMonitorEnvelopeGetBounds(TS ts,Vec *max,Vec *min)
7253: {
7254:   PetscInt i;

7257:   if (max) *max = NULL;
7258:   if (min) *min = NULL;
7259:   for (i=0; i<ts->numbermonitors; i++) {
7260:     if (ts->monitor[i] == TSMonitorEnvelope) {
7261:       TSMonitorEnvelopeCtx  ctx = (TSMonitorEnvelopeCtx) ts->monitorcontext[i];
7262:       if (max) *max = ctx->max;
7263:       if (min) *min = ctx->min;
7264:       break;
7265:     }
7266:   }
7267:   return(0);
7268: }

7270: /*@C
7271:    TSMonitorEnvelopeCtxDestroy - Destroys a context that was created  with TSMonitorEnvelopeCtxCreate().

7273:    Collective on TSMonitorEnvelopeCtx

7275:    Input Parameter:
7276: .  ctx - the monitor context

7278:    Level: intermediate

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

7282: .seealso: TSMonitorLGCtxCreate(),  TSMonitorSet(), TSMonitorLGTimeStep()
7283: @*/
7284: PetscErrorCode  TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx)
7285: {

7289:   VecDestroy(&(*ctx)->min);
7290:   VecDestroy(&(*ctx)->max);
7291:   PetscFree(*ctx);
7292:   return(0);
7293: }

7295: /*@
7296:    TSRestartStep - Flags the solver to restart the next step

7298:    Collective on TS

7300:    Input Parameter:
7301: .  ts - the TS context obtained from TSCreate()

7303:    Level: advanced

7305:    Notes:
7306:    Multistep methods like BDF or Runge-Kutta methods with FSAL property require restarting the solver in the event of
7307:    discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
7308:    vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
7309:    the sake of correctness and maximum safety, users are expected to call TSRestart() whenever they introduce
7310:    discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
7311:    discontinuous source terms).

7313: .keywords: TS, timestep, restart

7315: .seealso: TSSolve(), TSSetPreStep(), TSSetPostStep()
7316: @*/
7317: PetscErrorCode TSRestartStep(TS ts)
7318: {
7321:   ts->steprestart = PETSC_TRUE;
7322:   return(0);
7323: }

7325: /*@
7326:    TSRollBack - Rolls back one time step

7328:    Collective on TS

7330:    Input Parameter:
7331: .  ts - the TS context obtained from TSCreate()

7333:    Level: advanced

7335: .keywords: TS, timestep, rollback

7337: .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate()
7338: @*/
7339: PetscErrorCode  TSRollBack(TS ts)
7340: {

7345:   if (ts->steprollback) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"TSRollBack already called");
7346:   if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name);
7347:   (*ts->ops->rollback)(ts);
7348:   ts->time_step = ts->ptime - ts->ptime_prev;
7349:   ts->ptime = ts->ptime_prev;
7350:   ts->ptime_prev = ts->ptime_prev_rollback;
7351:   ts->steps--;
7352:   ts->steprollback = PETSC_TRUE;
7353:   return(0);
7354: }

7356: /*@
7357:    TSGetStages - Get the number of stages and stage values

7359:    Input Parameter:
7360: .  ts - the TS context obtained from TSCreate()

7362:    Output Parameters:
7363: +  ns - the number of stages
7364: -  Y - the current stage vectors

7366:    Level: advanced

7368:    Notes: Both ns and Y can be NULL.

7370: .keywords: TS, getstages

7372: .seealso: TSCreate()
7373: @*/
7374: PetscErrorCode  TSGetStages(TS ts,PetscInt *ns,Vec **Y)
7375: {

7382:   if (!ts->ops->getstages) {
7383:     if (ns) *ns = 0;
7384:     if (Y) *Y = NULL;
7385:   } else {
7386:     (*ts->ops->getstages)(ts,ns,Y);
7387:   }
7388:   return(0);
7389: }

7391: /*@C
7392:   TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.

7394:   Collective on SNES

7396:   Input Parameters:
7397: + ts - the TS context
7398: . t - current timestep
7399: . U - state vector
7400: . Udot - time derivative of state vector
7401: . shift - shift to apply, see note below
7402: - ctx - an optional user context

7404:   Output Parameters:
7405: + J - Jacobian matrix (not altered in this routine)
7406: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)

7408:   Level: intermediate

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

7413:   dF/dU + shift*dF/dUdot

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

7418:   This will first try to get the coloring from the DM.  If the DM type has no coloring
7419:   routine, then it will try to get the coloring from the matrix.  This requires that the
7420:   matrix have nonzero entries precomputed.

7422: .keywords: TS, finite differences, Jacobian, coloring, sparse
7423: .seealso: TSSetIJacobian(), MatFDColoringCreate(), MatFDColoringSetFunction()
7424: @*/
7425: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat J,Mat B,void *ctx)
7426: {
7427:   SNES           snes;
7428:   MatFDColoring  color;
7429:   PetscBool      hascolor, matcolor = PETSC_FALSE;

7433:   PetscOptionsGetBool(((PetscObject)ts)->options,((PetscObject) ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL);
7434:   PetscObjectQuery((PetscObject) B, "TSMatFDColoring", (PetscObject *) &color);
7435:   if (!color) {
7436:     DM         dm;
7437:     ISColoring iscoloring;

7439:     TSGetDM(ts, &dm);
7440:     DMHasColoring(dm, &hascolor);
7441:     if (hascolor && !matcolor) {
7442:       DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring);
7443:       MatFDColoringCreate(B, iscoloring, &color);
7444:       MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);
7445:       MatFDColoringSetFromOptions(color);
7446:       MatFDColoringSetUp(B, iscoloring, color);
7447:       ISColoringDestroy(&iscoloring);
7448:     } else {
7449:       MatColoring mc;

7451:       MatColoringCreate(B, &mc);
7452:       MatColoringSetDistance(mc, 2);
7453:       MatColoringSetType(mc, MATCOLORINGSL);
7454:       MatColoringSetFromOptions(mc);
7455:       MatColoringApply(mc, &iscoloring);
7456:       MatColoringDestroy(&mc);
7457:       MatFDColoringCreate(B, iscoloring, &color);
7458:       MatFDColoringSetFunction(color, (PetscErrorCode (*)(void)) SNESTSFormFunction, (void *) ts);
7459:       MatFDColoringSetFromOptions(color);
7460:       MatFDColoringSetUp(B, iscoloring, color);
7461:       ISColoringDestroy(&iscoloring);
7462:     }
7463:     PetscObjectCompose((PetscObject) B, "TSMatFDColoring", (PetscObject) color);
7464:     PetscObjectDereference((PetscObject) color);
7465:   }
7466:   TSGetSNES(ts, &snes);
7467:   MatFDColoringApply(B, color, U, snes);
7468:   if (J != B) {
7469:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
7470:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
7471:   }
7472:   return(0);
7473: }

7475: /*@
7476:     TSSetFunctionDomainError - Set the function testing if the current state vector is valid

7478:     Input Parameters:
7479:     ts - the TS context
7480:     func - function called within TSFunctionDomainError

7482:     Level: intermediate

7484: .keywords: TS, state, domain
7485: .seealso: TSAdaptCheckStage(), TSFunctionDomainError()
7486: @*/

7488: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS,PetscReal,Vec,PetscBool*))
7489: {
7492:   ts->functiondomainerror = func;
7493:   return(0);
7494: }

7496: /*@
7497:     TSFunctionDomainError - Check if the current state is valid

7499:     Input Parameters:
7500:     ts - the TS context
7501:     stagetime - time of the simulation
7502:     Y - state vector to check.

7504:     Output Parameter:
7505:     accept - Set to PETSC_FALSE if the current state vector is valid.

7507:     Note:
7508:     This function should be used to ensure the state is in a valid part of the space.
7509:     For example, one can ensure here all values are positive.

7511:     Level: advanced
7512: @*/
7513: PetscErrorCode TSFunctionDomainError(TS ts,PetscReal stagetime,Vec Y,PetscBool* accept)
7514: {
7517:   *accept = PETSC_TRUE;
7518:   if (ts->functiondomainerror) {
7519:     PetscStackCallStandard((*ts->functiondomainerror),(ts,stagetime,Y,accept));
7520:   }
7521:   return(0);
7522: }

7524: /*@C
7525:   TSClone - This function clones a time step object.

7527:   Collective on MPI_Comm

7529:   Input Parameter:
7530: . tsin    - The input TS

7532:   Output Parameter:
7533: . tsout   - The output TS (cloned)

7535:   Notes:
7536:   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.

7538:   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);

7540:   Level: developer

7542: .keywords: TS, clone
7543: .seealso: TSCreate(), TSSetType(), TSSetUp(), TSDestroy(), TSSetProblemType()
7544: @*/
7545: PetscErrorCode  TSClone(TS tsin, TS *tsout)
7546: {
7547:   TS             t;
7549:   SNES           snes_start;
7550:   DM             dm;
7551:   TSType         type;

7555:   *tsout = NULL;

7557:   PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView);

7559:   /* General TS description */
7560:   t->numbermonitors    = 0;
7561:   t->setupcalled       = 0;
7562:   t->ksp_its           = 0;
7563:   t->snes_its          = 0;
7564:   t->nwork             = 0;
7565:   t->rhsjacobian.time  = -1e20;
7566:   t->rhsjacobian.scale = 1.;
7567:   t->ijacobian.shift   = 1.;

7569:   TSGetSNES(tsin,&snes_start);
7570:   TSSetSNES(t,snes_start);

7572:   TSGetDM(tsin,&dm);
7573:   TSSetDM(t,dm);

7575:   t->adapt = tsin->adapt;
7576:   PetscObjectReference((PetscObject)t->adapt);

7578:   t->trajectory = tsin->trajectory;
7579:   PetscObjectReference((PetscObject)t->trajectory);

7581:   t->event = tsin->event;
7582:   if (t->event) t->event->refct++;

7584:   t->problem_type      = tsin->problem_type;
7585:   t->ptime             = tsin->ptime;
7586:   t->ptime_prev        = tsin->ptime_prev;
7587:   t->time_step         = tsin->time_step;
7588:   t->max_time          = tsin->max_time;
7589:   t->steps             = tsin->steps;
7590:   t->max_steps         = tsin->max_steps;
7591:   t->equation_type     = tsin->equation_type;
7592:   t->atol              = tsin->atol;
7593:   t->rtol              = tsin->rtol;
7594:   t->max_snes_failures = tsin->max_snes_failures;
7595:   t->max_reject        = tsin->max_reject;
7596:   t->errorifstepfailed = tsin->errorifstepfailed;

7598:   TSGetType(tsin,&type);
7599:   TSSetType(t,type);

7601:   t->vec_sol           = NULL;

7603:   t->cfltime          = tsin->cfltime;
7604:   t->cfltime_local    = tsin->cfltime_local;
7605:   t->exact_final_time = tsin->exact_final_time;

7607:   PetscMemcpy(t->ops,tsin->ops,sizeof(struct _TSOps));

7609:   if (((PetscObject)tsin)->fortran_func_pointers) {
7610:     PetscInt i;
7611:     PetscMalloc((10)*sizeof(void(*)(void)),&((PetscObject)t)->fortran_func_pointers);
7612:     for (i=0; i<10; i++) {
7613:       ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
7614:     }
7615:   }
7616:   *tsout = t;
7617:   return(0);
7618: }

7620: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void* ctx,Vec x,Vec y)
7621: {
7623:   TS             ts = (TS) ctx;

7626:   TSComputeRHSFunction(ts,0,x,y);
7627:   return(0);
7628: }

7630: /*@
7631:     TSRHSJacobianTest - Compares the multiply routine provided to the MATSHELL with differencing on the TS given RHS function.

7633:    Logically Collective on TS and Mat

7635:     Input Parameters:
7636:     TS - the time stepping routine

7638:    Output Parameter:
7639: .   flg - PETSC_TRUE if the multiply is likely correct

7641:    Options Database:
7642:  .   -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator

7644:    Level: advanced

7646:    Notes:
7647:     This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian

7649: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose(), TSRHSJacobianTestTranspose()
7650: @*/
7651: PetscErrorCode  TSRHSJacobianTest(TS ts,PetscBool *flg)
7652: {
7653:   Mat            J,B;
7655:   TSRHSJacobian  func;
7656:   void*          ctx;

7659:   TSGetRHSJacobian(ts,&J,&B,&func,&ctx);
7660:   (*func)(ts,0.0,ts->vec_sol,J,B,ctx);
7661:   MatShellTestMult(J,RHSWrapperFunction_TSRHSJacobianTest,ts->vec_sol,ts,flg);
7662:   return(0);
7663: }

7665: /*@C
7666:     TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the MATSHELL with differencing on the TS given RHS function.

7668:    Logically Collective on TS and Mat

7670:     Input Parameters:
7671:     TS - the time stepping routine

7673:    Output Parameter:
7674: .   flg - PETSC_TRUE if the multiply is likely correct

7676:    Options Database:
7677: .   -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator

7679:    Notes:
7680:     This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian

7682:    Level: advanced

7684: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose(), TSRHSJacobianTest()
7685: @*/
7686: PetscErrorCode  TSRHSJacobianTestTranspose(TS ts,PetscBool *flg)
7687: {
7688:   Mat            J,B;
7690:   void           *ctx;
7691:   TSRHSJacobian  func;

7694:   TSGetRHSJacobian(ts,&J,&B,&func,&ctx);
7695:   (*func)(ts,0.0,ts->vec_sol,J,B,ctx);
7696:   MatShellTestMultTranspose(J,RHSWrapperFunction_TSRHSJacobianTest,ts->vec_sol,ts,flg);
7697:   return(0);
7698: }