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

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

 33:    Options Database Keys:
 34: +  -ts_type <type> - EULER, BEULER, SUNDIALS, PSEUDO, CN, RK, THETA, ALPHA, GLLE,  SSP, GLEE, BSYMP, IRK
 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:    Level: beginner

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

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

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

 81: .seealso: [](chapter_ts), `TS`, `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

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

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

420:    Level: advanced

422:    Note:
423:    This routine should be called after all `TS` options have been set

425: .seealso: [](chapter_ts), `TSTrajectory`, `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

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

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

446:    Level: intermediate

448:    Notes:
449:    This routine should be called after all `TS` options have been set

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

454: .seealso: [](chapter_ts), `TSTrajectory`, `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

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

471:    Level: intermediate

473: .seealso: [](chapter_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
474: @*/
475: PetscErrorCode TSResetTrajectory(TS ts)
476: {
478:   if (ts->trajectory) {
479:     TSTrajectoryDestroy(&ts->trajectory);
480:     TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory);
481:   }
482:   return 0;
483: }

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

488:    Collective

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

493:    Level: intermediate

495: .seealso: [](chapter_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
496: @*/
497: PetscErrorCode TSRemoveTrajectory(TS ts)
498: {
500:   if (ts->trajectory) TSTrajectoryDestroy(&ts->trajectory);
501:   return 0;
502: }

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

508:    Collective

510:    Input Parameters:
511: +  ts - the `TS` context
512: .  t - current timestep
513: -  U - input vector

515:    Output Parameters:
516: +  A - Jacobian matrix
517: -  B - optional preconditioning matrix

519:    Level: developer

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

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

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

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

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

567: /*@
568:    TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS`

570:    Collective

572:    Input Parameters:
573: +  ts - the `TS` context
574: .  t - current time
575: -  U - state vector

577:    Output Parameter:
578: .  y - right hand side

580:    Level: developer

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

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

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


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

615: /*@
616:    TSComputeSolutionFunction - Evaluates the solution function.

618:    Collective

620:    Input Parameters:
621: +  ts - the `TS` context
622: -  t - current time

624:    Output Parameter:
625: .  U - the solution

627:    Level: developer

629: .seealso: [](chapter_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
630: @*/
631: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
632: {
633:   TSSolutionFunction solutionfunction;
634:   void              *ctx;
635:   DM                 dm;

639:   TSGetDM(ts, &dm);
640:   DMTSGetSolutionFunction(dm, &solutionfunction, &ctx);
641:   if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
642:   return 0;
643: }
644: /*@
645:    TSComputeForcingFunction - Evaluates the forcing function.

647:    Collective

649:    Input Parameters:
650: +  ts - the `TS` context
651: -  t - current time

653:    Output Parameter:
654: .  U - the function value

656:    Level: developer

658: .seealso: [](chapter_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
659: @*/
660: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
661: {
662:   void             *ctx;
663:   DM                dm;
664:   TSForcingFunction forcing;

668:   TSGetDM(ts, &dm);
669:   DMTSGetForcingFunction(dm, &forcing, &ctx);

671:   if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
672:   return 0;
673: }

675: static PetscErrorCode TSGetRHSVec_Private(TS ts, Vec *Frhs)
676: {
677:   Vec F;

679:   *Frhs = NULL;
680:   TSGetIFunction(ts, &F, NULL, NULL);
681:   if (!ts->Frhs) VecDuplicate(F, &ts->Frhs);
682:   *Frhs = ts->Frhs;
683:   return 0;
684: }

686: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
687: {
688:   Mat         A, B;
689:   TSIJacobian ijacobian;

691:   if (Arhs) *Arhs = NULL;
692:   if (Brhs) *Brhs = NULL;
693:   TSGetIJacobian(ts, &A, &B, &ijacobian, NULL);
694:   if (Arhs) {
695:     if (!ts->Arhs) {
696:       if (ijacobian) {
697:         MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs);
698:         TSSetMatStructure(ts, SAME_NONZERO_PATTERN);
699:       } else {
700:         ts->Arhs = A;
701:         PetscObjectReference((PetscObject)A);
702:       }
703:     } else {
704:       PetscBool flg;
705:       SNESGetUseMatrixFree(ts->snes, NULL, &flg);
706:       /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
707:       if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
708:         PetscObjectDereference((PetscObject)ts->Arhs);
709:         ts->Arhs = A;
710:         PetscObjectReference((PetscObject)A);
711:       }
712:     }
713:     *Arhs = ts->Arhs;
714:   }
715:   if (Brhs) {
716:     if (!ts->Brhs) {
717:       if (A != B) {
718:         if (ijacobian) {
719:           MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs);
720:         } else {
721:           ts->Brhs = B;
722:           PetscObjectReference((PetscObject)B);
723:         }
724:       } else {
725:         PetscObjectReference((PetscObject)ts->Arhs);
726:         ts->Brhs = ts->Arhs;
727:       }
728:     }
729:     *Brhs = ts->Brhs;
730:   }
731:   return 0;
732: }

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

737:    Collective

739:    Input Parameters:
740: +  ts - the `TS` context
741: .  t - current time
742: .  U - state vector
743: .  Udot - time derivative of state vector
744: -  imex - flag indicates if the method is `TSIMEX` so that the RHSFunction should be kept separate

746:    Output Parameter:
747: .  Y - right hand side

749:    Level: developer

751:    Note:
752:    Most users should not need to explicitly call this routine, as it
753:    is used internally within the nonlinear solvers.

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

758: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
759: @*/
760: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
761: {
762:   TSIFunction   ifunction;
763:   TSRHSFunction rhsfunction;
764:   void         *ctx;
765:   DM            dm;


772:   TSGetDM(ts, &dm);
773:   DMTSGetIFunction(dm, &ifunction, &ctx);
774:   DMTSGetRHSFunction(dm, &rhsfunction, NULL);


778:   PetscLogEventBegin(TS_FunctionEval, ts, U, Udot, Y);
779:   if (ifunction) {
780:     PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
781:     ts->ifuncs++;
782:   }
783:   if (imex) {
784:     if (!ifunction) VecCopy(Udot, Y);
785:   } else if (rhsfunction) {
786:     if (ifunction) {
787:       Vec Frhs;
788:       TSGetRHSVec_Private(ts, &Frhs);
789:       TSComputeRHSFunction(ts, t, U, Frhs);
790:       VecAXPY(Y, -1, Frhs);
791:     } else {
792:       TSComputeRHSFunction(ts, t, U, Y);
793:       VecAYPX(Y, -1, Udot);
794:     }
795:   }
796:   PetscLogEventEnd(TS_FunctionEval, ts, U, Udot, Y);
797:   return 0;
798: }

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

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

806: */
807: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
808: {

813:   if (ts->rhsjacobian.shift) MatShift(A, -ts->rhsjacobian.shift);
814:   if (ts->rhsjacobian.scale == -1.) MatScale(A, -1);
815:   if (B && B == ts->Brhs && A != B) {
816:     if (ts->rhsjacobian.shift) MatShift(B, -ts->rhsjacobian.shift);
817:     if (ts->rhsjacobian.scale == -1.) MatScale(B, -1);
818:   }
819:   ts->rhsjacobian.shift = 0;
820:   ts->rhsjacobian.scale = 1.;
821:   return 0;
822: }

824: /*@
825:    TSComputeIJacobian - Evaluates the Jacobian of the DAE

827:    Collective

829:    Input
830:       Input Parameters:
831: +  ts - the `TS` context
832: .  t - current timestep
833: .  U - state vector
834: .  Udot - time derivative of state vector
835: .  shift - shift to apply, see note below
836: -  imex - flag indicates if the method is `TSIMEX` so that the RHSJacobian should be kept separate

838:    Output Parameters:
839: +  A - Jacobian matrix
840: -  B - matrix from which the preconditioner is constructed; often the same as A

842:    Level: developer

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

847:    dF/dU + shift*dF/dUdot

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

852: .seealso: [](chapter_ts), `TS`, `TSSetIJacobian()`
853: @*/
854: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
855: {
856:   TSIJacobian   ijacobian;
857:   TSRHSJacobian rhsjacobian;
858:   DM            dm;
859:   void         *ctx;


869:   TSGetDM(ts, &dm);
870:   DMTSGetIJacobian(dm, &ijacobian, &ctx);
871:   DMTSGetRHSJacobian(dm, &rhsjacobian, NULL);


875:   PetscLogEventBegin(TS_JacobianEval, ts, U, A, B);
876:   if (ijacobian) {
877:     PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
878:     ts->ijacs++;
879:   }
880:   if (imex) {
881:     if (!ijacobian) { /* system was written as Udot = G(t,U) */
882:       PetscBool assembled;
883:       if (rhsjacobian) {
884:         Mat Arhs = NULL;
885:         TSGetRHSMats_Private(ts, &Arhs, NULL);
886:         if (A == Arhs) {
888:           ts->rhsjacobian.time = PETSC_MIN_REAL;
889:         }
890:       }
891:       MatZeroEntries(A);
892:       MatAssembled(A, &assembled);
893:       if (!assembled) {
894:         MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
895:         MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
896:       }
897:       MatShift(A, shift);
898:       if (A != B) {
899:         MatZeroEntries(B);
900:         MatAssembled(B, &assembled);
901:         if (!assembled) {
902:           MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
903:           MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
904:         }
905:         MatShift(B, shift);
906:       }
907:     }
908:   } else {
909:     Mat Arhs = NULL, Brhs = NULL;

911:     /* RHSJacobian needs to be converted to part of IJacobian if exists */
912:     if (rhsjacobian) TSGetRHSMats_Private(ts, &Arhs, &Brhs);
913:     if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
914:       PetscObjectState Ustate;
915:       PetscObjectId    Uid;
916:       TSRHSFunction    rhsfunction;

918:       DMTSGetRHSFunction(dm, &rhsfunction, NULL);
919:       PetscObjectStateGet((PetscObject)U, &Ustate);
920:       PetscObjectGetId((PetscObject)U, &Uid);
921:       if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
922:           ts->rhsjacobian.scale == -1.) {                      /* No need to recompute RHSJacobian */
923:         MatShift(A, shift - ts->rhsjacobian.shift); /* revert the old shift and add the new shift with a single call to MatShift */
924:         if (A != B) MatShift(B, shift - ts->rhsjacobian.shift);
925:       } else {
926:         PetscBool flg;

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

966: /*@C
967:     TSSetRHSFunction - Sets the routine for evaluating the function,
968:     where U_t = G(t,u).

970:     Logically Collective

972:     Input Parameters:
973: +   ts - the `TS` context obtained from `TSCreate()`
974: .   r - vector to put the computed right hand side (or NULL to have it created)
975: .   f - routine for evaluating the right-hand-side function
976: -   ctx - [optional] user-defined context for private data for the
977:           function evaluation routine (may be NULL)

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

982: +   ts - timestep context
983: .   t - current timestep
984: .   u - input vector
985: .   F - function vector
986: -   ctx - [optional] user-defined function context

988:     Level: beginner

990:     Note:
991:     You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.

993: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
994: @*/
995: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, PetscErrorCode (*f)(TS, PetscReal, Vec, Vec, void *), void *ctx)
996: {
997:   SNES snes;
998:   Vec  ralloc = NULL;
999:   DM   dm;


1004:   TSGetDM(ts, &dm);
1005:   DMTSSetRHSFunction(dm, f, ctx);
1006:   TSGetSNES(ts, &snes);
1007:   if (!r && !ts->dm && ts->vec_sol) {
1008:     VecDuplicate(ts->vec_sol, &ralloc);
1009:     r = ralloc;
1010:   }
1011:   SNESSetFunction(snes, r, SNESTSFormFunction, ts);
1012:   VecDestroy(&ralloc);
1013:   return 0;
1014: }

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

1019:     Logically Collective

1021:     Input Parameters:
1022: +   ts - the `TS` context obtained from `TSCreate()`
1023: .   f - routine for evaluating the solution
1024: -   ctx - [optional] user-defined context for private data for the
1025:           function evaluation routine (may be NULL)

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

1030: +   t - current timestep
1031: .   u - output vector
1032: -   ctx - [optional] user-defined function context

1034:     Options Database Keys:
1035: +  -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1036: -  -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`

1038:     Level: intermediate

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

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

1047: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1048: @*/
1049: PetscErrorCode TSSetSolutionFunction(TS ts, PetscErrorCode (*f)(TS, PetscReal, Vec, void *), void *ctx)
1050: {
1051:   DM dm;

1054:   TSGetDM(ts, &dm);
1055:   DMTSSetSolutionFunction(dm, f, ctx);
1056:   return 0;
1057: }

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

1062:     Logically Collective

1064:     Input Parameters:
1065: +   ts - the `TS` context obtained from `TSCreate()`
1066: .   func - routine for evaluating the forcing function
1067: -   ctx - [optional] user-defined context for private data for the
1068:           function evaluation routine (may be NULL)

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

1073: +   t - current timestep
1074: .   f - output vector
1075: -   ctx - [optional] user-defined function context

1077:     Level: intermediate

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

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

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

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

1091: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1092: @*/
1093: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFunction func, void *ctx)
1094: {
1095:   DM dm;

1098:   TSGetDM(ts, &dm);
1099:   DMTSSetForcingFunction(dm, func, ctx);
1100:   return 0;
1101: }

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

1107:    Logically Collective

1109:    Input Parameters:
1110: +  ts  - the `TS` context obtained from `TSCreate()`
1111: .  Amat - (approximate) Jacobian matrix
1112: .  Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1113: .  f   - the Jacobian evaluation routine
1114: -  ctx - [optional] user-defined context for private data for the
1115:          Jacobian evaluation routine (may be NULL)

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

1120: +  t - current timestep
1121: .  u - input vector
1122: .  Amat - (approximate) Jacobian matrix
1123: .  Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1124: -  ctx - [optional] user-defined context for matrix evaluation routine

1126:    Level: beginner

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

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

1134: .seealso: [](chapter_ts), `TS`, `SNESComputeJacobianDefaultColor()`, `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`
1135: @*/
1136: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobian f, void *ctx)
1137: {
1138:   SNES        snes;
1139:   DM          dm;
1140:   TSIJacobian ijacobian;


1148:   TSGetDM(ts, &dm);
1149:   DMTSSetRHSJacobian(dm, f, ctx);
1150:   DMTSGetIJacobian(dm, &ijacobian, NULL);
1151:   TSGetSNES(ts, &snes);
1152:   if (!ijacobian) SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts);
1153:   if (Amat) {
1154:     PetscObjectReference((PetscObject)Amat);
1155:     MatDestroy(&ts->Arhs);
1156:     ts->Arhs = Amat;
1157:   }
1158:   if (Pmat) {
1159:     PetscObjectReference((PetscObject)Pmat);
1160:     MatDestroy(&ts->Brhs);
1161:     ts->Brhs = Pmat;
1162:   }
1163:   return 0;
1164: }

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

1169:    Logically Collective

1171:    Input Parameters:
1172: +  ts  - the `TS` context obtained from `TSCreate()`
1173: .  r   - vector to hold the residual (or NULL to have it created internally)
1174: .  f   - the function evaluation routine
1175: -  ctx - user-defined context for private data for the function evaluation routine (may be NULL)

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

1180: +  t   - time at step/stage being solved
1181: .  u   - state vector
1182: .  u_t - time derivative of state vector
1183: .  F   - function vector
1184: -  ctx - [optional] user-defined context for matrix evaluation routine

1186:    Level: beginner

1188:    Note:
1189:    The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE.  When solving DAEs you must use this function.

1191: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`, `TSSetIJacobian()`
1192: @*/
1193: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunction f, void *ctx)
1194: {
1195:   SNES snes;
1196:   Vec  ralloc = NULL;
1197:   DM   dm;


1202:   TSGetDM(ts, &dm);
1203:   DMTSSetIFunction(dm, f, ctx);

1205:   TSGetSNES(ts, &snes);
1206:   if (!r && !ts->dm && ts->vec_sol) {
1207:     VecDuplicate(ts->vec_sol, &ralloc);
1208:     r = ralloc;
1209:   }
1210:   SNESSetFunction(snes, r, SNESTSFormFunction, ts);
1211:   VecDestroy(&ralloc);
1212:   return 0;
1213: }

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

1218:    Not Collective

1220:    Input Parameter:
1221: .  ts - the `TS` context

1223:    Output Parameters:
1224: +  r - vector to hold residual (or NULL)
1225: .  func - the function to compute residual (or NULL)
1226: -  ctx - the function context (or NULL)

1228:    Level: advanced

1230: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1231: @*/
1232: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunction *func, void **ctx)
1233: {
1234:   SNES snes;
1235:   DM   dm;

1238:   TSGetSNES(ts, &snes);
1239:   SNESGetFunction(snes, r, NULL, NULL);
1240:   TSGetDM(ts, &dm);
1241:   DMTSGetIFunction(dm, func, ctx);
1242:   return 0;
1243: }

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

1248:    Not Collective

1250:    Input Parameter:
1251: .  ts - the `TS` context

1253:    Output Parameters:
1254: +  r - vector to hold computed right hand side (or NULL)
1255: .  func - the function to compute right hand side (or NULL)
1256: -  ctx - the function context (or NULL)

