Actual source code: ts.c

  1: #include <petsc/private/tsimpl.h>
  2: #include <petscdmshell.h>
  3: #include <petscdmda.h>
  4: #include <petscviewer.h>
  5: #include <petscdraw.h>
  6: #include <petscconvest.h>

  8: #define SkipSmallValue(a, b, tol) \
  9:   if (PetscAbsScalar(a) < tol || PetscAbsScalar(b) < tol) continue;

 11: /* Logging support */
 12: PetscClassId  TS_CLASSID, DMTS_CLASSID;
 13: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;

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

 17: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt, TSAdaptType default_type)
 18: {
 21:   if (!((PetscObject)adapt)->type_name) TSAdaptSetType(adapt, default_type);
 22:   return 0;
 23: }

 25: /*@
 26:    TSSetFromOptions - Sets various TS parameters from user options.

 28:    Collective on TS

 30:    Input Parameter:
 31: .  ts - the TS context obtained from TSCreate()

 33:    Options Database Keys:
 34: +  -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSALPHA, TSGLLE, TSSSP, TSGLEE, TSBSYMP, TSIRK
 35: .  -ts_save_trajectory - checkpoint the solution at each time-step
 36: .  -ts_max_time <time> - maximum time to compute to
 37: .  -ts_time_span <t0,...tf> - sets the time span, solutions are computed and stored for each indicated time
 38: .  -ts_max_steps <steps> - maximum number of time-steps to take
 39: .  -ts_init_time <time> - initial time to start computation
 40: .  -ts_final_time <time> - final time to compute to (deprecated: use -ts_max_time)
 41: .  -ts_dt <dt> - initial time step
 42: .  -ts_exact_final_time <stepover,interpolate,matchstep> - whether to stop at the exact given final time and how to compute the solution at that time
 43: .  -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
 44: .  -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
 45: .  -ts_error_if_step_fails <true,false> - Error if no step succeeds
 46: .  -ts_rtol <rtol> - relative tolerance for local truncation error
 47: .  -ts_atol <atol> - Absolute tolerance for local truncation error
 48: .  -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function
 49: .  -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - test the Jacobian at each iteration against finite difference with RHS function
 50: .  -ts_adjoint_solve <yes,no> - After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
 51: .  -ts_fd_color - Use finite differences with coloring to compute IJacobian
 52: .  -ts_monitor - print information at each timestep
 53: .  -ts_monitor_cancel - Cancel all monitors
 54: .  -ts_monitor_lg_solution - Monitor solution graphically
 55: .  -ts_monitor_lg_error - Monitor error graphically
 56: .  -ts_monitor_error - Monitors norm of error
 57: .  -ts_monitor_lg_timestep - Monitor timestep size graphically
 58: .  -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
 59: .  -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
 60: .  -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
 61: .  -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
 62: .  -ts_monitor_draw_solution - Monitor solution graphically
 63: .  -ts_monitor_draw_solution_phase  <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
 64: .  -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
 65: .  -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
 66: .  -ts_monitor_solution_vtk <filename.vts,filename.vtu> - Save each time step to a binary file, use filename-%%03" PetscInt_FMT ".vts (filename-%%03" PetscInt_FMT ".vtu)
 67: -  -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time

 69:    Notes:
 70:      See SNESSetFromOptions() and KSPSetFromOptions() for how to control the nonlinear and linear solves used by the time-stepper.

 72:      Certain SNES options get reset for each new nonlinear solver, for example -snes_lag_jacobian <its> and -snes_lag_preconditioner <its>, in order
 73:      to retain them over the multiple nonlinear solves that TS uses you mush also provide -snes_lag_jacobian_persists true and
 74:      -snes_lag_preconditioner_persists true

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

 79:    Level: beginner

 81: .seealso: `TSGetType()`
 82: @*/
 83: PetscErrorCode TSSetFromOptions(TS ts)
 84: {
 85:   PetscBool              opt, flg, tflg;
 86:   char                   monfilename[PETSC_MAX_PATH_LEN];
 87:   PetscReal              time_step, tspan[100];
 88:   PetscInt               nt = PETSC_STATIC_ARRAY_LENGTH(tspan);
 89:   TSExactFinalTimeOption eftopt;
 90:   char                   dir[16];
 91:   TSIFunction            ifun;
 92:   const char            *defaultType;
 93:   char                   typeName[256];


 97:   TSRegisterAll();
 98:   TSGetIFunction(ts, NULL, &ifun, NULL);

100:   PetscObjectOptionsBegin((PetscObject)ts);
101:   if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
102:   else defaultType = ifun ? TSBEULER : TSEULER;
103:   PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, 256, &opt);
104:   if (opt) TSSetType(ts, typeName);
105:   else TSSetType(ts, defaultType);

107:   /* Handle generic TS options */
108:   PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL);
109:   PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL);
110:   PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", tspan, &nt, &flg);
111:   if (flg) TSSetTimeSpan(ts, nt, tspan);
112:   PetscOptionsInt("-ts_max_steps", "Maximum number of time steps", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL);
113:   PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL);
114:   PetscOptionsReal("-ts_dt", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg);
115:   if (flg) TSSetTimeStep(ts, time_step);
116:   PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg);
117:   if (flg) TSSetExactFinalTime(ts, eftopt);
118:   PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, NULL);
119:   PetscOptionsInt("-ts_max_reject", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, NULL);
120:   PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL);
121:   PetscOptionsReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL);
122:   PetscOptionsReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL);

124:   PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL);
125:   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);
126:   PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL);
127: #if defined(PETSC_HAVE_SAWS)
128:   {
129:     PetscBool set;
130:     flg = PETSC_FALSE;
131:     PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set);
132:     if (set) PetscObjectSAWsSetBlock((PetscObject)ts, flg);
133:   }
134: #endif

136:   /* Monitor options */
137:   PetscOptionsInt("-ts_monitor_frequency", "Number of time steps between monitor output", "TSMonitorSetFrequency", ts->monitorFrequency, &ts->monitorFrequency, NULL);
138:   TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL);
139:   TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL);
140:   TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, NULL);
141:   TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL);

143:   PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg);
144:   if (flg) PetscPythonMonitorSet((PetscObject)ts, monfilename);

146:   PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt);
147:   if (opt) {
148:     PetscInt  howoften = 1;
149:     DM        dm;
150:     PetscBool net;

152:     PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL);
153:     TSGetDM(ts, &dm);
154:     PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net);
155:     if (net) {
156:       TSMonitorLGCtxNetwork ctx;
157:       TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx);
158:       TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxNetworkDestroy);
159:       PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL);
160:     } else {
161:       TSMonitorLGCtx ctx;
162:       TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
163:       TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
164:     }
165:   }

167:   PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt);
168:   if (opt) {
169:     TSMonitorLGCtx ctx;
170:     PetscInt       howoften = 1;

172:     PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL);
173:     TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
174:     TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
175:   }
176:   TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL);

178:   PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt);
179:   if (opt) {
180:     TSMonitorLGCtx ctx;
181:     PetscInt       howoften = 1;

183:     PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL);
184:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
185:     TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
186:   }
187:   PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt);
188:   if (opt) {
189:     TSMonitorLGCtx ctx;
190:     PetscInt       howoften = 1;

192:     PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL);
193:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
194:     TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
195:     ctx->semilogy = PETSC_TRUE;
196:   }

198:   PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt);
199:   if (opt) {
200:     TSMonitorLGCtx ctx;
201:     PetscInt       howoften = 1;

203:     PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL);
204:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
205:     TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
206:   }
207:   PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt);
208:   if (opt) {
209:     TSMonitorLGCtx ctx;
210:     PetscInt       howoften = 1;

212:     PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL);
213:     TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
214:     TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
215:   }
216:   PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt);
217:   if (opt) {
218:     TSMonitorSPEigCtx ctx;
219:     PetscInt          howoften = 1;

221:     PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL);
222:     TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
223:     TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscErrorCode(*)(void **))TSMonitorSPEigCtxDestroy);
224:   }
225:   PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase from the DMSwarm", "TSMonitorSPSwarm", &opt);
226:   if (opt) {
227:     TSMonitorSPCtx ctx;
228:     PetscInt       howoften = 1, retain = 0;
229:     PetscBool      phase = PETSC_TRUE, create = PETSC_TRUE;

231:     for (PetscInt i = 0; i < ts->numbermonitors; ++i)
232:       if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
233:         create = PETSC_FALSE;
234:         break;
235:       }
236:     if (create) {
237:       PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL);
238:       PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL);
239:       PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL);
240:       TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, &ctx);
241:       TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorSPCtxDestroy);
242:     }
243:   }
244:   opt = PETSC_FALSE;
245:   PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt);
246:   if (opt) {
247:     TSMonitorDrawCtx ctx;
248:     PetscInt         howoften = 1;

250:     PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL);
251:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
252:     TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
253:   }
254:   opt = PETSC_FALSE;
255:   PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt);
256:   if (opt) {
257:     TSMonitorDrawCtx ctx;
258:     PetscReal        bounds[4];
259:     PetscInt         n = 4;
260:     PetscDraw        draw;
261:     PetscDrawAxis    axis;

263:     PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL);
265:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx);
266:     PetscViewerDrawGetDraw(ctx->viewer, 0, &draw);
267:     PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis);
268:     PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]);
269:     PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2");
270:     TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
271:   }
272:   opt = PETSC_FALSE;
273:   PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt);
274:   if (opt) {
275:     TSMonitorDrawCtx ctx;
276:     PetscInt         howoften = 1;

278:     PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL);
279:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
280:     TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
281:   }
282:   opt = PETSC_FALSE;
283:   PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt);
284:   if (opt) {
285:     TSMonitorDrawCtx ctx;
286:     PetscInt         howoften = 1;

288:     PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL);
289:     TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
290:     TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
291:   }

293:   opt = PETSC_FALSE;
294:   PetscOptionsString("-ts_monitor_solution_vtk", "Save each time step to a binary file, use filename-%%03" PetscInt_FMT ".vts", "TSMonitorSolutionVTK", NULL, monfilename, sizeof(monfilename), &flg);
295:   if (flg) {
296:     const char *ptr, *ptr2;
297:     char       *filetemplate;
299:     /* Do some cursory validation of the input. */
300:     PetscStrstr(monfilename, "%", (char **)&ptr);
302:     for (ptr++; ptr && *ptr; ptr++) {
303:       PetscStrchr("DdiouxX", *ptr, (char **)&ptr2);
305:       if (ptr2) break;
306:     }
307:     PetscStrallocpy(monfilename, &filetemplate);
308:     TSMonitorSet(ts, TSMonitorSolutionVTK, filetemplate, (PetscErrorCode(*)(void **))TSMonitorSolutionVTKDestroy);
309:   }

311:   PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg);
312:   if (flg) {
313:     TSMonitorDMDARayCtx *rayctx;
314:     int                  ray = 0;
315:     DMDirection          ddir;
316:     DM                   da;
317:     PetscMPIInt          rank;

320:     if (dir[0] == 'x') ddir = DM_X;
321:     else if (dir[0] == 'y') ddir = DM_Y;
322:     else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
323:     sscanf(dir + 2, "%d", &ray);

325:     PetscInfo(((PetscObject)ts), "Displaying DMDA ray %c = %d\n", dir[0], ray);
326:     PetscNew(&rayctx);
327:     TSGetDM(ts, &da);
328:     DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);
329:     MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank);
330:     if (rank == 0) PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer);
331:     rayctx->lgctx = NULL;
332:     TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy);
333:   }
334:   PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg);
335:   if (flg) {
336:     TSMonitorDMDARayCtx *rayctx;
337:     int                  ray = 0;
338:     DMDirection          ddir;
339:     DM                   da;
340:     PetscInt             howoften = 1;

343:     if (dir[0] == 'x') ddir = DM_X;
344:     else if (dir[0] == 'y') ddir = DM_Y;
345:     else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
346:     sscanf(dir + 2, "%d", &ray);

348:     PetscInfo(((PetscObject)ts), "Displaying LG DMDA ray %c = %d\n", dir[0], ray);
349:     PetscNew(&rayctx);
350:     TSGetDM(ts, &da);
351:     DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);
352:     TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx);
353:     TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy);
354:   }

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

360:     TSMonitorEnvelopeCtxCreate(ts, &ctx);
361:     TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscErrorCode(*)(void **))TSMonitorEnvelopeCtxDestroy);
362:   }
363:   flg = PETSC_FALSE;
364:   PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt);
365:   if (opt && flg) TSMonitorCancel(ts);

367:   flg = PETSC_FALSE;
368:   PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeJacobianDefaultColor", flg, &flg, NULL);
369:   if (flg) {
370:     DM dm;

372:     TSGetDM(ts, &dm);
373:     DMTSUnsetIJacobianContext_Internal(dm);
374:     TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL);
375:     PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n");
376:   }

378:   /* Handle specific TS options */
379:   PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);

381:   /* Handle TSAdapt options */
382:   TSGetAdapt(ts, &ts->adapt);
383:   TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type);
384:   TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject);

386:   /* TS trajectory must be set after TS, since it may use some TS options above */
387:   tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
388:   PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL);
389:   if (tflg) TSSetSaveTrajectory(ts);

391:   TSAdjointSetFromOptions(ts, PetscOptionsObject);

393:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
394:   PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject);
395:   PetscOptionsEnd();

397:   if (ts->trajectory) TSTrajectorySetFromOptions(ts->trajectory, ts);

399:   /* why do we have to do this here and not during TSSetUp? */
400:   TSGetSNES(ts, &ts->snes);
401:   if (ts->problem_type == TS_LINEAR) {
402:     PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, "");
403:     if (!flg) SNESSetType(ts->snes, SNESKSPONLY);
404:   }
405:   SNESSetFromOptions(ts->snes);
406:   return 0;
407: }

409: /*@
410:    TSGetTrajectory - Gets the trajectory from a TS if it exists

412:    Collective on TS

414:    Input Parameters:
415: .  ts - the TS context obtained from TSCreate()

417:    Output Parameters:
418: .  tr - the TSTrajectory object, if it exists

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

422:    Level: advanced

424: .seealso: `TSGetTrajectory()`, `TSAdjointSolve()`, `TSTrajectory`, `TSTrajectoryCreate()`

426: @*/
427: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
428: {
430:   *tr = ts->trajectory;
431:   return 0;
432: }

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

437:    Collective on TS

439:    Input Parameter:
440: .  ts - the TS context obtained from TSCreate()

442:    Options Database:
443: +  -ts_save_trajectory - saves the trajectory to a file
444: -  -ts_trajectory_type type - set trajectory type

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

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

451:    Level: intermediate

453: .seealso: `TSGetTrajectory()`, `TSAdjointSolve()`

455: @*/
456: PetscErrorCode TSSetSaveTrajectory(TS ts)
457: {
459:   if (!ts->trajectory) TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory);
460:   return 0;
461: }

463: /*@
464:    TSResetTrajectory - Destroys and recreates the internal TSTrajectory object

466:    Collective on TS

468:    Input Parameters:
469: .  ts - the TS context obtained from TSCreate()

471:    Level: intermediate

473: .seealso: `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`

475: @*/
476: PetscErrorCode TSResetTrajectory(TS ts)
477: {
479:   if (ts->trajectory) {
480:     TSTrajectoryDestroy(&ts->trajectory);
481:     TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory);
482:   }
483:   return 0;
484: }

486: /*@
487:    TSRemoveTrajectory - Destroys and removes the internal TSTrajectory object from TS

489:    Collective on TS

491:    Input Parameters:
492: .  ts - the TS context obtained from TSCreate()

494:    Level: intermediate

496: .seealso: `TSResetTrajectory()`, `TSAdjointSolve()`

498: @*/
499: PetscErrorCode TSRemoveTrajectory(TS ts)
500: {
502:   if (ts->trajectory) TSTrajectoryDestroy(&ts->trajectory);
503:   return 0;
504: }

506: /*@
507:    TSComputeRHSJacobian - Computes the Jacobian matrix that has been
508:       set with TSSetRHSJacobian().

510:    Collective on TS

512:    Input Parameters:
513: +  ts - the TS context
514: .  t - current timestep
515: -  U - input vector

517:    Output Parameters:
518: +  A - Jacobian matrix
519: -  B - optional preconditioning matrix

521:    Notes:
522:    Most users should not need to explicitly call this routine, as it
523:    is used internally within the nonlinear solvers.

525:    Level: developer

527: .seealso: `TSSetRHSJacobian()`, `KSPSetOperators()`
528: @*/
529: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
530: {
531:   PetscObjectState Ustate;
532:   PetscObjectId    Uid;
533:   DM               dm;
534:   DMTS             tsdm;
535:   TSRHSJacobian    rhsjacobianfunc;
536:   void            *ctx;
537:   TSRHSFunction    rhsfunction;

542:   TSGetDM(ts, &dm);
543:   DMGetDMTS(dm, &tsdm);
544:   DMTSGetRHSFunction(dm, &rhsfunction, NULL);
545:   DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx);
546:   PetscObjectStateGet((PetscObject)U, &Ustate);
547:   PetscObjectGetId((PetscObject)U, &Uid);

549:   if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) return 0;

552:   if (rhsjacobianfunc) {
553:     PetscLogEventBegin(TS_JacobianEval, ts, U, A, B);
554:     PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
555:     ts->rhsjacs++;
556:     PetscLogEventEnd(TS_JacobianEval, ts, U, A, B);
557:   } else {
558:     MatZeroEntries(A);
559:     if (B && A != B) MatZeroEntries(B);
560:   }
561:   ts->rhsjacobian.time  = t;
562:   ts->rhsjacobian.shift = 0;
563:   ts->rhsjacobian.scale = 1.;
564:   PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid);
565:   PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate);
566:   return 0;
567: }

569: /*@
570:    TSComputeRHSFunction - Evaluates the right-hand-side function.

572:    Collective on TS

574:    Input Parameters:
575: +  ts - the TS context
576: .  t - current time
577: -  U - state vector

579:    Output Parameter:
580: .  y - right hand side

582:    Note:
583:    Most users should not need to explicitly call this routine, as it
584:    is used internally within the nonlinear solvers.

586:    Level: developer

588: .seealso: `TSSetRHSFunction()`, `TSComputeIFunction()`
589: @*/
590: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
591: {
592:   TSRHSFunction rhsfunction;
593:   TSIFunction   ifunction;
594:   void         *ctx;
595:   DM            dm;

600:   TSGetDM(ts, &dm);
601:   DMTSGetRHSFunction(dm, &rhsfunction, &ctx);
602:   DMTSGetIFunction(dm, &ifunction, NULL);


606:   if (rhsfunction) {
607:     PetscLogEventBegin(TS_FunctionEval, ts, U, y, 0);
608:     VecLockReadPush(U);
609:     PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
610:     VecLockReadPop(U);
611:     ts->rhsfuncs++;
612:     PetscLogEventEnd(TS_FunctionEval, ts, U, y, 0);
613:   } else VecZeroEntries(y);
614:   return 0;
615: }