1258:    Level: advanced

1260: .seealso: [](chapter_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1261: @*/
1262: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunction *func, void **ctx)
1263: {
1264:   SNES snes;
1265:   DM   dm;

1268:   TSGetSNES(ts, &snes);
1269:   SNESGetFunction(snes, r, NULL, NULL);
1270:   TSGetDM(ts, &dm);
1271:   DMTSGetRHSFunction(dm, func, ctx);
1272:   return 0;
1273: }

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

1279:    Logically Collective

1281:    Input Parameters:
1282: +  ts  - the `TS` context obtained from `TSCreate()`
1283: .  Amat - (approximate) Jacobian matrix
1284: .  Pmat - matrix used to compute preconditioner (usually the same as Amat)
1285: .  f   - the Jacobian evaluation routine
1286: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)

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

1291: +  t    - time at step/stage being solved
1292: .  U    - state vector
1293: .  U_t  - time derivative of state vector
1294: .  a    - shift
1295: .  Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1296: .  Pmat - matrix used for constructing preconditioner, usually the same as Amat
1297: -  ctx  - [optional] user-defined context for matrix evaluation routine

1299:    Level: beginner

1301:    Notes:
1302:    The matrices Amat and Pmat are exactly the matrices that are used by `SNES` for the nonlinear solve.

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

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

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

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

1319: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `TSSetRHSJacobian()`, `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1320: @*/
1321: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobian f, void *ctx)
1322: {
1323:   SNES snes;
1324:   DM   dm;


1332:   TSGetDM(ts, &dm);
1333:   DMTSSetIJacobian(dm, f, ctx);

1335:   TSGetSNES(ts, &snes);
1336:   SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts);
1337:   return 0;
1338: }

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

1346:    Logically Collective

1348:    Input Parameters:
1349: +  ts - `TS` context obtained from `TSCreate()`
1350: -  reuse - `PETSC_TRUE` if the RHS Jacobian

1352:    Level: intermediate

1354: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1355: @*/
1356: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1357: {
1358:   ts->rhsjacobian.reuse = reuse;
1359:   return 0;
1360: }

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

1365:    Logically Collective

1367:    Input Parameters:
1368: +  ts  - the `TS` context obtained from `TSCreate()`
1369: .  F   - vector to hold the residual (or NULL to have it created internally)
1370: .  fun - the function evaluation routine
1371: -  ctx - user-defined context for private data for the function evaluation routine (may be NULL)

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

1376: +  t    - time at step/stage being solved
1377: .  U    - state vector
1378: .  U_t  - time derivative of state vector
1379: .  U_tt - second time derivative of state vector
1380: .  F    - function vector
1381: -  ctx  - [optional] user-defined context for matrix evaluation routine (may be NULL)

1383:    Level: beginner

1385: .seealso: [](chapter_ts), `TS`, `TSSetI2Jacobian()`, `TSSetIFunction()`, `TSCreate()`, `TSSetRHSFunction()`
1386: @*/
1387: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2Function fun, void *ctx)
1388: {
1389:   DM dm;

1393:   TSSetIFunction(ts, F, NULL, NULL);
1394:   TSGetDM(ts, &dm);
1395:   DMTSSetI2Function(dm, fun, ctx);
1396:   return 0;
1397: }

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

1402:   Not Collective

1404:   Input Parameter:
1405: . ts - the `TS` context

1407:   Output Parameters:
1408: + r - vector to hold residual (or NULL)
1409: . fun - the function to compute residual (or NULL)
1410: - ctx - the function context (or NULL)

1412:   Level: advanced

1414: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1415: @*/
1416: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2Function *fun, void **ctx)
1417: {
1418:   SNES snes;
1419:   DM   dm;

1422:   TSGetSNES(ts, &snes);
1423:   SNESGetFunction(snes, r, NULL, NULL);
1424:   TSGetDM(ts, &dm);
1425:   DMTSGetI2Function(dm, fun, ctx);
1426:   return 0;
1427: }

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

1433:    Logically Collective

1435:    Input Parameters:
1436: +  ts  - the `TS` context obtained from `TSCreate()`
1437: .  J   - Jacobian matrix
1438: .  P   - preconditioning matrix for J (may be same as J)
1439: .  jac - the Jacobian evaluation routine
1440: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)

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

1445: +  t    - time at step/stage being solved
1446: .  U    - state vector
1447: .  U_t  - time derivative of state vector
1448: .  U_tt - second time derivative of state vector
1449: .  v    - shift for U_t
1450: .  a    - shift for U_tt
1451: .  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
1452: .  P    - preconditioning matrix for J, may be same as J
1453: -  ctx  - [optional] user-defined context for matrix evaluation routine

1455:    Level: beginner

1457:    Notes:
1458:    The matrices J and P are exactly the matrices that are used by `SNES` for the nonlinear solve.

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

1465: .seealso: [](chapter_ts), `TS`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1466: @*/
1467: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2Jacobian jac, void *ctx)
1468: {
1469:   DM dm;

1474:   TSSetIJacobian(ts, J, P, NULL, NULL);
1475:   TSGetDM(ts, &dm);
1476:   DMTSSetI2Jacobian(dm, jac, ctx);
1477:   return 0;
1478: }

1480: /*@C
1481:   TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.

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

1485:   Input Parameter:
1486: . ts  - The `TS` context obtained from `TSCreate()`

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

1494:   Level: advanced

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

1499: .seealso: [](chapter_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1500: @*/
1501: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2Jacobian *jac, void **ctx)
1502: {
1503:   SNES snes;
1504:   DM   dm;

1506:   TSGetSNES(ts, &snes);
1507:   SNESSetUpMatrices(snes);
1508:   SNESGetJacobian(snes, J, P, NULL, NULL);
1509:   TSGetDM(ts, &dm);
1510:   DMTSGetI2Jacobian(dm, jac, ctx);
1511:   return 0;
1512: }

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

1517:   Collective

1519:   Input Parameters:
1520: + ts - the `TS` context
1521: . t - current time
1522: . U - state vector
1523: . V - time derivative of state vector (U_t)
1524: - A - second time derivative of state vector (U_tt)

1526:   Output Parameter:
1527: . F - the residual vector

1529:   Level: developer

1531:   Note:
1532:   Most users should not need to explicitly call this routine, as it
1533:   is used internally within the nonlinear solvers.

1535: .seealso: [](chapter_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1536: @*/
1537: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1538: {
1539:   DM            dm;
1540:   TSI2Function  I2Function;
1541:   void         *ctx;
1542:   TSRHSFunction rhsfunction;


1550:   TSGetDM(ts, &dm);
1551:   DMTSGetI2Function(dm, &I2Function, &ctx);
1552:   DMTSGetRHSFunction(dm, &rhsfunction, NULL);

1554:   if (!I2Function) {
1555:     TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE);
1556:     return 0;
1557:   }

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

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

1563:   if (rhsfunction) {
1564:     Vec Frhs;
1565:     TSGetRHSVec_Private(ts, &Frhs);
1566:     TSComputeRHSFunction(ts, t, U, Frhs);
1567:     VecAXPY(F, -1, Frhs);
1568:   }

1570:   PetscLogEventEnd(TS_FunctionEval, ts, U, V, F);
1571:   return 0;
1572: }

1574: /*@
1575:   TSComputeI2Jacobian - Evaluates the Jacobian of the DAE

1577:   Collective

1579:   Input Parameters:
1580: + ts - the `TS` context
1581: . t - current timestep
1582: . U - state vector
1583: . V - time derivative of state vector
1584: . A - second time derivative of state vector
1585: . shiftV - shift to apply, see note below
1586: - shiftA - shift to apply, see note below

1588:   Output Parameters:
1589: + J - Jacobian matrix
1590: - P - optional preconditioning matrix

1592:   Level: developer

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

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

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

1602: .seealso: [](chapter_ts), `TS`, `TSSetI2Jacobian()`
1603: @*/
1604: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1605: {
1606:   DM            dm;
1607:   TSI2Jacobian  I2Jacobian;
1608:   void         *ctx;
1609:   TSRHSJacobian rhsjacobian;


1618:   TSGetDM(ts, &dm);
1619:   DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx);
1620:   DMTSGetRHSJacobian(dm, &rhsjacobian, NULL);

1622:   if (!I2Jacobian) {
1623:     TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE);
1624:     return 0;
1625:   }

1627:   PetscLogEventBegin(TS_JacobianEval, ts, U, J, P);
1628:   PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1629:   if (rhsjacobian) {
1630:     Mat Jrhs, Prhs;
1631:     TSGetRHSMats_Private(ts, &Jrhs, &Prhs);
1632:     TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs);
1633:     MatAXPY(J, -1, Jrhs, ts->axpy_pattern);
1634:     if (P != J) MatAXPY(P, -1, Prhs, ts->axpy_pattern);
1635:   }

1637:   PetscLogEventEnd(TS_JacobianEval, ts, U, J, P);
1638:   return 0;
1639: }

1641: /*@C
1642:    TSSetTransientVariable - sets function to transform from state to transient variables

1644:    Logically Collective

1646:    Input Parameters:
1647: +  ts - time stepping context on which to change the transient variable
1648: .  tvar - a function that transforms to transient variables
1649: -  ctx - a context for tvar

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

1654: +   ts - timestep context
1655: .   p - input vector (primitive form)
1656: .   c - output vector, transient variables (conservative form)
1657: -   ctx - [optional] user-defined function context

1659:    Level: advanced

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

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

1670: .seealso: [](chapter_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1671: @*/
1672: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariable tvar, void *ctx)
1673: {
1674:   DM dm;

1677:   TSGetDM(ts, &dm);
1678:   DMTSSetTransientVariable(dm, tvar, ctx);
1679:   return 0;
1680: }

1682: /*@
1683:    TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables

1685:    Logically Collective

1687:    Input Parameters:
1688: +  ts - TS on which to compute
1689: -  U - state vector to be transformed to transient variables

1691:    Output Parameters:
1692: .  C - transient (conservative) variable

1694:    Level: developer

1696:    Developer Note:
1697:    If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = NULL is allowed.
1698:    This makes it safe to call without a guard.  One can use `TSHasTransientVariable()` to check if transient variables are
1699:    being used.

1701: .seealso: [](chapter_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1702: @*/
1703: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1704: {
1705:   DM   dm;
1706:   DMTS dmts;

1710:   TSGetDM(ts, &dm);
1711:   DMGetDMTS(dm, &dmts);
1712:   if (dmts->ops->transientvar) {
1714:     (*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx);
1715:   }
1716:   return 0;
1717: }

1719: /*@
1720:    TSHasTransientVariable - determine whether transient variables have been set

1722:    Logically Collective

1724:    Input Parameters:
1725: .  ts - TS on which to compute

1727:    Output Parameters:
1728: .  has - `PETSC_TRUE` if transient variables have been set

1730:    Level: developer

1732: .seealso: [](chapter_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1733: @*/
1734: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1735: {
1736:   DM   dm;
1737:   DMTS dmts;

1740:   TSGetDM(ts, &dm);
1741:   DMGetDMTS(dm, &dmts);
1742:   *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1743:   return 0;
1744: }

1746: /*@
1747:    TS2SetSolution - Sets the initial solution and time derivative vectors
1748:    for use by the `TS` routines handling second order equations.

1750:    Logically Collective

1752:    Input Parameters:
1753: +  ts - the `TS` context obtained from `TSCreate()`
1754: .  u - the solution vector
1755: -  v - the time derivative vector

1757:    Level: beginner

1759: .seealso: [](chapter_ts), `TS`
1760: @*/
1761: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1762: {
1766:   TSSetSolution(ts, u);
1767:   PetscObjectReference((PetscObject)v);
1768:   VecDestroy(&ts->vec_dot);
1769:   ts->vec_dot = v;
1770:   return 0;
1771: }

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

1779:    Not Collective, but Vec returned is parallel if `TS` is parallel

1781:    Input Parameter:
1782: .  ts - the `TS` context obtained from `TSCreate()`

1784:    Output Parameters:
1785: +  u - the vector containing the solution
1786: -  v - the vector containing the time derivative

1788:    Level: intermediate

1790: .seealso: [](chapter_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1791: @*/
1792: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1793: {
1797:   if (u) *u = ts->vec_sol;
1798:   if (v) *v = ts->vec_dot;
1799:   return 0;
1800: }

1802: /*@C
1803:   TSLoad - Loads a `TS` that has been stored in binary  with `TSView()`.

1805:   Collective

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

1812:    Level: intermediate

1814:   Note:
1815:  The type is determined by the data in the file, any type set into the `TS` before this call is ignored.

1817: .seealso: [](chapter_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1818: @*/
1819: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1820: {
1821:   PetscBool isbinary;
1822:   PetscInt  classid;
1823:   char      type[256];
1824:   DMTS      sdm;
1825:   DM        dm;

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

1832:   PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT);
1834:   PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR);
1835:   TSSetType(ts, type);
1836:   PetscTryTypeMethod(ts, load, viewer);
1837:   DMCreate(PetscObjectComm((PetscObject)ts), &dm);
1838:   DMLoad(dm, viewer);
1839:   TSSetDM(ts, dm);
1840:   DMCreateGlobalVector(ts->dm, &ts->vec_sol);
1841:   VecLoad(ts->vec_sol, viewer);
1842:   DMGetDMTS(ts->dm, &sdm);
1843:   DMTSLoad(sdm, viewer);
1844:   return 0;
1845: }

1847: #include <petscdraw.h>
1848: #if defined(PETSC_HAVE_SAWS)
1849: #include <petscviewersaws.h>
1850: #endif

1852: /*@C
1853:    TSViewFromOptions - View a `TS` based on values in the options database

1855:    Collective

1857:    Input Parameters:
1858: +  ts - the `TS` context
1859: .  obj - Optional object that provides the prefix for the options database keys
1860: -  name - command line option string to be passed by user

1862:    Level: intermediate

1864: .seealso: [](chapter_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1865: @*/
1866: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1867: {
1869:   PetscObjectViewFromOptions((PetscObject)ts, obj, name);
1870:   return 0;
1871: }

1873: /*@C
1874:     TSView - Prints the `TS` data structure.

1876:     Collective

1878:     Input Parameters:
1879: +   ts - the `TS` context obtained from `TSCreate()`
1880: -   viewer - visualization context

1882:     Options Database Key:
1883: .   -ts_view - calls `TSView()` at end of `TSStep()`

1885:     Level: beginner

1887:     Notes:
1888:     The available visualization contexts include
1889: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1890: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1891:          output where only the first processor opens
1892:          the file.  All other processors send their
1893:          data to the first processor to print.

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

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

1900: .seealso: [](chapter_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1901: @*/
1902: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1903: {
1904:   TSType    type;
1905:   PetscBool iascii, isstring, isundials, isbinary, isdraw;
1906:   DMTS      sdm;
1907: #if defined(PETSC_HAVE_SAWS)
1908:   PetscBool issaws;
1909: #endif

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

1916:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1917:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
1918:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
1919:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1920: #if defined(PETSC_HAVE_SAWS)
1921:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws);
1922: #endif
1923:   if (iascii) {
1924:     PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer);
1925:     if (ts->ops->view) {
1926:       PetscViewerASCIIPushTab(viewer);
1927:       PetscUseTypeMethod(ts, view, viewer);
1928:       PetscViewerASCIIPopTab(viewer);
1929:     }
1930:     if (ts->max_steps < PETSC_MAX_INT) PetscViewerASCIIPrintf(viewer, "  maximum steps=%" PetscInt_FMT "\n", ts->max_steps);
1931:     if (ts->max_time < PETSC_MAX_REAL) PetscViewerASCIIPrintf(viewer, "  maximum time=%g\n", (double)ts->max_time);
1932:     if (ts->ifuncs) PetscViewerASCIIPrintf(viewer, "  total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs);
1933:     if (ts->ijacs) PetscViewerASCIIPrintf(viewer, "  total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs);
1934:     if (ts->rhsfuncs) PetscViewerASCIIPrintf(viewer, "  total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs);
1935:     if (ts->rhsjacs) PetscViewerASCIIPrintf(viewer, "  total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs);
1936:     if (ts->usessnes) {
1937:       PetscBool lin;
1938:       if (ts->problem_type == TS_NONLINEAR) PetscViewerASCIIPrintf(viewer, "  total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its);
1939:       PetscViewerASCIIPrintf(viewer, "  total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its);
1940:       PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, "");
1941:       PetscViewerASCIIPrintf(viewer, "  total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures);
1942:     }
1943:     PetscViewerASCIIPrintf(viewer, "  total number of rejected steps=%" PetscInt_FMT "\n", ts->reject);
1944:     if (ts->vrtol) PetscViewerASCIIPrintf(viewer, "  using vector of relative error tolerances, ");
1945:     else PetscViewerASCIIPrintf(viewer, "  using relative error tolerance of %g, ", (double)ts->rtol);
1946:     if (ts->vatol) PetscViewerASCIIPrintf(viewer, "  using vector of absolute error tolerances\n");
1947:     else PetscViewerASCIIPrintf(viewer, "  using absolute error tolerance of %g\n", (double)ts->atol);
1948:     PetscViewerASCIIPushTab(viewer);
1949:     TSAdaptView(ts->adapt, viewer);
1950:     PetscViewerASCIIPopTab(viewer);
1951:   } else if (isstring) {
1952:     TSGetType(ts, &type);
1953:     PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type);
1954:     PetscTryTypeMethod(ts, view, viewer);
1955:   } else if (isbinary) {
1956:     PetscInt    classid = TS_FILE_CLASSID;
1957:     MPI_Comm    comm;
1958:     PetscMPIInt rank;
1959:     char        type[256];

1961:     PetscObjectGetComm((PetscObject)ts, &comm);
1962:     MPI_Comm_rank(comm, &rank);
1963:     if (rank == 0) {
1964:       PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT);
1965:       PetscStrncpy(type, ((PetscObject)ts)->type_name, 256);
1966:       PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR);
1967:     }
1968:     PetscTryTypeMethod(ts, view, viewer);
1969:     if (ts->adapt) TSAdaptView(ts->adapt, viewer);
1970:     DMView(ts->dm, viewer);
1971:     VecView(ts->vec_sol, viewer);
1972:     DMGetDMTS(ts->dm, &sdm);
1973:     DMTSView(sdm, viewer);
1974:   } else if (isdraw) {
1975:     PetscDraw draw;
1976:     char      str[36];
1977:     PetscReal x, y, bottom, h;

1979:     PetscViewerDrawGetDraw(viewer, 0, &draw);
1980:     PetscDrawGetCurrentPoint(draw, &x, &y);
1981:     PetscStrcpy(str, "TS: ");
1982:     PetscStrcat(str, ((PetscObject)ts)->type_name);
1983:     PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h);
1984:     bottom = y - h;
1985:     PetscDrawPushCurrentPoint(draw, x, bottom);
1986:     PetscTryTypeMethod(ts, view, viewer);
1987:     if (ts->adapt) TSAdaptView(ts->adapt, viewer);
1988:     if (ts->snes) SNESView(ts->snes, viewer);
1989:     PetscDrawPopCurrentPoint(draw);
1990: #if defined(PETSC_HAVE_SAWS)
1991:   } else if (issaws) {
1992:     PetscMPIInt rank;
1993:     const char *name;

1995:     PetscObjectGetName((PetscObject)ts, &name);
1996:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1997:     if (!((PetscObject)ts)->amsmem && rank == 0) {
1998:       char dir[1024];

2000:       PetscObjectViewSAWs((PetscObject)ts, viewer);
2001:       PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name);
2002:       SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT);
2003:       PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name);
2004:       SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE);
2005:     }
2006:     PetscTryTypeMethod(ts, view, viewer);
2007: #endif
2008:   }
2009:   if (ts->snes && ts->usessnes) {
2010:     PetscViewerASCIIPushTab(viewer);
2011:     SNESView(ts->snes, viewer);
2012:     PetscViewerASCIIPopTab(viewer);
2013:   }
2014:   DMGetDMTS(ts->dm, &sdm);
2015:   DMTSView(sdm, viewer);

2017:   PetscViewerASCIIPushTab(viewer);
2018:   PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials);
2019:   PetscViewerASCIIPopTab(viewer);
2020:   return 0;
2021: }

2023: /*@
2024:    TSSetApplicationContext - Sets an optional user-defined context for
2025:    the timesteppers.

2027:    Logically Collective

2029:    Input Parameters:
2030: +  ts - the `TS` context obtained from `TSCreate()`
2031: -  usrP -  user context

2033:    Level: intermediate

2035:    Fortran Note:
2036:     To use this from Fortran you must write a Fortran interface definition for this
2037:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

2039: .seealso: [](chapter_ts), `TS`, `TSGetApplicationContext()`
2040: @*/
2041: PetscErrorCode TSSetApplicationContext(TS ts, void *usrP)
2042: {
2044:   ts->user = usrP;
2045:   return 0;
2046: }

2048: /*@
2049:     TSGetApplicationContext - Gets the user-defined context for the
2050:     timestepper that was set with `TSSetApplicationContext()`

2052:     Not Collective

2054:     Input Parameter:
2055: .   ts - the `TS` context obtained from `TSCreate()`

2057:     Output Parameter:
2058: .   usrP - user context

2060:     Level: intermediate

2062:     Fortran Note:
2063:     To use this from Fortran you must write a Fortran interface definition for this
2064:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

2066: .seealso: [](chapter_ts), `TS`, `TSSetApplicationContext()`
2067: @*/
2068: PetscErrorCode TSGetApplicationContext(TS ts, void *usrP)
2069: {
2071:   *(void **)usrP = ts->user;
2072:   return 0;
2073: }

2075: /*@
2076:    TSGetStepNumber - Gets the number of time steps completed.

2078:    Not Collective

2080:    Input Parameter:
2081: .  ts - the `TS` context obtained from `TSCreate()`

2083:    Output Parameter:
2084: .  steps - number of steps completed so far

2086:    Level: intermediate

2088: .seealso: [](chapter_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2089: @*/
2090: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2091: {
2094:   *steps = ts->steps;
2095:   return 0;
2096: }

2098: /*@
2099:    TSSetStepNumber - Sets the number of steps completed.

2101:    Logically Collective

2103:    Input Parameters:
2104: +  ts - the `TS` context
2105: -  steps - number of steps completed so far

2107:    Level: developer

2109:    Note:
2110:    For most uses of the `TS` solvers the user need not explicitly call
2111:    `TSSetStepNumber()`, as the step counter is appropriately updated in
2112:    `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2113:    reinitialize timestepping by setting the step counter to zero (and time
2114:    to the initial time) to solve a similar problem with different initial
2115:    conditions or parameters. Other possible use case is to continue
2116:    timestepping from a previously interrupted run in such a way that TS
2117:    monitors will be called with a initial nonzero step counter.