617: /*@
618:    TSComputeSolutionFunction - Evaluates the solution function.

620:    Collective on TS

622:    Input Parameters:
623: +  ts - the TS context
624: -  t - current time

626:    Output Parameter:
627: .  U - the solution

629:    Note:
630:    Most users should not need to explicitly call this routine, as it
631:    is used internally within the nonlinear solvers.

633:    Level: developer

635: .seealso: `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
636: @*/
637: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
638: {
639:   TSSolutionFunction solutionfunction;
640:   void              *ctx;
641:   DM                 dm;

645:   TSGetDM(ts, &dm);
646:   DMTSGetSolutionFunction(dm, &solutionfunction, &ctx);

648:   if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
649:   return 0;
650: }
651: /*@
652:    TSComputeForcingFunction - Evaluates the forcing function.

654:    Collective on TS

656:    Input Parameters:
657: +  ts - the TS context
658: -  t - current time

660:    Output Parameter:
661: .  U - the function value

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

667:    Level: developer

669: .seealso: `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
670: @*/
671: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
672: {
673:   void             *ctx;
674:   DM                dm;
675:   TSForcingFunction forcing;

679:   TSGetDM(ts, &dm);
680:   DMTSGetForcingFunction(dm, &forcing, &ctx);

682:   if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
683:   return 0;
684: }

686: static PetscErrorCode TSGetRHSVec_Private(TS ts, Vec *Frhs)
687: {
688:   Vec F;

690:   *Frhs = NULL;
691:   TSGetIFunction(ts, &F, NULL, NULL);
692:   if (!ts->Frhs) VecDuplicate(F, &ts->Frhs);
693:   *Frhs = ts->Frhs;
694:   return 0;
695: }

697: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
698: {
699:   Mat         A, B;
700:   TSIJacobian ijacobian;

702:   if (Arhs) *Arhs = NULL;
703:   if (Brhs) *Brhs = NULL;
704:   TSGetIJacobian(ts, &A, &B, &ijacobian, NULL);
705:   if (Arhs) {
706:     if (!ts->Arhs) {
707:       if (ijacobian) {
708:         MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs);
709:         TSSetMatStructure(ts, SAME_NONZERO_PATTERN);
710:       } else {
711:         ts->Arhs = A;
712:         PetscObjectReference((PetscObject)A);
713:       }
714:     } else {
715:       PetscBool flg;
716:       SNESGetUseMatrixFree(ts->snes, NULL, &flg);
717:       /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
718:       if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
719:         PetscObjectDereference((PetscObject)ts->Arhs);
720:         ts->Arhs = A;
721:         PetscObjectReference((PetscObject)A);
722:       }
723:     }
724:     *Arhs = ts->Arhs;
725:   }
726:   if (Brhs) {
727:     if (!ts->Brhs) {
728:       if (A != B) {
729:         if (ijacobian) {
730:           MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs);
731:         } else {
732:           ts->Brhs = B;
733:           PetscObjectReference((PetscObject)B);
734:         }
735:       } else {
736:         PetscObjectReference((PetscObject)ts->Arhs);
737:         ts->Brhs = ts->Arhs;
738:       }
739:     }
740:     *Brhs = ts->Brhs;
741:   }
742:   return 0;
743: }

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

748:    Collective on TS

750:    Input Parameters:
751: +  ts - the TS context
752: .  t - current time
753: .  U - state vector
754: .  Udot - time derivative of state vector
755: -  imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate

757:    Output Parameter:
758: .  Y - right hand side

760:    Note:
761:    Most users should not need to explicitly call this routine, as it
762:    is used internally within the nonlinear solvers.

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

767:    Level: developer

769: .seealso: `TSSetIFunction()`, `TSComputeRHSFunction()`
770: @*/
771: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
772: {
773:   TSIFunction   ifunction;
774:   TSRHSFunction rhsfunction;
775:   void         *ctx;
776:   DM            dm;


783:   TSGetDM(ts, &dm);
784:   DMTSGetIFunction(dm, &ifunction, &ctx);
785:   DMTSGetRHSFunction(dm, &rhsfunction, NULL);


789:   PetscLogEventBegin(TS_FunctionEval, ts, U, Udot, Y);
790:   if (ifunction) {
791:     PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
792:     ts->ifuncs++;
793:   }
794:   if (imex) {
795:     if (!ifunction) VecCopy(Udot, Y);
796:   } else if (rhsfunction) {
797:     if (ifunction) {
798:       Vec Frhs;
799:       TSGetRHSVec_Private(ts, &Frhs);
800:       TSComputeRHSFunction(ts, t, U, Frhs);
801:       VecAXPY(Y, -1, Frhs);
802:     } else {
803:       TSComputeRHSFunction(ts, t, U, Y);
804:       VecAYPX(Y, -1, Udot);
805:     }
806:   }
807:   PetscLogEventEnd(TS_FunctionEval, ts, U, Udot, Y);
808:   return 0;
809: }

811: /*
812:    TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call TSComputeRHSJacobian() on it.

814:    Note:
815:    This routine is needed when one switches from TSComputeIJacobian() to TSComputeRHSJacobian() because the Jacobian matrix may be shifted or scaled in TSComputeIJacobian().

817: */
818: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
819: {

824:   if (ts->rhsjacobian.shift) MatShift(A, -ts->rhsjacobian.shift);
825:   if (ts->rhsjacobian.scale == -1.) MatScale(A, -1);
826:   if (B && B == ts->Brhs && A != B) {
827:     if (ts->rhsjacobian.shift) MatShift(B, -ts->rhsjacobian.shift);
828:     if (ts->rhsjacobian.scale == -1.) MatScale(B, -1);
829:   }
830:   ts->rhsjacobian.shift = 0;
831:   ts->rhsjacobian.scale = 1.;
832:   return 0;
833: }

835: /*@
836:    TSComputeIJacobian - Evaluates the Jacobian of the DAE

838:    Collective on TS

840:    Input
841:       Input Parameters:
842: +  ts - the TS context
843: .  t - current timestep
844: .  U - state vector
845: .  Udot - time derivative of state vector
846: .  shift - shift to apply, see note below
847: -  imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate

849:    Output Parameters:
850: +  A - Jacobian matrix
851: -  B - matrix from which the preconditioner is constructed; often the same as A

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

856:    dF/dU + shift*dF/dUdot

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

861:    Level: developer

863: .seealso: `TSSetIJacobian()`
864: @*/
865: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
866: {
867:   TSIJacobian   ijacobian;
868:   TSRHSJacobian rhsjacobian;
869:   DM            dm;
870:   void         *ctx;


880:   TSGetDM(ts, &dm);
881:   DMTSGetIJacobian(dm, &ijacobian, &ctx);
882:   DMTSGetRHSJacobian(dm, &rhsjacobian, NULL);


886:   PetscLogEventBegin(TS_JacobianEval, ts, U, A, B);
887:   if (ijacobian) {
888:     PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
889:     ts->ijacs++;
890:   }
891:   if (imex) {
892:     if (!ijacobian) { /* system was written as Udot = G(t,U) */
893:       PetscBool assembled;
894:       if (rhsjacobian) {
895:         Mat Arhs = NULL;
896:         TSGetRHSMats_Private(ts, &Arhs, NULL);
897:         if (A == Arhs) {
899:           ts->rhsjacobian.time = PETSC_MIN_REAL;
900:         }
901:       }
902:       MatZeroEntries(A);
903:       MatAssembled(A, &assembled);
904:       if (!assembled) {
905:         MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
906:         MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
907:       }
908:       MatShift(A, shift);
909:       if (A != B) {
910:         MatZeroEntries(B);
911:         MatAssembled(B, &assembled);
912:         if (!assembled) {
913:           MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
914:           MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
915:         }
916:         MatShift(B, shift);
917:       }
918:     }
919:   } else {
920:     Mat Arhs = NULL, Brhs = NULL;

922:     /* RHSJacobian needs to be converted to part of IJacobian if exists */
923:     if (rhsjacobian) TSGetRHSMats_Private(ts, &Arhs, &Brhs);
924:     if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
925:       PetscObjectState Ustate;
926:       PetscObjectId    Uid;
927:       TSRHSFunction    rhsfunction;

929:       DMTSGetRHSFunction(dm, &rhsfunction, NULL);
930:       PetscObjectStateGet((PetscObject)U, &Ustate);
931:       PetscObjectGetId((PetscObject)U, &Uid);
932:       if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
933:           ts->rhsjacobian.scale == -1.) {                      /* No need to recompute RHSJacobian */
934:         MatShift(A, shift - ts->rhsjacobian.shift); /* revert the old shift and add the new shift with a single call to MatShift */
935:         if (A != B) MatShift(B, shift - ts->rhsjacobian.shift);
936:       } else {
937:         PetscBool flg;

939:         if (ts->rhsjacobian.reuse) { /* Undo the damage */
940:           /* MatScale has a short path for this case.
941:              However, this code path is taken the first time TSComputeRHSJacobian is called
942:              and the matrices have not been assembled yet */
943:           TSRecoverRHSJacobian(ts, A, B);
944:         }
945:         TSComputeRHSJacobian(ts, t, U, A, B);
946:         SNESGetUseMatrixFree(ts->snes, NULL, &flg);
947:         /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
948:         if (!flg) {
949:           MatScale(A, -1);
950:           MatShift(A, shift);
951:         }
952:         if (A != B) {
953:           MatScale(B, -1);
954:           MatShift(B, shift);
955:         }
956:       }
957:       ts->rhsjacobian.scale = -1;
958:       ts->rhsjacobian.shift = shift;
959:     } else if (Arhs) {  /* Both IJacobian and RHSJacobian */
960:       if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
961:         MatZeroEntries(A);
962:         MatShift(A, shift);
963:         if (A != B) {
964:           MatZeroEntries(B);
965:           MatShift(B, shift);
966:         }
967:       }
968:       TSComputeRHSJacobian(ts, t, U, Arhs, Brhs);
969:       MatAXPY(A, -1, Arhs, ts->axpy_pattern);
970:       if (A != B) MatAXPY(B, -1, Brhs, ts->axpy_pattern);
971:     }
972:   }
973:   PetscLogEventEnd(TS_JacobianEval, ts, U, A, B);
974:   return 0;
975: }

977: /*@C
978:     TSSetRHSFunction - Sets the routine for evaluating the function,
979:     where U_t = G(t,u).

981:     Logically Collective on TS

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

990:     Calling sequence of f:
991: $     PetscErrorCode f(TS ts,PetscReal t,Vec u,Vec F,void *ctx);

993: +   ts - timestep context
994: .   t - current timestep
995: .   u - input vector
996: .   F - function vector
997: -   ctx - [optional] user-defined function context

999:     Level: beginner

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

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


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

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

1030:     Logically Collective on TS

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

1038:     Calling sequence of f:
1039: $     PetscErrorCode f(TS ts,PetscReal t,Vec u,void *ctx);

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

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

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

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

1056:     Level: beginner

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

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

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

1073:     Logically Collective on TS

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

1081:     Calling sequence of func:
1082: $     PetscErrorCode func (TS ts,PetscReal t,Vec f,void *ctx);

1084: +   t - current timestep
1085: .   f - output vector
1086: -   ctx - [optional] user-defined function context

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

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

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

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

1100:     Level: beginner

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

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

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

1118:    Logically Collective on TS

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

1128:    Calling sequence of f:
1129: $     PetscErrorCode f(TS ts,PetscReal t,Vec u,Mat A,Mat B,void *ctx);

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

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

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

1143:    Level: beginner

1145: .seealso: `SNESComputeJacobianDefaultColor()`, `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`

1147: @*/
1148: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobian f, void *ctx)
1149: {
1150:   SNES        snes;
1151:   DM          dm;
1152:   TSIJacobian ijacobian;


1160:   TSGetDM(ts, &dm);
1161:   DMTSSetRHSJacobian(dm, f, ctx);
1162:   DMTSGetIJacobian(dm, &ijacobian, NULL);
1163:   TSGetSNES(ts, &snes);
1164:   if (!ijacobian) SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts);
1165:   if (Amat) {
1166:     PetscObjectReference((PetscObject)Amat);
1167:     MatDestroy(&ts->Arhs);
1168:     ts->Arhs = Amat;
1169:   }
1170:   if (Pmat) {
1171:     PetscObjectReference((PetscObject)Pmat);
1172:     MatDestroy(&ts->Brhs);
1173:     ts->Brhs = Pmat;
1174:   }
1175:   return 0;
1176: }

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

1181:    Logically Collective on TS

1183:    Input Parameters:
1184: +  ts  - the TS context obtained from TSCreate()
1185: .  r   - vector to hold the residual (or NULL to have it created internally)
1186: .  f   - the function evaluation routine
1187: -  ctx - user-defined context for private data for the function evaluation routine (may be NULL)

1189:    Calling sequence of f:
1190: $     PetscErrorCode f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx);

1192: +  t   - time at step/stage being solved
1193: .  u   - state vector
1194: .  u_t - time derivative of state vector
1195: .  F   - function vector
1196: -  ctx - [optional] user-defined context for matrix evaluation routine

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

1201:    Level: beginner

1203: .seealso: `TSSetRHSJacobian()`, `TSSetRHSFunction()`, `TSSetIJacobian()`
1204: @*/
1205: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunction f, void *ctx)
1206: {
1207:   SNES snes;
1208:   Vec  ralloc = NULL;
1209:   DM   dm;


1214:   TSGetDM(ts, &dm);
1215:   DMTSSetIFunction(dm, f, ctx);

1217:   TSGetSNES(ts, &snes);
1218:   if (!r && !ts->dm && ts->vec_sol) {
1219:     VecDuplicate(ts->vec_sol, &ralloc);
1220:     r = ralloc;
1221:   }
1222:   SNESSetFunction(snes, r, SNESTSFormFunction, ts);
1223:   VecDestroy(&ralloc);
1224:   return 0;
1225: }

1227: /*@C
1228:    TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.

1230:    Not Collective

1232:    Input Parameter:
1233: .  ts - the TS context

1235:    Output Parameters:
1236: +  r - vector to hold residual (or NULL)
1237: .  func - the function to compute residual (or NULL)
1238: -  ctx - the function context (or NULL)

1240:    Level: advanced

1242: .seealso: `TSSetIFunction()`, `SNESGetFunction()`
1243: @*/
1244: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunction *func, void **ctx)
1245: {
1246:   SNES snes;
1247:   DM   dm;

1250:   TSGetSNES(ts, &snes);
1251:   SNESGetFunction(snes, r, NULL, NULL);
1252:   TSGetDM(ts, &dm);
1253:   DMTSGetIFunction(dm, func, ctx);
1254:   return 0;
1255: }

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

1260:    Not Collective

1262:    Input Parameter:
1263: .  ts - the TS context

1265:    Output Parameters:
1266: +  r - vector to hold computed right hand side (or NULL)
1267: .  func - the function to compute right hand side (or NULL)
1268: -  ctx - the function context (or NULL)

1270:    Level: advanced

1272: .seealso: `TSSetRHSFunction()`, `SNESGetFunction()`
1273: @*/
1274: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunction *func, void **ctx)
1275: {
1276:   SNES snes;
1277:   DM   dm;

1280:   TSGetSNES(ts, &snes);
1281:   SNESGetFunction(snes, r, NULL, NULL);
1282:   TSGetDM(ts, &dm);
1283:   DMTSGetRHSFunction(dm, func, ctx);
1284:   return 0;
1285: }

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

1291:    Logically Collective on TS

1293:    Input Parameters:
1294: +  ts  - the TS context obtained from TSCreate()
1295: .  Amat - (approximate) Jacobian matrix
1296: .  Pmat - matrix used to compute preconditioner (usually the same as Amat)
1297: .  f   - the Jacobian evaluation routine
1298: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)

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

1303: +  t    - time at step/stage being solved
1304: .  U    - state vector
1305: .  U_t  - time derivative of state vector
1306: .  a    - shift
1307: .  Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1308: .  Pmat - matrix used for constructing preconditioner, usually the same as Amat
1309: -  ctx  - [optional] user-defined context for matrix evaluation routine

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

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

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

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

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

1329:    Level: beginner

1331: .seealso: `TSSetIFunction()`, `TSSetRHSJacobian()`, `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`

1333: @*/
1334: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobian f, void *ctx)
1335: {
1336:   SNES snes;
1337:   DM   dm;


1345:   TSGetDM(ts, &dm);
1346:   DMTSSetIJacobian(dm, f, ctx);

1348:   TSGetSNES(ts, &snes);
1349:   SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts);
1350:   return 0;
1351: }

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

1359:    Logically Collective

1361:    Input Parameters:
1362: +  ts - TS context obtained from TSCreate()
1363: -  reuse - PETSC_TRUE if the RHS Jacobian

1365:    Level: intermediate

1367: .seealso: `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1368: @*/
1369: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1370: {
1371:   ts->rhsjacobian.reuse = reuse;
1372:   return 0;
1373: }

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

1378:    Logically Collective on TS

1380:    Input Parameters:
1381: +  ts  - the TS context obtained from TSCreate()
1382: .  F   - vector to hold the residual (or NULL to have it created internally)
1383: .  fun - the function evaluation routine
1384: -  ctx - user-defined context for private data for the function evaluation routine (may be NULL)

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

1389: +  t    - time at step/stage being solved
1390: .  U    - state vector
1391: .  U_t  - time derivative of state vector
1392: .  U_tt - second time derivative of state vector
1393: .  F    - function vector
1394: -  ctx  - [optional] user-defined context for matrix evaluation routine (may be NULL)

1396:    Level: beginner

1398: .seealso: `TSSetI2Jacobian()`, `TSSetIFunction()`, `TSCreate()`, `TSSetRHSFunction()`
1399: @*/
1400: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2Function fun, void *ctx)
1401: {
1402:   DM dm;

1406:   TSSetIFunction(ts, F, NULL, NULL);
1407:   TSGetDM(ts, &dm);
1408:   DMTSSetI2Function(dm, fun, ctx);
1409:   return 0;
1410: }

1412: /*@C
1413:   TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.

1415:   Not Collective

1417:   Input Parameter:
1418: . ts - the TS context

1420:   Output Parameters:
1421: + r - vector to hold residual (or NULL)
1422: . fun - the function to compute residual (or NULL)
1423: - ctx - the function context (or NULL)

1425:   Level: advanced

1427: .seealso: `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1428: @*/
1429: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2Function *fun, void **ctx)
1430: {
1431:   SNES snes;
1432:   DM   dm;

1435:   TSGetSNES(ts, &snes);
1436:   SNESGetFunction(snes, r, NULL, NULL);
1437:   TSGetDM(ts, &dm);
1438:   DMTSGetI2Function(dm, fun, ctx);
1439:   return 0;
1440: }

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

1446:    Logically Collective on TS

1448:    Input Parameters:
1449: +  ts  - the TS context obtained from TSCreate()
1450: .  J   - Jacobian matrix
1451: .  P   - preconditioning matrix for J (may be same as J)
1452: .  jac - the Jacobian evaluation routine
1453: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)

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

1458: +  t    - time at step/stage being solved
1459: .  U    - state vector
1460: .  U_t  - time derivative of state vector
1461: .  U_tt - second time derivative of state vector
1462: .  v    - shift for U_t
1463: .  a    - shift for U_tt
1464: .  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
1465: .  P    - preconditioning matrix for J, may be same as J
1466: -  ctx  - [optional] user-defined context for matrix evaluation routine

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

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

1476:    Level: beginner

1478: .seealso: `TSSetI2Function()`, `TSGetI2Jacobian()`
1479: @*/
1480: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2Jacobian jac, void *ctx)
1481: {
1482:   DM dm;

1487:   TSSetIJacobian(ts, J, P, NULL, NULL);
1488:   TSGetDM(ts, &dm);
1489:   DMTSSetI2Jacobian(dm, jac, ctx);
1490:   return 0;
1491: }