2119: .seealso: [](chapter_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2120: @*/
2121: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2122: {
2126:   ts->steps = steps;
2127:   return 0;
2128: }

2130: /*@
2131:    TSSetTimeStep - Allows one to reset the timestep at any time,
2132:    useful for simple pseudo-timestepping codes.

2134:    Logically Collective

2136:    Input Parameters:
2137: +  ts - the `TS` context obtained from `TSCreate()`
2138: -  time_step - the size of the timestep

2140:    Level: intermediate

2142: .seealso: [](chapter_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2143: @*/
2144: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2145: {
2148:   ts->time_step = time_step;
2149:   return 0;
2150: }

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

2157:   Logically Collective

2159:    Input Parameters:
2160: +   ts - the time-step context
2161: -   eftopt - exact final time option
2162: .vb
2163:   TS_EXACTFINALTIME_STEPOVER    - Don't do anything if final time is exceeded
2164:   TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2165:   TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2166: .ve

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

2171:    Level: beginner

2173:    Note:
2174:    If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2175:    then the final time you selected.

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

2187: /*@
2188:    TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`

2190:    Not Collective

2192:    Input Parameter:
2193: .  ts - the `TS` context

2195:    Output Parameter:
2196: .  eftopt - exact final time option

2198:    Level: beginner

2200: .seealso: [](chapter_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2201: @*/
2202: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2203: {
2206:   *eftopt = ts->exact_final_time;
2207:   return 0;
2208: }

2210: /*@
2211:    TSGetTimeStep - Gets the current timestep size.

2213:    Not Collective

2215:    Input Parameter:
2216: .  ts - the `TS` context obtained from `TSCreate()`

2218:    Output Parameter:
2219: .  dt - the current timestep size

2221:    Level: intermediate

2223: .seealso: [](chapter_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2224: @*/
2225: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2226: {
2229:   *dt = ts->time_step;
2230:   return 0;
2231: }

2233: /*@
2234:    TSGetSolution - Returns the solution at the present timestep. It
2235:    is valid to call this routine inside the function that you are evaluating
2236:    in order to move to the new timestep. This vector not changed until
2237:    the solution at the next timestep has been calculated.

2239:    Not Collective, but v returned is parallel if ts is parallel

2241:    Input Parameter:
2242: .  ts - the `TS` context obtained from `TSCreate()`

2244:    Output Parameter:
2245: .  v - the vector containing the solution

2247:    Level: intermediate

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

2253: .seealso: [](chapter_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2254: @*/
2255: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2256: {
2259:   *v = ts->vec_sol;
2260:   return 0;
2261: }

2263: /*@
2264:    TSGetSolutionComponents - Returns any solution components at the present
2265:    timestep, if available for the time integration method being used.
2266:    Solution components are quantities that share the same size and
2267:    structure as the solution vector.

2269:    Not Collective, but v returned is parallel if ts is parallel

2271:    Parameters :
2272: +  ts - the `TS` context obtained from `TSCreate()` (input parameter).
2273: .  n - If v is PETSC_NULL, then the number of solution components is
2274:        returned through n, else the n-th solution component is
2275:        returned in v.
2276: -  v - the vector containing the n-th solution component
2277:        (may be PETSC_NULL to use this function to find out
2278:         the number of solutions components).

2280:    Level: advanced

2282: .seealso: [](chapter_ts), `TS`, `TSGetSolution()`
2283: @*/
2284: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2285: {
2287:   if (!ts->ops->getsolutioncomponents) *n = 0;
2288:   else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2289:   return 0;
2290: }

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

2296:    Not Collective, but v returned is parallel if ts is parallel

2298:    Parameters :
2299: +  ts - the `TS` context obtained from `TSCreate()` (input parameter).
2300: -  v - the vector containing the auxiliary solution

2302:    Level: intermediate

2304: .seealso: [](chapter_ts), `TS`, `TSGetSolution()`
2305: @*/
2306: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2307: {
2309:   if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2310:   else VecZeroEntries(*v);
2311:   return 0;
2312: }

2314: /*@
2315:    TSGetTimeError - Returns the estimated error vector, if the chosen
2316:    `TSType` has an error estimation functionality and `TSSetTimeError()` was called

2318:    Not Collective, but v returned is parallel if ts is parallel

2320:    Parameters :
2321: +  ts - the `TS` context obtained from `TSCreate()` (input parameter).
2322: .  n - current estimate (n=0) or previous one (n=-1)
2323: -  v - the vector containing the error (same size as the solution).

2325:    Level: intermediate

2327:    Note:
2328:    MUST call after `TSSetUp()`

2330: .seealso: [](chapter_ts), `TSGetSolution()`, `TSSetTimeError()`
2331: @*/
2332: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2333: {
2335:   if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2336:   else VecZeroEntries(*v);
2337:   return 0;
2338: }

2340: /*@
2341:    TSSetTimeError - Sets the estimated error vector, if the chosen
2342:    `TSType` has an error estimation functionality. This can be used
2343:    to restart such a time integrator with a given error vector.

2345:    Not Collective, but v returned is parallel if ts is parallel

2347:    Parameters :
2348: +  ts - the `TS` context obtained from `TSCreate()` (input parameter).
2349: -  v - the vector containing the error (same size as the solution).

2351:    Level: intermediate

2353: .seealso: [](chapter_ts), `TS`, `TSSetSolution()`, `TSGetTimeError)`
2354: @*/
2355: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2356: {
2359:   PetscTryTypeMethod(ts, settimeerror, v);
2360:   return 0;
2361: }

2363: /* ----- Routines to initialize and destroy a timestepper ---- */
2364: /*@
2365:   TSSetProblemType - Sets the type of problem to be solved.

2367:   Not collective

2369:   Input Parameters:
2370: + ts   - The `TS`
2371: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2372: .vb
2373:          U_t - A U = 0      (linear)
2374:          U_t - A(t) U = 0   (linear)
2375:          F(t,U,U_t) = 0     (nonlinear)
2376: .ve

2378:    Level: beginner

2380: .seealso: [](chapter_ts), `TSSetUp()`, `TSProblemType`, `TS`
2381: @*/
2382: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2383: {
2385:   ts->problem_type = type;
2386:   if (type == TS_LINEAR) {
2387:     SNES snes;
2388:     TSGetSNES(ts, &snes);
2389:     SNESSetType(snes, SNESKSPONLY);
2390:   }
2391:   return 0;
2392: }

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

2397:   Not collective

2399:   Input Parameter:
2400: . ts   - The `TS`

2402:   Output Parameter:
2403: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2404: .vb
2405:          M U_t = A U
2406:          M(t) U_t = A(t) U
2407:          F(t,U,U_t)
2408: .ve

2410:    Level: beginner

2412: .seealso: [](chapter_ts), `TSSetUp()`, `TSProblemType`, `TS`
2413: @*/
2414: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2415: {
2418:   *type = ts->problem_type;
2419:   return 0;
2420: }

2422: /*
2423:     Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2424: */
2425: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2426: {
2427:   PetscBool isnone;

2429:   TSGetAdapt(ts, &ts->adapt);
2430:   TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type);

2432:   PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone);
2433:   if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2434:   else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2435:   return 0;
2436: }

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

2441:    Collective

2443:    Input Parameter:
2444: .  ts - the `TS` context obtained from `TSCreate()`

2446:    Level: advanced

2448:    Note:
2449:    For basic use of the `TS` solvers the user need not explicitly call
2450:    `TSSetUp()`, since these actions will automatically occur during
2451:    the call to `TSStep()` or `TSSolve()`.  However, if one wishes to control this
2452:    phase separately, `TSSetUp()` should be called after `TSCreate()`
2453:    and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.

2455: .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2456: @*/
2457: PetscErrorCode TSSetUp(TS ts)
2458: {
2459:   DM dm;
2460:   PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2461:   PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2462:   TSIFunction   ifun;
2463:   TSIJacobian   ijac;
2464:   TSI2Jacobian  i2jac;
2465:   TSRHSJacobian rhsjac;

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

2470:   if (!((PetscObject)ts)->type_name) {
2471:     TSGetIFunction(ts, NULL, &ifun, NULL);
2472:     TSSetType(ts, ifun ? TSBEULER : TSEULER);
2473:   }

2475:   if (!ts->vec_sol) {
2477:     DMCreateGlobalVector(ts->dm, &ts->vec_sol);
2478:   }

2480:   if (ts->tspan) {
2481:     if (!ts->tspan->vecs_sol) VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol);
2482:   }
2483:   if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2484:     PetscObjectReference((PetscObject)ts->Jacprhs);
2485:     ts->Jacp = ts->Jacprhs;
2486:   }

2488:   if (ts->quadraturets) {
2489:     TSSetUp(ts->quadraturets);
2490:     VecDestroy(&ts->vec_costintegrand);
2491:     VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand);
2492:   }

2494:   TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL);
2495:   if (rhsjac == TSComputeRHSJacobianConstant) {
2496:     Mat  Amat, Pmat;
2497:     SNES snes;
2498:     TSGetSNES(ts, &snes);
2499:     SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL);
2500:     /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2501:      * have displaced the RHS matrix */
2502:     if (Amat && Amat == ts->Arhs) {
2503:       /* 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 */
2504:       MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat);
2505:       SNESSetJacobian(snes, Amat, NULL, NULL, NULL);
2506:       MatDestroy(&Amat);
2507:     }
2508:     if (Pmat && Pmat == ts->Brhs) {
2509:       MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat);
2510:       SNESSetJacobian(snes, NULL, Pmat, NULL, NULL);
2511:       MatDestroy(&Pmat);
2512:     }
2513:   }

2515:   TSGetAdapt(ts, &ts->adapt);
2516:   TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type);

2518:   PetscTryTypeMethod(ts, setup);

2520:   TSSetExactFinalTimeDefault(ts);

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

2529:   /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2530:      Otherwise, the SNES will use coloring internally to form the Jacobian.
2531:    */
2532:   DMSNESGetJacobian(dm, &jac, NULL);
2533:   DMTSGetIJacobian(dm, &ijac, NULL);
2534:   DMTSGetI2Jacobian(dm, &i2jac, NULL);
2535:   DMTSGetRHSJacobian(dm, &rhsjac, NULL);
2536:   if (!jac && (ijac || i2jac || rhsjac)) DMSNESSetJacobian(dm, SNESTSFormJacobian, ts);

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

2541:   ts->setupcalled = PETSC_TRUE;
2542:   return 0;
2543: }

2545: /*@
2546:    TSReset - Resets a `TS` context and removes any allocated `Vec`s and `Mat`s.

2548:    Collective

2550:    Input Parameter:
2551: .  ts - the `TS` context obtained from `TSCreate()`

2553:    Level: beginner

2555: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()`
2556: @*/
2557: PetscErrorCode TSReset(TS ts)
2558: {
2559:   TS_RHSSplitLink ilink = ts->tsrhssplit, next;


2563:   PetscTryTypeMethod(ts, reset);
2564:   if (ts->snes) SNESReset(ts->snes);
2565:   if (ts->adapt) TSAdaptReset(ts->adapt);

2567:   MatDestroy(&ts->Arhs);
2568:   MatDestroy(&ts->Brhs);
2569:   VecDestroy(&ts->Frhs);
2570:   VecDestroy(&ts->vec_sol);
2571:   VecDestroy(&ts->vec_dot);
2572:   VecDestroy(&ts->vatol);
2573:   VecDestroy(&ts->vrtol);
2574:   VecDestroyVecs(ts->nwork, &ts->work);

2576:   MatDestroy(&ts->Jacprhs);
2577:   MatDestroy(&ts->Jacp);
2578:   if (ts->forward_solve) TSForwardReset(ts);
2579:   if (ts->quadraturets) {
2580:     TSReset(ts->quadraturets);
2581:     VecDestroy(&ts->vec_costintegrand);
2582:   }
2583:   while (ilink) {
2584:     next = ilink->next;
2585:     TSDestroy(&ilink->ts);
2586:     PetscFree(ilink->splitname);
2587:     ISDestroy(&ilink->is);
2588:     PetscFree(ilink);
2589:     ilink = next;
2590:   }
2591:   ts->tsrhssplit     = NULL;
2592:   ts->num_rhs_splits = 0;
2593:   if (ts->tspan) {
2594:     PetscFree(ts->tspan->span_times);
2595:     VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol);
2596:     PetscFree(ts->tspan);
2597:   }
2598:   ts->setupcalled = PETSC_FALSE;
2599:   return 0;
2600: }

2602: /*@C
2603:    TSDestroy - Destroys the timestepper context that was created
2604:    with `TSCreate()`.

2606:    Collective

2608:    Input Parameter:
2609: .  ts - the `TS` context obtained from `TSCreate()`

2611:    Level: beginner

2613: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2614: @*/
2615: PetscErrorCode TSDestroy(TS *ts)
2616: {
2617:   if (!*ts) return 0;
2619:   if (--((PetscObject)(*ts))->refct > 0) {
2620:     *ts = NULL;
2621:     return 0;
2622:   }

2624:   TSReset(*ts);
2625:   TSAdjointReset(*ts);
2626:   if ((*ts)->forward_solve) TSForwardReset(*ts);

2628:   /* if memory was published with SAWs then destroy it */
2629:   PetscObjectSAWsViewOff((PetscObject)*ts);
2630:   PetscTryTypeMethod((*ts), destroy);

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

2634:   TSAdaptDestroy(&(*ts)->adapt);
2635:   TSEventDestroy(&(*ts)->event);

2637:   SNESDestroy(&(*ts)->snes);
2638:   DMDestroy(&(*ts)->dm);
2639:   TSMonitorCancel((*ts));
2640:   TSAdjointMonitorCancel((*ts));

2642:   TSDestroy(&(*ts)->quadraturets);
2643:   PetscHeaderDestroy(ts);
2644:   return 0;
2645: }

2647: /*@
2648:    TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2649:    a `TS` (timestepper) context. Valid only for nonlinear problems.

2651:    Not Collective, but snes is parallel if ts is parallel

2653:    Input Parameter:
2654: .  ts - the `TS` context obtained from `TSCreate()`

2656:    Output Parameter:
2657: .  snes - the nonlinear solver context

2659:    Level: beginner

2661:    Notes:
2662:    The user can then directly manipulate the `SNES` context to set various
2663:    options, etc.  Likewise, the user can then extract and manipulate the
2664:    `KSP`, and `PC` contexts as well.

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

2669: .seealso: [](chapter_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2670: @*/
2671: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2672: {
2675:   if (!ts->snes) {
2676:     SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes);
2677:     PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options);
2678:     SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts);
2679:     PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1);
2680:     if (ts->dm) SNESSetDM(ts->snes, ts->dm);
2681:     if (ts->problem_type == TS_LINEAR) SNESSetType(ts->snes, SNESKSPONLY);
2682:   }
2683:   *snes = ts->snes;
2684:   return 0;
2685: }

2687: /*@
2688:    TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the timestepping context

2690:    Collective

2692:    Input Parameters:
2693: +  ts - the `TS` context obtained from `TSCreate()`
2694: -  snes - the nonlinear solver context

2696:    Level: developer

2698:    Note:
2699:    Most users should have the `TS` created by calling `TSGetSNES()`

2701: .seealso: [](chapter_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2702: @*/
2703: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2704: {
2705:   PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);

2709:   PetscObjectReference((PetscObject)snes);
2710:   SNESDestroy(&ts->snes);

2712:   ts->snes = snes;

2714:   SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts);
2715:   SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL);
2716:   if (func == SNESTSFormJacobian) SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts);
2717:   return 0;
2718: }

2720: /*@
2721:    TSGetKSP - Returns the `KSP` (linear solver) associated with
2722:    a `TS` (timestepper) context.

2724:    Not Collective, but ksp is parallel if ts is parallel

2726:    Input Parameter:
2727: .  ts - the `TS` context obtained from `TSCreate()`

2729:    Output Parameter:
2730: .  ksp - the nonlinear solver context

2732:    Level: beginner

2734:    Notes:
2735:    The user can then directly manipulate the `KSP` context to set various
2736:    options, etc.  Likewise, the user can then extract and manipulate the
2737:    `PC` context as well.

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

2742: .seealso: [](chapter_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2743: @*/
2744: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2745: {
2746:   SNES snes;

2752:   TSGetSNES(ts, &snes);
2753:   SNESGetKSP(snes, ksp);
2754:   return 0;
2755: }

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

2759: /*@
2760:    TSSetMaxSteps - Sets the maximum number of steps to use.

2762:    Logically Collective

2764:    Input Parameters:
2765: +  ts - the `TS` context obtained from `TSCreate()`
2766: -  maxsteps - maximum number of steps to use

2768:    Options Database Key:
2769: .  -ts_max_steps <maxsteps> - Sets maxsteps

2771:    Level: intermediate

2773:    Note:
2774:    The default maximum number of steps is 5000

2776: .seealso: [](chapter_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2777: @*/
2778: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2779: {
2783:   ts->max_steps = maxsteps;
2784:   return 0;
2785: }

2787: /*@
2788:    TSGetMaxSteps - Gets the maximum number of steps to use.

2790:    Not Collective

2792:    Input Parameters:
2793: .  ts - the `TS` context obtained from `TSCreate()`

2795:    Output Parameter:
2796: .  maxsteps - maximum number of steps to use

2798:    Level: advanced

2800: .seealso: [](chapter_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2801: @*/
2802: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2803: {
2806:   *maxsteps = ts->max_steps;
2807:   return 0;
2808: }

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

2813:    Logically Collective

2815:    Input Parameters:
2816: +  ts - the `TS` context obtained from `TSCreate()`
2817: -  maxtime - final time to step to

2819:    Options Database Key:
2820: .  -ts_max_time <maxtime> - Sets maxtime

2822:    Level: intermediate

2824:    Notes:
2825:    The default maximum time is 5.0

2827: .seealso: [](chapter_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2828: @*/
2829: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2830: {
2833:   ts->max_time = maxtime;
2834:   return 0;
2835: }

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

2840:    Not Collective

2842:    Input Parameters:
2843: .  ts - the `TS` context obtained from `TSCreate()`

2845:    Output Parameter:
2846: .  maxtime - final time to step to

2848:    Level: advanced

2850: .seealso: [](chapter_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2851: @*/
2852: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2853: {
2856:   *maxtime = ts->max_time;
2857:   return 0;
2858: }

2860: /*@
2861:    TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.

2863:    Level: deprecated

2865: @*/
2866: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2867: {
2869:   TSSetTime(ts, initial_time);
2870:   TSSetTimeStep(ts, time_step);
2871:   return 0;
2872: }

2874: /*@
2875:    TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.

2877:    Level: deprecated

2879: @*/
2880: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2881: {
2883:   if (maxsteps) {
2885:     *maxsteps = ts->max_steps;
2886:   }
2887:   if (maxtime) {
2889:     *maxtime = ts->max_time;
2890:   }
2891:   return 0;
2892: }

2894: /*@
2895:    TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.

2897:    Level: deprecated

2899: @*/
2900: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2901: {
2905:   if (maxsteps >= 0) ts->max_steps = maxsteps;
2906:   if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime;
2907:   return 0;
2908: }

2910: /*@
2911:    TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.

2913:    Level: deprecated

2915: @*/
2916: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
2917: {
2918:   return TSGetStepNumber(ts, steps);
2919: }

2921: /*@
2922:    TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.

2924:    Level: deprecated

2926: @*/
2927: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
2928: {
2929:   return TSGetStepNumber(ts, steps);
2930: }

2932: /*@
2933:    TSSetSolution - Sets the initial solution vector
2934:    for use by the `TS` routines.

2936:    Logically Collective

2938:    Input Parameters:
2939: +  ts - the `TS` context obtained from `TSCreate()`
2940: -  u - the solution vector

2942:    Level: beginner

2944: .seealso: [](chapter_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
2945: @*/
2946: PetscErrorCode TSSetSolution(TS ts, Vec u)
2947: {
2948:   DM dm;

2952:   PetscObjectReference((PetscObject)u);
2953:   VecDestroy(&ts->vec_sol);
2954:   ts->vec_sol = u;

2956:   TSGetDM(ts, &dm);
2957:   DMShellSetGlobalVector(dm, u);
2958:   return 0;
2959: }

2961: /*@C
2962:   TSSetPreStep - Sets the general-purpose function
2963:   called once at the beginning of each time step.

2965:   Logically Collective

2967:   Input Parameters:
2968: + ts   - The `TS` context obtained from `TSCreate()`
2969: - func - The function

2971:   Calling sequence of func:
2972: .vb
2973:   PetscErrorCode func (TS ts);
2974: .ve

2976:   Level: intermediate

2978: .seealso: [](chapter_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
2979: @*/
2980: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
2981: {
2983:   ts->prestep = func;
2984:   return 0;
2985: }

2987: /*@
2988:   TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`

2990:   Collective

2992:   Input Parameters:
2993: . ts   - The `TS` context obtained from `TSCreate()`

2995:   Level: developer

2997:   Note:
2998:   `TSPreStep()` is typically used within time stepping implementations,
2999:   so most users would not generally call this routine themselves.

3001: .seealso: [](chapter_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3002: @*/
3003: PetscErrorCode TSPreStep(TS ts)
3004: {
3006:   if (ts->prestep) {
3007:     Vec              U;
3008:     PetscObjectId    idprev;
3009:     PetscBool        sameObject;
3010:     PetscObjectState sprev, spost;

3012:     TSGetSolution(ts, &U);
3013:     PetscObjectGetId((PetscObject)U, &idprev);
3014:     PetscObjectStateGet((PetscObject)U, &sprev);
3015:     PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3016:     TSGetSolution(ts, &U);
3017:     PetscObjectCompareId((PetscObject)U, idprev, &sameObject);
3018:     PetscObjectStateGet((PetscObject)U, &spost);
3019:     if (!sameObject || sprev != spost) TSRestartStep(ts);
3020:   }
3021:   return 0;
3022: }

3024: /*@C
3025:   TSSetPreStage - Sets the general-purpose function
3026:   called once at the beginning of each stage.

3028:   Logically Collective

3030:   Input Parameters:
3031: + ts   - The `TS` context obtained from `TSCreate()`
3032: - func - The function

3034:   Calling sequence of func:
3035: .vb
3036:   PetscErrorCode func(TS ts, PetscReal stagetime);
3037: .ve

3039:   Level: intermediate

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

3046: .seealso: [](chapter_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3047: @*/
3048: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS, PetscReal))
3049: {
3051:   ts->prestage = func;
3052:   return 0;
3053: }

3055: /*@C
3056:   TSSetPostStage - Sets the general-purpose function, provided with `TSSetPostStep()`,
3057:   called once at the end of each stage.

3059:   Logically Collective

3061:   Input Parameters:
3062: + ts   - The `TS` context obtained from `TSCreate()`
3063: - func - The function

3065:   Calling sequence of func:
3066: .vb
3067:   PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y);
3068: .ve

3070:   Level: intermediate

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

3077: .seealso: [](chapter_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3078: @*/
3079: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS, PetscReal, PetscInt, Vec *))
3080: {
3082:   ts->poststage = func;
3083:   return 0;
3084: }

3086: /*@C
3087:   TSSetPostEvaluate - Sets the general-purpose function
3088:   called once at the end of each step evaluation.

3090:   Logically Collective

3092:   Input Parameters:
3093: + ts   - The `TS` context obtained from `TSCreate()`
3094: - func - The function

3096:   Calling sequence of func:
3097: .vb
3098:   PetscErrorCode func(TS ts);
3099: .ve

3101:   Level: intermediate

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

3110: .seealso: [](chapter_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3111: @*/
3112: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS))
3113: {
3115:   ts->postevaluate = func;
3116:   return 0;
3117: }

3119: /*@
3120:   TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`

3122:   Collective

3124:   Input Parameters:
3125: . ts          - The `TS` context obtained from `TSCreate()`
3126:   stagetime   - The absolute time of the current stage

3128:   Level: developer

3130:   Note:
3131:   `TSPreStage()` is typically used within time stepping implementations,
3132:   most users would not generally call this routine themselves.

3134: .seealso: [](chapter_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3135: @*/
3136: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3137: {
3139:   if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3140:   return 0;
3141: }

3143: /*@
3144:   TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`

3146:   Collective

3148:   Input Parameters:
3149: . ts          - The `TS` context obtained from `TSCreate()`
3150:   stagetime   - The absolute time of the current stage
3151:   stageindex  - Stage number
3152:   Y           - Array of vectors (of size = total number
3153:                 of stages) with the stage solutions

3155:   Level: developer

3157:   Note:
3158:   `TSPostStage()` is typically used within time stepping implementations,
3159:   most users would not generally call this routine themselves.

3161: .seealso: [](chapter_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3162: @*/
3163: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3164: {
3166:   if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3167:   return 0;
3168: }

3170: /*@
3171:   TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`

3173:   Collective

3175:   Input Parameters:
3176: . ts - The `TS` context obtained from `TSCreate()`

3178:   Level: developer

3180:   Note:
3181:   `TSPostEvaluate()` is typically used within time stepping implementations,
3182:   most users would not generally call this routine themselves.

3184: .seealso: [](chapter_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3185: @*/
3186: PetscErrorCode TSPostEvaluate(TS ts)
3187: {
3189:   if (ts->postevaluate) {
3190:     Vec              U;
3191:     PetscObjectState sprev, spost;

3193:     TSGetSolution(ts, &U);
3194:     PetscObjectStateGet((PetscObject)U, &sprev);
3195:     PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3196:     PetscObjectStateGet((PetscObject)U, &spost);
3197:     if (sprev != spost) TSRestartStep(ts);
3198:   }
3199:   return 0;
3200: }

3202: /*@C
3203:   TSSetPostStep - Sets the general-purpose function
3204:   called once at the end of each time step.

3206:   Logically Collective

3208:   Input Parameters:
3209: + ts   - The `TS` context obtained from `TSCreate()`
3210: - func - The function

3212:   Calling sequence of func:
3213: $ func (TS ts);

3215:   Level: intermediate

3217:   Note:
3218:   The function set by `TSSetPostStep()` is called after each successful step. The solution vector X
3219:   obtained by `TSGetSolution()` may be different than that computed at the step end if the event handler
3220:   locates an event and `TSPostEvent()` modifies it. Use `TSSetPostEvaluate()` if an unmodified solution is needed instead.

3222: .seealso: [](chapter_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3223: @*/
3224: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
3225: {
3227:   ts->poststep = func;
3228:   return 0;
3229: }

3231: /*@
3232:   TSPostStep - Runs the user-defined post-step function that was set with `TSSetPotsStep()`

3234:   Collective

3236:   Input Parameters:
3237: . ts   - The `TS` context obtained from `TSCreate()`

3239:   Note:
3240:   `TSPostStep()` is typically used within time stepping implementations,
3241:   so most users would not generally call this routine themselves.

3243:   Level: developer

3245: .seealso: [](chapter_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()`
3246: @*/
3247: PetscErrorCode TSPostStep(TS ts)
3248: {
3250:   if (ts->poststep) {
3251:     Vec              U;
3252:     PetscObjectId    idprev;
3253:     PetscBool        sameObject;
3254:     PetscObjectState sprev, spost;

3256:     TSGetSolution(ts, &U);
3257:     PetscObjectGetId((PetscObject)U, &idprev);
3258:     PetscObjectStateGet((PetscObject)U, &sprev);
3259:     PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3260:     TSGetSolution(ts, &U);
3261:     PetscObjectCompareId((PetscObject)U, idprev, &sameObject);
3262:     PetscObjectStateGet((PetscObject)U, &spost);
3263:     if (!sameObject || sprev != spost) TSRestartStep(ts);
3264:   }
3265:   return 0;
3266: }

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

3271:    Collective

3273:    Input Parameters:
3274: +  ts - time stepping context
3275: -  t - time to interpolate to

3277:    Output Parameter:
3278: .  U - state at given time

3280:    Level: intermediate