1493: /*@C
1494:   TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.

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

1498:   Input Parameter:
1499: . ts  - The TS context obtained from TSCreate()

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

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

1510:   Level: advanced

1512: .seealso: `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`

1514: @*/
1515: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2Jacobian *jac, void **ctx)
1516: {
1517:   SNES snes;
1518:   DM   dm;

1520:   TSGetSNES(ts, &snes);
1521:   SNESSetUpMatrices(snes);
1522:   SNESGetJacobian(snes, J, P, NULL, NULL);
1523:   TSGetDM(ts, &dm);
1524:   DMTSGetI2Jacobian(dm, jac, ctx);
1525:   return 0;
1526: }

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

1531:   Collective on TS

1533:   Input Parameters:
1534: + ts - the TS context
1535: . t - current time
1536: . U - state vector
1537: . V - time derivative of state vector (U_t)
1538: - A - second time derivative of state vector (U_tt)

1540:   Output Parameter:
1541: . F - the residual vector

1543:   Note:
1544:   Most users should not need to explicitly call this routine, as it
1545:   is used internally within the nonlinear solvers.

1547:   Level: developer

1549: .seealso: `TSSetI2Function()`, `TSGetI2Function()`
1550: @*/
1551: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1552: {
1553:   DM            dm;
1554:   TSI2Function  I2Function;
1555:   void         *ctx;
1556:   TSRHSFunction rhsfunction;


1564:   TSGetDM(ts, &dm);
1565:   DMTSGetI2Function(dm, &I2Function, &ctx);
1566:   DMTSGetRHSFunction(dm, &rhsfunction, NULL);

1568:   if (!I2Function) {
1569:     TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE);
1570:     return 0;
1571:   }

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

1575:   PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));

1577:   if (rhsfunction) {
1578:     Vec Frhs;
1579:     TSGetRHSVec_Private(ts, &Frhs);
1580:     TSComputeRHSFunction(ts, t, U, Frhs);
1581:     VecAXPY(F, -1, Frhs);
1582:   }

1584:   PetscLogEventEnd(TS_FunctionEval, ts, U, V, F);
1585:   return 0;
1586: }

1588: /*@
1589:   TSComputeI2Jacobian - Evaluates the Jacobian of the DAE

1591:   Collective on TS

1593:   Input Parameters:
1594: + ts - the TS context
1595: . t - current timestep
1596: . U - state vector
1597: . V - time derivative of state vector
1598: . A - second time derivative of state vector
1599: . shiftV - shift to apply, see note below
1600: - shiftA - shift to apply, see note below

1602:   Output Parameters:
1603: + J - Jacobian matrix
1604: - P - optional preconditioning matrix

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

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

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

1614:   Level: developer

1616: .seealso: `TSSetI2Jacobian()`
1617: @*/
1618: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1619: {
1620:   DM            dm;
1621:   TSI2Jacobian  I2Jacobian;
1622:   void         *ctx;
1623:   TSRHSJacobian rhsjacobian;


1632:   TSGetDM(ts, &dm);
1633:   DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx);
1634:   DMTSGetRHSJacobian(dm, &rhsjacobian, NULL);

1636:   if (!I2Jacobian) {
1637:     TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE);
1638:     return 0;
1639:   }

1641:   PetscLogEventBegin(TS_JacobianEval, ts, U, J, P);
1642:   PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1643:   if (rhsjacobian) {
1644:     Mat Jrhs, Prhs;
1645:     TSGetRHSMats_Private(ts, &Jrhs, &Prhs);
1646:     TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs);
1647:     MatAXPY(J, -1, Jrhs, ts->axpy_pattern);
1648:     if (P != J) MatAXPY(P, -1, Prhs, ts->axpy_pattern);
1649:   }

1651:   PetscLogEventEnd(TS_JacobianEval, ts, U, J, P);
1652:   return 0;
1653: }

1655: /*@C
1656:    TSSetTransientVariable - sets function to transform from state to transient variables

1658:    Logically Collective

1660:    Input Parameters:
1661: +  ts - time stepping context on which to change the transient variable
1662: .  tvar - a function that transforms to transient variables
1663: -  ctx - a context for tvar

1665:     Calling sequence of tvar:
1666: $     PetscErrorCode tvar(TS ts,Vec p,Vec c,void *ctx);

1668: +   ts - timestep context
1669: .   p - input vector (primitive form)
1670: .   c - output vector, transient variables (conservative form)
1671: -   ctx - [optional] user-defined function context

1673:    Level: advanced

1675:    Notes:
1676:    This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., TSBDF)
1677:    can be conservative.  In this context, primitive variables P are used to model the state (e.g., because they lead to
1678:    well-conditioned formulations even in limiting cases such as low-Mach or zero porosity).  The transient variable is
1679:    C(P), specified by calling this function.  An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1680:    evaluated via the chain rule, as in

1682:      dF/dP + shift * dF/dCdot dC/dP.

1684: .seealso: `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1685: @*/
1686: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariable tvar, void *ctx)
1687: {
1688:   DM dm;

1691:   TSGetDM(ts, &dm);
1692:   DMTSSetTransientVariable(dm, tvar, ctx);
1693:   return 0;
1694: }

1696: /*@
1697:    TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables

1699:    Logically Collective

1701:    Input Parameters:
1702: +  ts - TS on which to compute
1703: -  U - state vector to be transformed to transient variables

1705:    Output Parameters:
1706: .  C - transient (conservative) variable

1708:    Developer Notes:
1709:    If DMTSSetTransientVariable() has not been called, then C is not modified in this routine and C=NULL is allowed.
1710:    This makes it safe to call without a guard.  One can use TSHasTransientVariable() to check if transient variables are
1711:    being used.

1713:    Level: developer

1715: .seealso: `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1716: @*/
1717: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1718: {
1719:   DM   dm;
1720:   DMTS dmts;

1724:   TSGetDM(ts, &dm);
1725:   DMGetDMTS(dm, &dmts);
1726:   if (dmts->ops->transientvar) {
1728:     (*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx);
1729:   }
1730:   return 0;
1731: }

1733: /*@
1734:    TSHasTransientVariable - determine whether transient variables have been set

1736:    Logically Collective

1738:    Input Parameters:
1739: .  ts - TS on which to compute

1741:    Output Parameters:
1742: .  has - PETSC_TRUE if transient variables have been set

1744:    Level: developer

1746: .seealso: `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1747: @*/
1748: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1749: {
1750:   DM   dm;
1751:   DMTS dmts;

1754:   TSGetDM(ts, &dm);
1755:   DMGetDMTS(dm, &dmts);
1756:   *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1757:   return 0;
1758: }

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

1764:    Logically Collective on TS

1766:    Input Parameters:
1767: +  ts - the TS context obtained from TSCreate()
1768: .  u - the solution vector
1769: -  v - the time derivative vector

1771:    Level: beginner

1773: @*/
1774: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1775: {
1779:   TSSetSolution(ts, u);
1780:   PetscObjectReference((PetscObject)v);
1781:   VecDestroy(&ts->vec_dot);
1782:   ts->vec_dot = v;
1783:   return 0;
1784: }

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

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

1794:    Input Parameter:
1795: .  ts - the TS context obtained from TSCreate()

1797:    Output Parameters:
1798: +  u - the vector containing the solution
1799: -  v - the vector containing the time derivative

1801:    Level: intermediate

1803: .seealso: `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`

1805: @*/
1806: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1807: {
1811:   if (u) *u = ts->vec_sol;
1812:   if (v) *v = ts->vec_dot;
1813:   return 0;
1814: }

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

1819:   Collective on PetscViewer

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

1826:    Level: intermediate

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

1831:   Notes for advanced users:
1832:   Most users should not need to know the details of the binary storage
1833:   format, since TSLoad() and TSView() completely hide these details.
1834:   But for anyone who's interested, the standard binary matrix storage
1835:   format is
1836: .vb
1837:      has not yet been determined
1838: .ve

1840: .seealso: `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1841: @*/
1842: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1843: {
1844:   PetscBool isbinary;
1845:   PetscInt  classid;
1846:   char      type[256];
1847:   DMTS      sdm;
1848:   DM        dm;

1852:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);

1855:   PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT);
1857:   PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR);
1858:   TSSetType(ts, type);
1859:   PetscTryTypeMethod(ts, load, viewer);
1860:   DMCreate(PetscObjectComm((PetscObject)ts), &dm);
1861:   DMLoad(dm, viewer);
1862:   TSSetDM(ts, dm);
1863:   DMCreateGlobalVector(ts->dm, &ts->vec_sol);
1864:   VecLoad(ts->vec_sol, viewer);
1865:   DMGetDMTS(ts->dm, &sdm);
1866:   DMTSLoad(sdm, viewer);
1867:   return 0;
1868: }

1870: #include <petscdraw.h>
1871: #if defined(PETSC_HAVE_SAWS)
1872: #include <petscviewersaws.h>
1873: #endif

1875: /*@C
1876:    TSViewFromOptions - View from Options

1878:    Collective on TS

1880:    Input Parameters:
1881: +  A - the application ordering context
1882: .  obj - Optional object
1883: -  name - command line option

1885:    Level: intermediate
1886: .seealso: `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1887: @*/
1888: PetscErrorCode TSViewFromOptions(TS A, PetscObject obj, const char name[])
1889: {
1891:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
1892:   return 0;
1893: }

1895: /*@C
1896:     TSView - Prints the TS data structure.

1898:     Collective on TS

1900:     Input Parameters:
1901: +   ts - the TS context obtained from TSCreate()
1902: -   viewer - visualization context

1904:     Options Database Key:
1905: .   -ts_view - calls TSView() at end of TSStep()

1907:     Notes:
1908:     The available visualization contexts include
1909: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
1910: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1911:          output where only the first processor opens
1912:          the file.  All other processors send their
1913:          data to the first processor to print.

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

1918:     In the debugger you can do "call TSView(ts,0)" to display the TS solver. (The same holds for any PETSc object viewer).

1920:     Level: beginner

1922: .seealso: `PetscViewerASCIIOpen()`
1923: @*/
1924: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1925: {
1926:   TSType    type;
1927:   PetscBool iascii, isstring, isundials, isbinary, isdraw;
1928:   DMTS      sdm;
1929: #if defined(PETSC_HAVE_SAWS)
1930:   PetscBool issaws;
1931: #endif

1934:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer);

1938:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1939:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
1940:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
1941:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1942: #if defined(PETSC_HAVE_SAWS)
1943:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws);
1944: #endif
1945:   if (iascii) {
1946:     PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer);
1947:     if (ts->ops->view) {
1948:       PetscViewerASCIIPushTab(viewer);
1949:       PetscUseTypeMethod(ts, view, viewer);
1950:       PetscViewerASCIIPopTab(viewer);
1951:     }
1952:     if (ts->max_steps < PETSC_MAX_INT) PetscViewerASCIIPrintf(viewer, "  maximum steps=%" PetscInt_FMT "\n", ts->max_steps);
1953:     if (ts->max_time < PETSC_MAX_REAL) PetscViewerASCIIPrintf(viewer, "  maximum time=%g\n", (double)ts->max_time);
1954:     if (ts->ifuncs) PetscViewerASCIIPrintf(viewer, "  total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs);
1955:     if (ts->ijacs) PetscViewerASCIIPrintf(viewer, "  total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs);
1956:     if (ts->rhsfuncs) PetscViewerASCIIPrintf(viewer, "  total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs);
1957:     if (ts->rhsjacs) PetscViewerASCIIPrintf(viewer, "  total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs);
1958:     if (ts->usessnes) {
1959:       PetscBool lin;
1960:       if (ts->problem_type == TS_NONLINEAR) PetscViewerASCIIPrintf(viewer, "  total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its);
1961:       PetscViewerASCIIPrintf(viewer, "  total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its);
1962:       PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, "");
1963:       PetscViewerASCIIPrintf(viewer, "  total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures);
1964:     }
1965:     PetscViewerASCIIPrintf(viewer, "  total number of rejected steps=%" PetscInt_FMT "\n", ts->reject);
1966:     if (ts->vrtol) PetscViewerASCIIPrintf(viewer, "  using vector of relative error tolerances, ");
1967:     else PetscViewerASCIIPrintf(viewer, "  using relative error tolerance of %g, ", (double)ts->rtol);
1968:     if (ts->vatol) PetscViewerASCIIPrintf(viewer, "  using vector of absolute error tolerances\n");
1969:     else PetscViewerASCIIPrintf(viewer, "  using absolute error tolerance of %g\n", (double)ts->atol);
1970:     PetscViewerASCIIPushTab(viewer);
1971:     TSAdaptView(ts->adapt, viewer);
1972:     PetscViewerASCIIPopTab(viewer);
1973:   } else if (isstring) {
1974:     TSGetType(ts, &type);
1975:     PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type);
1976:     PetscTryTypeMethod(ts, view, viewer);
1977:   } else if (isbinary) {
1978:     PetscInt    classid = TS_FILE_CLASSID;
1979:     MPI_Comm    comm;
1980:     PetscMPIInt rank;
1981:     char        type[256];

1983:     PetscObjectGetComm((PetscObject)ts, &comm);
1984:     MPI_Comm_rank(comm, &rank);
1985:     if (rank == 0) {
1986:       PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT);
1987:       PetscStrncpy(type, ((PetscObject)ts)->type_name, 256);
1988:       PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR);
1989:     }
1990:     PetscTryTypeMethod(ts, view, viewer);
1991:     if (ts->adapt) TSAdaptView(ts->adapt, viewer);
1992:     DMView(ts->dm, viewer);
1993:     VecView(ts->vec_sol, viewer);
1994:     DMGetDMTS(ts->dm, &sdm);
1995:     DMTSView(sdm, viewer);
1996:   } else if (isdraw) {
1997:     PetscDraw draw;
1998:     char      str[36];
1999:     PetscReal x, y, bottom, h;

2001:     PetscViewerDrawGetDraw(viewer, 0, &draw);
2002:     PetscDrawGetCurrentPoint(draw, &x, &y);
2003:     PetscStrcpy(str, "TS: ");
2004:     PetscStrcat(str, ((PetscObject)ts)->type_name);
2005:     PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h);
2006:     bottom = y - h;
2007:     PetscDrawPushCurrentPoint(draw, x, bottom);
2008:     PetscTryTypeMethod(ts, view, viewer);
2009:     if (ts->adapt) TSAdaptView(ts->adapt, viewer);
2010:     if (ts->snes) SNESView(ts->snes, viewer);
2011:     PetscDrawPopCurrentPoint(draw);
2012: #if defined(PETSC_HAVE_SAWS)
2013:   } else if (issaws) {
2014:     PetscMPIInt rank;
2015:     const char *name;

2017:     PetscObjectGetName((PetscObject)ts, &name);
2018:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
2019:     if (!((PetscObject)ts)->amsmem && rank == 0) {
2020:       char dir[1024];

2022:       PetscObjectViewSAWs((PetscObject)ts, viewer);
2023:       PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name);
2024:       SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT);
2025:       PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name);
2026:       SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE);
2027:     }
2028:     PetscTryTypeMethod(ts, view, viewer);
2029: #endif
2030:   }
2031:   if (ts->snes && ts->usessnes) {
2032:     PetscViewerASCIIPushTab(viewer);
2033:     SNESView(ts->snes, viewer);
2034:     PetscViewerASCIIPopTab(viewer);
2035:   }
2036:   DMGetDMTS(ts->dm, &sdm);
2037:   DMTSView(sdm, viewer);

2039:   PetscViewerASCIIPushTab(viewer);
2040:   PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials);
2041:   PetscViewerASCIIPopTab(viewer);
2042:   return 0;
2043: }

2045: /*@
2046:    TSSetApplicationContext - Sets an optional user-defined context for
2047:    the timesteppers.

2049:    Logically Collective on TS

2051:    Input Parameters:
2052: +  ts - the TS context obtained from TSCreate()
2053: -  usrP - optional user context

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

2059:    Level: intermediate

2061: .seealso: `TSGetApplicationContext()`
2062: @*/
2063: PetscErrorCode TSSetApplicationContext(TS ts, void *usrP)
2064: {
2066:   ts->user = usrP;
2067:   return 0;
2068: }

2070: /*@
2071:     TSGetApplicationContext - Gets the user-defined context for the
2072:     timestepper.

2074:     Not Collective

2076:     Input Parameter:
2077: .   ts - the TS context obtained from TSCreate()

2079:     Output Parameter:
2080: .   usrP - user context

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

2086:     Level: intermediate

2088: .seealso: `TSSetApplicationContext()`
2089: @*/
2090: PetscErrorCode TSGetApplicationContext(TS ts, void *usrP)
2091: {
2093:   *(void **)usrP = ts->user;
2094:   return 0;
2095: }

2097: /*@
2098:    TSGetStepNumber - Gets the number of steps completed.

2100:    Not Collective

2102:    Input Parameter:
2103: .  ts - the TS context obtained from TSCreate()

2105:    Output Parameter:
2106: .  steps - number of steps completed so far

2108:    Level: intermediate

2110: .seealso: `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2111: @*/
2112: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2113: {
2116:   *steps = ts->steps;
2117:   return 0;
2118: }

2120: /*@
2121:    TSSetStepNumber - Sets the number of steps completed.

2123:    Logically Collective on TS

2125:    Input Parameters:
2126: +  ts - the TS context
2127: -  steps - number of steps completed so far

2129:    Notes:
2130:    For most uses of the TS solvers the user need not explicitly call
2131:    TSSetStepNumber(), as the step counter is appropriately updated in
2132:    TSSolve()/TSStep()/TSRollBack(). Power users may call this routine to
2133:    reinitialize timestepping by setting the step counter to zero (and time
2134:    to the initial time) to solve a similar problem with different initial
2135:    conditions or parameters. Other possible use case is to continue
2136:    timestepping from a previously interrupted run in such a way that TS
2137:    monitors will be called with a initial nonzero step counter.

2139:    Level: advanced

2141: .seealso: `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2142: @*/
2143: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2144: {
2148:   ts->steps = steps;
2149:   return 0;
2150: }

2152: /*@
2153:    TSSetTimeStep - Allows one to reset the timestep at any time,
2154:    useful for simple pseudo-timestepping codes.

2156:    Logically Collective on TS

2158:    Input Parameters:
2159: +  ts - the TS context obtained from TSCreate()
2160: -  time_step - the size of the timestep

2162:    Level: intermediate

2164: .seealso: `TSGetTimeStep()`, `TSSetTime()`

2166: @*/
2167: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2168: {
2171:   ts->time_step = time_step;
2172:   return 0;
2173: }

2175: /*@
2176:    TSSetExactFinalTime - Determines whether to adapt the final time step to
2177:      match the exact final time, interpolate solution to the exact final time,
2178:      or just return at the final time TS computed.

2180:   Logically Collective on TS

2182:    Input Parameters:
2183: +   ts - the time-step context
2184: -   eftopt - exact final time option

2186: $  TS_EXACTFINALTIME_STEPOVER    - Don't do anything if final time is exceeded
2187: $  TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2188: $  TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time

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

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

2196:    Level: beginner

2198: .seealso: `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2199: @*/
2200: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2201: {
2204:   ts->exact_final_time = eftopt;
2205:   return 0;
2206: }

2208: /*@
2209:    TSGetExactFinalTime - Gets the exact final time option.

2211:    Not Collective

2213:    Input Parameter:
2214: .  ts - the TS context

2216:    Output Parameter:
2217: .  eftopt - exact final time option

2219:    Level: beginner

2221: .seealso: `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2222: @*/
2223: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2224: {
2227:   *eftopt = ts->exact_final_time;
2228:   return 0;
2229: }

2231: /*@
2232:    TSGetTimeStep - Gets the current timestep size.

2234:    Not Collective

2236:    Input Parameter:
2237: .  ts - the TS context obtained from TSCreate()

2239:    Output Parameter:
2240: .  dt - the current timestep size

2242:    Level: intermediate

2244: .seealso: `TSSetTimeStep()`, `TSGetTime()`

2246: @*/
2247: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2248: {
2251:   *dt = ts->time_step;
2252:   return 0;
2253: }

2255: /*@
2256:    TSGetSolution - Returns the solution at the present timestep. It
2257:    is valid to call this routine inside the function that you are evaluating
2258:    in order to move to the new timestep. This vector not changed until
2259:    the solution at the next timestep has been calculated.

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

2263:    Input Parameter:
2264: .  ts - the TS context obtained from TSCreate()

2266:    Output Parameter:
2267: .  v - the vector containing the solution

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

2272:    Level: intermediate

2274: .seealso: `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`

2276: @*/
2277: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2278: {
2281:   *v = ts->vec_sol;
2282:   return 0;
2283: }

2285: /*@
2286:    TSGetSolutionComponents - Returns any solution components at the present
2287:    timestep, if available for the time integration method being used.
2288:    Solution components are quantities that share the same size and
2289:    structure as the solution vector.

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

2293:    Parameters :
2294: +  ts - the TS context obtained from TSCreate() (input parameter).
2295: .  n - If v is PETSC_NULL, then the number of solution components is
2296:        returned through n, else the n-th solution component is
2297:        returned in v.
2298: -  v - the vector containing the n-th solution component
2299:        (may be PETSC_NULL to use this function to find out
2300:         the number of solutions components).

2302:    Level: advanced

2304: .seealso: `TSGetSolution()`

2306: @*/
2307: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2308: {
2310:   if (!ts->ops->getsolutioncomponents) *n = 0;
2311:   else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2312:   return 0;
2313: }

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

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

2321:    Parameters :
2322: +  ts - the TS context obtained from TSCreate() (input parameter).
2323: -  v - the vector containing the auxiliary solution

2325:    Level: intermediate

2327: .seealso: `TSGetSolution()`

2329: @*/
2330: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2331: {
2333:   if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2334:   else VecZeroEntries(*v);
2335:   return 0;
2336: }

2338: /*@
2339:    TSGetTimeError - Returns the estimated error vector, if the chosen
2340:    TSType has an error estimation functionality.

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

2344:    Note: MUST call after TSSetUp()

2346:    Parameters :
2347: +  ts - the TS context obtained from TSCreate() (input parameter).
2348: .  n - current estimate (n=0) or previous one (n=-1)
2349: -  v - the vector containing the error (same size as the solution).

2351:    Level: intermediate

2353: .seealso: `TSGetSolution()`, `TSSetTimeError()`

2355: @*/
2356: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2357: {
2359:   if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2360:   else VecZeroEntries(*v);
2361:   return 0;
2362: }

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

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

2371:    Parameters :
2372: +  ts - the TS context obtained from TSCreate() (input parameter).
2373: -  v - the vector containing the error (same size as the solution).

2375:    Level: intermediate

2377: .seealso: `TSSetSolution()`, `TSGetTimeError)`

2379: @*/
2380: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2381: {
2384:   PetscTryTypeMethod(ts, settimeerror, v);
2385:   return 0;
2386: }

2388: /* ----- Routines to initialize and destroy a timestepper ---- */
2389: /*@
2390:   TSSetProblemType - Sets the type of problem to be solved.

2392:   Not collective

2394:   Input Parameters:
2395: + ts   - The TS
2396: - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2397: .vb
2398:          U_t - A U = 0      (linear)
2399:          U_t - A(t) U = 0   (linear)
2400:          F(t,U,U_t) = 0     (nonlinear)
2401: .ve

2403:    Level: beginner

2405: .seealso: `TSSetUp()`, `TSProblemType`, `TS`
2406: @*/
2407: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2408: {
2410:   ts->problem_type = type;
2411:   if (type == TS_LINEAR) {
2412:     SNES snes;
2413:     TSGetSNES(ts, &snes);
2414:     SNESSetType(snes, SNESKSPONLY);
2415:   }
2416:   return 0;
2417: }

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

2422:   Not collective

2424:   Input Parameter:
2425: . ts   - The TS

2427:   Output Parameter:
2428: . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms
2429: .vb
2430:          M U_t = A U
2431:          M(t) U_t = A(t) U
2432:          F(t,U,U_t)
2433: .ve

2435:    Level: beginner

2437: .seealso: `TSSetUp()`, `TSProblemType`, `TS`
2438: @*/
2439: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2440: {
2443:   *type = ts->problem_type;
2444:   return 0;
2445: }

2447: /*
2448:     Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2449: */
2450: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2451: {
2452:   PetscBool isnone;

2454:   TSGetAdapt(ts, &ts->adapt);
2455:   TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type);

2457:   PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone);
2458:   if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2459:   else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2460:   return 0;
2461: }

2463: /*@
2464:    TSSetUp - Sets up the internal data structures for the later use of a timestepper.

2466:    Collective on TS

2468:    Input Parameter:
2469: .  ts - the TS context obtained from TSCreate()

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

2478:    Level: advanced

2480: .seealso: `TSCreate()`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2481: @*/
2482: PetscErrorCode TSSetUp(TS ts)
2483: {
2484:   DM dm;
2485:   PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2486:   PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2487:   TSIFunction   ifun;
2488:   TSIJacobian   ijac;
2489:   TSI2Jacobian  i2jac;
2490:   TSRHSJacobian rhsjac;

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

2495:   if (!((PetscObject)ts)->type_name) {
2496:     TSGetIFunction(ts, NULL, &ifun, NULL);
2497:     TSSetType(ts, ifun ? TSBEULER : TSEULER);
2498:   }

2500:   if (!ts->vec_sol) {
2502:     DMCreateGlobalVector(ts->dm, &ts->vec_sol);
2503:   }

2505:   if (ts->tspan) {
2506:     if (!ts->tspan->vecs_sol) VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol);
2507:   }
2508:   if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2509:     PetscObjectReference((PetscObject)ts->Jacprhs);
2510:     ts->Jacp = ts->Jacprhs;
2511:   }

2513:   if (ts->quadraturets) {
2514:     TSSetUp(ts->quadraturets);
2515:     VecDestroy(&ts->vec_costintegrand);
2516:     VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand);
2517:   }

2519:   TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL);
2520:   if (rhsjac == TSComputeRHSJacobianConstant) {
2521:     Mat  Amat, Pmat;
2522:     SNES snes;
2523:     TSGetSNES(ts, &snes);
2524:     SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL);
2525:     /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2526:      * have displaced the RHS matrix */
2527:     if (Amat && Amat == ts->Arhs) {
2528:       /* 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 */
2529:       MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat);
2530:       SNESSetJacobian(snes, Amat, NULL, NULL, NULL);
2531:       MatDestroy(&Amat);
2532:     }
2533:     if (Pmat && Pmat == ts->Brhs) {
2534:       MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat);
2535:       SNESSetJacobian(snes, NULL, Pmat, NULL, NULL);
2536:       MatDestroy(&Pmat);
2537:     }
2538:   }

2540:   TSGetAdapt(ts, &ts->adapt);
2541:   TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type);

2543:   PetscTryTypeMethod(ts, setup);

2545:   TSSetExactFinalTimeDefault(ts);

2547:   /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2548:      to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2549:    */
2550:   TSGetDM(ts, &dm);
2551:   DMSNESGetFunction(dm, &func, NULL);
2552:   if (!func) DMSNESSetFunction(dm, SNESTSFormFunction, ts);

2554:   /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2555:      Otherwise, the SNES will use coloring internally to form the Jacobian.
2556:    */
2557:   DMSNESGetJacobian(dm, &jac, NULL);
2558:   DMTSGetIJacobian(dm, &ijac, NULL);
2559:   DMTSGetI2Jacobian(dm, &i2jac, NULL);
2560:   DMTSGetRHSJacobian(dm, &rhsjac, NULL);
2561:   if (!jac && (ijac || i2jac || rhsjac)) DMSNESSetJacobian(dm, SNESTSFormJacobian, ts);

2563:   /* if time integration scheme has a starting method, call it */
2564:   PetscTryTypeMethod(ts, startingmethod);

2566:   ts->setupcalled = PETSC_TRUE;
2567:   return 0;
2568: }

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

2573:    Collective on TS

2575:    Input Parameter:
2576: .  ts - the TS context obtained from TSCreate()

2578:    Level: beginner

2580: .seealso: `TSCreate()`, `TSSetup()`, `TSDestroy()`
2581: @*/
2582: PetscErrorCode TSReset(TS ts)
2583: {
2584:   TS_RHSSplitLink ilink = ts->tsrhssplit, next;


2588:   PetscTryTypeMethod(ts, reset);
2589:   if (ts->snes) SNESReset(ts->snes);
2590:   if (ts->adapt) TSAdaptReset(ts->adapt);

2592:   MatDestroy(&ts->Arhs);
2593:   MatDestroy(&ts->Brhs);
2594:   VecDestroy(&ts->Frhs);
2595:   VecDestroy(&ts->vec_sol);
2596:   VecDestroy(&ts->vec_dot);
2597:   VecDestroy(&ts->vatol);
2598:   VecDestroy(&ts->vrtol);
2599:   VecDestroyVecs(ts->nwork, &ts->work);

2601:   MatDestroy(&ts->Jacprhs);
2602:   MatDestroy(&ts->Jacp);
2603:   if (ts->forward_solve) TSForwardReset(ts);
2604:   if (ts->quadraturets) {
2605:     TSReset(ts->quadraturets);
2606:     VecDestroy(&ts->vec_costintegrand);
2607:   }
2608:   while (ilink) {
2609:     next = ilink->next;
2610:     TSDestroy(&ilink->ts);
2611:     PetscFree(ilink->splitname);
2612:     ISDestroy(&ilink->is);
2613:     PetscFree(ilink);
2614:     ilink = next;
2615:   }
2616:   ts->tsrhssplit     = NULL;
2617:   ts->num_rhs_splits = 0;
2618:   if (ts->tspan) {
2619:     PetscFree(ts->tspan->span_times);
2620:     VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol);
2621:     PetscFree(ts->tspan);
2622:   }
2623:   ts->setupcalled = PETSC_FALSE;
2624:   return 0;
2625: }

2627: /*@C
2628:    TSDestroy - Destroys the timestepper context that was created
2629:    with TSCreate().

2631:    Collective on TS

2633:    Input Parameter:
2634: .  ts - the TS context obtained from TSCreate()

2636:    Level: beginner

2638: .seealso: `TSCreate()`, `TSSetUp()`, `TSSolve()`
2639: @*/
2640: PetscErrorCode TSDestroy(TS *ts)
2641: {
2642:   if (!*ts) return 0;
2644:   if (--((PetscObject)(*ts))->refct > 0) {
2645:     *ts = NULL;
2646:     return 0;
2647:   }

2649:   TSReset(*ts);
2650:   TSAdjointReset(*ts);
2651:   if ((*ts)->forward_solve) TSForwardReset(*ts);

2653:   /* if memory was published with SAWs then destroy it */
2654:   PetscObjectSAWsViewOff((PetscObject)*ts);
2655:   PetscTryTypeMethod((*ts), destroy);

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

2659:   TSAdaptDestroy(&(*ts)->adapt);
2660:   TSEventDestroy(&(*ts)->event);

2662:   SNESDestroy(&(*ts)->snes);
2663:   DMDestroy(&(*ts)->dm);
2664:   TSMonitorCancel((*ts));
2665:   TSAdjointMonitorCancel((*ts));

2667:   TSDestroy(&(*ts)->quadraturets);
2668:   PetscHeaderDestroy(ts);
2669:   return 0;
2670: }

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

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

2678:    Input Parameter:
2679: .  ts - the TS context obtained from TSCreate()

2681:    Output Parameter:
2682: .  snes - the nonlinear solver context

2684:    Notes:
2685:    The user can then directly manipulate the SNES context to set various
2686:    options, etc.  Likewise, the user can then extract and manipulate the
2687:    KSP, KSP, and PC contexts as well.

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

2692:    Level: beginner

2694: @*/
2695: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2696: {
2699:   if (!ts->snes) {
2700:     SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes);
2701:     PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options);
2702:     SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts);
2703:     PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1);
2704:     if (ts->dm) SNESSetDM(ts->snes, ts->dm);
2705:     if (ts->problem_type == TS_LINEAR) SNESSetType(ts->snes, SNESKSPONLY);
2706:   }
2707:   *snes = ts->snes;
2708:   return 0;
2709: }

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

2714:    Collective

2716:    Input Parameters:
2717: +  ts - the TS context obtained from TSCreate()
2718: -  snes - the nonlinear solver context

2720:    Notes:
2721:    Most users should have the TS created by calling TSGetSNES()

2723:    Level: developer

2725: @*/
2726: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2727: {
2728:   PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);

2732:   PetscObjectReference((PetscObject)snes);
2733:   SNESDestroy(&ts->snes);

2735:   ts->snes = snes;

2737:   SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts);
2738:   SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL);
2739:   if (func == SNESTSFormJacobian) SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts);
2740:   return 0;
2741: }

2743: /*@
2744:    TSGetKSP - Returns the KSP (linear solver) associated with
2745:    a TS (timestepper) context.

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

2749:    Input Parameter:
2750: .  ts - the TS context obtained from TSCreate()

2752:    Output Parameter:
2753: .  ksp - the nonlinear solver context

2755:    Notes:
2756:    The user can then directly manipulate the KSP context to set various
2757:    options, etc.  Likewise, the user can then extract and manipulate the
2758:    KSP and PC contexts as well.

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

2763:    Level: beginner

2765: @*/
2766: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2767: {
2768:   SNES snes;

2774:   TSGetSNES(ts, &snes);
2775:   SNESGetKSP(snes, ksp);
2776:   return 0;
2777: }

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

2781: /*@
2782:    TSSetMaxSteps - Sets the maximum number of steps to use.

2784:    Logically Collective on TS

2786:    Input Parameters:
2787: +  ts - the TS context obtained from TSCreate()
2788: -  maxsteps - maximum number of steps to use

2790:    Options Database Keys:
2791: .  -ts_max_steps <maxsteps> - Sets maxsteps

2793:    Notes:
2794:    The default maximum number of steps is 5000

2796:    Level: intermediate

2798: .seealso: `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2799: @*/
2800: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2801: {
2805:   ts->max_steps = maxsteps;
2806:   return 0;
2807: }

2809: /*@
2810:    TSGetMaxSteps - Gets the maximum number of steps to use.

2812:    Not Collective

2814:    Input Parameters:
2815: .  ts - the TS context obtained from TSCreate()

2817:    Output Parameter:
2818: .  maxsteps - maximum number of steps to use

2820:    Level: advanced

2822: .seealso: `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2823: @*/
2824: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2825: {
2828:   *maxsteps = ts->max_steps;
2829:   return 0;
2830: }

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

2835:    Logically Collective on TS

2837:    Input Parameters:
2838: +  ts - the TS context obtained from TSCreate()
2839: -  maxtime - final time to step to

2841:    Options Database Keys:
2842: .  -ts_max_time <maxtime> - Sets maxtime

2844:    Notes:
2845:    The default maximum time is 5.0

2847:    Level: intermediate

2849: .seealso: `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2850: @*/
2851: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2852: {
2855:   ts->max_time = maxtime;
2856:   return 0;
2857: }

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

2862:    Not Collective

2864:    Input Parameters:
2865: .  ts - the TS context obtained from TSCreate()

2867:    Output Parameter:
2868: .  maxtime - final time to step to

2870:    Level: advanced

2872: .seealso: `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2873: @*/
2874: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2875: {
2878:   *maxtime = ts->max_time;
2879:   return 0;
2880: }

2882: /*@
2883:    TSSetInitialTimeStep - Deprecated, use TSSetTime() and TSSetTimeStep().

2885:    Level: deprecated

2887: @*/
2888: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2889: {
2891:   TSSetTime(ts, initial_time);
2892:   TSSetTimeStep(ts, time_step);
2893:   return 0;
2894: }

2896: /*@
2897:    TSGetDuration - Deprecated, use TSGetMaxSteps() and TSGetMaxTime().

2899:    Level: deprecated

2901: @*/
2902: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2903: {
2905:   if (maxsteps) {
2907:     *maxsteps = ts->max_steps;
2908:   }
2909:   if (maxtime) {
2911:     *maxtime = ts->max_time;
2912:   }
2913:   return 0;
2914: }

2916: /*@
2917:    TSSetDuration - Deprecated, use TSSetMaxSteps() and TSSetMaxTime().

2919:    Level: deprecated

2921: @*/
2922: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2923: {
2927:   if (maxsteps >= 0) ts->max_steps = maxsteps;
2928:   if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime;
2929:   return 0;
2930: }

2932: /*@
2933:    TSGetTimeStepNumber - Deprecated, use TSGetStepNumber().

2935:    Level: deprecated

2937: @*/
2938: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
2939: {
2940:   return TSGetStepNumber(ts, steps);
2941: }

2943: /*@
2944:    TSGetTotalSteps - Deprecated, use TSGetStepNumber().

2946:    Level: deprecated

2948: @*/
2949: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
2950: {
2951:   return TSGetStepNumber(ts, steps);
2952: }

2954: /*@
2955:    TSSetSolution - Sets the initial solution vector
2956:    for use by the TS routines.

2958:    Logically Collective on TS

2960:    Input Parameters:
2961: +  ts - the TS context obtained from TSCreate()
2962: -  u - the solution vector

2964:    Level: beginner

2966: .seealso: `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
2967: @*/
2968: PetscErrorCode TSSetSolution(TS ts, Vec u)
2969: {
2970:   DM dm;

2974:   PetscObjectReference((PetscObject)u);
2975:   VecDestroy(&ts->vec_sol);
2976:   ts->vec_sol = u;

2978:   TSGetDM(ts, &dm);
2979:   DMShellSetGlobalVector(dm, u);
2980:   return 0;
2981: }

2983: /*@C
2984:   TSSetPreStep - Sets the general-purpose function
2985:   called once at the beginning of each time step.

2987:   Logically Collective on TS

2989:   Input Parameters:
2990: + ts   - The TS context obtained from TSCreate()
2991: - func - The function

2993:   Calling sequence of func:
2994: .vb
2995:   PetscErrorCode func (TS ts);
2996: .ve

2998:   Level: intermediate

3000: .seealso: `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3001: @*/
3002: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
3003: {
3005:   ts->prestep = func;
3006:   return 0;
3007: }

3009: /*@
3010:   TSPreStep - Runs the user-defined pre-step function.

3012:   Collective on TS

3014:   Input Parameters:
3015: . ts   - The TS context obtained from TSCreate()

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

3021:   Level: developer

3023: .seealso: `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3024: @*/
3025: PetscErrorCode TSPreStep(TS ts)
3026: {
3028:   if (ts->prestep) {
3029:     Vec              U;
3030:     PetscObjectId    idprev;
3031:     PetscBool        sameObject;
3032:     PetscObjectState sprev, spost;

3034:     TSGetSolution(ts, &U);
3035:     PetscObjectGetId((PetscObject)U, &idprev);
3036:     PetscObjectStateGet((PetscObject)U, &sprev);
3037:     PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3038:     TSGetSolution(ts, &U);
3039:     PetscObjectCompareId((PetscObject)U, idprev, &sameObject);
3040:     PetscObjectStateGet((PetscObject)U, &spost);
3041:     if (!sameObject || sprev != spost) TSRestartStep(ts);
3042:   }
3043:   return 0;
3044: }