3282:    Developer Note:
3283:    `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.

3285: .seealso: [](chapter_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3286: @*/
3287: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3288: {
3292:   PetscUseTypeMethod(ts, interpolate, t, U);
3293:   return 0;
3294: }

3296: /*@
3297:    TSStep - Steps one time step

3299:    Collective

3301:    Input Parameter:
3302: .  ts - the `TS` context obtained from `TSCreate()`

3304:    Level: developer

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

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

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

3315: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3316: @*/
3317: PetscErrorCode TSStep(TS ts)
3318: {
3319:   static PetscBool cite = PETSC_FALSE;
3320:   PetscReal        ptime;

3323:   PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3324:                                    "  title         = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3325:                                    "  author        = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3326:                                    "  journal       = {arXiv e-preprints},\n"
3327:                                    "  eprint        = {1806.01437},\n"
3328:                                    "  archivePrefix = {arXiv},\n"
3329:                                    "  year          = {2018}\n}\n",
3330:                                    &cite));
3331:   TSSetUp(ts);
3332:   TSTrajectorySetUp(ts->trajectory, ts);


3338:   if (!ts->steps) ts->ptime_prev = ts->ptime;
3339:   ptime                   = ts->ptime;
3340:   ts->ptime_prev_rollback = ts->ptime_prev;
3341:   ts->reason              = TS_CONVERGED_ITERATING;

3343:   PetscLogEventBegin(TS_Step, ts, 0, 0, 0);
3344:   PetscUseTypeMethod(ts, step);
3345:   PetscLogEventEnd(TS_Step, ts, 0, 0, 0);

3347:   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)
3348:     VecCopy(ts->vec_sol, ts->tspan->vecs_sol[ts->tspan->spanctr++]);
3349:   if (ts->reason >= 0) {
3350:     ts->ptime_prev = ptime;
3351:     ts->steps++;
3352:     ts->steprollback = PETSC_FALSE;
3353:     ts->steprestart  = PETSC_FALSE;
3354:   }
3355:   if (!ts->reason) {
3356:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3357:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3358:   }

3362:   return 0;
3363: }

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

3369:    Collective

3371:    Input Parameters:
3372: +  ts - time stepping context
3373: -  wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

3375:    Input/Output Parameter:
3376: .  order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3377:            on output, the actual order of the error evaluation

3379:    Output Parameter:
3380: .  wlte - the weighted local truncation error norm

3382:    Level: advanced

3384:    Note:
3385:    If the timestepper cannot evaluate the error in a particular step
3386:    (eg. in the first step or restart steps after event handling),
3387:    this routine returns wlte=-1.0 .

3389: .seealso: [](chapter_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3390: @*/
3391: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3392: {
3400:   PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3401:   return 0;
3402: }

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

3407:    Collective

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

3414:    Output Parameter:
3415: .  U - state at the end of the current step

3417:    Level: advanced

3419:    Notes:
3420:    This function cannot be called until all stages have been evaluated.

3422:    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.

3424: .seealso: [](chapter_ts), `TS`, `TSStep()`, `TSAdapt`
3425: @*/
3426: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3427: {
3431:   PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3432:   return 0;
3433: }

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

3438:   Not collective

3440:   Input Parameter:
3441: . ts - time stepping context

3443:   Output Parameter:
3444: . initConditions - The function which computes an initial condition

3446:   The calling sequence for the function is
3447: .vb
3448:  initCondition(TS ts, Vec u)
3449:  ts - The timestepping context
3450:  u  - The input vector in which the initial condition is stored
3451: .ve

3453:    Level: advanced

3455: .seealso: [](chapter_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3456: @*/
3457: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS, Vec))
3458: {
3461:   *initCondition = ts->ops->initcondition;
3462:   return 0;
3463: }

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

3468:   Logically collective

3470:   Input Parameters:
3471: + ts  - time stepping context
3472: - initCondition - The function which computes an initial condition

3474:   Calling sequence for initCondition:
3475: $ PetscErrorCode initCondition(TS ts, Vec u)
3476: + ts - The timestepping context
3477: - u  - The input vector in which the initial condition is to be stored

3479:   Level: advanced

3481: .seealso: [](chapter_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3482: @*/
3483: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS, Vec))
3484: {
3487:   ts->ops->initcondition = initCondition;
3488:   return 0;
3489: }

3491: /*@
3492:   TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`

3494:   Collective

3496:   Input Parameters:
3497: + ts - time stepping context
3498: - u  - The `Vec` to store the condition in which will be used in `TSSolve()`

3500:   Level: advanced

3502: .seealso: [](chapter_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3503: @*/
3504: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3505: {
3508:   PetscTryTypeMethod(ts, initcondition, u);
3509:   return 0;
3510: }

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

3515:   Not collective

3517:   Input Parameter:
3518: . ts - time stepping context

3520:   Output Parameter:
3521: . exactError - The function which computes the solution error

3523:   Calling sequence for exactError:
3524: $ PetscErrorCode exactError(TS ts, Vec u)
3525: + ts - The timestepping context
3526: . u  - The approximate solution vector
3527: - e  - The input vector in which the error is stored

3529:   Level: advanced

3531: .seealso: [](chapter_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3532: @*/
3533: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS, Vec, Vec))
3534: {
3537:   *exactError = ts->ops->exacterror;
3538:   return 0;
3539: }

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

3544:   Logically collective

3546:   Input Parameters:
3547: + ts - time stepping context
3548: - exactError - The function which computes the solution error

3550:   Calling sequence for exactError:
3551: $ PetscErrorCode exactError(TS ts, Vec u)
3552: + ts - The timestepping context
3553: . u  - The approximate solution vector
3554: - e  - The input vector in which the error is stored

3556:   Level: advanced

3558: .seealso: [](chapter_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3559: @*/
3560: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS, Vec, Vec))
3561: {
3564:   ts->ops->exacterror = exactError;
3565:   return 0;
3566: }

3568: /*@
3569:   TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`

3571:   Collective

3573:   Input Parameters:
3574: + ts - time stepping context
3575: . u  - The approximate solution
3576: - e  - The `Vec` used to store the error

3578:   Level: advanced

3580: .seealso: [](chapter_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3581: @*/
3582: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3583: {
3587:   PetscTryTypeMethod(ts, exacterror, u, e);
3588:   return 0;
3589: }

3591: /*@
3592:    TSSolve - Steps the requested number of timesteps.

3594:    Collective

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

3601:    Level: beginner

3603:    Notes:
3604:    The final time returned by this function may be different from the time of the internally
3605:    held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
3606:    stepped over the final time.

3608: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3609: @*/
3610: PetscErrorCode TSSolve(TS ts, Vec u)
3611: {
3612:   Vec solution;


3617:   TSSetExactFinalTimeDefault(ts);
3618:   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 */
3619:     if (!ts->vec_sol || u == ts->vec_sol) {
3620:       VecDuplicate(u, &solution);
3621:       TSSetSolution(ts, solution);
3622:       VecDestroy(&solution); /* grant ownership */
3623:     }
3624:     VecCopy(u, ts->vec_sol);
3626:   } else if (u) TSSetSolution(ts, u);
3627:   TSSetUp(ts);
3628:   TSTrajectorySetUp(ts->trajectory, ts);


3635:   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 */
3636:     VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]);
3637:     ts->tspan->spanctr = 1;
3638:   }

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

3642:   /* reset number of steps only when the step is not restarted. ARKIMEX
3643:      restarts the step after an event. Resetting these counters in such case causes
3644:      TSTrajectory to incorrectly save the output files
3645:   */
3646:   /* reset time step and iteration counters */
3647:   if (!ts->steps) {
3648:     ts->ksp_its           = 0;
3649:     ts->snes_its          = 0;
3650:     ts->num_snes_failures = 0;
3651:     ts->reject            = 0;
3652:     ts->steprestart       = PETSC_TRUE;
3653:     ts->steprollback      = PETSC_FALSE;
3654:     ts->rhsjacobian.time  = PETSC_MIN_REAL;
3655:   }

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

3662:     if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
3663:     else maxdt = ts->max_time - ts->ptime;
3664:     ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
3665:   }
3666:   ts->reason = TS_CONVERGED_ITERATING;

3668:   {
3669:     PetscViewer       viewer;
3670:     PetscViewerFormat format;
3671:     PetscBool         flg;
3672:     static PetscBool  incall = PETSC_FALSE;

3674:     if (!incall) {
3675:       /* Estimate the convergence rate of the time discretization */
3676:       PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg);
3677:       if (flg) {
3678:         PetscConvEst conv;
3679:         DM           dm;
3680:         PetscReal   *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
3681:         PetscInt     Nf;
3682:         PetscBool    checkTemporal = PETSC_TRUE;

3684:         incall = PETSC_TRUE;
3685:         PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg);
3686:         TSGetDM(ts, &dm);
3687:         DMGetNumFields(dm, &Nf);
3688:         PetscCalloc1(PetscMax(Nf, 1), &alpha);
3689:         PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv);
3690:         PetscConvEstUseTS(conv, checkTemporal);
3691:         PetscConvEstSetSolver(conv, (PetscObject)ts);
3692:         PetscConvEstSetFromOptions(conv);
3693:         PetscConvEstSetUp(conv);
3694:         PetscConvEstGetConvRate(conv, alpha);
3695:         PetscViewerPushFormat(viewer, format);
3696:         PetscConvEstRateView(conv, alpha, viewer);
3697:         PetscViewerPopFormat(viewer);
3698:         PetscViewerDestroy(&viewer);
3699:         PetscConvEstDestroy(&conv);
3700:         PetscFree(alpha);
3701:         incall = PETSC_FALSE;
3702:       }
3703:     }
3704:   }

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

3708:   if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
3709:     PetscUseTypeMethod(ts, solve);
3710:     if (u) VecCopy(ts->vec_sol, u);
3711:     ts->solvetime = ts->ptime;
3712:     solution      = ts->vec_sol;
3713:   } else { /* Step the requested number of timesteps. */
3714:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3715:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;

3717:     if (!ts->steps) {
3718:       TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol);
3719:       TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol);
3720:     }

3722:     while (!ts->reason) {
3723:       TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);
3724:       if (!ts->steprollback) TSPreStep(ts);
3725:       TSStep(ts);
3726:       if (ts->testjacobian) TSRHSJacobianTest(ts, NULL);
3727:       if (ts->testjacobiantranspose) TSRHSJacobianTestTranspose(ts, NULL);
3728:       if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
3729:         if (ts->reason >= 0) ts->steps--;            /* Revert the step number changed by TSStep() */
3730:         TSForwardCostIntegral(ts);
3731:         if (ts->reason >= 0) ts->steps++;
3732:       }
3733:       if (ts->forward_solve) {            /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
3734:         if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
3735:         TSForwardStep(ts);
3736:         if (ts->reason >= 0) ts->steps++;
3737:       }
3738:       TSPostEvaluate(ts);
3739:       TSEventHandler(ts); /* The right-hand side may be changed due to event. Be careful with Any computation using the RHS information after this point. */
3740:       if (ts->steprollback) TSPostEvaluate(ts);
3741:       if (!ts->steprollback) {
3742:         TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol);
3743:         TSPostStep(ts);
3744:       }
3745:     }
3746:     TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);

3748:     if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
3749:       TSInterpolate(ts, ts->max_time, u);
3750:       ts->solvetime = ts->max_time;
3751:       solution      = u;
3752:       TSMonitor(ts, -1, ts->solvetime, solution);
3753:     } else {
3754:       if (u) VecCopy(ts->vec_sol, u);
3755:       ts->solvetime = ts->ptime;
3756:       solution      = ts->vec_sol;
3757:     }
3758:   }

3760:   TSViewFromOptions(ts, NULL, "-ts_view");
3761:   VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution");
3762:   PetscObjectSAWsBlock((PetscObject)ts);
3763:   if (ts->adjoint_solve) TSAdjointSolve(ts);
3764:   return 0;
3765: }

3767: /*@
3768:    TSGetTime - Gets the time of the most recently completed step.

3770:    Not Collective

3772:    Input Parameter:
3773: .  ts - the `TS` context obtained from `TSCreate()`

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

3778:    Level: beginner

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

3784: .seealso: [](chapter_ts), TS`, ``TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
3785: @*/
3786: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
3787: {
3790:   *t = ts->ptime;
3791:   return 0;
3792: }

3794: /*@
3795:    TSGetPrevTime - Gets the starting time of the previously completed step.

3797:    Not Collective

3799:    Input Parameter:
3800: .  ts - the `TS` context obtained from `TSCreate()`

3802:    Output Parameter:
3803: .  t  - the previous time

3805:    Level: beginner

3807: .seealso: [](chapter_ts), TS`, ``TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
3808: @*/
3809: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
3810: {
3813:   *t = ts->ptime_prev;
3814:   return 0;
3815: }

3817: /*@
3818:    TSSetTime - Allows one to reset the time.

3820:    Logically Collective

3822:    Input Parameters:
3823: +  ts - the `TS` context obtained from `TSCreate()`
3824: -  time - the time

3826:    Level: intermediate

3828: .seealso: [](chapter_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
3829: @*/
3830: PetscErrorCode TSSetTime(TS ts, PetscReal t)
3831: {
3834:   ts->ptime = t;
3835:   return 0;
3836: }

3838: /*@C
3839:    TSSetOptionsPrefix - Sets the prefix used for searching for all
3840:    TS options in the database.

3842:    Logically Collective

3844:    Input Parameters:
3845: +  ts     - The `TS` context
3846: -  prefix - The prefix to prepend to all option names

3848:    Level: advanced

3850:    Note:
3851:    A hyphen (-) must NOT be given at the beginning of the prefix name.
3852:    The first character of all runtime options is AUTOMATICALLY the
3853:    hyphen.

3855: .seealso: [](chapter_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
3856: @*/
3857: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
3858: {
3859:   SNES snes;

3862:   PetscObjectSetOptionsPrefix((PetscObject)ts, prefix);
3863:   TSGetSNES(ts, &snes);
3864:   SNESSetOptionsPrefix(snes, prefix);
3865:   return 0;
3866: }

3868: /*@C
3869:    TSAppendOptionsPrefix - Appends to the prefix used for searching for all
3870:    TS options in the database.

3872:    Logically Collective

3874:    Input Parameters:
3875: +  ts     - The `TS` context
3876: -  prefix - The prefix to prepend to all option names

3878:    Level: advanced

3880:    Note:
3881:    A hyphen (-) must NOT be given at the beginning of the prefix name.
3882:    The first character of all runtime options is AUTOMATICALLY the
3883:    hyphen.

3885: .seealso: [](chapter_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
3886: @*/
3887: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
3888: {
3889:   SNES snes;

3892:   PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix);
3893:   TSGetSNES(ts, &snes);
3894:   SNESAppendOptionsPrefix(snes, prefix);
3895:   return 0;
3896: }

3898: /*@C
3899:    TSGetOptionsPrefix - Sets the prefix used for searching for all
3900:    `TS` options in the database.

3902:    Not Collective

3904:    Input Parameter:
3905: .  ts - The `TS` context

3907:    Output Parameter:
3908: .  prefix - A pointer to the prefix string used

3910:    Level: intermediate

3912:    Fortran Note:
3913:    On the fortran side, the user should pass in a string 'prefix' of
3914:    sufficient length to hold the prefix.

3916: .seealso: [](chapter_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
3917: @*/
3918: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
3919: {
3922:   PetscObjectGetOptionsPrefix((PetscObject)ts, prefix);
3923:   return 0;
3924: }

3926: /*@C
3927:    TSGetRHSJacobian - Returns the Jacobian J at the present timestep.

3929:    Not Collective, but parallel objects are returned if ts is parallel

3931:    Input Parameter:
3932: .  ts  - The `TS` context obtained from `TSCreate()`

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

3940:    Level: intermediate

3942:    Note:
3943:     You can pass in NULL for any return argument you do not need.

3945: .seealso: [](chapter_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`

3947: @*/
3948: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobian *func, void **ctx)
3949: {
3950:   DM dm;

3952:   if (Amat || Pmat) {
3953:     SNES snes;
3954:     TSGetSNES(ts, &snes);
3955:     SNESSetUpMatrices(snes);
3956:     SNESGetJacobian(snes, Amat, Pmat, NULL, NULL);
3957:   }
3958:   TSGetDM(ts, &dm);
3959:   DMTSGetRHSJacobian(dm, func, ctx);
3960:   return 0;
3961: }

3963: /*@C
3964:    TSGetIJacobian - Returns the implicit Jacobian at the present timestep.

3966:    Not Collective, but parallel objects are returned if ts is parallel

3968:    Input Parameter:
3969: .  ts  - The `TS` context obtained from `TSCreate()`

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

3977:    Level: advanced

3979:    Note:
3980:     You can pass in NULL for any return argument you do not need.

3982: .seealso: [](chapter_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`

3984: @*/
3985: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobian *f, void **ctx)
3986: {
3987:   DM dm;

3989:   if (Amat || Pmat) {
3990:     SNES snes;
3991:     TSGetSNES(ts, &snes);
3992:     SNESSetUpMatrices(snes);
3993:     SNESGetJacobian(snes, Amat, Pmat, NULL, NULL);
3994:   }
3995:   TSGetDM(ts, &dm);
3996:   DMTSGetIJacobian(dm, f, ctx);
3997:   return 0;
3998: }

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

4004:    Logically Collective

4006:    Input Parameters:
4007: +  ts - the `TS` integrator object
4008: -  dm - the dm, cannot be NULL

4010:    Level: intermediate

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