3046: /*@C
3047:   TSSetPreStage - Sets the general-purpose function
3048:   called once at the beginning of each stage.

3050:   Logically Collective on TS

3052:   Input Parameters:
3053: + ts   - The TS context obtained from TSCreate()
3054: - func - The function

3056:   Calling sequence of func:
3057: .vb
3058:   PetscErrorCode func(TS ts, PetscReal stagetime);
3059: .ve

3061:   Level: intermediate

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

3068: .seealso: `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3069: @*/
3070: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS, PetscReal))
3071: {
3073:   ts->prestage = func;
3074:   return 0;
3075: }

3077: /*@C
3078:   TSSetPostStage - Sets the general-purpose function
3079:   called once at the end of each stage.

3081:   Logically Collective on TS

3083:   Input Parameters:
3084: + ts   - The TS context obtained from TSCreate()
3085: - func - The function

3087:   Calling sequence of func:
3088: .vb
3089:   PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y);
3090: .ve

3092:   Level: intermediate

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

3099: .seealso: `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3100: @*/
3101: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS, PetscReal, PetscInt, Vec *))
3102: {
3104:   ts->poststage = func;
3105:   return 0;
3106: }

3108: /*@C
3109:   TSSetPostEvaluate - Sets the general-purpose function
3110:   called once at the end of each step evaluation.

3112:   Logically Collective on TS

3114:   Input Parameters:
3115: + ts   - The TS context obtained from TSCreate()
3116: - func - The function

3118:   Calling sequence of func:
3119: .vb
3120:   PetscErrorCode func(TS ts);
3121: .ve

3123:   Level: intermediate

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

3132: .seealso: `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3133: @*/
3134: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS))
3135: {
3137:   ts->postevaluate = func;
3138:   return 0;
3139: }

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

3144:   Collective on TS

3146:   Input Parameters:
3147: . ts          - The TS context obtained from TSCreate()
3148:   stagetime   - The absolute time of the current stage

3150:   Notes:
3151:   TSPreStage() is typically used within time stepping implementations,
3152:   most users would not generally call this routine themselves.

3154:   Level: developer

3156: .seealso: `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3157: @*/
3158: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3159: {
3161:   if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3162:   return 0;
3163: }

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

3168:   Collective on TS

3170:   Input Parameters:
3171: . ts          - The TS context obtained from TSCreate()
3172:   stagetime   - The absolute time of the current stage
3173:   stageindex  - Stage number
3174:   Y           - Array of vectors (of size = total number
3175:                 of stages) with the stage solutions

3177:   Notes:
3178:   TSPostStage() is typically used within time stepping implementations,
3179:   most users would not generally call this routine themselves.

3181:   Level: developer

3183: .seealso: `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3184: @*/
3185: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3186: {
3188:   if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3189:   return 0;
3190: }

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

3195:   Collective on TS

3197:   Input Parameters:
3198: . ts          - The TS context obtained from TSCreate()

3200:   Notes:
3201:   TSPostEvaluate() is typically used within time stepping implementations,
3202:   most users would not generally call this routine themselves.

3204:   Level: developer

3206: .seealso: `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3207: @*/
3208: PetscErrorCode TSPostEvaluate(TS ts)
3209: {
3211:   if (ts->postevaluate) {
3212:     Vec              U;
3213:     PetscObjectState sprev, spost;

3215:     TSGetSolution(ts, &U);
3216:     PetscObjectStateGet((PetscObject)U, &sprev);
3217:     PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3218:     PetscObjectStateGet((PetscObject)U, &spost);
3219:     if (sprev != spost) TSRestartStep(ts);
3220:   }
3221:   return 0;
3222: }

3224: /*@C
3225:   TSSetPostStep - Sets the general-purpose function
3226:   called once at the end of each time step.

3228:   Logically Collective on TS

3230:   Input Parameters:
3231: + ts   - The TS context obtained from TSCreate()
3232: - func - The function

3234:   Calling sequence of func:
3235: $ func (TS ts);

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

3242:   Level: intermediate

3244: .seealso: `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3245: @*/
3246: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
3247: {
3249:   ts->poststep = func;
3250:   return 0;
3251: }

3253: /*@
3254:   TSPostStep - Runs the user-defined post-step function.

3256:   Collective on TS

3258:   Input Parameters:
3259: . ts   - The TS context obtained from TSCreate()

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

3265:   Level: developer

3267: @*/
3268: PetscErrorCode TSPostStep(TS ts)
3269: {
3271:   if (ts->poststep) {
3272:     Vec              U;
3273:     PetscObjectId    idprev;
3274:     PetscBool        sameObject;
3275:     PetscObjectState sprev, spost;

3277:     TSGetSolution(ts, &U);
3278:     PetscObjectGetId((PetscObject)U, &idprev);
3279:     PetscObjectStateGet((PetscObject)U, &sprev);
3280:     PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3281:     TSGetSolution(ts, &U);
3282:     PetscObjectCompareId((PetscObject)U, idprev, &sameObject);
3283:     PetscObjectStateGet((PetscObject)U, &spost);
3284:     if (!sameObject || sprev != spost) TSRestartStep(ts);
3285:   }
3286:   return 0;
3287: }

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

3292:    Collective on TS

3294:    Input Parameters:
3295: +  ts - time stepping context
3296: -  t - time to interpolate to

3298:    Output Parameter:
3299: .  U - state at given time

3301:    Level: intermediate

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

3306: .seealso: `TSSetExactFinalTime()`, `TSSolve()`
3307: @*/
3308: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3309: {
3313:   PetscUseTypeMethod(ts, interpolate, t, U);
3314:   return 0;
3315: }

3317: /*@
3318:    TSStep - Steps one time step

3320:    Collective on TS

3322:    Input Parameter:
3323: .  ts - the TS context obtained from TSCreate()

3325:    Level: developer

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

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

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

3336: .seealso: `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3337: @*/
3338: PetscErrorCode TSStep(TS ts)
3339: {
3340:   static PetscBool cite = PETSC_FALSE;
3341:   PetscReal        ptime;

3344:   PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3345:                                    "  title         = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3346:                                    "  author        = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3347:                                    "  journal       = {arXiv e-preprints},\n"
3348:                                    "  eprint        = {1806.01437},\n"
3349:                                    "  archivePrefix = {arXiv},\n"
3350:                                    "  year          = {2018}\n}\n",
3351:                                    &cite));
3352:   TSSetUp(ts);
3353:   TSTrajectorySetUp(ts->trajectory, ts);


3359:   if (!ts->steps) ts->ptime_prev = ts->ptime;
3360:   ptime                   = ts->ptime;
3361:   ts->ptime_prev_rollback = ts->ptime_prev;
3362:   ts->reason              = TS_CONVERGED_ITERATING;

3364:   PetscLogEventBegin(TS_Step, ts, 0, 0, 0);
3365:   PetscUseTypeMethod(ts, step);
3366:   PetscLogEventEnd(TS_Step, ts, 0, 0, 0);

3368:   if (ts->tspan && PetscIsCloseAtTol(ts->ptime, ts->tspan->span_times[ts->tspan->spanctr], ts->tspan->reltol * ts->time_step + ts->tspan->abstol, 0) && ts->tspan->spanctr < ts->tspan->num_span_times)
3369:     VecCopy(ts->vec_sol, ts->tspan->vecs_sol[ts->tspan->spanctr++]);
3370:   if (ts->reason >= 0) {
3371:     ts->ptime_prev = ptime;
3372:     ts->steps++;
3373:     ts->steprollback = PETSC_FALSE;
3374:     ts->steprestart  = PETSC_FALSE;
3375:   }
3376:   if (!ts->reason) {
3377:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3378:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3379:   }

3383:   return 0;
3384: }

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

3390:    Collective on TS

3392:    Input Parameters:
3393: +  ts - time stepping context
3394: -  wnormtype - norm type, either NORM_2 or NORM_INFINITY

3396:    Input/Output Parameter:
3397: .  order - optional, desired order for the error evaluation or PETSC_DECIDE;
3398:            on output, the actual order of the error evaluation

3400:    Output Parameter:
3401: .  wlte - the weighted local truncation error norm

3403:    Level: advanced

3405:    Notes:
3406:    If the timestepper cannot evaluate the error in a particular step
3407:    (eg. in the first step or restart steps after event handling),
3408:    this routine returns wlte=-1.0 .

3410: .seealso: `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3411: @*/
3412: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3413: {
3421:   PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3422:   return 0;
3423: }

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

3428:    Collective on TS

3430:    Input Parameters:
3431: +  ts - time stepping context
3432: .  order - desired order of accuracy
3433: -  done - whether the step was evaluated at this order (pass NULL to generate an error if not available)

3435:    Output Parameter:
3436: .  U - state at the end of the current step

3438:    Level: advanced

3440:    Notes:
3441:    This function cannot be called until all stages have been evaluated.
3442:    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.

3444: .seealso: `TSStep()`, `TSAdapt`
3445: @*/
3446: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3447: {
3451:   PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3452:   return 0;
3453: }

3455: /*@C
3456:   TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.

3458:   Not collective

3460:   Input Parameter:
3461: . ts        - time stepping context

3463:   Output Parameter:
3464: . initConditions - The function which computes an initial condition

3466:    Level: advanced

3468:    Notes:
3469:    The calling sequence for the function is
3470: $ initCondition(TS ts, Vec u)
3471: $ ts - The timestepping context
3472: $ u  - The input vector in which the initial condition is stored

3474: .seealso: `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3475: @*/
3476: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS, Vec))
3477: {
3480:   *initCondition = ts->ops->initcondition;
3481:   return 0;
3482: }

3484: /*@C
3485:   TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.

3487:   Logically collective on ts

3489:   Input Parameters:
3490: + ts        - time stepping context
3491: - initCondition - The function which computes an initial condition

3493:   Level: advanced

3495:   Calling sequence for initCondition:
3496: $ PetscErrorCode initCondition(TS ts, Vec u)

3498: + ts - The timestepping context
3499: - u  - The input vector in which the initial condition is to be stored

3501: .seealso: `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3502: @*/
3503: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS, Vec))
3504: {
3507:   ts->ops->initcondition = initCondition;
3508:   return 0;
3509: }

3511: /*@
3512:   TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set.

3514:   Collective on ts

3516:   Input Parameters:
3517: + ts - time stepping context
3518: - u  - The Vec to store the condition in which will be used in TSSolve()

3520:   Level: advanced

3522: .seealso: `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3523: @*/
3524: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3525: {
3528:   PetscTryTypeMethod(ts, initcondition, u);
3529:   return 0;
3530: }

3532: /*@C
3533:   TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.

3535:   Not collective

3537:   Input Parameter:
3538: . ts         - time stepping context

3540:   Output Parameter:
3541: . exactError - The function which computes the solution error

3543:   Level: advanced

3545:   Calling sequence for exactError:
3546: $ PetscErrorCode exactError(TS ts, Vec u)

3548: + ts - The timestepping context
3549: . u  - The approximate solution vector
3550: - e  - The input vector in which the error is stored

3552: .seealso: `TSGetComputeExactError()`, `TSComputeExactError()`
3553: @*/
3554: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS, Vec, Vec))
3555: {
3558:   *exactError = ts->ops->exacterror;
3559:   return 0;
3560: }

3562: /*@C
3563:   TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.

3565:   Logically collective on ts

3567:   Input Parameters:
3568: + ts         - time stepping context
3569: - exactError - The function which computes the solution error

3571:   Level: advanced

3573:   Calling sequence for exactError:
3574: $ PetscErrorCode exactError(TS ts, Vec u)

3576: + ts - The timestepping context
3577: . u  - The approximate solution vector
3578: - e  - The input vector in which the error is stored

3580: .seealso: `TSGetComputeExactError()`, `TSComputeExactError()`
3581: @*/
3582: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS, Vec, Vec))
3583: {
3586:   ts->ops->exacterror = exactError;
3587:   return 0;
3588: }

3590: /*@
3591:   TSComputeExactError - Compute the solution error for the timestepping using the function previously set.

3593:   Collective on ts

3595:   Input Parameters:
3596: + ts - time stepping context
3597: . u  - The approximate solution
3598: - e  - The Vec used to store the error

3600:   Level: advanced

3602: .seealso: `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3603: @*/
3604: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3605: {
3609:   PetscTryTypeMethod(ts, exacterror, u, e);
3610:   return 0;
3611: }

3613: /*@
3614:    TSSolve - Steps the requested number of timesteps.

3616:    Collective on TS

3618:    Input Parameters:
3619: +  ts - the TS context obtained from TSCreate()
3620: -  u - the solution vector  (can be null if TSSetSolution() was used and TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP) was not used,
3621:                              otherwise must contain the initial conditions and will contain the solution at the final requested time

3623:    Level: beginner

3625:    Notes:
3626:    The final time returned by this function may be different from the time of the internally
3627:    held state accessible by TSGetSolution() and TSGetTime() because the method may have
3628:    stepped over the final time.

3630: .seealso: `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3631: @*/
3632: PetscErrorCode TSSolve(TS ts, Vec u)
3633: {
3634:   Vec solution;


3639:   TSSetExactFinalTimeDefault(ts);
3640:   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 */
3641:     if (!ts->vec_sol || u == ts->vec_sol) {
3642:       VecDuplicate(u, &solution);
3643:       TSSetSolution(ts, solution);
3644:       VecDestroy(&solution); /* grant ownership */
3645:     }
3646:     VecCopy(u, ts->vec_sol);
3648:   } else if (u) TSSetSolution(ts, u);
3649:   TSSetUp(ts);
3650:   TSTrajectorySetUp(ts->trajectory, ts);


3657:   if (ts->tspan && PetscIsCloseAtTol(ts->ptime, ts->tspan->span_times[0], ts->tspan->reltol * ts->time_step + ts->tspan->abstol, 0)) { /* starting point in time span */
3658:     VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]);
3659:     ts->tspan->spanctr = 1;
3660:   }

3662:   if (ts->forward_solve) TSForwardSetUp(ts);

3664:   /* reset number of steps only when the step is not restarted. ARKIMEX
3665:      restarts the step after an event. Resetting these counters in such case causes
3666:      TSTrajectory to incorrectly save the output files
3667:   */
3668:   /* reset time step and iteration counters */
3669:   if (!ts->steps) {
3670:     ts->ksp_its           = 0;
3671:     ts->snes_its          = 0;
3672:     ts->num_snes_failures = 0;
3673:     ts->reject            = 0;
3674:     ts->steprestart       = PETSC_TRUE;
3675:     ts->steprollback      = PETSC_FALSE;
3676:     ts->rhsjacobian.time  = PETSC_MIN_REAL;
3677:   }

3679:   /* make sure initial time step does not overshoot final time or the next point in tspan */
3680:   if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
3681:     PetscReal maxdt;
3682:     PetscReal dt = ts->time_step;

3684:     if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
3685:     else maxdt = ts->max_time - ts->ptime;
3686:     ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
3687:   }
3688:   ts->reason = TS_CONVERGED_ITERATING;

3690:   {
3691:     PetscViewer       viewer;
3692:     PetscViewerFormat format;
3693:     PetscBool         flg;
3694:     static PetscBool  incall = PETSC_FALSE;

3696:     if (!incall) {
3697:       /* Estimate the convergence rate of the time discretization */
3698:       PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg);
3699:       if (flg) {
3700:         PetscConvEst conv;
3701:         DM           dm;
3702:         PetscReal   *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
3703:         PetscInt     Nf;
3704:         PetscBool    checkTemporal = PETSC_TRUE;

3706:         incall = PETSC_TRUE;
3707:         PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg);
3708:         TSGetDM(ts, &dm);
3709:         DMGetNumFields(dm, &Nf);
3710:         PetscCalloc1(PetscMax(Nf, 1), &alpha);
3711:         PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv);
3712:         PetscConvEstUseTS(conv, checkTemporal);
3713:         PetscConvEstSetSolver(conv, (PetscObject)ts);
3714:         PetscConvEstSetFromOptions(conv);
3715:         PetscConvEstSetUp(conv);
3716:         PetscConvEstGetConvRate(conv, alpha);
3717:         PetscViewerPushFormat(viewer, format);
3718:         PetscConvEstRateView(conv, alpha, viewer);
3719:         PetscViewerPopFormat(viewer);
3720:         PetscViewerDestroy(&viewer);
3721:         PetscConvEstDestroy(&conv);
3722:         PetscFree(alpha);
3723:         incall = PETSC_FALSE;
3724:       }
3725:     }
3726:   }

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

3730:   if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
3731:     PetscUseTypeMethod(ts, solve);
3732:     if (u) VecCopy(ts->vec_sol, u);
3733:     ts->solvetime = ts->ptime;
3734:     solution      = ts->vec_sol;
3735:   } else { /* Step the requested number of timesteps. */
3736:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3737:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;

3739:     if (!ts->steps) {
3740:       TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol);
3741:       TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol);
3742:     }

3744:     while (!ts->reason) {
3745:       TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);
3746:       if (!ts->steprollback) TSPreStep(ts);
3747:       TSStep(ts);
3748:       if (ts->testjacobian) TSRHSJacobianTest(ts, NULL);
3749:       if (ts->testjacobiantranspose) TSRHSJacobianTestTranspose(ts, NULL);
3750:       if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
3751:         if (ts->reason >= 0) ts->steps--;            /* Revert the step number changed by TSStep() */
3752:         TSForwardCostIntegral(ts);
3753:         if (ts->reason >= 0) ts->steps++;
3754:       }
3755:       if (ts->forward_solve) {            /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
3756:         if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
3757:         TSForwardStep(ts);
3758:         if (ts->reason >= 0) ts->steps++;
3759:       }
3760:       TSPostEvaluate(ts);
3761:       TSEventHandler(ts); /* The right-hand side may be changed due to event. Be careful with Any computation using the RHS information after this point. */
3762:       if (ts->steprollback) TSPostEvaluate(ts);
3763:       if (!ts->steprollback) {
3764:         TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol);
3765:         TSPostStep(ts);
3766:       }
3767:     }
3768:     TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);

3770:     if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
3771:       TSInterpolate(ts, ts->max_time, u);
3772:       ts->solvetime = ts->max_time;
3773:       solution      = u;
3774:       TSMonitor(ts, -1, ts->solvetime, solution);
3775:     } else {
3776:       if (u) VecCopy(ts->vec_sol, u);
3777:       ts->solvetime = ts->ptime;
3778:       solution      = ts->vec_sol;
3779:     }
3780:   }

3782:   TSViewFromOptions(ts, NULL, "-ts_view");
3783:   VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution");
3784:   PetscObjectSAWsBlock((PetscObject)ts);
3785:   if (ts->adjoint_solve) TSAdjointSolve(ts);
3786:   return 0;
3787: }

3789: /*@
3790:    TSGetTime - Gets the time of the most recently completed step.

3792:    Not Collective

3794:    Input Parameter:
3795: .  ts - the TS context obtained from TSCreate()

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

3800:    Level: beginner

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

3806: .seealso: `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`

3808: @*/
3809: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
3810: {
3813:   *t = ts->ptime;
3814:   return 0;
3815: }

3817: /*@
3818:    TSGetPrevTime - Gets the starting time of the previously completed step.

3820:    Not Collective

3822:    Input Parameter:
3823: .  ts - the TS context obtained from TSCreate()

3825:    Output Parameter:
3826: .  t  - the previous time

3828:    Level: beginner

3830: .seealso: `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`

3832: @*/
3833: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
3834: {
3837:   *t = ts->ptime_prev;
3838:   return 0;
3839: }

3841: /*@
3842:    TSSetTime - Allows one to reset the time.

3844:    Logically Collective on TS

3846:    Input Parameters:
3847: +  ts - the TS context obtained from TSCreate()
3848: -  time - the time

3850:    Level: intermediate

3852: .seealso: `TSGetTime()`, `TSSetMaxSteps()`

3854: @*/
3855: PetscErrorCode TSSetTime(TS ts, PetscReal t)
3856: {
3859:   ts->ptime = t;
3860:   return 0;
3861: }

3863: /*@C
3864:    TSSetOptionsPrefix - Sets the prefix used for searching for all
3865:    TS options in the database.

3867:    Logically Collective on TS

3869:    Input Parameters:
3870: +  ts     - The TS context
3871: -  prefix - The prefix to prepend to all option names

3873:    Notes:
3874:    A hyphen (-) must NOT be given at the beginning of the prefix name.
3875:    The first character of all runtime options is AUTOMATICALLY the
3876:    hyphen.

3878:    Level: advanced

3880: .seealso: `TSSetFromOptions()`

3882: @*/
3883: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
3884: {
3885:   SNES snes;

3888:   PetscObjectSetOptionsPrefix((PetscObject)ts, prefix);
3889:   TSGetSNES(ts, &snes);
3890:   SNESSetOptionsPrefix(snes, prefix);
3891:   return 0;
3892: }

3894: /*@C
3895:    TSAppendOptionsPrefix - Appends to the prefix used for searching for all
3896:    TS options in the database.

3898:    Logically Collective on TS

3900:    Input Parameters:
3901: +  ts     - The TS context
3902: -  prefix - The prefix to prepend to all option names

3904:    Notes:
3905:    A hyphen (-) must NOT be given at the beginning of the prefix name.
3906:    The first character of all runtime options is AUTOMATICALLY the
3907:    hyphen.

3909:    Level: advanced

3911: .seealso: `TSGetOptionsPrefix()`

3913: @*/
3914: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
3915: {
3916:   SNES snes;

3919:   PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix);
3920:   TSGetSNES(ts, &snes);
3921:   SNESAppendOptionsPrefix(snes, prefix);
3922:   return 0;
3923: }

3925: /*@C
3926:    TSGetOptionsPrefix - Sets the prefix used for searching for all
3927:    TS options in the database.

3929:    Not Collective

3931:    Input Parameter:
3932: .  ts - The TS context

3934:    Output Parameter:
3935: .  prefix - A pointer to the prefix string used

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

3941:    Level: intermediate

3943: .seealso: `TSAppendOptionsPrefix()`
3944: @*/
3945: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
3946: {
3949:   PetscObjectGetOptionsPrefix((PetscObject)ts, prefix);
3950:   return 0;
3951: }

3953: /*@C
3954:    TSGetRHSJacobian - Returns the Jacobian J at the present timestep.

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

3958:    Input Parameter:
3959: .  ts  - The TS context obtained from TSCreate()

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

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

3970:    Level: intermediate

3972: .seealso: `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`

3974: @*/
3975: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobian *func, void **ctx)
3976: {
3977:   DM dm;

3979:   if (Amat || Pmat) {
3980:     SNES snes;
3981:     TSGetSNES(ts, &snes);
3982:     SNESSetUpMatrices(snes);
3983:     SNESGetJacobian(snes, Amat, Pmat, NULL, NULL);
3984:   }
3985:   TSGetDM(ts, &dm);
3986:   DMTSGetRHSJacobian(dm, func, ctx);
3987:   return 0;
3988: }

3990: /*@C
3991:    TSGetIJacobian - Returns the implicit Jacobian at the present timestep.

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

3995:    Input Parameter:
3996: .  ts  - The TS context obtained from TSCreate()

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

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

4007:    Level: advanced

4009: .seealso: `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`

4011: @*/
4012: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobian *f, void **ctx)
4013: {
4014:   DM dm;

4016:   if (Amat || Pmat) {
4017:     SNES snes;
4018:     TSGetSNES(ts, &snes);
4019:     SNESSetUpMatrices(snes);
4020:     SNESGetJacobian(snes, Amat, Pmat, NULL, NULL);
4021:   }
4022:   TSGetDM(ts, &dm);
4023:   DMTSGetIJacobian(dm, f, ctx);
4024:   return 0;
4025: }

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

4031:    Logically Collective on ts

4033:    Input Parameters:
4034: +  ts - the ODE integrator object
4035: -  dm - the dm, cannot be NULL

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

4042:    Level: intermediate

4044: .seealso: `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4045: @*/
4046: PetscErrorCode TSSetDM(TS ts, DM dm)
4047: {
4048:   SNES snes;
4049:   DMTS tsdm;

4053:   PetscObjectReference((PetscObject)dm);
4054:   if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4055:     if (ts->dm->dmts && !dm->dmts) {
4056:       DMCopyDMTS(ts->dm, dm);
4057:       DMGetDMTS(ts->dm, &tsdm);
4058:       /* Grant write privileges to the replacement DM */
4059:       if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4060:     }
4061:     DMDestroy(&ts->dm);
4062:   }
4063:   ts->dm = dm;

4065:   TSGetSNES(ts, &snes);
4066:   SNESSetDM(snes, dm);
4067:   return 0;
4068: }

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

4073:    Not Collective

4075:    Input Parameter:
4076: . ts - the preconditioner context

4078:    Output Parameter:
4079: .  dm - the dm

4081:    Level: intermediate

4083: .seealso: `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4084: @*/
4085: PetscErrorCode TSGetDM(TS ts, DM *dm)
4086: {
4088:   if (!ts->dm) {
4089:     DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm);
4090:     if (ts->snes) SNESSetDM(ts->snes, ts->dm);
4091:   }
4092:   *dm = ts->dm;
4093:   return 0;
4094: }

4096: /*@
4097:    SNESTSFormFunction - Function to evaluate nonlinear residual

4099:    Logically Collective on SNES

4101:    Input Parameters:
4102: + snes - nonlinear solver
4103: . U - the current state at which to evaluate the residual
4104: - ctx - user context, must be a TS

4106:    Output Parameter:
4107: . F - the nonlinear residual

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

4113:    Level: advanced

4115: .seealso: `SNESSetFunction()`, `MatFDColoringSetFunction()`
4116: @*/
4117: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4118: {
4119:   TS ts = (TS)ctx;

4125:   (ts->ops->snesfunction)(snes, U, F, ts);
4126:   return 0;
4127: }

4129: /*@
4130:    SNESTSFormJacobian - Function to evaluate the Jacobian

4132:    Collective on SNES

4134:    Input Parameters:
4135: + snes - nonlinear solver
4136: . U - the current state at which to evaluate the residual
4137: - ctx - user context, must be a TS

4139:    Output Parameters:
4140: + A - the Jacobian
4141: - B - the preconditioning matrix (may be the same as A)

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

4146:    Level: developer

4148: .seealso: `SNESSetJacobian()`
4149: @*/
4150: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4151: {
4152:   TS ts = (TS)ctx;

4161:   (ts->ops->snesjacobian)(snes, U, A, B, ts);
4162:   return 0;
4163: }

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

4168:    Collective on TS

4170:    Input Parameters:
4171: +  ts - time stepping context
4172: .  t - time at which to evaluate
4173: .  U - state at which to evaluate
4174: -  ctx - context

4176:    Output Parameter:
4177: .  F - right hand side

4179:    Level: intermediate

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

4185: .seealso: `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4186: @*/
4187: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4188: {
4189:   Mat Arhs, Brhs;

4191:   TSGetRHSMats_Private(ts, &Arhs, &Brhs);
4192:   /* undo the damage caused by shifting */
4193:   TSRecoverRHSJacobian(ts, Arhs, Brhs);
4194:   TSComputeRHSJacobian(ts, t, U, Arhs, Brhs);
4195:   MatMult(Arhs, U, F);
4196:   return 0;
4197: }

4199: /*@C
4200:    TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.

4202:    Collective on TS

4204:    Input Parameters:
4205: +  ts - time stepping context
4206: .  t - time at which to evaluate
4207: .  U - state at which to evaluate
4208: -  ctx - context

4210:    Output Parameters:
4211: +  A - pointer to operator
4212: -  B - pointer to preconditioning matrix

4214:    Level: intermediate

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

4219: .seealso: `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4220: @*/
4221: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4222: {
4223:   return 0;
4224: }

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

4229:    Collective on TS

4231:    Input Parameters:
4232: +  ts - time stepping context
4233: .  t - time at which to evaluate
4234: .  U - state at which to evaluate
4235: .  Udot - time derivative of state vector
4236: -  ctx - context

4238:    Output Parameter:
4239: .  F - left hand side

4241:    Level: intermediate

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

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

4251: .seealso: `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4252: @*/
4253: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4254: {
4255:   Mat A, B;

4257:   TSGetIJacobian(ts, &A, &B, NULL, NULL);
4258:   TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE);
4259:   MatMult(A, Udot, F);
4260:   return 0;
4261: }

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

4266:    Collective on TS

4268:    Input Parameters:
4269: +  ts - time stepping context
4270: .  t - time at which to evaluate
4271: .  U - state at which to evaluate
4272: .  Udot - time derivative of state vector
4273: .  shift - shift to apply
4274: -  ctx - context

4276:    Output Parameters:
4277: +  A - pointer to operator
4278: -  B - pointer to preconditioning matrix

4280:    Level: advanced

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

4285:    It is only appropriate for problems of the form

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

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

4293: $    shift*M + J

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

4298: .seealso: `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4299: @*/
4300: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4301: {
4302:   MatScale(A, shift / ts->ijacobian.shift);
4303:   ts->ijacobian.shift = shift;
4304:   return 0;
4305: }

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

4310:    Not Collective

4312:    Input Parameter:
4313: .  ts - the TS context

4315:    Output Parameter:
4316: .  equation_type - see TSEquationType

4318:    Level: beginner

4320: .seealso: `TSSetEquationType()`, `TSEquationType`
4321: @*/
4322: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4323: {
4326:   *equation_type = ts->equation_type;
4327:   return 0;
4328: }

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

4333:    Not Collective

4335:    Input Parameters:
4336: +  ts - the TS context
4337: -  equation_type - see TSEquationType

4339:    Level: advanced

4341: .seealso: `TSGetEquationType()`, `TSEquationType`
4342: @*/
4343: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4344: {
4346:   ts->equation_type = equation_type;
4347:   return 0;
4348: }

4350: /*@
4351:    TSGetConvergedReason - Gets the reason the TS iteration was stopped.

4353:    Not Collective

4355:    Input Parameter:
4356: .  ts - the TS context

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

4362:    Level: beginner

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

4367: .seealso: `TSSetConvergenceTest()`, `TSConvergedReason`
4368: @*/
4369: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4370: {
4373:   *reason = ts->reason;
4374:   return 0;
4375: }

4377: /*@
4378:    TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve.

4380:    Logically Collective; reason must contain common value

4382:    Input Parameters:
4383: +  ts - the TS context
4384: -  reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the
4385:             manual pages for the individual convergence tests for complete lists

4387:    Level: advanced

4389:    Notes:
4390:    Can only be called while TSSolve() is active.

4392: .seealso: `TSConvergedReason`
4393: @*/
4394: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4395: {
4397:   ts->reason = reason;
4398:   return 0;
4399: }

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

4404:    Not Collective

4406:    Input Parameter:
4407: .  ts - the TS context

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

4412:    Level: beginner

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

4417: .seealso: `TSSetConvergenceTest()`, `TSConvergedReason`
4418: @*/
4419: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4420: {
4423:   *ftime = ts->solvetime;
4424:   return 0;
4425: }

4427: /*@
4428:    TSGetSNESIterations - Gets the total number of nonlinear iterations
4429:    used by the time integrator.

4431:    Not Collective

4433:    Input Parameter:
4434: .  ts - TS context

4436:    Output Parameter:
4437: .  nits - number of nonlinear iterations

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

4442:    Level: intermediate

4444: .seealso: `TSGetKSPIterations()`
4445: @*/
4446: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4447: {
4450:   *nits = ts->snes_its;
4451:   return 0;
4452: }

4454: /*@
4455:    TSGetKSPIterations - Gets the total number of linear iterations
4456:    used by the time integrator.

4458:    Not Collective

4460:    Input Parameter:
4461: .  ts - TS context

4463:    Output Parameter:
4464: .  lits - number of linear iterations

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

4469:    Level: intermediate

4471: .seealso: `TSGetSNESIterations()`, `SNESGetKSPIterations()`
4472: @*/
4473: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4474: {
4477:   *lits = ts->ksp_its;
4478:   return 0;
4479: }

4481: /*@
4482:    TSGetStepRejections - Gets the total number of rejected steps.

4484:    Not Collective

4486:    Input Parameter:
4487: .  ts - TS context

4489:    Output Parameter:
4490: .  rejects - number of steps rejected

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

4495:    Level: intermediate

4497: .seealso: `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4498: @*/
4499: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4500: {
4503:   *rejects = ts->reject;
4504:   return 0;
4505: }

4507: /*@
4508:    TSGetSNESFailures - Gets the total number of failed SNES solves

4510:    Not Collective

4512:    Input Parameter:
4513: .  ts - TS context

4515:    Output Parameter:
4516: .  fails - number of failed nonlinear solves

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

4521:    Level: intermediate

4523: .seealso: `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4524: @*/
4525: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4526: {
4529:   *fails = ts->num_snes_failures;
4530:   return 0;
4531: }

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

4536:    Not Collective

4538:    Input Parameters:
4539: +  ts - TS context
4540: -  rejects - maximum number of rejected steps, pass -1 for unlimited

4542:    Notes:
4543:    The counter is reset to zero for each step

4545:    Options Database Key:
4546: .  -ts_max_reject - Maximum number of step rejections before a step fails

4548:    Level: intermediate

4550: .seealso: `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4551: @*/
4552: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4553: {
4555:   ts->max_reject = rejects;
4556:   return 0;
4557: }

4559: /*@
4560:    TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves

4562:    Not Collective

4564:    Input Parameters:
4565: +  ts - TS context
4566: -  fails - maximum number of failed nonlinear solves, pass -1 for unlimited

4568:    Notes:
4569:    The counter is reset to zero for each successive call to TSSolve().

4571:    Options Database Key:
4572: .  -ts_max_snes_failures - Maximum number of nonlinear solve failures

4574:    Level: intermediate

4576: .seealso: `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4577: @*/
4578: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4579: {
4581:   ts->max_snes_failures = fails;
4582:   return 0;
4583: }

4585: /*@
4586:    TSSetErrorIfStepFails - Error if no step succeeds

4588:    Not Collective

4590:    Input Parameters:
4591: +  ts - TS context
4592: -  err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure

4594:    Options Database Key:
4595: .  -ts_error_if_step_fails - Error if no step succeeds

4597:    Level: intermediate

4599: .seealso: `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4600: @*/
4601: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
4602: {
4604:   ts->errorifstepfailed = err;
4605:   return 0;
4606: }

4608: /*@
4609:    TSGetAdapt - Get the adaptive controller context for the current method

4611:    Collective on TS if controller has not been created yet

4613:    Input Parameter:
4614: .  ts - time stepping context

4616:    Output Parameter:
4617: .  adapt - adaptive controller

4619:    Level: intermediate

4621: .seealso: `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
4622: @*/
4623: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
4624: {
4627:   if (!ts->adapt) {
4628:     TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt);
4629:     PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1);
4630:   }
4631:   *adapt = ts->adapt;
4632:   return 0;
4633: }

4635: /*@
4636:    TSSetTolerances - Set tolerances for local truncation error when using adaptive controller

4638:    Logically Collective

4640:    Input Parameters:
4641: +  ts - time integration context
4642: .  atol - scalar absolute tolerances, PETSC_DECIDE to leave current value
4643: .  vatol - vector of absolute tolerances or NULL, used in preference to atol if present
4644: .  rtol - scalar relative tolerances, PETSC_DECIDE to leave current value
4645: -  vrtol - vector of relative tolerances or NULL, used in preference to atol if present

4647:    Options Database keys:
4648: +  -ts_rtol <rtol> - relative tolerance for local truncation error
4649: -  -ts_atol <atol> - Absolute tolerance for local truncation error

4651:    Notes:
4652:    With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
4653:    (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
4654:    computed only for the differential or the algebraic part then this can be done using the vector of
4655:    tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
4656:    differential part and infinity for the algebraic part, the LTE calculation will include only the
4657:    differential variables.

4659:    Level: beginner

4661: .seealso: `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
4662: @*/
4663: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
4664: {
4665:   if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol;
4666:   if (vatol) {
4667:     PetscObjectReference((PetscObject)vatol);
4668:     VecDestroy(&ts->vatol);
4669:     ts->vatol = vatol;
4670:   }
4671:   if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol;
4672:   if (vrtol) {
4673:     PetscObjectReference((PetscObject)vrtol);
4674:     VecDestroy(&ts->vrtol);
4675:     ts->vrtol = vrtol;
4676:   }
4677:   return 0;
4678: }

4680: /*@
4681:    TSGetTolerances - Get tolerances for local truncation error when using adaptive controller

4683:    Logically Collective

4685:    Input Parameter:
4686: .  ts - time integration context

4688:    Output Parameters:
4689: +  atol - scalar absolute tolerances, NULL to ignore
4690: .  vatol - vector of absolute tolerances, NULL to ignore
4691: .  rtol - scalar relative tolerances, NULL to ignore
4692: -  vrtol - vector of relative tolerances, NULL to ignore

4694:    Level: beginner

4696: .seealso: `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
4697: @*/
4698: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
4699: {
4700:   if (atol) *atol = ts->atol;
4701:   if (vatol) *vatol = ts->vatol;
4702:   if (rtol) *rtol = ts->rtol;
4703:   if (vrtol) *vrtol = ts->vrtol;
4704:   return 0;
4705: }