4017: .seealso: [](chapter_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4018: @*/
4019: PetscErrorCode TSSetDM(TS ts, DM dm)
4020: {
4021:   SNES snes;
4022:   DMTS tsdm;

4026:   PetscObjectReference((PetscObject)dm);
4027:   if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4028:     if (ts->dm->dmts && !dm->dmts) {
4029:       DMCopyDMTS(ts->dm, dm);
4030:       DMGetDMTS(ts->dm, &tsdm);
4031:       /* Grant write privileges to the replacement DM */
4032:       if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4033:     }
4034:     DMDestroy(&ts->dm);
4035:   }
4036:   ts->dm = dm;

4038:   TSGetSNES(ts, &snes);
4039:   SNESSetDM(snes, dm);
4040:   return 0;
4041: }

4043: /*@
4044:    TSGetDM - Gets the `DM` that may be used by some preconditioners

4046:    Not Collective

4048:    Input Parameter:
4049: . ts - the `TS`

4051:    Output Parameter:
4052: .  dm - the dm

4054:    Level: intermediate

4056: .seealso: [](chapter_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4057: @*/
4058: PetscErrorCode TSGetDM(TS ts, DM *dm)
4059: {
4061:   if (!ts->dm) {
4062:     DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm);
4063:     if (ts->snes) SNESSetDM(ts->snes, ts->dm);
4064:   }
4065:   *dm = ts->dm;
4066:   return 0;
4067: }

4069: /*@
4070:    SNESTSFormFunction - Function to evaluate nonlinear residual

4072:    Logically Collective

4074:    Input Parameters:
4075: + snes - nonlinear solver
4076: . U - the current state at which to evaluate the residual
4077: - ctx - user context, must be a TS

4079:    Output Parameter:
4080: . F - the nonlinear residual

4082:    Level: advanced

4084:    Note:
4085:    This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4086:    It is most frequently passed to `MatFDColoringSetFunction()`.

4088: .seealso: [](chapter_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4089: @*/
4090: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4091: {
4092:   TS ts = (TS)ctx;

4098:   (ts->ops->snesfunction)(snes, U, F, ts);
4099:   return 0;
4100: }

4102: /*@
4103:    SNESTSFormJacobian - Function to evaluate the Jacobian

4105:    Collective

4107:    Input Parameters:
4108: + snes - nonlinear solver
4109: . U - the current state at which to evaluate the residual
4110: - ctx - user context, must be a `TS`

4112:    Output Parameters:
4113: + A - the Jacobian
4114: - B - the preconditioning matrix (may be the same as A)

4116:    Level: developer

4118:    Note:
4119:    This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.

4121: .seealso: [](chapter_ts), `SNESSetJacobian()`
4122: @*/
4123: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4124: {
4125:   TS ts = (TS)ctx;

4134:   (ts->ops->snesjacobian)(snes, U, A, B, ts);
4135:   return 0;
4136: }

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

4141:    Collective

4143:    Input Parameters:
4144: +  ts - time stepping context
4145: .  t - time at which to evaluate
4146: .  U - state at which to evaluate
4147: -  ctx - context

4149:    Output Parameter:
4150: .  F - right hand side

4152:    Level: intermediate

4154:    Note:
4155:    This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right hand side for linear problems.
4156:    The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.

4158: .seealso: [](chapter_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4159: @*/
4160: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4161: {
4162:   Mat Arhs, Brhs;

4164:   TSGetRHSMats_Private(ts, &Arhs, &Brhs);
4165:   /* undo the damage caused by shifting */
4166:   TSRecoverRHSJacobian(ts, Arhs, Brhs);
4167:   TSComputeRHSJacobian(ts, t, U, Arhs, Brhs);
4168:   MatMult(Arhs, U, F);
4169:   return 0;
4170: }

4172: /*@C
4173:    TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.

4175:    Collective

4177:    Input Parameters:
4178: +  ts - time stepping context
4179: .  t - time at which to evaluate
4180: .  U - state at which to evaluate
4181: -  ctx - context

4183:    Output Parameters:
4184: +  A - pointer to operator
4185: -  B - pointer to preconditioning matrix

4187:    Level: intermediate

4189:    Note:
4190:    This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.

4192: .seealso: [](chapter_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4193: @*/
4194: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4195: {
4196:   return 0;
4197: }

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

4202:    Collective

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

4211:    Output Parameter:
4212: .  F - left hand side

4214:    Level: intermediate

4216:    Notes:
4217:    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
4218:    user is required to write their own `TSComputeIFunction()`.
4219:    This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4220:    The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.

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

4224: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4225: @*/
4226: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4227: {
4228:   Mat A, B;

4230:   TSGetIJacobian(ts, &A, &B, NULL, NULL);
4231:   TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE);
4232:   MatMult(A, Udot, F);
4233:   return 0;
4234: }

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

4239:    Collective

4241:    Input Parameters:
4242: +  ts - time stepping context
4243: .  t - time at which to evaluate
4244: .  U - state at which to evaluate
4245: .  Udot - time derivative of state vector
4246: .  shift - shift to apply
4247: -  ctx - context

4249:    Output Parameters:
4250: +  A - pointer to operator
4251: -  B - pointer to preconditioning matrix

4253:    Level: advanced

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

4258:    It is only appropriate for problems of the form

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

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

4266: $    shift*M + J

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

4271: .seealso: [](chapter_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4272: @*/
4273: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4274: {
4275:   MatScale(A, shift / ts->ijacobian.shift);
4276:   ts->ijacobian.shift = shift;
4277:   return 0;
4278: }

4280: /*@
4281:    TSGetEquationType - Gets the type of the equation that `TS` is solving.

4283:    Not Collective

4285:    Input Parameter:
4286: .  ts - the `TS` context

4288:    Output Parameter:
4289: .  equation_type - see `TSEquationType`

4291:    Level: beginner

4293: .seealso: [](chapter_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4294: @*/
4295: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4296: {
4299:   *equation_type = ts->equation_type;
4300:   return 0;
4301: }

4303: /*@
4304:    TSSetEquationType - Sets the type of the equation that `TS` is solving.

4306:    Not Collective

4308:    Input Parameters:
4309: +  ts - the `TS` context
4310: -  equation_type - see `TSEquationType`

4312:    Level: advanced

4314: .seealso: [](chapter_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4315: @*/
4316: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4317: {
4319:   ts->equation_type = equation_type;
4320:   return 0;
4321: }

4323: /*@
4324:    TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.

4326:    Not Collective

4328:    Input Parameter:
4329: .  ts - the `TS` context

4331:    Output Parameter:
4332: .  reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4333:             manual pages for the individual convergence tests for complete lists

4335:    Level: beginner

4337:    Note:
4338:    Can only be called after the call to `TSSolve()` is complete.

4340: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4341: @*/
4342: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4343: {
4346:   *reason = ts->reason;
4347:   return 0;
4348: }

4350: /*@
4351:    TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.

4353:    Logically Collective; reason must contain common value

4355:    Input Parameters:
4356: +  ts - the `TS` context
4357: -  reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4358:             manual pages for the individual convergence tests for complete lists

4360:    Level: advanced

4362:    Note:
4363:    Can only be called while `TSSolve()` is active.

4365: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4366: @*/
4367: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4368: {
4370:   ts->reason = reason;
4371:   return 0;
4372: }

4374: /*@
4375:    TSGetSolveTime - Gets the time after a call to `TSSolve()`

4377:    Not Collective

4379:    Input Parameter:
4380: .  ts - the `TS` context

4382:    Output Parameter:
4383: .  ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`

4385:    Level: beginner

4387:    Note:
4388:    Can only be called after the call to `TSSolve()` is complete.

4390: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4391: @*/
4392: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4393: {
4396:   *ftime = ts->solvetime;
4397:   return 0;
4398: }

4400: /*@
4401:    TSGetSNESIterations - Gets the total number of nonlinear iterations
4402:    used by the time integrator.

4404:    Not Collective

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

4409:    Output Parameter:
4410: .  nits - number of nonlinear iterations

4412:    Level: intermediate

4414:    Notes:
4415:    This counter is reset to zero for each successive call to `TSSolve()`.

4417: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4418: @*/
4419: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4420: {
4423:   *nits = ts->snes_its;
4424:   return 0;
4425: }

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

4431:    Not Collective

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

4436:    Output Parameter:
4437: .  lits - number of linear iterations

4439:    Level: intermediate

4441:    Note:
4442:    This counter is reset to zero for each successive call to `TSSolve()`.

4444: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()`
4445: @*/
4446: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4447: {
4450:   *lits = ts->ksp_its;
4451:   return 0;
4452: }

4454: /*@
4455:    TSGetStepRejections - Gets the total number of rejected steps.

4457:    Not Collective

4459:    Input Parameter:
4460: .  ts - `TS` context

4462:    Output Parameter:
4463: .  rejects - number of steps rejected

4465:    Level: intermediate

4467:    Note:
4468:    This counter is reset to zero for each successive call to `TSSolve()`.

4470: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4471: @*/
4472: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4473: {
4476:   *rejects = ts->reject;
4477:   return 0;
4478: }

4480: /*@
4481:    TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`

4483:    Not Collective

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

4488:    Output Parameter:
4489: .  fails - number of failed nonlinear solves

4491:    Level: intermediate

4493:    Note:
4494:    This counter is reset to zero for each successive call to `TSSolve()`.

4496: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4497: @*/
4498: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4499: {
4502:   *fails = ts->num_snes_failures;
4503:   return 0;
4504: }

4506: /*@
4507:    TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails

4509:    Not Collective

4511:    Input Parameters:
4512: +  ts - `TS` context
4513: -  rejects - maximum number of rejected steps, pass -1 for unlimited

4515:    Options Database Key:
4516: .  -ts_max_reject - Maximum number of step rejections before a step fails

4518:    Level: intermediate

4520: .seealso: [](chapter_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4521: @*/
4522: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4523: {
4525:   ts->max_reject = rejects;
4526:   return 0;
4527: }

4529: /*@
4530:    TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves

4532:    Not Collective

4534:    Input Parameters:
4535: +  ts - `TS` context
4536: -  fails - maximum number of failed nonlinear solves, pass -1 for unlimited

4538:    Options Database Key:
4539: .  -ts_max_snes_failures - Maximum number of nonlinear solve failures

4541:    Level: intermediate

4543: .seealso: [](chapter_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4544: @*/
4545: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4546: {
4548:   ts->max_snes_failures = fails;
4549:   return 0;
4550: }

4552: /*@
4553:    TSSetErrorIfStepFails - Immediately error if no step succeeds

4555:    Not Collective

4557:    Input Parameters:
4558: +  ts - `TS` context
4559: -  err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure

4561:    Options Database Key:
4562: .  -ts_error_if_step_fails - Error if no step succeeds

4564:    Level: intermediate

4566: .seealso: [](chapter_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4567: @*/
4568: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
4569: {
4571:   ts->errorifstepfailed = err;
4572:   return 0;
4573: }

4575: /*@
4576:    TSGetAdapt - Get the adaptive controller context for the current method

4578:    Collective on ts if controller has not been created yet

4580:    Input Parameter:
4581: .  ts - time stepping context

4583:    Output Parameter:
4584: .  adapt - adaptive controller

4586:    Level: intermediate

4588: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
4589: @*/
4590: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
4591: {
4594:   if (!ts->adapt) {
4595:     TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt);
4596:     PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1);
4597:   }
4598:   *adapt = ts->adapt;
4599:   return 0;
4600: }

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

4605:    Logically Collective

4607:    Input Parameters:
4608: +  ts - time integration context
4609: .  atol - scalar absolute tolerances, `PETSC_DECIDE` to leave current value
4610: .  vatol - vector of absolute tolerances or NULL, used in preference to atol if present
4611: .  rtol - scalar relative tolerances, `PETSC_DECIDE` to leave current value
4612: -  vrtol - vector of relative tolerances or NULL, used in preference to atol if present

4614:    Options Database keys:
4615: +  -ts_rtol <rtol> - relative tolerance for local truncation error
4616: -  -ts_atol <atol> - Absolute tolerance for local truncation error

4618:    Level: beginner

4620:    Notes:
4621:    With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
4622:    (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
4623:    computed only for the differential or the algebraic part then this can be done using the vector of
4624:    tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
4625:    differential part and infinity for the algebraic part, the LTE calculation will include only the
4626:    differential variables.

4628: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
4629: @*/
4630: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
4631: {
4632:   if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol;
4633:   if (vatol) {
4634:     PetscObjectReference((PetscObject)vatol);
4635:     VecDestroy(&ts->vatol);
4636:     ts->vatol = vatol;
4637:   }
4638:   if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol;
4639:   if (vrtol) {
4640:     PetscObjectReference((PetscObject)vrtol);
4641:     VecDestroy(&ts->vrtol);
4642:     ts->vrtol = vrtol;
4643:   }
4644:   return 0;
4645: }

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

4650:    Logically Collective

4652:    Input Parameter:
4653: .  ts - time integration context

4655:    Output Parameters:
4656: +  atol - scalar absolute tolerances, NULL to ignore
4657: .  vatol - vector of absolute tolerances, NULL to ignore
4658: .  rtol - scalar relative tolerances, NULL to ignore
4659: -  vrtol - vector of relative tolerances, NULL to ignore

4661:    Level: beginner

4663: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
4664: @*/
4665: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
4666: {
4667:   if (atol) *atol = ts->atol;
4668:   if (vatol) *vatol = ts->vatol;
4669:   if (rtol) *rtol = ts->rtol;
4670:   if (vrtol) *vrtol = ts->vrtol;
4671:   return 0;
4672: }

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

4677:    Collective

4679:    Input Parameters:
4680: +  ts - time stepping context
4681: .  U - state vector, usually ts->vec_sol
4682: -  Y - state vector to be compared to U

4684:    Output Parameters:
4685: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
4686: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
4687: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

4689:    Level: developer

4691: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedNorm()`, `TSErrorWeightedNormInfinity()`
4692: @*/
4693: PetscErrorCode TSErrorWeightedNorm2(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4694: {
4695:   PetscInt           i, n, N, rstart;
4696:   PetscInt           n_loc, na_loc, nr_loc;
4697:   PetscReal          n_glb, na_glb, nr_glb;
4698:   const PetscScalar *u, *y;
4699:   PetscReal          sum, suma, sumr, gsum, gsuma, gsumr, diff;
4700:   PetscReal          tol, tola, tolr;
4701:   PetscReal          err_loc[6], err_glb[6];


4714:   VecGetSize(U, &N);
4715:   VecGetLocalSize(U, &n);
4716:   VecGetOwnershipRange(U, &rstart, NULL);
4717:   VecGetArrayRead(U, &u);
4718:   VecGetArrayRead(Y, &y);
4719:   sum    = 0.;
4720:   n_loc  = 0;
4721:   suma   = 0.;
4722:   na_loc = 0;
4723:   sumr   = 0.;
4724:   nr_loc = 0;
4725:   if (ts->vatol && ts->vrtol) {
4726:     const PetscScalar *atol, *rtol;
4727:     VecGetArrayRead(ts->vatol, &atol);
4728:     VecGetArrayRead(ts->vrtol, &rtol);
4729:     for (i = 0; i < n; i++) {
4730:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4731:       diff = PetscAbsScalar(y[i] - u[i]);
4732:       tola = PetscRealPart(atol[i]);
4733:       if (tola > 0.) {
4734:         suma += PetscSqr(diff / tola);
4735:         na_loc++;
4736:       }
4737:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4738:       if (tolr > 0.) {
4739:         sumr += PetscSqr(diff / tolr);
4740:         nr_loc++;
4741:       }
4742:       tol = tola + tolr;
4743:       if (tol > 0.) {
4744:         sum += PetscSqr(diff / tol);
4745:         n_loc++;
4746:       }
4747:     }
4748:     VecRestoreArrayRead(ts->vatol, &atol);
4749:     VecRestoreArrayRead(ts->vrtol, &rtol);
4750:   } else if (ts->vatol) { /* vector atol, scalar rtol */
4751:     const PetscScalar *atol;
4752:     VecGetArrayRead(ts->vatol, &atol);
4753:     for (i = 0; i < n; i++) {
4754:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4755:       diff = PetscAbsScalar(y[i] - u[i]);
4756:       tola = PetscRealPart(atol[i]);
4757:       if (tola > 0.) {
4758:         suma += PetscSqr(diff / tola);
4759:         na_loc++;
4760:       }
4761:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4762:       if (tolr > 0.) {
4763:         sumr += PetscSqr(diff / tolr);
4764:         nr_loc++;
4765:       }
4766:       tol = tola + tolr;
4767:       if (tol > 0.) {
4768:         sum += PetscSqr(diff / tol);
4769:         n_loc++;
4770:       }
4771:     }
4772:     VecRestoreArrayRead(ts->vatol, &atol);
4773:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
4774:     const PetscScalar *rtol;
4775:     VecGetArrayRead(ts->vrtol, &rtol);
4776:     for (i = 0; i < n; i++) {
4777:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4778:       diff = PetscAbsScalar(y[i] - u[i]);
4779:       tola = ts->atol;
4780:       if (tola > 0.) {
4781:         suma += PetscSqr(diff / tola);
4782:         na_loc++;
4783:       }
4784:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4785:       if (tolr > 0.) {
4786:         sumr += PetscSqr(diff / tolr);
4787:         nr_loc++;
4788:       }
4789:       tol = tola + tolr;
4790:       if (tol > 0.) {
4791:         sum += PetscSqr(diff / tol);
4792:         n_loc++;
4793:       }
4794:     }
4795:     VecRestoreArrayRead(ts->vrtol, &rtol);
4796:   } else { /* scalar atol, scalar rtol */
4797:     for (i = 0; i < n; i++) {
4798:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4799:       diff = PetscAbsScalar(y[i] - u[i]);
4800:       tola = ts->atol;
4801:       if (tola > 0.) {
4802:         suma += PetscSqr(diff / tola);
4803:         na_loc++;
4804:       }
4805:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4806:       if (tolr > 0.) {
4807:         sumr += PetscSqr(diff / tolr);
4808:         nr_loc++;
4809:       }
4810:       tol = tola + tolr;
4811:       if (tol > 0.) {
4812:         sum += PetscSqr(diff / tol);
4813:         n_loc++;
4814:       }
4815:     }
4816:   }
4817:   VecRestoreArrayRead(U, &u);
4818:   VecRestoreArrayRead(Y, &y);

4820:   err_loc[0] = sum;
4821:   err_loc[1] = suma;
4822:   err_loc[2] = sumr;
4823:   err_loc[3] = (PetscReal)n_loc;
4824:   err_loc[4] = (PetscReal)na_loc;
4825:   err_loc[5] = (PetscReal)nr_loc;

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

4829:   gsum   = err_glb[0];
4830:   gsuma  = err_glb[1];
4831:   gsumr  = err_glb[2];
4832:   n_glb  = err_glb[3];
4833:   na_glb = err_glb[4];
4834:   nr_glb = err_glb[5];

4836:   *norm = 0.;
4837:   if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
4838:   *norma = 0.;
4839:   if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
4840:   *normr = 0.;
4841:   if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);

4846:   return 0;
4847: }

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

4852:    Collective

4854:    Input Parameters:
4855: +  ts - time stepping context
4856: .  U - state vector, usually ts->vec_sol
4857: -  Y - state vector to be compared to U

4859:    Output Parameters:
4860: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
4861: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
4862: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

4864:    Level: developer

4866: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedNorm()`, `TSErrorWeightedNorm2()`
4867: @*/
4868: PetscErrorCode TSErrorWeightedNormInfinity(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4869: {
4870:   PetscInt           i, n, N, rstart;
4871:   const PetscScalar *u, *y;
4872:   PetscReal          max, gmax, maxa, gmaxa, maxr, gmaxr;
4873:   PetscReal          tol, tola, tolr, diff;
4874:   PetscReal          err_loc[3], err_glb[3];


4887:   VecGetSize(U, &N);
4888:   VecGetLocalSize(U, &n);
4889:   VecGetOwnershipRange(U, &rstart, NULL);
4890:   VecGetArrayRead(U, &u);
4891:   VecGetArrayRead(Y, &y);

4893:   max  = 0.;
4894:   maxa = 0.;
4895:   maxr = 0.;

4897:   if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
4898:     const PetscScalar *atol, *rtol;
4899:     VecGetArrayRead(ts->vatol, &atol);
4900:     VecGetArrayRead(ts->vrtol, &rtol);

4902:     for (i = 0; i < n; i++) {
4903:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4904:       diff = PetscAbsScalar(y[i] - u[i]);
4905:       tola = PetscRealPart(atol[i]);
4906:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4907:       tol  = tola + tolr;
4908:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4909:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4910:       if (tol > 0.) max = PetscMax(max, diff / tol);
4911:     }
4912:     VecRestoreArrayRead(ts->vatol, &atol);
4913:     VecRestoreArrayRead(ts->vrtol, &rtol);
4914:   } else if (ts->vatol) { /* vector atol, scalar rtol */
4915:     const PetscScalar *atol;
4916:     VecGetArrayRead(ts->vatol, &atol);
4917:     for (i = 0; i < n; i++) {
4918:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4919:       diff = PetscAbsScalar(y[i] - u[i]);
4920:       tola = PetscRealPart(atol[i]);
4921:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4922:       tol  = tola + tolr;
4923:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4924:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4925:       if (tol > 0.) max = PetscMax(max, diff / tol);
4926:     }
4927:     VecRestoreArrayRead(ts->vatol, &atol);
4928:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
4929:     const PetscScalar *rtol;
4930:     VecGetArrayRead(ts->vrtol, &rtol);

4932:     for (i = 0; i < n; i++) {
4933:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4934:       diff = PetscAbsScalar(y[i] - u[i]);
4935:       tola = ts->atol;
4936:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4937:       tol  = tola + tolr;
4938:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4939:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4940:       if (tol > 0.) max = PetscMax(max, diff / tol);
4941:     }
4942:     VecRestoreArrayRead(ts->vrtol, &rtol);
4943:   } else { /* scalar atol, scalar rtol */

4945:     for (i = 0; i < n; i++) {
4946:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4947:       diff = PetscAbsScalar(y[i] - u[i]);
4948:       tola = ts->atol;
4949:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4950:       tol  = tola + tolr;
4951:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4952:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4953:       if (tol > 0.) max = PetscMax(max, diff / tol);
4954:     }
4955:   }
4956:   VecRestoreArrayRead(U, &u);
4957:   VecRestoreArrayRead(Y, &y);
4958:   err_loc[0] = max;
4959:   err_loc[1] = maxa;
4960:   err_loc[2] = maxr;
4961:   MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
4962:   gmax  = err_glb[0];
4963:   gmaxa = err_glb[1];
4964:   gmaxr = err_glb[2];

4966:   *norm  = gmax;
4967:   *norma = gmaxa;
4968:   *normr = gmaxr;
4972:   return 0;
4973: }

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

4978:    Collective

4980:    Input Parameters:
4981: +  ts - time stepping context
4982: .  U - state vector, usually ts->vec_sol
4983: .  Y - state vector to be compared to U
4984: -  wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

4986:    Output Parameters:
4987: +  norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
4988: .  norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
4989: -  normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

4991:    Options Database Key:
4992: .  -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

4994:    Level: developer

4996: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`, `TSErrorWeightedENorm`
4997: @*/
4998: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4999: {
5000:   if (wnormtype == NORM_2) TSErrorWeightedNorm2(ts, U, Y, norm, norma, normr);
5001:   else if (wnormtype == NORM_INFINITY) TSErrorWeightedNormInfinity(ts, U, Y, norm, norma, normr);
5002:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5003:   return 0;
5004: }

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

5009:    Collective

5011:    Input Parameters:
5012: +  ts - time stepping context
5013: .  E - error vector
5014: .  U - state vector, usually ts->vec_sol
5015: -  Y - state vector, previous time step

5017:    Output Parameters:
5018: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5019: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5020: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5022:    Level: developer

5024: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedENorm()`, `TSErrorWeightedENormInfinity()`
5025: @*/
5026: PetscErrorCode TSErrorWeightedENorm2(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5027: {
5028:   PetscInt           i, n, N, rstart;
5029:   PetscInt           n_loc, na_loc, nr_loc;
5030:   PetscReal          n_glb, na_glb, nr_glb;
5031:   const PetscScalar *e, *u, *y;
5032:   PetscReal          err, sum, suma, sumr, gsum, gsuma, gsumr;
5033:   PetscReal          tol, tola, tolr;
5034:   PetscReal          err_loc[6], err_glb[6];


5049:   VecGetSize(E, &N);
5050:   VecGetLocalSize(E, &n);
5051:   VecGetOwnershipRange(E, &rstart, NULL);
5052:   VecGetArrayRead(E, &e);
5053:   VecGetArrayRead(U, &u);
5054:   VecGetArrayRead(Y, &y);
5055:   sum    = 0.;
5056:   n_loc  = 0;
5057:   suma   = 0.;
5058:   na_loc = 0;
5059:   sumr   = 0.;
5060:   nr_loc = 0;
5061:   if (ts->vatol && ts->vrtol) {
5062:     const PetscScalar *atol, *rtol;
5063:     VecGetArrayRead(ts->vatol, &atol);
5064:     VecGetArrayRead(ts->vrtol, &rtol);
5065:     for (i = 0; i < n; i++) {
5066:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5067:       err  = PetscAbsScalar(e[i]);
5068:       tola = PetscRealPart(atol[i]);
5069:       if (tola > 0.) {
5070:         suma += PetscSqr(err / tola);
5071:         na_loc++;
5072:       }
5073:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5074:       if (tolr > 0.) {
5075:         sumr += PetscSqr(err / tolr);
5076:         nr_loc++;
5077:       }
5078:       tol = tola + tolr;
5079:       if (tol > 0.) {
5080:         sum += PetscSqr(err / tol);
5081:         n_loc++;
5082:       }
5083:     }
5084:     VecRestoreArrayRead(ts->vatol, &atol);
5085:     VecRestoreArrayRead(ts->vrtol, &rtol);
5086:   } else if (ts->vatol) { /* vector atol, scalar rtol */
5087:     const PetscScalar *atol;
5088:     VecGetArrayRead(ts->vatol, &atol);
5089:     for (i = 0; i < n; i++) {
5090:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5091:       err  = PetscAbsScalar(e[i]);
5092:       tola = PetscRealPart(atol[i]);
5093:       if (tola > 0.) {
5094:         suma += PetscSqr(err / tola);
5095:         na_loc++;
5096:       }
5097:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5098:       if (tolr > 0.) {
5099:         sumr += PetscSqr(err / tolr);
5100:         nr_loc++;
5101:       }
5102:       tol = tola + tolr;
5103:       if (tol > 0.) {
5104:         sum += PetscSqr(err / tol);
5105:         n_loc++;
5106:       }
5107:     }
5108:     VecRestoreArrayRead(ts->vatol, &atol);
5109:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
5110:     const PetscScalar *rtol;
5111:     VecGetArrayRead(ts->vrtol, &rtol);
5112:     for (i = 0; i < n; i++) {
5113:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5114:       err  = PetscAbsScalar(e[i]);
5115:       tola = ts->atol;
5116:       if (tola > 0.) {
5117:         suma += PetscSqr(err / tola);
5118:         na_loc++;
5119:       }
5120:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5121:       if (tolr > 0.) {
5122:         sumr += PetscSqr(err / tolr);
5123:         nr_loc++;
5124:       }
5125:       tol = tola + tolr;
5126:       if (tol > 0.) {
5127:         sum += PetscSqr(err / tol);
5128:         n_loc++;
5129:       }
5130:     }
5131:     VecRestoreArrayRead(ts->vrtol, &rtol);
5132:   } else { /* scalar atol, scalar rtol */
5133:     for (i = 0; i < n; i++) {
5134:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5135:       err  = PetscAbsScalar(e[i]);
5136:       tola = ts->atol;
5137:       if (tola > 0.) {
5138:         suma += PetscSqr(err / tola);
5139:         na_loc++;
5140:       }
5141:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5142:       if (tolr > 0.) {
5143:         sumr += PetscSqr(err / tolr);
5144:         nr_loc++;
5145:       }
5146:       tol = tola + tolr;
5147:       if (tol > 0.) {
5148:         sum += PetscSqr(err / tol);
5149:         n_loc++;
5150:       }
5151:     }
5152:   }
5153:   VecRestoreArrayRead(E, &e);
5154:   VecRestoreArrayRead(U, &u);
5155:   VecRestoreArrayRead(Y, &y);

5157:   err_loc[0] = sum;
5158:   err_loc[1] = suma;
5159:   err_loc[2] = sumr;
5160:   err_loc[3] = (PetscReal)n_loc;
5161:   err_loc[4] = (PetscReal)na_loc;
5162:   err_loc[5] = (PetscReal)nr_loc;

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

5166:   gsum   = err_glb[0];
5167:   gsuma  = err_glb[1];
5168:   gsumr  = err_glb[2];
5169:   n_glb  = err_glb[3];
5170:   na_glb = err_glb[4];
5171:   nr_glb = err_glb[5];

5173:   *norm = 0.;
5174:   if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
5175:   *norma = 0.;
5176:   if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
5177:   *normr = 0.;
5178:   if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);

5183:   return 0;
5184: }

5186: /*@
5187:    TSErrorWeightedENormInfinity - compute a weighted infinity error norm based on supplied absolute and relative tolerances

5189:    Collective

5191:    Input Parameters:
5192: +  ts - time stepping context
5193: .  E - error vector
5194: .  U - state vector, usually ts->vec_sol
5195: -  Y - state vector, previous time step

5197:    Output Parameters:
5198: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5199: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5200: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5202:    Level: developer

5204: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedENorm()`, `TSErrorWeightedENorm2()`
5205: @*/
5206: PetscErrorCode TSErrorWeightedENormInfinity(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5207: {
5208:   PetscInt           i, n, N, rstart;
5209:   const PetscScalar *e, *u, *y;
5210:   PetscReal          err, max, gmax, maxa, gmaxa, maxr, gmaxr;
5211:   PetscReal          tol, tola, tolr;
5212:   PetscReal          err_loc[3], err_glb[3];


5227:   VecGetSize(E, &N);
5228:   VecGetLocalSize(E, &n);
5229:   VecGetOwnershipRange(E, &rstart, NULL);
5230:   VecGetArrayRead(E, &e);
5231:   VecGetArrayRead(U, &u);
5232:   VecGetArrayRead(Y, &y);

5234:   max  = 0.;
5235:   maxa = 0.;
5236:   maxr = 0.;

5238:   if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
5239:     const PetscScalar *atol, *rtol;
5240:     VecGetArrayRead(ts->vatol, &atol);
5241:     VecGetArrayRead(ts->vrtol, &rtol);

5243:     for (i = 0; i < n; i++) {
5244:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5245:       err  = PetscAbsScalar(e[i]);
5246:       tola = PetscRealPart(atol[i]);
5247:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5248:       tol  = tola + tolr;
5249:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5250:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5251:       if (tol > 0.) max = PetscMax(max, err / tol);
5252:     }
5253:     VecRestoreArrayRead(ts->vatol, &atol);
5254:     VecRestoreArrayRead(ts->vrtol, &rtol);
5255:   } else if (ts->vatol) { /* vector atol, scalar rtol */
5256:     const PetscScalar *atol;
5257:     VecGetArrayRead(ts->vatol, &atol);
5258:     for (i = 0; i < n; i++) {
5259:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5260:       err  = PetscAbsScalar(e[i]);
5261:       tola = PetscRealPart(atol[i]);
5262:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5263:       tol  = tola + tolr;
5264:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5265:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5266:       if (tol > 0.) max = PetscMax(max, err / tol);
5267:     }
5268:     VecRestoreArrayRead(ts->vatol, &atol);
5269:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
5270:     const PetscScalar *rtol;
5271:     VecGetArrayRead(ts->vrtol, &rtol);

5273:     for (i = 0; i < n; i++) {
5274:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5275:       err  = PetscAbsScalar(e[i]);
5276:       tola = ts->atol;
5277:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5278:       tol  = tola + tolr;
5279:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5280:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5281:       if (tol > 0.) max = PetscMax(max, err / tol);
5282:     }
5283:     VecRestoreArrayRead(ts->vrtol, &rtol);
5284:   } else { /* scalar atol, scalar rtol */

5286:     for (i = 0; i < n; i++) {
5287:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5288:       err  = PetscAbsScalar(e[i]);
5289:       tola = ts->atol;
5290:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5291:       tol  = tola + tolr;
5292:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5293:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5294:       if (tol > 0.) max = PetscMax(max, err / tol);
5295:     }
5296:   }
5297:   VecRestoreArrayRead(E, &e);
5298:   VecRestoreArrayRead(U, &u);
5299:   VecRestoreArrayRead(Y, &y);
5300:   err_loc[0] = max;
5301:   err_loc[1] = maxa;
5302:   err_loc[2] = maxr;
5303:   MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
5304:   gmax  = err_glb[0];
5305:   gmaxa = err_glb[1];
5306:   gmaxr = err_glb[2];

5308:   *norm  = gmax;
5309:   *norma = gmaxa;
5310:   *normr = gmaxr;
5314:   return 0;
5315: }

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

5320:    Collective

5322:    Input Parameters:
5323: +  ts - time stepping context
5324: .  E - error vector
5325: .  U - state vector, usually ts->vec_sol
5326: .  Y - state vector, previous time step
5327: -  wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

5329:    Output Parameters:
5330: +  norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5331: .  norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5332: -  normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

5334:    Options Database Key:
5335: .  -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5337:    Level: developer

5339: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedENormInfinity()`, `TSErrorWeightedENorm2()`, `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`
5340: @*/
5341: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5342: {
5343:   if (wnormtype == NORM_2) TSErrorWeightedENorm2(ts, E, U, Y, norm, norma, normr);
5344:   else if (wnormtype == NORM_INFINITY) TSErrorWeightedENormInfinity(ts, E, U, Y, norm, norma, normr);
5345:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5346:   return 0;
5347: }

5349: /*@
5350:    TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler

5352:    Logically Collective

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

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

5361:    Level: intermediate

5363: .seealso: [](chapter_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5364: @*/
5365: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5366: {
5368:   ts->cfltime_local = cfltime;
5369:   ts->cfltime       = -1.;
5370:   return 0;
5371: }

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

5376:    Collective

5378:    Input Parameter:
5379: .  ts - time stepping context

5381:    Output Parameter:
5382: .  cfltime - maximum stable time step for forward Euler

5384:    Level: advanced

5386: .seealso: [](chapter_ts), `TSSetCFLTimeLocal()`
5387: @*/
5388: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5389: {
5390:   if (ts->cfltime < 0) MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts));
5391:   *cfltime = ts->cfltime;
5392:   return 0;
5393: }

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

5398:    Input Parameters:
5399: +  ts   - the `TS` context.
5400: .  xl   - lower bound.
5401: -  xu   - upper bound.

5403:    Level: advanced

5405:    Note:
5406:    If this routine is not called then the lower and upper bounds are set to
5407:    `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.

5409: .seealso: [](chapter_ts), `TS`
5410: @*/
5411: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5412: {
5413:   SNES snes;

5415:   TSGetSNES(ts, &snes);
5416:   SNESVISetVariableBounds(snes, xl, xu);
5417:   return 0;
5418: }

5420: /*@
5421:    TSComputeLinearStability - computes the linear stability function at a point

5423:    Collective

5425:    Input Parameters:
5426: +  ts - the `TS` context
5427: -  xr,xi - real and imaginary part of input arguments

5429:    Output Parameters:
5430: .  yr,yi - real and imaginary part of function value

5432:    Level: developer

5434: .seealso: [](chapter_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5435: @*/
5436: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5437: {
5439:   PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5440:   return 0;
5441: }

5443: /*@
5444:    TSRestartStep - Flags the solver to restart the next step

5446:    Collective

5448:    Input Parameter:
5449: .  ts - the `TS` context obtained from `TSCreate()`

5451:    Level: advanced

5453:    Notes:
5454:    Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5455:    discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5456:    vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5457:    the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5458:    discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5459:    discontinuous source terms).

5461: .seealso: [](chapter_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5462: @*/
5463: PetscErrorCode TSRestartStep(TS ts)
5464: {
5466:   ts->steprestart = PETSC_TRUE;
5467:   return 0;
5468: }

5470: /*@
5471:    TSRollBack - Rolls back one time step

5473:    Collective

5475:    Input Parameter:
5476: .  ts - the `TS` context obtained from `TSCreate()`

5478:    Level: advanced

5480: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5481: @*/
5482: PetscErrorCode TSRollBack(TS ts)
5483: {
5486:   PetscUseTypeMethod(ts, rollback);
5487:   ts->time_step  = ts->ptime - ts->ptime_prev;
5488:   ts->ptime      = ts->ptime_prev;
5489:   ts->ptime_prev = ts->ptime_prev_rollback;
5490:   ts->steps--;
5491:   ts->steprollback = PETSC_TRUE;
5492:   return 0;
5493: }

5495: /*@
5496:    TSGetStages - Get the number of stages and stage values

5498:    Input Parameter:
5499: .  ts - the `TS` context obtained from `TSCreate()`

5501:    Output Parameters:
5502: +  ns - the number of stages
5503: -  Y - the current stage vectors

5505:    Level: advanced

5507:    Note:
5508:    Both ns and Y can be NULL.

5510: .seealso: [](chapter_ts), `TS`, `TSCreate()`
5511: @*/
5512: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5513: {
5517:   if (!ts->ops->getstages) {
5518:     if (ns) *ns = 0;
5519:     if (Y) *Y = NULL;
5520:   } else PetscUseTypeMethod(ts, getstages, ns, Y);
5521:   return 0;
5522: }

5524: /*@C
5525:   TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.

5527:   Collective

5529:   Input Parameters:
5530: + ts - the `TS` context
5531: . t - current timestep
5532: . U - state vector
5533: . Udot - time derivative of state vector
5534: . shift - shift to apply, see note below
5535: - ctx - an optional user context

5537:   Output Parameters:
5538: + J - Jacobian matrix (not altered in this routine)
5539: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)

5541:   Level: intermediate

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

5546:   dF/dU + shift*dF/dUdot

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

5551:   This will first try to get the coloring from the `DM`.  If the `DM` type has no coloring
5552:   routine, then it will try to get the coloring from the matrix.  This requires that the
5553:   matrix have nonzero entries precomputed.

5555: .seealso: [](chapter_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5556: @*/
5557: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5558: {
5559:   SNES          snes;
5560:   MatFDColoring color;
5561:   PetscBool     hascolor, matcolor = PETSC_FALSE;

5563:   PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL);
5564:   PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color);
5565:   if (!color) {
5566:     DM         dm;
5567:     ISColoring iscoloring;

5569:     TSGetDM(ts, &dm);
5570:     DMHasColoring(dm, &hascolor);
5571:     if (hascolor && !matcolor) {
5572:       DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring);
5573:       MatFDColoringCreate(B, iscoloring, &color);
5574:       MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts);
5575:       MatFDColoringSetFromOptions(color);
5576:       MatFDColoringSetUp(B, iscoloring, color);
5577:       ISColoringDestroy(&iscoloring);
5578:     } else {
5579:       MatColoring mc;

5581:       MatColoringCreate(B, &mc);
5582:       MatColoringSetDistance(mc, 2);
5583:       MatColoringSetType(mc, MATCOLORINGSL);
5584:       MatColoringSetFromOptions(mc);
5585:       MatColoringApply(mc, &iscoloring);
5586:       MatColoringDestroy(&mc);
5587:       MatFDColoringCreate(B, iscoloring, &color);
5588:       MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts);
5589:       MatFDColoringSetFromOptions(color);
5590:       MatFDColoringSetUp(B, iscoloring, color);
5591:       ISColoringDestroy(&iscoloring);
5592:     }
5593:     PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color);
5594:     PetscObjectDereference((PetscObject)color);
5595:   }
5596:   TSGetSNES(ts, &snes);
5597:   MatFDColoringApply(B, color, U, snes);
5598:   if (J != B) {
5599:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
5600:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
5601:   }
5602:   return 0;
5603: }