4707: /*@
4708:    TSErrorWeightedNorm2 - compute a weighted 2-norm of the difference between two state vectors

4710:    Collective on TS

4712:    Input Parameters:
4713: +  ts - time stepping context
4714: .  U - state vector, usually ts->vec_sol
4715: -  Y - state vector to be compared to U

4717:    Output Parameters:
4718: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
4719: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
4720: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

4722:    Level: developer

4724: .seealso: `TSErrorWeightedNorm()`, `TSErrorWeightedNormInfinity()`
4725: @*/
4726: PetscErrorCode TSErrorWeightedNorm2(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4727: {
4728:   PetscInt           i, n, N, rstart;
4729:   PetscInt           n_loc, na_loc, nr_loc;
4730:   PetscReal          n_glb, na_glb, nr_glb;
4731:   const PetscScalar *u, *y;
4732:   PetscReal          sum, suma, sumr, gsum, gsuma, gsumr, diff;
4733:   PetscReal          tol, tola, tolr;
4734:   PetscReal          err_loc[6], err_glb[6];


4747:   VecGetSize(U, &N);
4748:   VecGetLocalSize(U, &n);
4749:   VecGetOwnershipRange(U, &rstart, NULL);
4750:   VecGetArrayRead(U, &u);
4751:   VecGetArrayRead(Y, &y);
4752:   sum    = 0.;
4753:   n_loc  = 0;
4754:   suma   = 0.;
4755:   na_loc = 0;
4756:   sumr   = 0.;
4757:   nr_loc = 0;
4758:   if (ts->vatol && ts->vrtol) {
4759:     const PetscScalar *atol, *rtol;
4760:     VecGetArrayRead(ts->vatol, &atol);
4761:     VecGetArrayRead(ts->vrtol, &rtol);
4762:     for (i = 0; i < n; i++) {
4763:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4764:       diff = PetscAbsScalar(y[i] - u[i]);
4765:       tola = PetscRealPart(atol[i]);
4766:       if (tola > 0.) {
4767:         suma += PetscSqr(diff / tola);
4768:         na_loc++;
4769:       }
4770:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4771:       if (tolr > 0.) {
4772:         sumr += PetscSqr(diff / tolr);
4773:         nr_loc++;
4774:       }
4775:       tol = tola + tolr;
4776:       if (tol > 0.) {
4777:         sum += PetscSqr(diff / tol);
4778:         n_loc++;
4779:       }
4780:     }
4781:     VecRestoreArrayRead(ts->vatol, &atol);
4782:     VecRestoreArrayRead(ts->vrtol, &rtol);
4783:   } else if (ts->vatol) { /* vector atol, scalar rtol */
4784:     const PetscScalar *atol;
4785:     VecGetArrayRead(ts->vatol, &atol);
4786:     for (i = 0; i < n; i++) {
4787:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4788:       diff = PetscAbsScalar(y[i] - u[i]);
4789:       tola = PetscRealPart(atol[i]);
4790:       if (tola > 0.) {
4791:         suma += PetscSqr(diff / tola);
4792:         na_loc++;
4793:       }
4794:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4795:       if (tolr > 0.) {
4796:         sumr += PetscSqr(diff / tolr);
4797:         nr_loc++;
4798:       }
4799:       tol = tola + tolr;
4800:       if (tol > 0.) {
4801:         sum += PetscSqr(diff / tol);
4802:         n_loc++;
4803:       }
4804:     }
4805:     VecRestoreArrayRead(ts->vatol, &atol);
4806:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
4807:     const PetscScalar *rtol;
4808:     VecGetArrayRead(ts->vrtol, &rtol);
4809:     for (i = 0; i < n; i++) {
4810:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4811:       diff = PetscAbsScalar(y[i] - u[i]);
4812:       tola = ts->atol;
4813:       if (tola > 0.) {
4814:         suma += PetscSqr(diff / tola);
4815:         na_loc++;
4816:       }
4817:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4818:       if (tolr > 0.) {
4819:         sumr += PetscSqr(diff / tolr);
4820:         nr_loc++;
4821:       }
4822:       tol = tola + tolr;
4823:       if (tol > 0.) {
4824:         sum += PetscSqr(diff / tol);
4825:         n_loc++;
4826:       }
4827:     }
4828:     VecRestoreArrayRead(ts->vrtol, &rtol);
4829:   } else { /* scalar atol, scalar rtol */
4830:     for (i = 0; i < n; i++) {
4831:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4832:       diff = PetscAbsScalar(y[i] - u[i]);
4833:       tola = ts->atol;
4834:       if (tola > 0.) {
4835:         suma += PetscSqr(diff / tola);
4836:         na_loc++;
4837:       }
4838:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4839:       if (tolr > 0.) {
4840:         sumr += PetscSqr(diff / tolr);
4841:         nr_loc++;
4842:       }
4843:       tol = tola + tolr;
4844:       if (tol > 0.) {
4845:         sum += PetscSqr(diff / tol);
4846:         n_loc++;
4847:       }
4848:     }
4849:   }
4850:   VecRestoreArrayRead(U, &u);
4851:   VecRestoreArrayRead(Y, &y);

4853:   err_loc[0] = sum;
4854:   err_loc[1] = suma;
4855:   err_loc[2] = sumr;
4856:   err_loc[3] = (PetscReal)n_loc;
4857:   err_loc[4] = (PetscReal)na_loc;
4858:   err_loc[5] = (PetscReal)nr_loc;

4860:   MPIU_Allreduce(err_loc, err_glb, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts));

4862:   gsum   = err_glb[0];
4863:   gsuma  = err_glb[1];
4864:   gsumr  = err_glb[2];
4865:   n_glb  = err_glb[3];
4866:   na_glb = err_glb[4];
4867:   nr_glb = err_glb[5];

4869:   *norm = 0.;
4870:   if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
4871:   *norma = 0.;
4872:   if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
4873:   *normr = 0.;
4874:   if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);

4879:   return 0;
4880: }

4882: /*@
4883:    TSErrorWeightedNormInfinity - compute a weighted infinity-norm of the difference between two state vectors

4885:    Collective on TS

4887:    Input Parameters:
4888: +  ts - time stepping context
4889: .  U - state vector, usually ts->vec_sol
4890: -  Y - state vector to be compared to U

4892:    Output Parameters:
4893: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
4894: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
4895: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

4897:    Level: developer

4899: .seealso: `TSErrorWeightedNorm()`, `TSErrorWeightedNorm2()`
4900: @*/
4901: PetscErrorCode TSErrorWeightedNormInfinity(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4902: {
4903:   PetscInt           i, n, N, rstart;
4904:   const PetscScalar *u, *y;
4905:   PetscReal          max, gmax, maxa, gmaxa, maxr, gmaxr;
4906:   PetscReal          tol, tola, tolr, diff;
4907:   PetscReal          err_loc[3], err_glb[3];


4920:   VecGetSize(U, &N);
4921:   VecGetLocalSize(U, &n);
4922:   VecGetOwnershipRange(U, &rstart, NULL);
4923:   VecGetArrayRead(U, &u);
4924:   VecGetArrayRead(Y, &y);

4926:   max  = 0.;
4927:   maxa = 0.;
4928:   maxr = 0.;

4930:   if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
4931:     const PetscScalar *atol, *rtol;
4932:     VecGetArrayRead(ts->vatol, &atol);
4933:     VecGetArrayRead(ts->vrtol, &rtol);

4935:     for (i = 0; i < n; i++) {
4936:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4937:       diff = PetscAbsScalar(y[i] - u[i]);
4938:       tola = PetscRealPart(atol[i]);
4939:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4940:       tol  = tola + tolr;
4941:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4942:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4943:       if (tol > 0.) max = PetscMax(max, diff / tol);
4944:     }
4945:     VecRestoreArrayRead(ts->vatol, &atol);
4946:     VecRestoreArrayRead(ts->vrtol, &rtol);
4947:   } else if (ts->vatol) { /* vector atol, scalar rtol */
4948:     const PetscScalar *atol;
4949:     VecGetArrayRead(ts->vatol, &atol);
4950:     for (i = 0; i < n; i++) {
4951:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4952:       diff = PetscAbsScalar(y[i] - u[i]);
4953:       tola = PetscRealPart(atol[i]);
4954:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4955:       tol  = tola + tolr;
4956:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4957:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4958:       if (tol > 0.) max = PetscMax(max, diff / tol);
4959:     }
4960:     VecRestoreArrayRead(ts->vatol, &atol);
4961:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
4962:     const PetscScalar *rtol;
4963:     VecGetArrayRead(ts->vrtol, &rtol);

4965:     for (i = 0; i < n; i++) {
4966:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4967:       diff = PetscAbsScalar(y[i] - u[i]);
4968:       tola = ts->atol;
4969:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4970:       tol  = tola + tolr;
4971:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4972:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4973:       if (tol > 0.) max = PetscMax(max, diff / tol);
4974:     }
4975:     VecRestoreArrayRead(ts->vrtol, &rtol);
4976:   } else { /* scalar atol, scalar rtol */

4978:     for (i = 0; i < n; i++) {
4979:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4980:       diff = PetscAbsScalar(y[i] - u[i]);
4981:       tola = ts->atol;
4982:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4983:       tol  = tola + tolr;
4984:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4985:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4986:       if (tol > 0.) max = PetscMax(max, diff / tol);
4987:     }
4988:   }
4989:   VecRestoreArrayRead(U, &u);
4990:   VecRestoreArrayRead(Y, &y);
4991:   err_loc[0] = max;
4992:   err_loc[1] = maxa;
4993:   err_loc[2] = maxr;
4994:   MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
4995:   gmax  = err_glb[0];
4996:   gmaxa = err_glb[1];
4997:   gmaxr = err_glb[2];

4999:   *norm  = gmax;
5000:   *norma = gmaxa;
5001:   *normr = gmaxr;
5005:   return 0;
5006: }

5008: /*@
5009:    TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances

5011:    Collective on TS

5013:    Input Parameters:
5014: +  ts - time stepping context
5015: .  U - state vector, usually ts->vec_sol
5016: .  Y - state vector to be compared to U
5017: -  wnormtype - norm type, either NORM_2 or NORM_INFINITY

5019:    Output Parameters:
5020: +  norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5021: .  norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5022: -  normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

5024:    Options Database Keys:
5025: .  -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5027:    Level: developer

5029: .seealso: `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`, `TSErrorWeightedENorm`
5030: @*/
5031: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5032: {
5033:   if (wnormtype == NORM_2) TSErrorWeightedNorm2(ts, U, Y, norm, norma, normr);
5034:   else if (wnormtype == NORM_INFINITY) TSErrorWeightedNormInfinity(ts, U, Y, norm, norma, normr);
5035:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5036:   return 0;
5037: }

5039: /*@
5040:    TSErrorWeightedENorm2 - compute a weighted 2 error norm based on supplied absolute and relative tolerances

5042:    Collective on TS

5044:    Input Parameters:
5045: +  ts - time stepping context
5046: .  E - error vector
5047: .  U - state vector, usually ts->vec_sol
5048: -  Y - state vector, previous time step

5050:    Output Parameters:
5051: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5052: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5053: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5055:    Level: developer

5057: .seealso: `TSErrorWeightedENorm()`, `TSErrorWeightedENormInfinity()`
5058: @*/
5059: PetscErrorCode TSErrorWeightedENorm2(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5060: {
5061:   PetscInt           i, n, N, rstart;
5062:   PetscInt           n_loc, na_loc, nr_loc;
5063:   PetscReal          n_glb, na_glb, nr_glb;
5064:   const PetscScalar *e, *u, *y;
5065:   PetscReal          err, sum, suma, sumr, gsum, gsuma, gsumr;
5066:   PetscReal          tol, tola, tolr;
5067:   PetscReal          err_loc[6], err_glb[6];


5082:   VecGetSize(E, &N);
5083:   VecGetLocalSize(E, &n);
5084:   VecGetOwnershipRange(E, &rstart, NULL);
5085:   VecGetArrayRead(E, &e);
5086:   VecGetArrayRead(U, &u);
5087:   VecGetArrayRead(Y, &y);
5088:   sum    = 0.;
5089:   n_loc  = 0;
5090:   suma   = 0.;
5091:   na_loc = 0;
5092:   sumr   = 0.;
5093:   nr_loc = 0;
5094:   if (ts->vatol && ts->vrtol) {
5095:     const PetscScalar *atol, *rtol;
5096:     VecGetArrayRead(ts->vatol, &atol);
5097:     VecGetArrayRead(ts->vrtol, &rtol);
5098:     for (i = 0; i < n; i++) {
5099:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5100:       err  = PetscAbsScalar(e[i]);
5101:       tola = PetscRealPart(atol[i]);
5102:       if (tola > 0.) {
5103:         suma += PetscSqr(err / tola);
5104:         na_loc++;
5105:       }
5106:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5107:       if (tolr > 0.) {
5108:         sumr += PetscSqr(err / tolr);
5109:         nr_loc++;
5110:       }
5111:       tol = tola + tolr;
5112:       if (tol > 0.) {
5113:         sum += PetscSqr(err / tol);
5114:         n_loc++;
5115:       }
5116:     }
5117:     VecRestoreArrayRead(ts->vatol, &atol);
5118:     VecRestoreArrayRead(ts->vrtol, &rtol);
5119:   } else if (ts->vatol) { /* vector atol, scalar rtol */
5120:     const PetscScalar *atol;
5121:     VecGetArrayRead(ts->vatol, &atol);
5122:     for (i = 0; i < n; i++) {
5123:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5124:       err  = PetscAbsScalar(e[i]);
5125:       tola = PetscRealPart(atol[i]);
5126:       if (tola > 0.) {
5127:         suma += PetscSqr(err / tola);
5128:         na_loc++;
5129:       }
5130:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5131:       if (tolr > 0.) {
5132:         sumr += PetscSqr(err / tolr);
5133:         nr_loc++;
5134:       }
5135:       tol = tola + tolr;
5136:       if (tol > 0.) {
5137:         sum += PetscSqr(err / tol);
5138:         n_loc++;
5139:       }
5140:     }
5141:     VecRestoreArrayRead(ts->vatol, &atol);
5142:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
5143:     const PetscScalar *rtol;
5144:     VecGetArrayRead(ts->vrtol, &rtol);
5145:     for (i = 0; i < n; i++) {
5146:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5147:       err  = PetscAbsScalar(e[i]);
5148:       tola = ts->atol;
5149:       if (tola > 0.) {
5150:         suma += PetscSqr(err / tola);
5151:         na_loc++;
5152:       }
5153:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5154:       if (tolr > 0.) {
5155:         sumr += PetscSqr(err / tolr);
5156:         nr_loc++;
5157:       }
5158:       tol = tola + tolr;
5159:       if (tol > 0.) {
5160:         sum += PetscSqr(err / tol);
5161:         n_loc++;
5162:       }
5163:     }
5164:     VecRestoreArrayRead(ts->vrtol, &rtol);
5165:   } else { /* scalar atol, scalar rtol */
5166:     for (i = 0; i < n; i++) {
5167:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5168:       err  = PetscAbsScalar(e[i]);
5169:       tola = ts->atol;
5170:       if (tola > 0.) {
5171:         suma += PetscSqr(err / tola);
5172:         na_loc++;
5173:       }
5174:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5175:       if (tolr > 0.) {
5176:         sumr += PetscSqr(err / tolr);
5177:         nr_loc++;
5178:       }
5179:       tol = tola + tolr;
5180:       if (tol > 0.) {
5181:         sum += PetscSqr(err / tol);
5182:         n_loc++;
5183:       }
5184:     }
5185:   }
5186:   VecRestoreArrayRead(E, &e);
5187:   VecRestoreArrayRead(U, &u);
5188:   VecRestoreArrayRead(Y, &y);

5190:   err_loc[0] = sum;
5191:   err_loc[1] = suma;
5192:   err_loc[2] = sumr;
5193:   err_loc[3] = (PetscReal)n_loc;
5194:   err_loc[4] = (PetscReal)na_loc;
5195:   err_loc[5] = (PetscReal)nr_loc;

5197:   MPIU_Allreduce(err_loc, err_glb, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts));

5199:   gsum   = err_glb[0];
5200:   gsuma  = err_glb[1];
5201:   gsumr  = err_glb[2];
5202:   n_glb  = err_glb[3];
5203:   na_glb = err_glb[4];
5204:   nr_glb = err_glb[5];

5206:   *norm = 0.;
5207:   if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
5208:   *norma = 0.;
5209:   if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
5210:   *normr = 0.;
5211:   if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);

5216:   return 0;
5217: }

5219: /*@
5220:    TSErrorWeightedENormInfinity - compute a weighted infinity error norm based on supplied absolute and relative tolerances
5221:    Collective on TS

5223:    Input Parameters:
5224: +  ts - time stepping context
5225: .  E - error vector
5226: .  U - state vector, usually ts->vec_sol
5227: -  Y - state vector, previous time step

5229:    Output Parameters:
5230: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5231: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5232: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5234:    Level: developer

5236: .seealso: `TSErrorWeightedENorm()`, `TSErrorWeightedENorm2()`
5237: @*/
5238: PetscErrorCode TSErrorWeightedENormInfinity(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5239: {
5240:   PetscInt           i, n, N, rstart;
5241:   const PetscScalar *e, *u, *y;
5242:   PetscReal          err, max, gmax, maxa, gmaxa, maxr, gmaxr;
5243:   PetscReal          tol, tola, tolr;
5244:   PetscReal          err_loc[3], err_glb[3];


5259:   VecGetSize(E, &N);
5260:   VecGetLocalSize(E, &n);
5261:   VecGetOwnershipRange(E, &rstart, NULL);
5262:   VecGetArrayRead(E, &e);
5263:   VecGetArrayRead(U, &u);
5264:   VecGetArrayRead(Y, &y);

5266:   max  = 0.;
5267:   maxa = 0.;
5268:   maxr = 0.;

5270:   if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
5271:     const PetscScalar *atol, *rtol;
5272:     VecGetArrayRead(ts->vatol, &atol);
5273:     VecGetArrayRead(ts->vrtol, &rtol);

5275:     for (i = 0; i < n; i++) {
5276:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5277:       err  = PetscAbsScalar(e[i]);
5278:       tola = PetscRealPart(atol[i]);
5279:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5280:       tol  = tola + tolr;
5281:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5282:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5283:       if (tol > 0.) max = PetscMax(max, err / tol);
5284:     }
5285:     VecRestoreArrayRead(ts->vatol, &atol);
5286:     VecRestoreArrayRead(ts->vrtol, &rtol);
5287:   } else if (ts->vatol) { /* vector atol, scalar rtol */
5288:     const PetscScalar *atol;
5289:     VecGetArrayRead(ts->vatol, &atol);
5290:     for (i = 0; i < n; i++) {
5291:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5292:       err  = PetscAbsScalar(e[i]);
5293:       tola = PetscRealPart(atol[i]);
5294:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5295:       tol  = tola + tolr;
5296:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5297:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5298:       if (tol > 0.) max = PetscMax(max, err / tol);
5299:     }
5300:     VecRestoreArrayRead(ts->vatol, &atol);
5301:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
5302:     const PetscScalar *rtol;
5303:     VecGetArrayRead(ts->vrtol, &rtol);

5305:     for (i = 0; i < n; i++) {
5306:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5307:       err  = PetscAbsScalar(e[i]);
5308:       tola = ts->atol;
5309:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5310:       tol  = tola + tolr;
5311:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5312:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5313:       if (tol > 0.) max = PetscMax(max, err / tol);
5314:     }
5315:     VecRestoreArrayRead(ts->vrtol, &rtol);
5316:   } else { /* scalar atol, scalar rtol */

5318:     for (i = 0; i < n; i++) {
5319:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5320:       err  = PetscAbsScalar(e[i]);
5321:       tola = ts->atol;
5322:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5323:       tol  = tola + tolr;
5324:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5325:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5326:       if (tol > 0.) max = PetscMax(max, err / tol);
5327:     }
5328:   }
5329:   VecRestoreArrayRead(E, &e);
5330:   VecRestoreArrayRead(U, &u);
5331:   VecRestoreArrayRead(Y, &y);
5332:   err_loc[0] = max;
5333:   err_loc[1] = maxa;
5334:   err_loc[2] = maxr;
5335:   MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
5336:   gmax  = err_glb[0];
5337:   gmaxa = err_glb[1];
5338:   gmaxr = err_glb[2];

5340:   *norm  = gmax;
5341:   *norma = gmaxa;
5342:   *normr = gmaxr;
5346:   return 0;
5347: }