5605: /*@
5606:     TSSetFunctionDomainError - Set a function that tests if the current state vector is valid

5608:     Input Parameters:
5609: +    ts - the `TS` context
5610: -    func - function called within `TSFunctionDomainError()`

5612:     Calling sequence of func:
5613: $     PetscErrorCode func(TS ts,PetscReal time,Vec state,PetscBool reject)

5615: +   ts - the TS context
5616: .   time - the current time (of the stage)
5617: .   state - the state to check if it is valid
5618: -   reject - (output parameter) PETSC_FALSE if the state is acceptable, PETSC_TRUE if not acceptable

5620:     Level: intermediate

5622:     Notes:
5623:       If an implicit ODE solver is being used then, in addition to providing this routine, the
5624:       user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5625:       function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5626:       Use `TSGetSNES()` to obtain the `SNES` object

5628:     Developer Note:
5629:       The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5630:       since one takes a function pointer and the other does not.

5632: .seealso: [](chapter_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5633: @*/

5635: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS, PetscReal, Vec, PetscBool *))
5636: {
5638:   ts->functiondomainerror = func;
5639:   return 0;
5640: }

5642: /*@
5643:     TSFunctionDomainError - Checks if the current state is valid

5645:     Input Parameters:
5646: +    ts - the `TS` context
5647: .    stagetime - time of the simulation
5648: -    Y - state vector to check.

5650:     Output Parameter:
5651: .    accept - Set to `PETSC_FALSE` if the current state vector is valid.

5653:     Level: developer

5655:     Note:
5656:     This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5657:     to check if the current state is valid.

5659: .seealso: [](chapter_ts), `TS`, `TSSetFunctionDomainError()`
5660: @*/
5661: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5662: {
5664:   *accept = PETSC_TRUE;
5665:   if (ts->functiondomainerror) (*ts->functiondomainerror)(ts, stagetime, Y, accept);
5666:   return 0;
5667: }

5669: /*@C
5670:   TSClone - This function clones a time step object.

5672:   Collective

5674:   Input Parameter:
5675: . tsin    - The input `TS`

5677:   Output Parameter:
5678: . tsout   - The output `TS` (cloned)

5680:   Level: developer

5682:   Notes:
5683:   This function is used to create a clone of a `TS` object. It is used in `TSARKIMEX` for initializing the slope for first stage explicit methods.
5684:   It will likely be replaced in the future with a mechanism of switching methods on the fly.

5686:   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);

5688: .seealso: [](chapter_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5689: @*/
5690: PetscErrorCode TSClone(TS tsin, TS *tsout)
5691: {
5692:   TS     t;
5693:   SNES   snes_start;
5694:   DM     dm;
5695:   TSType type;

5698:   *tsout = NULL;

5700:   PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView);

5702:   /* General TS description */
5703:   t->numbermonitors    = 0;
5704:   t->monitorFrequency  = 1;
5705:   t->setupcalled       = 0;
5706:   t->ksp_its           = 0;
5707:   t->snes_its          = 0;
5708:   t->nwork             = 0;
5709:   t->rhsjacobian.time  = PETSC_MIN_REAL;
5710:   t->rhsjacobian.scale = 1.;
5711:   t->ijacobian.shift   = 1.;

5713:   TSGetSNES(tsin, &snes_start);
5714:   TSSetSNES(t, snes_start);

5716:   TSGetDM(tsin, &dm);
5717:   TSSetDM(t, dm);

5719:   t->adapt = tsin->adapt;
5720:   PetscObjectReference((PetscObject)t->adapt);

5722:   t->trajectory = tsin->trajectory;
5723:   PetscObjectReference((PetscObject)t->trajectory);

5725:   t->event = tsin->event;
5726:   if (t->event) t->event->refct++;

5728:   t->problem_type      = tsin->problem_type;
5729:   t->ptime             = tsin->ptime;
5730:   t->ptime_prev        = tsin->ptime_prev;
5731:   t->time_step         = tsin->time_step;
5732:   t->max_time          = tsin->max_time;
5733:   t->steps             = tsin->steps;
5734:   t->max_steps         = tsin->max_steps;
5735:   t->equation_type     = tsin->equation_type;
5736:   t->atol              = tsin->atol;
5737:   t->rtol              = tsin->rtol;
5738:   t->max_snes_failures = tsin->max_snes_failures;
5739:   t->max_reject        = tsin->max_reject;
5740:   t->errorifstepfailed = tsin->errorifstepfailed;

5742:   TSGetType(tsin, &type);
5743:   TSSetType(t, type);

5745:   t->vec_sol = NULL;

5747:   t->cfltime          = tsin->cfltime;
5748:   t->cfltime_local    = tsin->cfltime_local;
5749:   t->exact_final_time = tsin->exact_final_time;

5751:   PetscMemcpy(t->ops, tsin->ops, sizeof(struct _TSOps));

5753:   if (((PetscObject)tsin)->fortran_func_pointers) {
5754:     PetscInt i;
5755:     PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers);
5756:     for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5757:   }
5758:   *tsout = t;
5759:   return 0;
5760: }

5762: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5763: {
5764:   TS ts = (TS)ctx;

5766:   TSComputeRHSFunction(ts, 0, x, y);
5767:   return 0;
5768: }

5770: /*@
5771:     TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.

5773:    Logically Collective

5775:     Input Parameters:
5776:     TS - the time stepping routine

5778:    Output Parameter:
5779: .   flg - `PETSC_TRUE` if the multiply is likely correct

5781:    Options Database Key:
5782:  .   -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator

5784:    Level: advanced

5786:    Note:
5787:     This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian

5789: .seealso: [](chapter_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5790: @*/
5791: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5792: {
5793:   Mat           J, B;
5794:   TSRHSJacobian func;
5795:   void         *ctx;

5797:   TSGetRHSJacobian(ts, &J, &B, &func, &ctx);
5798:   (*func)(ts, 0.0, ts->vec_sol, J, B, ctx);
5799:   MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg);
5800:   return 0;
5801: }

5803: /*@C
5804:     TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.

5806:    Logically Collective

5808:     Input Parameters:
5809:     TS - the time stepping routine

5811:    Output Parameter:
5812: .   flg - `PETSC_TRUE` if the multiply is likely correct

5814:    Options Database Key:
5815: .   -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator

5817:    Level: advanced

5819:    Notes:
5820:     This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian

5822: .seealso: [](chapter_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5823: @*/
5824: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5825: {
5826:   Mat           J, B;
5827:   void         *ctx;
5828:   TSRHSJacobian func;

5830:   TSGetRHSJacobian(ts, &J, &B, &func, &ctx);
5831:   (*func)(ts, 0.0, ts->vec_sol, J, B, ctx);
5832:   MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg);
5833:   return 0;
5834: }

5836: /*@
5837:   TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.

5839:   Logically collective

5841:   Input Parameters:
5842: +  ts - timestepping context
5843: -  use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

5845:   Options Database Key:
5846: .   -ts_use_splitrhsfunction - <true,false>

5848:   Level: intermediate

5850:   Note:
5851:     This is only useful for multirate methods

5853: .seealso: [](chapter_ts), `TS`, `TSGetUseSplitRHSFunction()`
5854: @*/
5855: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5856: {
5858:   ts->use_splitrhsfunction = use_splitrhsfunction;
5859:   return 0;
5860: }

5862: /*@
5863:   TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.

5865:   Not collective

5867:   Input Parameter:
5868: .  ts - timestepping context

5870:   Output Parameter:
5871: .  use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

5873:   Level: intermediate

5875: .seealso: [](chapter_ts), `TS`, `TSSetUseSplitRHSFunction()`
5876: @*/
5877: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5878: {
5880:   *use_splitrhsfunction = ts->use_splitrhsfunction;
5881:   return 0;
5882: }

5884: /*@
5885:     TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.

5887:    Logically  Collective

5889:    Input Parameters:
5890: +  ts - the time-stepper
5891: -  str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)

5893:    Level: intermediate

5895:    Note:
5896:      When the relationship between the nonzero structures is known and supplied the solution process can be much faster

5898: .seealso: [](chapter_ts), `TS`, `MatAXPY()`, `MatStructure`
5899:  @*/
5900: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5901: {
5903:   ts->axpy_pattern = str;
5904:   return 0;
5905: }

5907: /*@
5908:   TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span

5910:   Collective

5912:   Input Parameters:
5913: + ts - the time-stepper
5914: . n - number of the time points (>=2)
5915: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.

5917:   Options Database Key:
5918: . -ts_time_span <t0,...tf> - Sets the time span

5920:   Level: intermediate

5922:   Notes:
5923:   The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
5924:   `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
5925:   The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may
5926:   pressure the memory system when using a large number of span points.

5928: .seealso: [](chapter_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()`
5929:  @*/
5930: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
5931: {
5934:   if (ts->tspan && n != ts->tspan->num_span_times) {
5935:     PetscFree(ts->tspan->span_times);
5936:     VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol);
5937:     PetscMalloc1(n, &ts->tspan->span_times);
5938:   }
5939:   if (!ts->tspan) {
5940:     TSTimeSpan tspan;
5941:     PetscNew(&tspan);
5942:     PetscMalloc1(n, &tspan->span_times);
5943:     tspan->reltol = 1e-6;
5944:     tspan->abstol = 10 * PETSC_MACHINE_EPSILON;
5945:     ts->tspan     = tspan;
5946:   }
5947:   ts->tspan->num_span_times = n;
5948:   PetscArraycpy(ts->tspan->span_times, span_times, n);
5949:   TSSetTime(ts, ts->tspan->span_times[0]);
5950:   TSSetMaxTime(ts, ts->tspan->span_times[n - 1]);
5951:   return 0;
5952: }

5954: /*@C
5955:   TSGetTimeSpan - gets the time span.

5957:   Not Collective

5959:   Input Parameter:
5960: . ts - the time-stepper

5962:   Output Parameters:
5963: + n - number of the time points (>=2)
5964: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5965:   The values are valid until the `TS` object is destroyed.

5967:   Level: beginner

5969:   Note:
5970:   Both n and span_times can be NULL.

5972: .seealso: [](chapter_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()`
5973:  @*/
5974: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal **span_times)
5975: {
5979:   if (!ts->tspan) {
5980:     if (n) *n = 0;
5981:     if (span_times) *span_times = NULL;
5982:   } else {
5983:     if (n) *n = ts->tspan->num_span_times;
5984:     if (span_times) *span_times = ts->tspan->span_times;
5985:   }
5986:   return 0;
5987: }

5989: /*@
5990:    TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.

5992:    Input Parameter:
5993: .  ts - the `TS` context obtained from `TSCreate()`

5995:    Output Parameters:
5996: +  nsol - the number of solutions
5997: -  Sols - the solution vectors

5999:    Level: intermediate

6001:    Notes:
6002:     Both nsol and Sols can be NULL.

6004:     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()`.
6005:     For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span.

6007: .seealso: [](chapter_ts), `TS`, `TSSetTimeSpan()`
6008: @*/
6009: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
6010: {
6014:   if (!ts->tspan) {
6015:     if (nsol) *nsol = 0;
6016:     if (Sols) *Sols = NULL;
6017:   } else {
6018:     if (nsol) *nsol = ts->tspan->spanctr;
6019:     if (Sols) *Sols = ts->tspan->vecs_sol;
6020:   }
6021:   return 0;
6022: }