5349: /*@
5350:    TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances

5352:    Collective on TS

5354:    Input Parameters:
5355: +  ts - time stepping context
5356: .  E - error vector
5357: .  U - state vector, usually ts->vec_sol
5358: .  Y - state vector, previous time step
5359: -  wnormtype - norm type, either NORM_2 or NORM_INFINITY

5361:    Output Parameters:
5362: +  norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5363: .  norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5364: -  normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

5366:    Options Database Keys:
5367: .  -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5369:    Level: developer

5371: .seealso: `TSErrorWeightedENormInfinity()`, `TSErrorWeightedENorm2()`, `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`
5372: @*/
5373: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5374: {
5375:   if (wnormtype == NORM_2) TSErrorWeightedENorm2(ts, E, U, Y, norm, norma, normr);
5376:   else if (wnormtype == NORM_INFINITY) TSErrorWeightedENormInfinity(ts, E, U, Y, norm, norma, normr);
5377:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5378:   return 0;
5379: }

5381: /*@
5382:    TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler

5384:    Logically Collective on TS

5386:    Input Parameters:
5387: +  ts - time stepping context
5388: -  cfltime - maximum stable time step if using forward Euler (value can be different on each process)

5390:    Note:
5391:    After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()

5393:    Level: intermediate

5395: .seealso: `TSGetCFLTime()`, `TSADAPTCFL`
5396: @*/
5397: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5398: {
5400:   ts->cfltime_local = cfltime;
5401:   ts->cfltime       = -1.;
5402:   return 0;
5403: }

5405: /*@
5406:    TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler

5408:    Collective on TS

5410:    Input Parameter:
5411: .  ts - time stepping context

5413:    Output Parameter:
5414: .  cfltime - maximum stable time step for forward Euler

5416:    Level: advanced

5418: .seealso: `TSSetCFLTimeLocal()`
5419: @*/
5420: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5421: {
5422:   if (ts->cfltime < 0) MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts));
5423:   *cfltime = ts->cfltime;
5424:   return 0;
5425: }

5427: /*@
5428:    TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu

5430:    Input Parameters:
5431: +  ts   - the TS context.
5432: .  xl   - lower bound.
5433: -  xu   - upper bound.

5435:    Notes:
5436:    If this routine is not called then the lower and upper bounds are set to
5437:    PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().

5439:    Level: advanced

5441: @*/
5442: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5443: {
5444:   SNES snes;

5446:   TSGetSNES(ts, &snes);
5447:   SNESVISetVariableBounds(snes, xl, xu);
5448:   return 0;
5449: }

5451: /*@
5452:    TSComputeLinearStability - computes the linear stability function at a point

5454:    Collective on TS

5456:    Input Parameters:
5457: +  ts - the TS context
5458: -  xr,xi - real and imaginary part of input arguments

5460:    Output Parameters:
5461: .  yr,yi - real and imaginary part of function value

5463:    Level: developer

5465: .seealso: `TSSetRHSFunction()`, `TSComputeIFunction()`
5466: @*/
5467: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5468: {
5470:   PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5471:   return 0;
5472: }

5474: /*@
5475:    TSRestartStep - Flags the solver to restart the next step

5477:    Collective on TS

5479:    Input Parameter:
5480: .  ts - the TS context obtained from TSCreate()

5482:    Level: advanced

5484:    Notes:
5485:    Multistep methods like BDF or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5486:    discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5487:    vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5488:    the sake of correctness and maximum safety, users are expected to call TSRestart() whenever they introduce
5489:    discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5490:    discontinuous source terms).

5492: .seealso: `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5493: @*/
5494: PetscErrorCode TSRestartStep(TS ts)
5495: {
5497:   ts->steprestart = PETSC_TRUE;
5498:   return 0;
5499: }

5501: /*@
5502:    TSRollBack - Rolls back one time step

5504:    Collective on TS

5506:    Input Parameter:
5507: .  ts - the TS context obtained from TSCreate()

5509:    Level: advanced

5511: .seealso: `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5512: @*/
5513: PetscErrorCode TSRollBack(TS ts)
5514: {
5517:   PetscUseTypeMethod(ts, rollback);
5518:   ts->time_step  = ts->ptime - ts->ptime_prev;
5519:   ts->ptime      = ts->ptime_prev;
5520:   ts->ptime_prev = ts->ptime_prev_rollback;
5521:   ts->steps--;
5522:   ts->steprollback = PETSC_TRUE;
5523:   return 0;
5524: }

5526: /*@
5527:    TSGetStages - Get the number of stages and stage values

5529:    Input Parameter:
5530: .  ts - the TS context obtained from TSCreate()

5532:    Output Parameters:
5533: +  ns - the number of stages
5534: -  Y - the current stage vectors

5536:    Level: advanced

5538:    Notes: Both ns and Y can be NULL.

5540: .seealso: `TSCreate()`
5541: @*/
5542: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5543: {
5547:   if (!ts->ops->getstages) {
5548:     if (ns) *ns = 0;
5549:     if (Y) *Y = NULL;
5550:   } else PetscUseTypeMethod(ts, getstages, ns, Y);
5551:   return 0;
5552: }

5554: /*@C
5555:   TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.

5557:   Collective on SNES

5559:   Input Parameters:
5560: + ts - the TS context
5561: . t - current timestep
5562: . U - state vector
5563: . Udot - time derivative of state vector
5564: . shift - shift to apply, see note below
5565: - ctx - an optional user context

5567:   Output Parameters:
5568: + J - Jacobian matrix (not altered in this routine)
5569: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)

5571:   Level: intermediate

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

5576:   dF/dU + shift*dF/dUdot

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

5581:   This will first try to get the coloring from the DM.  If the DM type has no coloring
5582:   routine, then it will try to get the coloring from the matrix.  This requires that the
5583:   matrix have nonzero entries precomputed.

5585: .seealso: `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5586: @*/
5587: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5588: {
5589:   SNES          snes;
5590:   MatFDColoring color;
5591:   PetscBool     hascolor, matcolor = PETSC_FALSE;

5593:   PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL);
5594:   PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color);
5595:   if (!color) {
5596:     DM         dm;
5597:     ISColoring iscoloring;

5599:     TSGetDM(ts, &dm);
5600:     DMHasColoring(dm, &hascolor);
5601:     if (hascolor && !matcolor) {
5602:       DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring);
5603:       MatFDColoringCreate(B, iscoloring, &color);
5604:       MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts);
5605:       MatFDColoringSetFromOptions(color);
5606:       MatFDColoringSetUp(B, iscoloring, color);
5607:       ISColoringDestroy(&iscoloring);
5608:     } else {
5609:       MatColoring mc;

5611:       MatColoringCreate(B, &mc);
5612:       MatColoringSetDistance(mc, 2);
5613:       MatColoringSetType(mc, MATCOLORINGSL);
5614:       MatColoringSetFromOptions(mc);
5615:       MatColoringApply(mc, &iscoloring);
5616:       MatColoringDestroy(&mc);
5617:       MatFDColoringCreate(B, iscoloring, &color);
5618:       MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts);
5619:       MatFDColoringSetFromOptions(color);
5620:       MatFDColoringSetUp(B, iscoloring, color);
5621:       ISColoringDestroy(&iscoloring);
5622:     }
5623:     PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color);
5624:     PetscObjectDereference((PetscObject)color);
5625:   }
5626:   TSGetSNES(ts, &snes);
5627:   MatFDColoringApply(B, color, U, snes);
5628:   if (J != B) {
5629:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
5630:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
5631:   }
5632:   return 0;
5633: }

5635: /*@
5636:     TSSetFunctionDomainError - Set a function that tests if the current state vector is valid

5638:     Input Parameters:
5639: +    ts - the TS context
5640: -    func - function called within TSFunctionDomainError

5642:     Calling sequence of func:
5643: $     PetscErrorCode func(TS ts,PetscReal time,Vec state,PetscBool reject)

5645: +   ts - the TS context
5646: .   time - the current time (of the stage)
5647: .   state - the state to check if it is valid
5648: -   reject - (output parameter) PETSC_FALSE if the state is acceptable, PETSC_TRUE if not acceptable

5650:     Level: intermediate

5652:     Notes:
5653:       If an implicit ODE solver is being used then, in addition to providing this routine, the
5654:       user's code should call SNESSetFunctionDomainError() when domain errors occur during
5655:       function evaluations where the functions are provided by TSSetIFunction() or TSSetRHSFunction().
5656:       Use TSGetSNES() to obtain the SNES object

5658:     Developer Notes:
5659:       The naming of this function is inconsistent with the SNESSetFunctionDomainError()
5660:       since one takes a function pointer and the other does not.

5662: .seealso: `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5663: @*/

5665: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS, PetscReal, Vec, PetscBool *))
5666: {
5668:   ts->functiondomainerror = func;
5669:   return 0;
5670: }

5672: /*@
5673:     TSFunctionDomainError - Checks if the current state is valid

5675:     Input Parameters:
5676: +    ts - the TS context
5677: .    stagetime - time of the simulation
5678: -    Y - state vector to check.

5680:     Output Parameter:
5681: .    accept - Set to PETSC_FALSE if the current state vector is valid.

5683:     Note:
5684:     This function is called by the TS integration routines and calls the user provided function (set with TSSetFunctionDomainError())
5685:     to check if the current state is valid.

5687:     Level: developer

5689: .seealso: `TSSetFunctionDomainError()`
5690: @*/
5691: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5692: {
5694:   *accept = PETSC_TRUE;
5695:   if (ts->functiondomainerror) (*ts->functiondomainerror)(ts, stagetime, Y, accept);
5696:   return 0;
5697: }

5699: /*@C
5700:   TSClone - This function clones a time step object.

5702:   Collective

5704:   Input Parameter:
5705: . tsin    - The input TS

5707:   Output Parameter:
5708: . tsout   - The output TS (cloned)

5710:   Notes:
5711:   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.

5713:   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);

5715:   Level: developer

5717: .seealso: `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5718: @*/
5719: PetscErrorCode TSClone(TS tsin, TS *tsout)
5720: {
5721:   TS     t;
5722:   SNES   snes_start;
5723:   DM     dm;
5724:   TSType type;

5727:   *tsout = NULL;

5729:   PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView);

5731:   /* General TS description */
5732:   t->numbermonitors    = 0;
5733:   t->monitorFrequency  = 1;
5734:   t->setupcalled       = 0;
5735:   t->ksp_its           = 0;
5736:   t->snes_its          = 0;
5737:   t->nwork             = 0;
5738:   t->rhsjacobian.time  = PETSC_MIN_REAL;
5739:   t->rhsjacobian.scale = 1.;
5740:   t->ijacobian.shift   = 1.;

5742:   TSGetSNES(tsin, &snes_start);
5743:   TSSetSNES(t, snes_start);

5745:   TSGetDM(tsin, &dm);
5746:   TSSetDM(t, dm);

5748:   t->adapt = tsin->adapt;
5749:   PetscObjectReference((PetscObject)t->adapt);

5751:   t->trajectory = tsin->trajectory;
5752:   PetscObjectReference((PetscObject)t->trajectory);

5754:   t->event = tsin->event;
5755:   if (t->event) t->event->refct++;

5757:   t->problem_type      = tsin->problem_type;
5758:   t->ptime             = tsin->ptime;
5759:   t->ptime_prev        = tsin->ptime_prev;
5760:   t->time_step         = tsin->time_step;
5761:   t->max_time          = tsin->max_time;
5762:   t->steps             = tsin->steps;
5763:   t->max_steps         = tsin->max_steps;
5764:   t->equation_type     = tsin->equation_type;
5765:   t->atol              = tsin->atol;
5766:   t->rtol              = tsin->rtol;
5767:   t->max_snes_failures = tsin->max_snes_failures;
5768:   t->max_reject        = tsin->max_reject;
5769:   t->errorifstepfailed = tsin->errorifstepfailed;

5771:   TSGetType(tsin, &type);
5772:   TSSetType(t, type);

5774:   t->vec_sol = NULL;

5776:   t->cfltime          = tsin->cfltime;
5777:   t->cfltime_local    = tsin->cfltime_local;
5778:   t->exact_final_time = tsin->exact_final_time;

5780:   PetscMemcpy(t->ops, tsin->ops, sizeof(struct _TSOps));

5782:   if (((PetscObject)tsin)->fortran_func_pointers) {
5783:     PetscInt i;
5784:     PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers);
5785:     for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5786:   }
5787:   *tsout = t;
5788:   return 0;
5789: }

5791: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5792: {
5793:   TS ts = (TS)ctx;

5795:   TSComputeRHSFunction(ts, 0, x, y);
5796:   return 0;
5797: }

5799: /*@
5800:     TSRHSJacobianTest - Compares the multiply routine provided to the MATSHELL with differencing on the TS given RHS function.

5802:    Logically Collective on TS

5804:     Input Parameters:
5805:     TS - the time stepping routine

5807:    Output Parameter:
5808: .   flg - PETSC_TRUE if the multiply is likely correct

5810:    Options Database:
5811:  .   -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator

5813:    Level: advanced

5815:    Notes:
5816:     This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian

5818: .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5819: @*/
5820: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5821: {
5822:   Mat           J, B;
5823:   TSRHSJacobian func;
5824:   void         *ctx;

5826:   TSGetRHSJacobian(ts, &J, &B, &func, &ctx);
5827:   (*func)(ts, 0.0, ts->vec_sol, J, B, ctx);
5828:   MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg);
5829:   return 0;
5830: }

5832: /*@C
5833:     TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the MATSHELL with differencing on the TS given RHS function.

5835:    Logically Collective on TS

5837:     Input Parameters:
5838:     TS - the time stepping routine

5840:    Output Parameter:
5841: .   flg - PETSC_TRUE if the multiply is likely correct

5843:    Options Database:
5844: .   -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator

5846:    Notes:
5847:     This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian

5849:    Level: advanced

5851: .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5852: @*/
5853: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5854: {
5855:   Mat           J, B;
5856:   void         *ctx;
5857:   TSRHSJacobian func;

5859:   TSGetRHSJacobian(ts, &J, &B, &func, &ctx);
5860:   (*func)(ts, 0.0, ts->vec_sol, J, B, ctx);
5861:   MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg);
5862:   return 0;
5863: }

5865: /*@
5866:   TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.

5868:   Logically collective

5870:   Input Parameters:
5871: +  ts - timestepping context
5872: -  use_splitrhsfunction - PETSC_TRUE indicates that the split RHSFunction will be used

5874:   Options Database:
5875: .   -ts_use_splitrhsfunction - <true,false>

5877:   Notes:
5878:     This is only useful for multirate methods

5880:   Level: intermediate

5882: .seealso: `TSGetUseSplitRHSFunction()`
5883: @*/
5884: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5885: {
5887:   ts->use_splitrhsfunction = use_splitrhsfunction;
5888:   return 0;
5889: }

5891: /*@
5892:   TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.

5894:   Not collective

5896:   Input Parameter:
5897: .  ts - timestepping context

5899:   Output Parameter:
5900: .  use_splitrhsfunction - PETSC_TRUE indicates that the split RHSFunction will be used

5902:   Level: intermediate

5904: .seealso: `TSSetUseSplitRHSFunction()`
5905: @*/
5906: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5907: {
5909:   *use_splitrhsfunction = ts->use_splitrhsfunction;
5910:   return 0;
5911: }

5913: /*@
5914:     TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.

5916:    Logically  Collective on ts

5918:    Input Parameters:
5919: +  ts - the time-stepper
5920: -  str - the structure (the default is UNKNOWN_NONZERO_PATTERN)

5922:    Level: intermediate

5924:    Notes:
5925:      When the relationship between the nonzero structures is known and supplied the solution process can be much faster

5927: .seealso: `MatAXPY()`, `MatStructure`
5928:  @*/
5929: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5930: {
5932:   ts->axpy_pattern = str;
5933:   return 0;
5934: }

5936: /*@
5937:   TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested.

5939:   Collective on ts

5941:   Input Parameters:
5942: + ts - the time-stepper
5943: . n - number of the time points (>=2)
5944: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.

5946:   Options Database Keys:
5947: . -ts_time_span <t0,...tf> - Sets the time span

5949:   Level: beginner

5951:   Notes:
5952:   The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
5953:   TS_EXACTFINALTIME_MATCHSTEP must be used to make the last time step in each sub-interval match the intermediate points specified.
5954:   The intermediate solutions are saved in a vector array that can be accessed with TSGetSolutions(). Thus using time span may
5955:   pressure the memory system when using a large number of span points.

5957: .seealso: `TSGetTimeSpan()`, `TSGetSolutions()`
5958:  @*/
5959: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
5960: {
5963:   if (ts->tspan && n != ts->tspan->num_span_times) {
5964:     PetscFree(ts->tspan->span_times);
5965:     VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol);
5966:     PetscMalloc1(n, &ts->tspan->span_times);
5967:   }
5968:   if (!ts->tspan) {
5969:     TSTimeSpan tspan;
5970:     PetscNew(&tspan);
5971:     PetscMalloc1(n, &tspan->span_times);
5972:     tspan->reltol = 1e-6;
5973:     tspan->abstol = 10 * PETSC_MACHINE_EPSILON;
5974:     ts->tspan     = tspan;
5975:   }
5976:   ts->tspan->num_span_times = n;
5977:   PetscArraycpy(ts->tspan->span_times, span_times, n);
5978:   TSSetTime(ts, ts->tspan->span_times[0]);
5979:   TSSetMaxTime(ts, ts->tspan->span_times[n - 1]);
5980:   return 0;
5981: }

5983: /*@C
5984:   TSGetTimeSpan - gets the time span.

5986:   Not Collective

5988:   Input Parameter:
5989: . ts - the time-stepper

5991:   Output Parameters:
5992: + n - number of the time points (>=2)
5993: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively. The values are valid until the TS object is destroyed.

5995:   Level: beginner
5996:   Notes: Both n and span_times can be NULL.

5998: .seealso: `TSSetTimeSpan()`, `TSGetSolutions()`
5999:  @*/
6000: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal **span_times)
6001: {
6005:   if (!ts->tspan) {
6006:     if (n) *n = 0;
6007:     if (span_times) *span_times = NULL;
6008:   } else {
6009:     if (n) *n = ts->tspan->num_span_times;
6010:     if (span_times) *span_times = ts->tspan->span_times;
6011:   }
6012:   return 0;
6013: }

6015: /*@
6016:    TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.

6018:    Input Parameter:
6019: .  ts - the TS context obtained from TSCreate()

6021:    Output Parameters:
6022: +  nsol - the number of solutions
6023: -  Sols - the solution vectors

6025:    Level: beginner

6027:    Notes:
6028:     Both nsol and Sols can be NULL.
6029:     Some time points in the time span may be skipped by TS so that nsol is less than the number of points specified by TSSetTimeSpan(). For example, manipulating the step size, especially with a reduced precision, may cause TS to step over certain points in the span.

6031: .seealso: `TSSetTimeSpan()`
6032: @*/
6033: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
6034: {
6038:   if (!ts->tspan) {
6039:     if (nsol) *nsol = 0;
6040:     if (Sols) *Sols = NULL;
6041:   } else {
6042:     if (nsol) *nsol = ts->tspan->spanctr;
6043:     if (Sols) *Sols = ts->tspan->vecs_sol;
6044:   }
6045:   return 0;
6046: }