Actual source code: ts.c

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

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

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

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

 19: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt, TSAdaptType default_type)
 20: {
 21:   PetscFunctionBegin;
 24:   if (!((PetscObject)adapt)->type_name) PetscCall(TSAdaptSetType(adapt, default_type));
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

 28: /*@
 29:    TSSetFromOptions - Sets various `TS` parameters from the options database

 31:    Collective

 33:    Input Parameter:
 34: .  ts - the `TS` context obtained from `TSCreate()`

 36:    Options Database Keys:
 37: +  -ts_type <type> - EULER, BEULER, SUNDIALS, PSEUDO, CN, RK, THETA, ALPHA, GLLE,  SSP, GLEE, BSYMP, IRK, see `TSType`
 38: .  -ts_save_trajectory - checkpoint the solution at each time-step
 39: .  -ts_max_time <time> - maximum time to compute to
 40: .  -ts_time_span <t0,...tf> - sets the time span, solutions are computed and stored for each indicated time
 41: .  -ts_max_steps <steps> - maximum number of time-steps to take
 42: .  -ts_init_time <time> - initial time to start computation
 43: .  -ts_final_time <time> - final time to compute to (deprecated: use `-ts_max_time`)
 44: .  -ts_dt <dt> - initial time step
 45: .  -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
 46: .  -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
 47: .  -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
 48: .  -ts_error_if_step_fails <true,false> - Error if no step succeeds
 49: .  -ts_rtol <rtol> - relative tolerance for local truncation error
 50: .  -ts_atol <atol> - Absolute tolerance for local truncation error
 51: .  -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function
 52: .  -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - test the Jacobian at each iteration against finite difference with RHS function
 53: .  -ts_adjoint_solve <yes,no> - After solving the ODE/DAE solve the adjoint problem (requires `-ts_save_trajectory`)
 54: .  -ts_fd_color - Use finite differences with coloring to compute IJacobian
 55: .  -ts_monitor - print information at each timestep
 56: .  -ts_monitor_cancel - Cancel all monitors
 57: .  -ts_monitor_lg_solution - Monitor solution graphically
 58: .  -ts_monitor_lg_error - Monitor error graphically
 59: .  -ts_monitor_error - Monitors norm of error
 60: .  -ts_monitor_lg_timestep - Monitor timestep size graphically
 61: .  -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
 62: .  -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
 63: .  -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
 64: .  -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
 65: .  -ts_monitor_draw_solution - Monitor solution graphically
 66: .  -ts_monitor_draw_solution_phase  <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
 67: .  -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
 68: .  -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
 69:    -ts_monitor_solution_interval <interval> - output once every interval (default=1) time steps; used with -ts_monitor_solution
 70: .  -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)
 71: -  -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time

 73:    Level: beginner

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

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

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

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

 99:   PetscFunctionBegin;

102:   PetscCall(TSRegisterAll());
103:   PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));

105:   PetscObjectOptionsBegin((PetscObject)ts);
106:   if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
107:   else defaultType = ifun ? TSBEULER : TSEULER;
108:   PetscCall(PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, 256, &opt));
109:   if (opt) PetscCall(TSSetType(ts, typeName));
110:   else PetscCall(TSSetType(ts, defaultType));

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

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

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

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

151:   PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt));
152:   if (opt) {
153:     PetscInt  howoften = 1;
154:     DM        dm;
155:     PetscBool net;

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

172:   PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt));
173:   if (opt) {
174:     TSMonitorLGCtx ctx;
175:     PetscInt       howoften = 1;

177:     PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL));
178:     PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
179:     PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
180:   }
181:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL));

183:   PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt));
184:   if (opt) {
185:     TSMonitorLGCtx ctx;
186:     PetscInt       howoften = 1;

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

197:     PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
198:     PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
199:     PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
200:     ctx->semilogy = PETSC_TRUE;
201:   }

203:   PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt));
204:   if (opt) {
205:     TSMonitorLGCtx ctx;
206:     PetscInt       howoften = 1;

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

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

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

236:     for (PetscInt i = 0; i < ts->numbermonitors; ++i)
237:       if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
238:         create = PETSC_FALSE;
239:         break;
240:       }
241:     if (create) {
242:       PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase space from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL));
243:       PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL));
244:       PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL));
245:       PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_multi_species", "Color particles by particle species", "TSMonitorSPSwarm", multispecies, &multispecies, NULL));
246:       PetscCall(TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, multispecies, &ctx));
247:       PetscCall(TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorSPCtxDestroy));
248:     }
249:   }
250:   PetscCall(PetscOptionsName("-ts_monitor_hg_swarm", "Display particle histogram from the DMSwarm", "TSMonitorHGSwarm", &opt));
251:   if (opt) {
252:     TSMonitorHGCtx ctx;
253:     PetscInt       howoften = 1, Ns = 1;
254:     PetscBool      velocity = PETSC_FALSE, create = PETSC_TRUE;

256:     for (PetscInt i = 0; i < ts->numbermonitors; ++i)
257:       if (ts->monitor[i] == TSMonitorHGSwarmSolution) {
258:         create = PETSC_FALSE;
259:         break;
260:       }
261:     if (create) {
262:       DM       sw, dm;
263:       PetscInt Nc, Nb;

265:       PetscCall(TSGetDM(ts, &sw));
266:       PetscCall(DMSwarmGetCellDM(sw, &dm));
267:       PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc));
268:       Nb = PetscMin(20, PetscMax(10, Nc));
269:       PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL));
270:       PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL));
271:       PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL));
272:       PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL));
273:       PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx));
274:       PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorHGCtxDestroy));
275:     }
276:   }
277:   opt = PETSC_FALSE;
278:   PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt));
279:   if (opt) {
280:     TSMonitorDrawCtx ctx;
281:     PetscInt         howoften = 1;

283:     PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL));
284:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
285:     PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
286:   }
287:   opt = PETSC_FALSE;
288:   PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt));
289:   if (opt) {
290:     TSMonitorDrawCtx ctx;
291:     PetscReal        bounds[4];
292:     PetscInt         n = 4;
293:     PetscDraw        draw;
294:     PetscDrawAxis    axis;

296:     PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL));
297:     PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field");
298:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx));
299:     PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw));
300:     PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis));
301:     PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]));
302:     PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2"));
303:     PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
304:   }
305:   opt = PETSC_FALSE;
306:   PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt));
307:   if (opt) {
308:     TSMonitorDrawCtx ctx;
309:     PetscInt         howoften = 1;

311:     PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL));
312:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
313:     PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
314:   }
315:   opt = PETSC_FALSE;
316:   PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt));
317:   if (opt) {
318:     TSMonitorDrawCtx ctx;
319:     PetscInt         howoften = 1;

321:     PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL));
322:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
323:     PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
324:   }

326:   opt = PETSC_FALSE;
327:   PetscCall(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));
328:   if (flg) {
329:     const char *ptr = NULL, *ptr2 = NULL;
330:     char       *filetemplate;
331:     PetscCheck(monfilename[0], PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03" PetscInt_FMT ".vts");
332:     /* Do some cursory validation of the input. */
333:     PetscCall(PetscStrstr(monfilename, "%", (char **)&ptr));
334:     PetscCheck(ptr, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03" PetscInt_FMT ".vts");
335:     for (ptr++; ptr && *ptr; ptr++) {
336:       PetscCall(PetscStrchr("DdiouxX", *ptr, (char **)&ptr2));
337:       PetscCheck(ptr2 || (*ptr >= '0' && *ptr <= '9'), PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Invalid file template argument to -ts_monitor_solution_vtk, should look like filename-%%03" PetscInt_FMT ".vts");
338:       if (ptr2) break;
339:     }
340:     PetscCall(PetscStrallocpy(monfilename, &filetemplate));
341:     PetscCall(TSMonitorSet(ts, TSMonitorSolutionVTK, filetemplate, (PetscErrorCode(*)(void **))TSMonitorSolutionVTKDestroy));
342:   }

344:   PetscCall(PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg));
345:   if (flg) {
346:     TSMonitorDMDARayCtx *rayctx;
347:     int                  ray = 0;
348:     DMDirection          ddir;
349:     DM                   da;
350:     PetscMPIInt          rank;

352:     PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
353:     if (dir[0] == 'x') ddir = DM_X;
354:     else if (dir[0] == 'y') ddir = DM_Y;
355:     else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
356:     sscanf(dir + 2, "%d", &ray);

358:     PetscCall(PetscInfo(((PetscObject)ts), "Displaying DMDA ray %c = %d\n", dir[0], ray));
359:     PetscCall(PetscNew(&rayctx));
360:     PetscCall(TSGetDM(ts, &da));
361:     PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
362:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank));
363:     if (rank == 0) PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer));
364:     rayctx->lgctx = NULL;
365:     PetscCall(TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy));
366:   }
367:   PetscCall(PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg));
368:   if (flg) {
369:     TSMonitorDMDARayCtx *rayctx;
370:     int                  ray = 0;
371:     DMDirection          ddir;
372:     DM                   da;
373:     PetscInt             howoften = 1;

375:     PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
376:     if (dir[0] == 'x') ddir = DM_X;
377:     else if (dir[0] == 'y') ddir = DM_Y;
378:     else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
379:     sscanf(dir + 2, "%d", &ray);

381:     PetscCall(PetscInfo(((PetscObject)ts), "Displaying LG DMDA ray %c = %d\n", dir[0], ray));
382:     PetscCall(PetscNew(&rayctx));
383:     PetscCall(TSGetDM(ts, &da));
384:     PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
385:     PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx));
386:     PetscCall(TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy));
387:   }

389:   PetscCall(PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt));
390:   if (opt) {
391:     TSMonitorEnvelopeCtx ctx;

393:     PetscCall(TSMonitorEnvelopeCtxCreate(ts, &ctx));
394:     PetscCall(TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscErrorCode(*)(void **))TSMonitorEnvelopeCtxDestroy));
395:   }
396:   flg = PETSC_FALSE;
397:   PetscCall(PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt));
398:   if (opt && flg) PetscCall(TSMonitorCancel(ts));

400:   flg = PETSC_FALSE;
401:   PetscCall(PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeIJacobianDefaultColor", flg, &flg, NULL));
402:   if (flg) {
403:     DM dm;

405:     PetscCall(TSGetDM(ts, &dm));
406:     PetscCall(DMTSUnsetIJacobianContext_Internal(dm));
407:     PetscCall(TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL));
408:     PetscCall(PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n"));
409:   }

411:   /* Handle specific TS options */
412:   PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);

414:   /* Handle TSAdapt options */
415:   PetscCall(TSGetAdapt(ts, &ts->adapt));
416:   PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
417:   PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject));

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

424:   PetscCall(TSAdjointSetFromOptions(ts, PetscOptionsObject));

426:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
427:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject));
428:   PetscOptionsEnd();

430:   if (ts->trajectory) PetscCall(TSTrajectorySetFromOptions(ts->trajectory, ts));

432:   /* why do we have to do this here and not during TSSetUp? */
433:   PetscCall(TSGetSNES(ts, &ts->snes));
434:   if (ts->problem_type == TS_LINEAR) {
435:     PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
436:     if (!flg) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
437:   }
438:   PetscCall(SNESSetFromOptions(ts->snes));
439:   PetscFunctionReturn(PETSC_SUCCESS);
440: }

442: /*@
443:    TSGetTrajectory - Gets the trajectory from a `TS` if it exists

445:    Collective

447:    Input Parameter:
448: .  ts - the `TS` context obtained from `TSCreate()`

450:    Output Parameter:
451: .  tr - the `TSTrajectory` object, if it exists

453:    Level: advanced

455:    Note:
456:    This routine should be called after all `TS` options have been set

458: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectory`, `TSTrajectoryCreate()`
459: @*/
460: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
461: {
462:   PetscFunctionBegin;
464:   *tr = ts->trajectory;
465:   PetscFunctionReturn(PETSC_SUCCESS);
466: }

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

471:    Collective

473:    Input Parameter:
474: .  ts - the `TS` context obtained from `TSCreate()`

476:    Options Database Keys:
477: +  -ts_save_trajectory - saves the trajectory to a file
478: -  -ts_trajectory_type type - set trajectory type

480:    Level: intermediate

482:    Notes:
483:    This routine should be called after all `TS` options have been set

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

488: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`
489: @*/
490: PetscErrorCode TSSetSaveTrajectory(TS ts)
491: {
492:   PetscFunctionBegin;
494:   if (!ts->trajectory) PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
495:   PetscFunctionReturn(PETSC_SUCCESS);
496: }

498: /*@
499:    TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object

501:    Collective

503:    Input Parameter:
504: .  ts - the `TS` context obtained from `TSCreate()`

506:    Level: intermediate

508: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
509: @*/
510: PetscErrorCode TSResetTrajectory(TS ts)
511: {
512:   PetscFunctionBegin;
514:   if (ts->trajectory) {
515:     PetscCall(TSTrajectoryDestroy(&ts->trajectory));
516:     PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
517:   }
518:   PetscFunctionReturn(PETSC_SUCCESS);
519: }

521: /*@
522:    TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from a `TS`

524:    Collective

526:    Input Parameter:
527: .  ts - the `TS` context obtained from `TSCreate()`

529:    Level: intermediate

531: .seealso: [](ch_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
532: @*/
533: PetscErrorCode TSRemoveTrajectory(TS ts)
534: {
535:   PetscFunctionBegin;
537:   if (ts->trajectory) PetscCall(TSTrajectoryDestroy(&ts->trajectory));
538:   PetscFunctionReturn(PETSC_SUCCESS);
539: }

541: /*@
542:    TSComputeRHSJacobian - Computes the Jacobian matrix that has been
543:       set with `TSSetRHSJacobian()`.

545:    Collective

547:    Input Parameters:
548: +  ts - the `TS` context
549: .  t - current timestep
550: -  U - input vector

552:    Output Parameters:
553: +  A - Jacobian matrix
554: -  B - optional preconditioning matrix

556:    Level: developer

558:    Note:
559:    Most users should not need to explicitly call this routine, as it
560:    is used internally within the nonlinear solvers.

562: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()`
563: @*/
564: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
565: {
566:   PetscObjectState Ustate;
567:   PetscObjectId    Uid;
568:   DM               dm;
569:   DMTS             tsdm;
570:   TSRHSJacobian    rhsjacobianfunc;
571:   void            *ctx;
572:   TSRHSFunction    rhsfunction;

574:   PetscFunctionBegin;
577:   PetscCheckSameComm(ts, 1, U, 3);
578:   PetscCall(TSGetDM(ts, &dm));
579:   PetscCall(DMGetDMTS(dm, &tsdm));
580:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
581:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx));
582:   PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
583:   PetscCall(PetscObjectGetId((PetscObject)U, &Uid));

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

587:   PetscCheck(ts->rhsjacobian.shift == 0.0 || !ts->rhsjacobian.reuse, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Should not call TSComputeRHSJacobian() on a shifted matrix (shift=%lf) when RHSJacobian is reusable.", (double)ts->rhsjacobian.shift);
588:   if (rhsjacobianfunc) {
589:     PetscCall(PetscLogEventBegin(TS_JacobianEval, ts, U, A, B));
590:     PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
591:     ts->rhsjacs++;
592:     PetscCall(PetscLogEventEnd(TS_JacobianEval, ts, U, A, B));
593:   } else {
594:     PetscCall(MatZeroEntries(A));
595:     if (B && A != B) PetscCall(MatZeroEntries(B));
596:   }
597:   ts->rhsjacobian.time  = t;
598:   ts->rhsjacobian.shift = 0;
599:   ts->rhsjacobian.scale = 1.;
600:   PetscCall(PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid));
601:   PetscCall(PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate));
602:   PetscFunctionReturn(PETSC_SUCCESS);
603: }

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

608:    Collective

610:    Input Parameters:
611: +  ts - the `TS` context
612: .  t - current time
613: -  U - state vector

615:    Output Parameter:
616: .  y - right hand side

618:    Level: developer

620:    Note:
621:    Most users should not need to explicitly call this routine, as it
622:    is used internally within the nonlinear solvers.

624: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
625: @*/
626: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
627: {
628:   TSRHSFunction rhsfunction;
629:   TSIFunction   ifunction;
630:   void         *ctx;
631:   DM            dm;

633:   PetscFunctionBegin;
637:   PetscCall(TSGetDM(ts, &dm));
638:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx));
639:   PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));

641:   PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");

643:   if (rhsfunction) {
644:     PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, y, 0));
645:     PetscCall(VecLockReadPush(U));
646:     PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
647:     PetscCall(VecLockReadPop(U));
648:     ts->rhsfuncs++;
649:     PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, y, 0));
650:   } else PetscCall(VecZeroEntries(y));
651:   PetscFunctionReturn(PETSC_SUCCESS);
652: }

654: /*@
655:    TSComputeSolutionFunction - Evaluates the solution function.

657:    Collective

659:    Input Parameters:
660: +  ts - the `TS` context
661: -  t - current time

663:    Output Parameter:
664: .  U - the solution

666:    Level: developer

668: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
669: @*/
670: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
671: {
672:   TSSolutionFunction solutionfunction;
673:   void              *ctx;
674:   DM                 dm;

676:   PetscFunctionBegin;
679:   PetscCall(TSGetDM(ts, &dm));
680:   PetscCall(DMTSGetSolutionFunction(dm, &solutionfunction, &ctx));
681:   if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
682:   PetscFunctionReturn(PETSC_SUCCESS);
683: }
684: /*@
685:    TSComputeForcingFunction - Evaluates the forcing function.

687:    Collective

689:    Input Parameters:
690: +  ts - the `TS` context
691: -  t - current time

693:    Output Parameter:
694: .  U - the function value

696:    Level: developer

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

706:   PetscFunctionBegin;
709:   PetscCall(TSGetDM(ts, &dm));
710:   PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));

712:   if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
713:   PetscFunctionReturn(PETSC_SUCCESS);
714: }

716: static PetscErrorCode TSGetRHSVec_Private(TS ts, Vec *Frhs)
717: {
718:   Vec F;

720:   PetscFunctionBegin;
721:   *Frhs = NULL;
722:   PetscCall(TSGetIFunction(ts, &F, NULL, NULL));
723:   if (!ts->Frhs) PetscCall(VecDuplicate(F, &ts->Frhs));
724:   *Frhs = ts->Frhs;
725:   PetscFunctionReturn(PETSC_SUCCESS);
726: }

728: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
729: {
730:   Mat         A, B;
731:   TSIJacobian ijacobian;

733:   PetscFunctionBegin;
734:   if (Arhs) *Arhs = NULL;
735:   if (Brhs) *Brhs = NULL;
736:   PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL));
737:   if (Arhs) {
738:     if (!ts->Arhs) {
739:       if (ijacobian) {
740:         PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs));
741:         PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN));
742:       } else {
743:         ts->Arhs = A;
744:         PetscCall(PetscObjectReference((PetscObject)A));
745:       }
746:     } else {
747:       PetscBool flg;
748:       PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
749:       /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
750:       if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
751:         PetscCall(PetscObjectDereference((PetscObject)ts->Arhs));
752:         ts->Arhs = A;
753:         PetscCall(PetscObjectReference((PetscObject)A));
754:       }
755:     }
756:     *Arhs = ts->Arhs;
757:   }
758:   if (Brhs) {
759:     if (!ts->Brhs) {
760:       if (A != B) {
761:         if (ijacobian) {
762:           PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs));
763:         } else {
764:           ts->Brhs = B;
765:           PetscCall(PetscObjectReference((PetscObject)B));
766:         }
767:       } else {
768:         PetscCall(PetscObjectReference((PetscObject)ts->Arhs));
769:         ts->Brhs = ts->Arhs;
770:       }
771:     }
772:     *Brhs = ts->Brhs;
773:   }
774:   PetscFunctionReturn(PETSC_SUCCESS);
775: }

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

780:    Collective

782:    Input Parameters:
783: +  ts - the `TS` context
784: .  t - current time
785: .  U - state vector
786: .  Udot - time derivative of state vector
787: -  imex - flag indicates if the method is `TSIMEX` so that the RHSFunction should be kept separate

789:    Output Parameter:
790: .  Y - right hand side

792:    Level: developer

794:    Note:
795:    Most users should not need to explicitly call this routine, as it
796:    is used internally within the nonlinear solvers.

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

801: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
802: @*/
803: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
804: {
805:   TSIFunction   ifunction;
806:   TSRHSFunction rhsfunction;
807:   void         *ctx;
808:   DM            dm;

810:   PetscFunctionBegin;

816:   PetscCall(TSGetDM(ts, &dm));
817:   PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
818:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));

820:   PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");

822:   PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, Udot, Y));
823:   if (ifunction) {
824:     PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
825:     ts->ifuncs++;
826:   }
827:   if (imex) {
828:     if (!ifunction) PetscCall(VecCopy(Udot, Y));
829:   } else if (rhsfunction) {
830:     if (ifunction) {
831:       Vec Frhs;
832:       PetscCall(TSGetRHSVec_Private(ts, &Frhs));
833:       PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
834:       PetscCall(VecAXPY(Y, -1, Frhs));
835:     } else {
836:       PetscCall(TSComputeRHSFunction(ts, t, U, Y));
837:       PetscCall(VecAYPX(Y, -1, Udot));
838:     }
839:   }
840:   PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, Udot, Y));
841:   PetscFunctionReturn(PETSC_SUCCESS);
842: }

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

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

850: */
851: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
852: {
853:   PetscFunctionBegin;
855:   PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
856:   PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");

858:   if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
859:   if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
860:   if (B && B == ts->Brhs && A != B) {
861:     if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
862:     if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
863:   }
864:   ts->rhsjacobian.shift = 0;
865:   ts->rhsjacobian.scale = 1.;
866:   PetscFunctionReturn(PETSC_SUCCESS);
867: }

869: /*@
870:    TSComputeIJacobian - Evaluates the Jacobian of the DAE

872:    Collective

874:    Input
875:       Input Parameters:
876: +  ts - the `TS` context
877: .  t - current timestep
878: .  U - state vector
879: .  Udot - time derivative of state vector
880: .  shift - shift to apply, see note below
881: -  imex - flag indicates if the method is `TSIMEX` so that the RHSJacobian should be kept separate

883:    Output Parameters:
884: +  A - Jacobian matrix
885: -  B - matrix from which the preconditioner is constructed; often the same as `A`

887:    Level: developer

889:    Notes:
890:    If F(t,U,Udot)=0 is the DAE, the required Jacobian is
891: .vb
892:    dF/dU + shift*dF/dUdot
893: .ve
894:    Most users should not need to explicitly call this routine, as it
895:    is used internally within the nonlinear solvers.

897: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
898: @*/
899: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
900: {
901:   TSIJacobian   ijacobian;
902:   TSRHSJacobian rhsjacobian;
903:   DM            dm;
904:   void         *ctx;

906:   PetscFunctionBegin;

915:   PetscCall(TSGetDM(ts, &dm));
916:   PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
917:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));

919:   PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");

921:   PetscCall(PetscLogEventBegin(TS_JacobianEval, ts, U, A, B));
922:   if (ijacobian) {
923:     PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
924:     ts->ijacs++;
925:   }
926:   if (imex) {
927:     if (!ijacobian) { /* system was written as Udot = G(t,U) */
928:       PetscBool assembled;
929:       if (rhsjacobian) {
930:         Mat Arhs = NULL;
931:         PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
932:         if (A == Arhs) {
933:           PetscCheck(rhsjacobian != TSComputeRHSJacobianConstant, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Unsupported operation! cannot use TSComputeRHSJacobianConstant"); /* there is no way to reconstruct shift*M-J since J cannot be reevaluated */
934:           ts->rhsjacobian.time = PETSC_MIN_REAL;
935:         }
936:       }
937:       PetscCall(MatZeroEntries(A));
938:       PetscCall(MatAssembled(A, &assembled));
939:       if (!assembled) {
940:         PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
941:         PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
942:       }
943:       PetscCall(MatShift(A, shift));
944:       if (A != B) {
945:         PetscCall(MatZeroEntries(B));
946:         PetscCall(MatAssembled(B, &assembled));
947:         if (!assembled) {
948:           PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
949:           PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
950:         }
951:         PetscCall(MatShift(B, shift));
952:       }
953:     }
954:   } else {
955:     Mat Arhs = NULL, Brhs = NULL;

957:     /* RHSJacobian needs to be converted to part of IJacobian if exists */
958:     if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
959:     if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
960:       PetscObjectState Ustate;
961:       PetscObjectId    Uid;
962:       TSRHSFunction    rhsfunction;

964:       PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
965:       PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
966:       PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
967:       if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
968:           ts->rhsjacobian.scale == -1.) {                      /* No need to recompute RHSJacobian */
969:         PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
970:         if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
971:       } else {
972:         PetscBool flg;

974:         if (ts->rhsjacobian.reuse) { /* Undo the damage */
975:           /* MatScale has a short path for this case.
976:              However, this code path is taken the first time TSComputeRHSJacobian is called
977:              and the matrices have not been assembled yet */
978:           PetscCall(TSRecoverRHSJacobian(ts, A, B));
979:         }
980:         PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
981:         PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
982:         /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
983:         if (!flg) {
984:           PetscCall(MatScale(A, -1));
985:           PetscCall(MatShift(A, shift));
986:         }
987:         if (A != B) {
988:           PetscCall(MatScale(B, -1));
989:           PetscCall(MatShift(B, shift));
990:         }
991:       }
992:       ts->rhsjacobian.scale = -1;
993:       ts->rhsjacobian.shift = shift;
994:     } else if (Arhs) {  /* Both IJacobian and RHSJacobian */
995:       if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
996:         PetscCall(MatZeroEntries(A));
997:         PetscCall(MatShift(A, shift));
998:         if (A != B) {
999:           PetscCall(MatZeroEntries(B));
1000:           PetscCall(MatShift(B, shift));
1001:         }
1002:       }
1003:       PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
1004:       PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
1005:       if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
1006:     }
1007:   }
1008:   PetscCall(PetscLogEventEnd(TS_JacobianEval, ts, U, A, B));
1009:   PetscFunctionReturn(PETSC_SUCCESS);
1010: }

1012: /*@C
1013:     TSSetRHSFunction - Sets the routine for evaluating the function,
1014:     where U_t = G(t,u).

1016:     Logically Collective

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

1024:     Calling sequence of `f`:
1025: $     PetscErrorCode f(TS ts, PetscReal t, Vec u, Vec F, void *ctx)
1026: +   ts - timestep context
1027: .   t - current timestep
1028: .   u - input vector
1029: .   F - function vector
1030: -   ctx - [optional] user-defined function context

1032:     Level: beginner

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

1037: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1038: @*/
1039: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, PetscErrorCode (*f)(TS, PetscReal, Vec, Vec, void *), void *ctx)
1040: {
1041:   SNES snes;
1042:   Vec  ralloc = NULL;
1043:   DM   dm;

1045:   PetscFunctionBegin;

1049:   PetscCall(TSGetDM(ts, &dm));
1050:   PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1051:   PetscCall(TSGetSNES(ts, &snes));
1052:   if (!r && !ts->dm && ts->vec_sol) {
1053:     PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1054:     r = ralloc;
1055:   }
1056:   PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1057:   PetscCall(VecDestroy(&ralloc));
1058:   PetscFunctionReturn(PETSC_SUCCESS);
1059: }

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

1064:     Logically Collective

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

1072:     Calling sequence of `f`:
1073: $   PetscErrorCode f(TS ts, PetscReal t, Vec u, void *ctx)
1074: +   ts - the integrator context
1075: .   t - current timestep
1076: .   u - output vector
1077: -   ctx - [optional] user-defined function context

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

1083:     Level: intermediate

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

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

1092: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1093: @*/
1094: PetscErrorCode TSSetSolutionFunction(TS ts, PetscErrorCode (*f)(TS, PetscReal, Vec, void *), void *ctx)
1095: {
1096:   DM dm;

1098:   PetscFunctionBegin;
1100:   PetscCall(TSGetDM(ts, &dm));
1101:   PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1102:   PetscFunctionReturn(PETSC_SUCCESS);
1103: }

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

1108:     Logically Collective

1110:     Input Parameters:
1111: +   ts - the `TS` context obtained from `TSCreate()`
1112: .   func - routine for evaluating the forcing function
1113: -   ctx - [optional] user-defined context for private data for the
1114:           function evaluation routine (may be `NULL`)

1116:     Calling sequence of func:
1117: $   PetscErrorCode func(TS ts, PetscReal t, Vec f, void *ctx)
1118: +   ts - the integrator context
1119: .   t - current timestep
1120: .   f - output vector
1121: -   ctx - [optional] user-defined function context

1123:     Level: intermediate

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

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

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

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

1137: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1138: @*/
1139: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFunction func, void *ctx)
1140: {
1141:   DM dm;

1143:   PetscFunctionBegin;
1145:   PetscCall(TSGetDM(ts, &dm));
1146:   PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1147:   PetscFunctionReturn(PETSC_SUCCESS);
1148: }

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

1154:    Logically Collective

1156:    Input Parameters:
1157: +  ts  - the `TS` context obtained from `TSCreate()`
1158: .  Amat - (approximate) location to store Jacobian matrix entries computed by `f`
1159: .  Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1160: .  f   - the Jacobian evaluation routine
1161: -  ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)

1163:    Calling sequence of `f`:
1164: $  PetscErrorCode f(TS ts, PetscReal t, Vec u, Mat A, Mat B, void *ctx)
1165: +  ts  - the `TS` context obtained from `TSCreate()`
1166: .  t - current timestep
1167: .  u - input vector
1168: .  Amat - (approximate) Jacobian matrix
1169: .  Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1170: -  ctx - [optional] user-defined context for matrix evaluation routine

1172:    Level: beginner

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

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

1180: .seealso: [](ch_ts), `TS`, `SNESComputeJacobianDefaultColor()`, `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`
1181: @*/
1182: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobian f, void *ctx)
1183: {
1184:   SNES        snes;
1185:   DM          dm;
1186:   TSIJacobian ijacobian;

1188:   PetscFunctionBegin;
1192:   if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1193:   if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);

1195:   PetscCall(TSGetDM(ts, &dm));
1196:   PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1197:   PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1198:   PetscCall(TSGetSNES(ts, &snes));
1199:   if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1200:   if (Amat) {
1201:     PetscCall(PetscObjectReference((PetscObject)Amat));
1202:     PetscCall(MatDestroy(&ts->Arhs));
1203:     ts->Arhs = Amat;
1204:   }
1205:   if (Pmat) {
1206:     PetscCall(PetscObjectReference((PetscObject)Pmat));
1207:     PetscCall(MatDestroy(&ts->Brhs));
1208:     ts->Brhs = Pmat;
1209:   }
1210:   PetscFunctionReturn(PETSC_SUCCESS);
1211: }

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

1216:    Logically Collective

1218:    Input Parameters:
1219: +  ts  - the `TS` context obtained from `TSCreate()`
1220: .  r   - vector to hold the residual (or `NULL` to have it created internally)
1221: .  f   - the function evaluation routine
1222: -  ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)

1224:    Calling sequence of `f`:
1225: $   PetscErrorCode f(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
1226: +  ts  - the `TS` context obtained from `TSCreate()`
1227: .  t   - time at step/stage being solved
1228: .  u   - state vector
1229: .  u_t - time derivative of state vector
1230: .  F   - function vector
1231: -  ctx - [optional] user-defined context for matrix evaluation routine

1233:    Level: beginner

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

1238: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`, `TSSetIJacobian()`
1239: @*/
1240: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunction f, void *ctx)
1241: {
1242:   SNES snes;
1243:   Vec  ralloc = NULL;
1244:   DM   dm;

1246:   PetscFunctionBegin;

1250:   PetscCall(TSGetDM(ts, &dm));
1251:   PetscCall(DMTSSetIFunction(dm, f, ctx));

1253:   PetscCall(TSGetSNES(ts, &snes));
1254:   if (!r && !ts->dm && ts->vec_sol) {
1255:     PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1256:     r = ralloc;
1257:   }
1258:   PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1259:   PetscCall(VecDestroy(&ralloc));
1260:   PetscFunctionReturn(PETSC_SUCCESS);
1261: }

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

1266:    Not Collective

1268:    Input Parameter:
1269: .  ts - the `TS` context

1271:    Output Parameters:
1272: +  r - vector to hold residual (or `NULL`)
1273: .  func - the function to compute residual (or `NULL`)
1274: -  ctx - the function context (or `NULL`)

1276:    Level: advanced

1278: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1279: @*/
1280: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunction *func, void **ctx)
1281: {
1282:   SNES snes;
1283:   DM   dm;

1285:   PetscFunctionBegin;
1287:   PetscCall(TSGetSNES(ts, &snes));
1288:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1289:   PetscCall(TSGetDM(ts, &dm));
1290:   PetscCall(DMTSGetIFunction(dm, func, ctx));
1291:   PetscFunctionReturn(PETSC_SUCCESS);
1292: }

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

1297:    Not Collective

1299:    Input Parameter:
1300: .  ts - the `TS` context

1302:    Output Parameters:
1303: +  r - vector to hold computed right hand side (or `NULL`)
1304: .  func - the function to compute right hand side (or `NULL`)
1305: -  ctx - the function context (or `NULL`)

1307:    Level: advanced

1309: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1310: @*/
1311: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunction *func, void **ctx)
1312: {
1313:   SNES snes;
1314:   DM   dm;

1316:   PetscFunctionBegin;
1318:   PetscCall(TSGetSNES(ts, &snes));
1319:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1320:   PetscCall(TSGetDM(ts, &dm));
1321:   PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1322:   PetscFunctionReturn(PETSC_SUCCESS);
1323: }

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

1329:    Logically Collective

1331:    Input Parameters:
1332: +  ts  - the `TS` context obtained from `TSCreate()`
1333: .  Amat - (approximate) matrix to store Jacobian entries computed by `f`
1334: .  Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1335: .  f   - the Jacobian evaluation routine
1336: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)

1338:    Calling sequence of `f`:
1339: $    PetscErrorCode f(TS ts, PetscReal t, Vec U, Vec U_t, PetscReal a, Mat Amat, Mat Pmat, void *ctx)
1340: +  ts  - the `TS` context obtained from `TSCreate()`
1341: .  t    - time at step/stage being solved
1342: .  U    - state vector
1343: .  U_t  - time derivative of state vector
1344: .  a    - shift
1345: .  Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1346: .  Pmat - matrix used for constructing preconditioner, usually the same as `Amat`
1347: -  ctx  - [optional] user-defined context for matrix evaluation routine

1349:    Level: beginner

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

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

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

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

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

1369: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetRHSJacobian()`, `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1370: @*/
1371: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobian f, void *ctx)
1372: {
1373:   SNES snes;
1374:   DM   dm;

1376:   PetscFunctionBegin;
1380:   if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1381:   if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);

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

1386:   PetscCall(TSGetSNES(ts, &snes));
1387:   PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1388:   PetscFunctionReturn(PETSC_SUCCESS);
1389: }

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

1397:    Logically Collective

1399:    Input Parameters:
1400: +  ts - `TS` context obtained from `TSCreate()`
1401: -  reuse - `PETSC_TRUE` if the RHS Jacobian

1403:    Level: intermediate

1405: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1406: @*/
1407: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1408: {
1409:   PetscFunctionBegin;
1410:   ts->rhsjacobian.reuse = reuse;
1411:   PetscFunctionReturn(PETSC_SUCCESS);
1412: }

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

1417:    Logically Collective

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

1425:    Calling sequence of `fun`:
1426: $   PetscErrorCode fun(TS ts, PetscReal t, Vec U, Vec U_t, Vec U_tt, Vec F,void *ctx);
1427: +  ts  - the `TS` context obtained from `TSCreate()`
1428: .  t    - time at step/stage being solved
1429: .  U    - state vector
1430: .  U_t  - time derivative of state vector
1431: .  U_tt - second time derivative of state vector
1432: .  F    - function vector
1433: -  ctx  - [optional] user-defined context for matrix evaluation routine (may be `NULL`)

1435:    Level: beginner

1437: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`, `TSSetIFunction()`, `TSCreate()`, `TSSetRHSFunction()`
1438: @*/
1439: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2Function fun, void *ctx)
1440: {
1441:   DM dm;

1443:   PetscFunctionBegin;
1446:   PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1447:   PetscCall(TSGetDM(ts, &dm));
1448:   PetscCall(DMTSSetI2Function(dm, fun, ctx));
1449:   PetscFunctionReturn(PETSC_SUCCESS);
1450: }

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

1455:   Not Collective

1457:   Input Parameter:
1458: . ts - the `TS` context

1460:   Output Parameters:
1461: + r - vector to hold residual (or `NULL`)
1462: . fun - the function to compute residual (or `NULL`)
1463: - ctx - the function context (or `NULL`)

1465:   Level: advanced

1467: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1468: @*/
1469: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2Function *fun, void **ctx)
1470: {
1471:   SNES snes;
1472:   DM   dm;

1474:   PetscFunctionBegin;
1476:   PetscCall(TSGetSNES(ts, &snes));
1477:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1478:   PetscCall(TSGetDM(ts, &dm));
1479:   PetscCall(DMTSGetI2Function(dm, fun, ctx));
1480:   PetscFunctionReturn(PETSC_SUCCESS);
1481: }

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

1487:    Logically Collective

1489:    Input Parameters:
1490: +  ts  - the `TS` context obtained from `TSCreate()`
1491: .  J   - matrix to hold the Jacobian values
1492: .  P   - matrix for constructing the preconditioner (may be same as `J`)
1493: .  jac - the Jacobian evaluation routine
1494: -  ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)

1496:    Calling sequence of `jac`:
1497: $    PetscErrorCode jac(TS ts, PetscReal t, Vec U, Vec U_t, Vec U_tt, PetscReal v, PetscReal a, Mat J, Mat P, void *ctx)
1498: +  ts  - the `TS` context obtained from `TSCreate()`
1499: .  t    - time at step/stage being solved
1500: .  U    - state vector
1501: .  U_t  - time derivative of state vector
1502: .  U_tt - second time derivative of state vector
1503: .  v    - shift for U_t
1504: .  a    - shift for U_tt
1505: .  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
1506: .  P    - preconditioning matrix for J, may be same as J
1507: -  ctx  - [optional] user-defined context for matrix evaluation routine

1509:    Level: beginner

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

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

1519: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1520: @*/
1521: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2Jacobian jac, void *ctx)
1522: {
1523:   DM dm;

1525:   PetscFunctionBegin;
1529:   PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1530:   PetscCall(TSGetDM(ts, &dm));
1531:   PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1532:   PetscFunctionReturn(PETSC_SUCCESS);
1533: }

1535: /*@C
1536:   TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.

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

1540:   Input Parameter:
1541: . ts  - The `TS` context obtained from `TSCreate()`

1543:   Output Parameters:
1544: + J  - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1545: . P - The matrix from which the preconditioner is constructed, often the same as `J`
1546: . jac - The function to compute the Jacobian matrices
1547: - ctx - User-defined context for Jacobian evaluation routine

1549:   Level: advanced

1551:   Note:
1552:     You can pass in `NULL` for any return argument you do not need.

1554: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1555: @*/
1556: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2Jacobian *jac, void **ctx)
1557: {
1558:   SNES snes;
1559:   DM   dm;

1561:   PetscFunctionBegin;
1562:   PetscCall(TSGetSNES(ts, &snes));
1563:   PetscCall(SNESSetUpMatrices(snes));
1564:   PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1565:   PetscCall(TSGetDM(ts, &dm));
1566:   PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1567:   PetscFunctionReturn(PETSC_SUCCESS);
1568: }

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

1573:   Collective

1575:   Input Parameters:
1576: + ts - the `TS` context
1577: . t - current time
1578: . U - state vector
1579: . V - time derivative of state vector (U_t)
1580: - A - second time derivative of state vector (U_tt)

1582:   Output Parameter:
1583: . F - the residual vector

1585:   Level: developer

1587:   Note:
1588:   Most users should not need to explicitly call this routine, as it
1589:   is used internally within the nonlinear solvers.

1591: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1592: @*/
1593: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1594: {
1595:   DM            dm;
1596:   TSI2Function  I2Function;
1597:   void         *ctx;
1598:   TSRHSFunction rhsfunction;

1600:   PetscFunctionBegin;

1607:   PetscCall(TSGetDM(ts, &dm));
1608:   PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1609:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));

1611:   if (!I2Function) {
1612:     PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1613:     PetscFunctionReturn(PETSC_SUCCESS);
1614:   }

1616:   PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, V, F));

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

1620:   if (rhsfunction) {
1621:     Vec Frhs;
1622:     PetscCall(TSGetRHSVec_Private(ts, &Frhs));
1623:     PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1624:     PetscCall(VecAXPY(F, -1, Frhs));
1625:   }

1627:   PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, V, F));
1628:   PetscFunctionReturn(PETSC_SUCCESS);
1629: }

1631: /*@
1632:   TSComputeI2Jacobian - Evaluates the Jacobian of the DAE

1634:   Collective

1636:   Input Parameters:
1637: + ts - the `TS` context
1638: . t - current timestep
1639: . U - state vector
1640: . V - time derivative of state vector
1641: . A - second time derivative of state vector
1642: . shiftV - shift to apply, see note below
1643: - shiftA - shift to apply, see note below

1645:   Output Parameters:
1646: + J - Jacobian matrix
1647: - P - optional preconditioning matrix

1649:   Level: developer

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

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

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

1659: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1660: @*/
1661: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1662: {
1663:   DM            dm;
1664:   TSI2Jacobian  I2Jacobian;
1665:   void         *ctx;
1666:   TSRHSJacobian rhsjacobian;

1668:   PetscFunctionBegin;

1676:   PetscCall(TSGetDM(ts, &dm));
1677:   PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1678:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));

1680:   if (!I2Jacobian) {
1681:     PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1682:     PetscFunctionReturn(PETSC_SUCCESS);
1683:   }

1685:   PetscCall(PetscLogEventBegin(TS_JacobianEval, ts, U, J, P));
1686:   PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1687:   if (rhsjacobian) {
1688:     Mat Jrhs, Prhs;
1689:     PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1690:     PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1691:     PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1692:     if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1693:   }

1695:   PetscCall(PetscLogEventEnd(TS_JacobianEval, ts, U, J, P));
1696:   PetscFunctionReturn(PETSC_SUCCESS);
1697: }

1699: /*@C
1700:    TSSetTransientVariable - sets function to transform from state to transient variables

1702:    Logically Collective

1704:    Input Parameters:
1705: +  ts - time stepping context on which to change the transient variable
1706: .  tvar - a function that transforms to transient variables
1707: -  ctx - a context for tvar

1709:     Calling sequence of `tvar`:
1710: $   PetscErrorCode tvar(TS ts, Vec p, Vec c, void *ctx)
1711: +   ts - timestep context
1712: .   p - input vector (primitive form)
1713: .   c - output vector, transient variables (conservative form)
1714: -   ctx - [optional] user-defined function context

1716:    Level: advanced

1718:    Notes:
1719:    This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1720:    can be conservative.  In this context, primitive variables P are used to model the state (e.g., because they lead to
1721:    well-conditioned formulations even in limiting cases such as low-Mach or zero porosity).  The transient variable is
1722:    C(P), specified by calling this function.  An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1723:    evaluated via the chain rule, as in
1724: .vb
1725:      dF/dP + shift * dF/dCdot dC/dP.
1726: .ve

1728: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1729: @*/
1730: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariable tvar, void *ctx)
1731: {
1732:   DM dm;

1734:   PetscFunctionBegin;
1736:   PetscCall(TSGetDM(ts, &dm));
1737:   PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1738:   PetscFunctionReturn(PETSC_SUCCESS);
1739: }

1741: /*@
1742:    TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables

1744:    Logically Collective

1746:    Input Parameters:
1747: +  ts - TS on which to compute
1748: -  U - state vector to be transformed to transient variables

1750:    Output Parameter:
1751: .  C - transient (conservative) variable

1753:    Level: developer

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

1760: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1761: @*/
1762: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1763: {
1764:   DM   dm;
1765:   DMTS dmts;

1767:   PetscFunctionBegin;
1770:   PetscCall(TSGetDM(ts, &dm));
1771:   PetscCall(DMGetDMTS(dm, &dmts));
1772:   if (dmts->ops->transientvar) {
1774:     PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1775:   }
1776:   PetscFunctionReturn(PETSC_SUCCESS);
1777: }

1779: /*@
1780:    TSHasTransientVariable - determine whether transient variables have been set

1782:    Logically Collective

1784:    Input Parameter:
1785: .  ts - `TS` on which to compute

1787:    Output Parameter:
1788: .  has - `PETSC_TRUE` if transient variables have been set

1790:    Level: developer

1792: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1793: @*/
1794: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1795: {
1796:   DM   dm;
1797:   DMTS dmts;

1799:   PetscFunctionBegin;
1801:   PetscCall(TSGetDM(ts, &dm));
1802:   PetscCall(DMGetDMTS(dm, &dmts));
1803:   *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1804:   PetscFunctionReturn(PETSC_SUCCESS);
1805: }

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

1811:    Logically Collective

1813:    Input Parameters:
1814: +  ts - the `TS` context obtained from `TSCreate()`
1815: .  u - the solution vector
1816: -  v - the time derivative vector

1818:    Level: beginner

1820: .seealso: [](ch_ts), `TS`
1821: @*/
1822: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1823: {
1824:   PetscFunctionBegin;
1828:   PetscCall(TSSetSolution(ts, u));
1829:   PetscCall(PetscObjectReference((PetscObject)v));
1830:   PetscCall(VecDestroy(&ts->vec_dot));
1831:   ts->vec_dot = v;
1832:   PetscFunctionReturn(PETSC_SUCCESS);
1833: }

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

1841:    Not Collective, but `u` returned is parallel if `TS` is parallel

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

1846:    Output Parameters:
1847: +  u - the vector containing the solution
1848: -  v - the vector containing the time derivative

1850:    Level: intermediate

1852: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1853: @*/
1854: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1855: {
1856:   PetscFunctionBegin;
1860:   if (u) *u = ts->vec_sol;
1861:   if (v) *v = ts->vec_dot;
1862:   PetscFunctionReturn(PETSC_SUCCESS);
1863: }

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

1868:   Collective

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

1875:    Level: intermediate

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

1880: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1881: @*/
1882: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1883: {
1884:   PetscBool isbinary;
1885:   PetscInt  classid;
1886:   char      type[256];
1887:   DMTS      sdm;
1888:   DM        dm;

1890:   PetscFunctionBegin;
1893:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1894:   PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");

1896:   PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1897:   PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1898:   PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1899:   PetscCall(TSSetType(ts, type));
1900:   PetscTryTypeMethod(ts, load, viewer);
1901:   PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1902:   PetscCall(DMLoad(dm, viewer));
1903:   PetscCall(TSSetDM(ts, dm));
1904:   PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1905:   PetscCall(VecLoad(ts->vec_sol, viewer));
1906:   PetscCall(DMGetDMTS(ts->dm, &sdm));
1907:   PetscCall(DMTSLoad(sdm, viewer));
1908:   PetscFunctionReturn(PETSC_SUCCESS);
1909: }

1911: #include <petscdraw.h>
1912: #if defined(PETSC_HAVE_SAWS)
1913: #include <petscviewersaws.h>
1914: #endif

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

1919:    Collective

1921:    Input Parameters:
1922: +  ts - the `TS` context
1923: .  obj - Optional object that provides the prefix for the options database keys
1924: -  name - command line option string to be passed by user

1926:    Level: intermediate

1928: .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1929: @*/
1930: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1931: {
1932:   PetscFunctionBegin;
1934:   PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1935:   PetscFunctionReturn(PETSC_SUCCESS);
1936: }

1938: /*@C
1939:     TSView - Prints the `TS` data structure.

1941:     Collective

1943:     Input Parameters:
1944: +   ts - the `TS` context obtained from `TSCreate()`
1945: -   viewer - visualization context

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

1950:     Level: beginner

1952:     Notes:
1953:     The available visualization contexts include
1954: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1955: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1956:          output where only the first processor opens
1957:          the file.  All other processors send their
1958:          data to the first processor to print.

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

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

1965: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1966: @*/
1967: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1968: {
1969:   TSType    type;
1970:   PetscBool iascii, isstring, isundials, isbinary, isdraw;
1971:   DMTS      sdm;
1972: #if defined(PETSC_HAVE_SAWS)
1973:   PetscBool issaws;
1974: #endif

1976:   PetscFunctionBegin;
1978:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1980:   PetscCheckSameComm(ts, 1, viewer, 2);

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

2027:     PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
2028:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
2029:     if (rank == 0) {
2030:       PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
2031:       PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
2032:       PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
2033:     }
2034:     PetscTryTypeMethod(ts, view, viewer);
2035:     if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
2036:     PetscCall(DMView(ts->dm, viewer));
2037:     PetscCall(VecView(ts->vec_sol, viewer));
2038:     PetscCall(DMGetDMTS(ts->dm, &sdm));
2039:     PetscCall(DMTSView(sdm, viewer));
2040:   } else if (isdraw) {
2041:     PetscDraw draw;
2042:     char      str[36];
2043:     PetscReal x, y, bottom, h;

2045:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2046:     PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
2047:     PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
2048:     PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
2049:     PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
2050:     bottom = y - h;
2051:     PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
2052:     PetscTryTypeMethod(ts, view, viewer);
2053:     if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
2054:     if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
2055:     PetscCall(PetscDrawPopCurrentPoint(draw));
2056: #if defined(PETSC_HAVE_SAWS)
2057:   } else if (issaws) {
2058:     PetscMPIInt rank;
2059:     const char *name;

2061:     PetscCall(PetscObjectGetName((PetscObject)ts, &name));
2062:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
2063:     if (!((PetscObject)ts)->amsmem && rank == 0) {
2064:       char dir[1024];

2066:       PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
2067:       PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
2068:       PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
2069:       PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
2070:       PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
2071:     }
2072:     PetscTryTypeMethod(ts, view, viewer);
2073: #endif
2074:   }
2075:   if (ts->snes && ts->usessnes) {
2076:     PetscCall(PetscViewerASCIIPushTab(viewer));
2077:     PetscCall(SNESView(ts->snes, viewer));
2078:     PetscCall(PetscViewerASCIIPopTab(viewer));
2079:   }
2080:   PetscCall(DMGetDMTS(ts->dm, &sdm));
2081:   PetscCall(DMTSView(sdm, viewer));

2083:   PetscCall(PetscViewerASCIIPushTab(viewer));
2084:   PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials));
2085:   PetscCall(PetscViewerASCIIPopTab(viewer));
2086:   PetscFunctionReturn(PETSC_SUCCESS);
2087: }

2089: /*@
2090:    TSSetApplicationContext - Sets an optional user-defined context for
2091:    the timesteppers.

2093:    Logically Collective

2095:    Input Parameters:
2096: +  ts - the `TS` context obtained from `TSCreate()`
2097: -  usrP -  user context

2099:    Level: intermediate

2101:    Fortran Note:
2102:     You must write a Fortran interface definition for this
2103:     function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.

2105: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2106: @*/
2107: PetscErrorCode TSSetApplicationContext(TS ts, void *usrP)
2108: {
2109:   PetscFunctionBegin;
2111:   ts->user = usrP;
2112:   PetscFunctionReturn(PETSC_SUCCESS);
2113: }

2115: /*@
2116:     TSGetApplicationContext - Gets the user-defined context for the
2117:     timestepper that was set with `TSSetApplicationContext()`

2119:     Not Collective

2121:     Input Parameter:
2122: .   ts - the `TS` context obtained from `TSCreate()`

2124:     Output Parameter:
2125: .   usrP - user context

2127:     Level: intermediate

2129:     Fortran Note:
2130:     You must write a Fortran interface definition for this
2131:     function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.

2133: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2134: @*/
2135: PetscErrorCode TSGetApplicationContext(TS ts, void *usrP)
2136: {
2137:   PetscFunctionBegin;
2139:   *(void **)usrP = ts->user;
2140:   PetscFunctionReturn(PETSC_SUCCESS);
2141: }

2143: /*@
2144:    TSGetStepNumber - Gets the number of time steps completed.

2146:    Not Collective

2148:    Input Parameter:
2149: .  ts - the `TS` context obtained from `TSCreate()`

2151:    Output Parameter:
2152: .  steps - number of steps completed so far

2154:    Level: intermediate

2156: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2157: @*/
2158: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2159: {
2160:   PetscFunctionBegin;
2163:   *steps = ts->steps;
2164:   PetscFunctionReturn(PETSC_SUCCESS);
2165: }

2167: /*@
2168:    TSSetStepNumber - Sets the number of steps completed.

2170:    Logically Collective

2172:    Input Parameters:
2173: +  ts - the `TS` context
2174: -  steps - number of steps completed so far

2176:    Level: developer

2178:    Note:
2179:    For most uses of the `TS` solvers the user need not explicitly call
2180:    `TSSetStepNumber()`, as the step counter is appropriately updated in
2181:    `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2182:    reinitialize timestepping by setting the step counter to zero (and time
2183:    to the initial time) to solve a similar problem with different initial
2184:    conditions or parameters. Other possible use case is to continue
2185:    timestepping from a previously interrupted run in such a way that `TS`
2186:    monitors will be called with a initial nonzero step counter.

2188: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2189: @*/
2190: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2191: {
2192:   PetscFunctionBegin;
2195:   PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2196:   ts->steps = steps;
2197:   PetscFunctionReturn(PETSC_SUCCESS);
2198: }

2200: /*@
2201:    TSSetTimeStep - Allows one to reset the timestep at any time,
2202:    useful for simple pseudo-timestepping codes.

2204:    Logically Collective

2206:    Input Parameters:
2207: +  ts - the `TS` context obtained from `TSCreate()`
2208: -  time_step - the size of the timestep

2210:    Level: intermediate

2212: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2213: @*/
2214: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2215: {
2216:   PetscFunctionBegin;
2219:   ts->time_step = time_step;
2220:   PetscFunctionReturn(PETSC_SUCCESS);
2221: }

2223: /*@
2224:    TSSetExactFinalTime - Determines whether to adapt the final time step to
2225:      match the exact final time, interpolate solution to the exact final time,
2226:      or just return at the final time `TS` computed.

2228:   Logically Collective

2230:    Input Parameters:
2231: +   ts - the time-step context
2232: -   eftopt - exact final time option
2233: .vb
2234:   TS_EXACTFINALTIME_STEPOVER    - Don't do anything if final time is exceeded
2235:   TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2236:   TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2237: .ve

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

2242:    Level: beginner

2244:    Note:
2245:    If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2246:    then the final time you selected.

2248: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2249: @*/
2250: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2251: {
2252:   PetscFunctionBegin;
2255:   ts->exact_final_time = eftopt;
2256:   PetscFunctionReturn(PETSC_SUCCESS);
2257: }

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

2262:    Not Collective

2264:    Input Parameter:
2265: .  ts - the `TS` context

2267:    Output Parameter:
2268: .  eftopt - exact final time option

2270:    Level: beginner

2272: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2273: @*/
2274: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2275: {
2276:   PetscFunctionBegin;
2279:   *eftopt = ts->exact_final_time;
2280:   PetscFunctionReturn(PETSC_SUCCESS);
2281: }

2283: /*@
2284:    TSGetTimeStep - Gets the current timestep size.

2286:    Not Collective

2288:    Input Parameter:
2289: .  ts - the `TS` context obtained from `TSCreate()`

2291:    Output Parameter:
2292: .  dt - the current timestep size

2294:    Level: intermediate

2296: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2297: @*/
2298: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2299: {
2300:   PetscFunctionBegin;
2303:   *dt = ts->time_step;
2304:   PetscFunctionReturn(PETSC_SUCCESS);
2305: }

2307: /*@
2308:    TSGetSolution - Returns the solution at the present timestep. It
2309:    is valid to call this routine inside the function that you are evaluating
2310:    in order to move to the new timestep. This vector not changed until
2311:    the solution at the next timestep has been calculated.

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

2315:    Input Parameter:
2316: .  ts - the `TS` context obtained from `TSCreate()`

2318:    Output Parameter:
2319: .  v - the vector containing the solution

2321:    Level: intermediate

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

2327: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2328: @*/
2329: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2330: {
2331:   PetscFunctionBegin;
2334:   *v = ts->vec_sol;
2335:   PetscFunctionReturn(PETSC_SUCCESS);
2336: }

2338: /*@
2339:    TSGetSolutionComponents - Returns any solution components at the present
2340:    timestep, if available for the time integration method being used.
2341:    Solution components are quantities that share the same size and
2342:    structure as the solution vector.

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

2346:    Parameters :
2347: +  ts - the `TS` context obtained from `TSCreate()` (input parameter).
2348: .  n - If v is `NULL`, then the number of solution components is
2349:        returned through n, else the n-th solution component is
2350:        returned in v.
2351: -  v - the vector containing the n-th solution component
2352:        (may be `NULL` to use this function to find out
2353:         the number of solutions components).

2355:    Level: advanced

2357: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2358: @*/
2359: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2360: {
2361:   PetscFunctionBegin;
2363:   if (!ts->ops->getsolutioncomponents) *n = 0;
2364:   else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2365:   PetscFunctionReturn(PETSC_SUCCESS);
2366: }

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

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

2374:    Parameters :
2375: +  ts - the `TS` context obtained from `TSCreate()` (input parameter).
2376: -  v - the vector containing the auxiliary solution

2378:    Level: intermediate

2380: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2381: @*/
2382: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2383: {
2384:   PetscFunctionBegin;
2386:   if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2387:   else PetscCall(VecZeroEntries(*v));
2388:   PetscFunctionReturn(PETSC_SUCCESS);
2389: }

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

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

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

2402:    Level: intermediate

2404:    Note:
2405:    MUST call after `TSSetUp()`

2407: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2408: @*/
2409: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2410: {
2411:   PetscFunctionBegin;
2413:   if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2414:   else PetscCall(VecZeroEntries(*v));
2415:   PetscFunctionReturn(PETSC_SUCCESS);
2416: }

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

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

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

2429:    Level: intermediate

2431: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError)`
2432: @*/
2433: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2434: {
2435:   PetscFunctionBegin;
2437:   PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2438:   PetscTryTypeMethod(ts, settimeerror, v);
2439:   PetscFunctionReturn(PETSC_SUCCESS);
2440: }

2442: /* ----- Routines to initialize and destroy a timestepper ---- */
2443: /*@
2444:   TSSetProblemType - Sets the type of problem to be solved.

2446:   Not collective

2448:   Input Parameters:
2449: + ts   - The `TS`
2450: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2451: .vb
2452:          U_t - A U = 0      (linear)
2453:          U_t - A(t) U = 0   (linear)
2454:          F(t,U,U_t) = 0     (nonlinear)
2455: .ve

2457:    Level: beginner

2459: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2460: @*/
2461: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2462: {
2463:   PetscFunctionBegin;
2465:   ts->problem_type = type;
2466:   if (type == TS_LINEAR) {
2467:     SNES snes;
2468:     PetscCall(TSGetSNES(ts, &snes));
2469:     PetscCall(SNESSetType(snes, SNESKSPONLY));
2470:   }
2471:   PetscFunctionReturn(PETSC_SUCCESS);
2472: }

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

2477:   Not collective

2479:   Input Parameter:
2480: . ts   - The `TS`

2482:   Output Parameter:
2483: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2484: .vb
2485:          M U_t = A U
2486:          M(t) U_t = A(t) U
2487:          F(t,U,U_t)
2488: .ve

2490:    Level: beginner

2492: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2493: @*/
2494: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2495: {
2496:   PetscFunctionBegin;
2499:   *type = ts->problem_type;
2500:   PetscFunctionReturn(PETSC_SUCCESS);
2501: }

2503: /*
2504:     Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2505: */
2506: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2507: {
2508:   PetscBool isnone;

2510:   PetscFunctionBegin;
2511:   PetscCall(TSGetAdapt(ts, &ts->adapt));
2512:   PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));

2514:   PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2515:   if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2516:   else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2517:   PetscFunctionReturn(PETSC_SUCCESS);
2518: }

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

2523:    Collective

2525:    Input Parameter:
2526: .  ts - the `TS` context obtained from `TSCreate()`

2528:    Level: advanced

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

2537: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2538: @*/
2539: PetscErrorCode TSSetUp(TS ts)
2540: {
2541:   DM dm;
2542:   PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2543:   PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2544:   TSIFunction   ifun;
2545:   TSIJacobian   ijac;
2546:   TSI2Jacobian  i2jac;
2547:   TSRHSJacobian rhsjac;

2549:   PetscFunctionBegin;
2551:   if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);

2553:   if (!((PetscObject)ts)->type_name) {
2554:     PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2555:     PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2556:   }

2558:   if (!ts->vec_sol) {
2559:     PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2560:     PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2561:   }

2563:   if (ts->tspan) {
2564:     if (!ts->tspan->vecs_sol) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2565:   }
2566:   if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2567:     PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2568:     ts->Jacp = ts->Jacprhs;
2569:   }

2571:   if (ts->quadraturets) {
2572:     PetscCall(TSSetUp(ts->quadraturets));
2573:     PetscCall(VecDestroy(&ts->vec_costintegrand));
2574:     PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2575:   }

2577:   PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
2578:   if (rhsjac == TSComputeRHSJacobianConstant) {
2579:     Mat  Amat, Pmat;
2580:     SNES snes;
2581:     PetscCall(TSGetSNES(ts, &snes));
2582:     PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
2583:     /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2584:      * have displaced the RHS matrix */
2585:     if (Amat && Amat == ts->Arhs) {
2586:       /* 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 */
2587:       PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
2588:       PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
2589:       PetscCall(MatDestroy(&Amat));
2590:     }
2591:     if (Pmat && Pmat == ts->Brhs) {
2592:       PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
2593:       PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
2594:       PetscCall(MatDestroy(&Pmat));
2595:     }
2596:   }

2598:   PetscCall(TSGetAdapt(ts, &ts->adapt));
2599:   PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));

2601:   PetscTryTypeMethod(ts, setup);

2603:   PetscCall(TSSetExactFinalTimeDefault(ts));

2605:   /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2606:      to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2607:    */
2608:   PetscCall(TSGetDM(ts, &dm));
2609:   PetscCall(DMSNESGetFunction(dm, &func, NULL));
2610:   if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));

2612:   /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2613:      Otherwise, the SNES will use coloring internally to form the Jacobian.
2614:    */
2615:   PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2616:   PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2617:   PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2618:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2619:   if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));

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

2624:   ts->setupcalled = PETSC_TRUE;
2625:   PetscFunctionReturn(PETSC_SUCCESS);
2626: }

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

2631:    Collective

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

2636:    Level: beginner

2638: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()`
2639: @*/
2640: PetscErrorCode TSReset(TS ts)
2641: {
2642:   TS_RHSSplitLink ilink = ts->tsrhssplit, next;

2644:   PetscFunctionBegin;

2647:   PetscTryTypeMethod(ts, reset);
2648:   if (ts->snes) PetscCall(SNESReset(ts->snes));
2649:   if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));

2651:   PetscCall(MatDestroy(&ts->Arhs));
2652:   PetscCall(MatDestroy(&ts->Brhs));
2653:   PetscCall(VecDestroy(&ts->Frhs));
2654:   PetscCall(VecDestroy(&ts->vec_sol));
2655:   PetscCall(VecDestroy(&ts->vec_dot));
2656:   PetscCall(VecDestroy(&ts->vatol));
2657:   PetscCall(VecDestroy(&ts->vrtol));
2658:   PetscCall(VecDestroyVecs(ts->nwork, &ts->work));

2660:   PetscCall(MatDestroy(&ts->Jacprhs));
2661:   PetscCall(MatDestroy(&ts->Jacp));
2662:   if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2663:   if (ts->quadraturets) {
2664:     PetscCall(TSReset(ts->quadraturets));
2665:     PetscCall(VecDestroy(&ts->vec_costintegrand));
2666:   }
2667:   while (ilink) {
2668:     next = ilink->next;
2669:     PetscCall(TSDestroy(&ilink->ts));
2670:     PetscCall(PetscFree(ilink->splitname));
2671:     PetscCall(ISDestroy(&ilink->is));
2672:     PetscCall(PetscFree(ilink));
2673:     ilink = next;
2674:   }
2675:   ts->tsrhssplit     = NULL;
2676:   ts->num_rhs_splits = 0;
2677:   if (ts->tspan) {
2678:     PetscCall(PetscFree(ts->tspan->span_times));
2679:     PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2680:     PetscCall(PetscFree(ts->tspan));
2681:   }
2682:   ts->setupcalled = PETSC_FALSE;
2683:   PetscFunctionReturn(PETSC_SUCCESS);
2684: }

2686: /*@C
2687:    TSDestroy - Destroys the timestepper context that was created
2688:    with `TSCreate()`.

2690:    Collective

2692:    Input Parameter:
2693: .  ts - the `TS` context obtained from `TSCreate()`

2695:    Level: beginner

2697: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2698: @*/
2699: PetscErrorCode TSDestroy(TS *ts)
2700: {
2701:   PetscFunctionBegin;
2702:   if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2704:   if (--((PetscObject)(*ts))->refct > 0) {
2705:     *ts = NULL;
2706:     PetscFunctionReturn(PETSC_SUCCESS);
2707:   }

2709:   PetscCall(TSReset(*ts));
2710:   PetscCall(TSAdjointReset(*ts));
2711:   if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));

2713:   /* if memory was published with SAWs then destroy it */
2714:   PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2715:   PetscTryTypeMethod((*ts), destroy);

2717:   PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));

2719:   PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2720:   PetscCall(TSEventDestroy(&(*ts)->event));

2722:   PetscCall(SNESDestroy(&(*ts)->snes));
2723:   PetscCall(DMDestroy(&(*ts)->dm));
2724:   PetscCall(TSMonitorCancel((*ts)));
2725:   PetscCall(TSAdjointMonitorCancel((*ts)));

2727:   PetscCall(TSDestroy(&(*ts)->quadraturets));
2728:   PetscCall(PetscHeaderDestroy(ts));
2729:   PetscFunctionReturn(PETSC_SUCCESS);
2730: }

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

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

2738:    Input Parameter:
2739: .  ts - the `TS` context obtained from `TSCreate()`

2741:    Output Parameter:
2742: .  snes - the nonlinear solver context

2744:    Level: beginner

2746:    Notes:
2747:    The user can then directly manipulate the `SNES` context to set various
2748:    options, etc.  Likewise, the user can then extract and manipulate the
2749:    `KSP`, and `PC` contexts as well.

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

2754: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2755: @*/
2756: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2757: {
2758:   PetscFunctionBegin;
2761:   if (!ts->snes) {
2762:     PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2763:     PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2764:     PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2765:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2766:     if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2767:     if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2768:   }
2769:   *snes = ts->snes;
2770:   PetscFunctionReturn(PETSC_SUCCESS);
2771: }

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

2776:    Collective

2778:    Input Parameters:
2779: +  ts - the `TS` context obtained from `TSCreate()`
2780: -  snes - the nonlinear solver context

2782:    Level: developer

2784:    Note:
2785:    Most users should have the `TS` created by calling `TSGetSNES()`

2787: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2788: @*/
2789: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2790: {
2791:   PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);

2793:   PetscFunctionBegin;
2796:   PetscCall(PetscObjectReference((PetscObject)snes));
2797:   PetscCall(SNESDestroy(&ts->snes));

2799:   ts->snes = snes;

2801:   PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2802:   PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2803:   if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2804:   PetscFunctionReturn(PETSC_SUCCESS);
2805: }

2807: /*@
2808:    TSGetKSP - Returns the `KSP` (linear solver) associated with
2809:    a `TS` (timestepper) context.

2811:    Not Collective, but `ksp` is parallel if `ts` is parallel

2813:    Input Parameter:
2814: .  ts - the `TS` context obtained from `TSCreate()`

2816:    Output Parameter:
2817: .  ksp - the nonlinear solver context

2819:    Level: beginner

2821:    Notes:
2822:    The user can then directly manipulate the `KSP` context to set various
2823:    options, etc.  Likewise, the user can then extract and manipulate the
2824:    `PC` context as well.

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

2829: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2830: @*/
2831: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2832: {
2833:   SNES snes;

2835:   PetscFunctionBegin;
2838:   PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2839:   PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
2840:   PetscCall(TSGetSNES(ts, &snes));
2841:   PetscCall(SNESGetKSP(snes, ksp));
2842:   PetscFunctionReturn(PETSC_SUCCESS);
2843: }

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

2847: /*@
2848:    TSSetMaxSteps - Sets the maximum number of steps to use.

2850:    Logically Collective

2852:    Input Parameters:
2853: +  ts - the `TS` context obtained from `TSCreate()`
2854: -  maxsteps - maximum number of steps to use

2856:    Options Database Key:
2857: .  -ts_max_steps <maxsteps> - Sets maxsteps

2859:    Level: intermediate

2861:    Note:
2862:    The default maximum number of steps is 5000

2864: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2865: @*/
2866: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2867: {
2868:   PetscFunctionBegin;
2871:   PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2872:   ts->max_steps = maxsteps;
2873:   PetscFunctionReturn(PETSC_SUCCESS);
2874: }

2876: /*@
2877:    TSGetMaxSteps - Gets the maximum number of steps to use.

2879:    Not Collective

2881:    Input Parameter:
2882: .  ts - the `TS` context obtained from `TSCreate()`

2884:    Output Parameter:
2885: .  maxsteps - maximum number of steps to use

2887:    Level: advanced

2889: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2890: @*/
2891: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2892: {
2893:   PetscFunctionBegin;
2896:   *maxsteps = ts->max_steps;
2897:   PetscFunctionReturn(PETSC_SUCCESS);
2898: }

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

2903:    Logically Collective

2905:    Input Parameters:
2906: +  ts - the `TS` context obtained from `TSCreate()`
2907: -  maxtime - final time to step to

2909:    Options Database Key:
2910: .  -ts_max_time <maxtime> - Sets maxtime

2912:    Level: intermediate

2914:    Notes:
2915:    The default maximum time is 5.0

2917: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2918: @*/
2919: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2920: {
2921:   PetscFunctionBegin;
2924:   ts->max_time = maxtime;
2925:   PetscFunctionReturn(PETSC_SUCCESS);
2926: }

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

2931:    Not Collective

2933:    Input Parameter:
2934: .  ts - the `TS` context obtained from `TSCreate()`

2936:    Output Parameter:
2937: .  maxtime - final time to step to

2939:    Level: advanced

2941: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2942: @*/
2943: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2944: {
2945:   PetscFunctionBegin;
2948:   *maxtime = ts->max_time;
2949:   PetscFunctionReturn(PETSC_SUCCESS);
2950: }

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

2955:    Level: deprecated

2957: @*/
2958: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2959: {
2960:   PetscFunctionBegin;
2962:   PetscCall(TSSetTime(ts, initial_time));
2963:   PetscCall(TSSetTimeStep(ts, time_step));
2964:   PetscFunctionReturn(PETSC_SUCCESS);
2965: }

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

2970:    Level: deprecated

2972: @*/
2973: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2974: {
2975:   PetscFunctionBegin;
2977:   if (maxsteps) {
2979:     *maxsteps = ts->max_steps;
2980:   }
2981:   if (maxtime) {
2983:     *maxtime = ts->max_time;
2984:   }
2985:   PetscFunctionReturn(PETSC_SUCCESS);
2986: }

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

2991:    Level: deprecated

2993: @*/
2994: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2995: {
2996:   PetscFunctionBegin;
3000:   if (maxsteps >= 0) ts->max_steps = maxsteps;
3001:   if (maxtime != (PetscReal)PETSC_DEFAULT) ts->max_time = maxtime;
3002:   PetscFunctionReturn(PETSC_SUCCESS);
3003: }

3005: /*@
3006:    TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.

3008:    Level: deprecated

3010: @*/
3011: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
3012: {
3013:   return TSGetStepNumber(ts, steps);
3014: }

3016: /*@
3017:    TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.

3019:    Level: deprecated

3021: @*/
3022: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
3023: {
3024:   return TSGetStepNumber(ts, steps);
3025: }

3027: /*@
3028:    TSSetSolution - Sets the initial solution vector
3029:    for use by the `TS` routines.

3031:    Logically Collective

3033:    Input Parameters:
3034: +  ts - the `TS` context obtained from `TSCreate()`
3035: -  u - the solution vector

3037:    Level: beginner

3039: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
3040: @*/
3041: PetscErrorCode TSSetSolution(TS ts, Vec u)
3042: {
3043:   DM dm;

3045:   PetscFunctionBegin;
3048:   PetscCall(PetscObjectReference((PetscObject)u));
3049:   PetscCall(VecDestroy(&ts->vec_sol));
3050:   ts->vec_sol = u;

3052:   PetscCall(TSGetDM(ts, &dm));
3053:   PetscCall(DMShellSetGlobalVector(dm, u));
3054:   PetscFunctionReturn(PETSC_SUCCESS);
3055: }

3057: /*@C
3058:   TSSetPreStep - Sets the general-purpose function
3059:   called once at the beginning of each time step.

3061:   Logically Collective

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

3067:   Calling sequence of `func`:
3068: .vb
3069:   PetscErrorCode func (TS ts)
3070: .ve

3072:   Level: intermediate

3074: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3075: @*/
3076: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
3077: {
3078:   PetscFunctionBegin;
3080:   ts->prestep = func;
3081:   PetscFunctionReturn(PETSC_SUCCESS);
3082: }

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

3087:   Collective

3089:   Input Parameter:
3090: . ts   - The `TS` context obtained from `TSCreate()`

3092:   Level: developer

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

3098: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3099: @*/
3100: PetscErrorCode TSPreStep(TS ts)
3101: {
3102:   PetscFunctionBegin;
3104:   if (ts->prestep) {
3105:     Vec              U;
3106:     PetscObjectId    idprev;
3107:     PetscBool        sameObject;
3108:     PetscObjectState sprev, spost;

3110:     PetscCall(TSGetSolution(ts, &U));
3111:     PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3112:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3113:     PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3114:     PetscCall(TSGetSolution(ts, &U));
3115:     PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3116:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3117:     if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3118:   }
3119:   PetscFunctionReturn(PETSC_SUCCESS);
3120: }

3122: /*@C
3123:   TSSetPreStage - Sets the general-purpose function
3124:   called once at the beginning of each stage.

3126:   Logically Collective

3128:   Input Parameters:
3129: + ts   - The `TS` context obtained from `TSCreate()`
3130: - func - The function

3132:   Calling sequence of `func`:
3133: .vb
3134:   PetscErrorCode func(TS ts, PetscReal stagetime)
3135: .ve

3137:   Level: intermediate

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

3144: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3145: @*/
3146: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS, PetscReal))
3147: {
3148:   PetscFunctionBegin;
3150:   ts->prestage = func;
3151:   PetscFunctionReturn(PETSC_SUCCESS);
3152: }

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

3158:   Logically Collective

3160:   Input Parameters:
3161: + ts   - The `TS` context obtained from `TSCreate()`
3162: - func - The function

3164:   Calling sequence of `func`:
3165: .vb
3166:   PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y)
3167: .ve

3169:   Level: intermediate

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

3176: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3177: @*/
3178: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS, PetscReal, PetscInt, Vec *))
3179: {
3180:   PetscFunctionBegin;
3182:   ts->poststage = func;
3183:   PetscFunctionReturn(PETSC_SUCCESS);
3184: }

3186: /*@C
3187:   TSSetPostEvaluate - Sets the general-purpose function
3188:   called once at the end of each step evaluation.

3190:   Logically Collective

3192:   Input Parameters:
3193: + ts   - The `TS` context obtained from `TSCreate()`
3194: - func - The function

3196:   Calling sequence of `func`:
3197: .vb
3198:   PetscErrorCode func(TS ts)
3199: .ve

3201:   Level: intermediate

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

3210: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3211: @*/
3212: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS))
3213: {
3214:   PetscFunctionBegin;
3216:   ts->postevaluate = func;
3217:   PetscFunctionReturn(PETSC_SUCCESS);
3218: }

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

3223:   Collective

3225:   Input Parameters:
3226: . ts          - The `TS` context obtained from `TSCreate()`
3227:   stagetime   - The absolute time of the current stage

3229:   Level: developer

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

3235: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3236: @*/
3237: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3238: {
3239:   PetscFunctionBegin;
3241:   if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3242:   PetscFunctionReturn(PETSC_SUCCESS);
3243: }

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

3248:   Collective

3250:   Input Parameters:
3251: . ts          - The `TS` context obtained from `TSCreate()`
3252:   stagetime   - The absolute time of the current stage
3253:   stageindex  - Stage number
3254:   Y           - Array of vectors (of size = total number
3255:                 of stages) with the stage solutions

3257:   Level: developer

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

3263: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3264: @*/
3265: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3266: {
3267:   PetscFunctionBegin;
3269:   if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3270:   PetscFunctionReturn(PETSC_SUCCESS);
3271: }

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

3276:   Collective

3278:   Input Parameter:
3279: . ts - The `TS` context obtained from `TSCreate()`

3281:   Level: developer

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

3287: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3288: @*/
3289: PetscErrorCode TSPostEvaluate(TS ts)
3290: {
3291:   PetscFunctionBegin;
3293:   if (ts->postevaluate) {
3294:     Vec              U;
3295:     PetscObjectState sprev, spost;

3297:     PetscCall(TSGetSolution(ts, &U));
3298:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3299:     PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3300:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3301:     if (sprev != spost) PetscCall(TSRestartStep(ts));
3302:   }
3303:   PetscFunctionReturn(PETSC_SUCCESS);
3304: }

3306: /*@C
3307:   TSSetPostStep - Sets the general-purpose function
3308:   called once at the end of each time step.

3310:   Logically Collective

3312:   Input Parameters:
3313: + ts   - The `TS` context obtained from `TSCreate()`
3314: - func - The function

3316:   Calling sequence of `func`:
3317: $ PetscErrorCode func(TS ts)

3319:   Level: intermediate

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

3326: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3327: @*/
3328: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
3329: {
3330:   PetscFunctionBegin;
3332:   ts->poststep = func;
3333:   PetscFunctionReturn(PETSC_SUCCESS);
3334: }

3336: /*@
3337:   TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`

3339:   Collective

3341:   Input Parameter:
3342: . ts   - The `TS` context obtained from `TSCreate()`

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

3348:   Level: developer

3350: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()`
3351: @*/
3352: PetscErrorCode TSPostStep(TS ts)
3353: {
3354:   PetscFunctionBegin;
3356:   if (ts->poststep) {
3357:     Vec              U;
3358:     PetscObjectId    idprev;
3359:     PetscBool        sameObject;
3360:     PetscObjectState sprev, spost;

3362:     PetscCall(TSGetSolution(ts, &U));
3363:     PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3364:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3365:     PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3366:     PetscCall(TSGetSolution(ts, &U));
3367:     PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3368:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3369:     if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3370:   }
3371:   PetscFunctionReturn(PETSC_SUCCESS);
3372: }

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

3377:    Collective

3379:    Input Parameters:
3380: +  ts - time stepping context
3381: -  t - time to interpolate to

3383:    Output Parameter:
3384: .  U - state at given time

3386:    Level: intermediate

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

3391: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3392: @*/
3393: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3394: {
3395:   PetscFunctionBegin;
3398:   PetscCheck(t >= ts->ptime_prev && t <= ts->ptime, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Requested time %g not in last time steps [%g,%g]", (double)t, (double)ts->ptime_prev, (double)ts->ptime);
3399:   PetscUseTypeMethod(ts, interpolate, t, U);
3400:   PetscFunctionReturn(PETSC_SUCCESS);
3401: }

3403: /*@
3404:    TSStep - Steps one time step

3406:    Collective

3408:    Input Parameter:
3409: .  ts - the `TS` context obtained from `TSCreate()`

3411:    Level: developer

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

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

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

3422: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3423: @*/
3424: PetscErrorCode TSStep(TS ts)
3425: {
3426:   static PetscBool cite = PETSC_FALSE;
3427:   PetscReal        ptime;

3429:   PetscFunctionBegin;
3431:   PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3432:                                    "  title         = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3433:                                    "  author        = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3434:                                    "  journal       = {arXiv e-preprints},\n"
3435:                                    "  eprint        = {1806.01437},\n"
3436:                                    "  archivePrefix = {arXiv},\n"
3437:                                    "  year          = {2018}\n}\n",
3438:                                    &cite));
3439:   PetscCall(TSSetUp(ts));
3440:   PetscCall(TSTrajectorySetUp(ts->trajectory, ts));

3442:   PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_MAX_INT, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
3443:   PetscCheck(ts->exact_final_time != TS_EXACTFINALTIME_UNSPECIFIED, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetExactFinalTime() or use -ts_exact_final_time <stepover,interpolate,matchstep> before calling TSStep()");
3444:   PetscCheck(ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP || ts->adapt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Since TS is not adaptive you cannot use TS_EXACTFINALTIME_MATCHSTEP, suggest TS_EXACTFINALTIME_INTERPOLATE");

3446:   if (!ts->steps) ts->ptime_prev = ts->ptime;
3447:   ptime                   = ts->ptime;
3448:   ts->ptime_prev_rollback = ts->ptime_prev;
3449:   ts->reason              = TS_CONVERGED_ITERATING;

3451:   PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3452:   PetscUseTypeMethod(ts, step);
3453:   PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));

3455:   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)
3456:     PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[ts->tspan->spanctr++]));
3457:   if (ts->reason >= 0) {
3458:     ts->ptime_prev = ptime;
3459:     ts->steps++;
3460:     ts->steprollback = PETSC_FALSE;
3461:     ts->steprestart  = PETSC_FALSE;
3462:   }
3463:   if (!ts->reason) {
3464:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3465:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3466:   }

3468:   if (ts->reason < 0 && ts->errorifstepfailed) {
3469:     PetscCall(TSMonitorCancel(ts));
3470:     PetscCheck(ts->reason != TS_DIVERGED_NONLINEAR_SOLVE, PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s, increase -ts_max_snes_failures or make negative to attempt recovery", TSConvergedReasons[ts->reason]);
3471:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3472:   }
3473:   PetscFunctionReturn(PETSC_SUCCESS);
3474: }

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

3480:    Collective

3482:    Input Parameters:
3483: +  ts - time stepping context
3484: -  wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

3486:    Input/Output Parameter:
3487: .  order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3488:            on output, the actual order of the error evaluation

3490:    Output Parameter:
3491: .  wlte - the weighted local truncation error norm

3493:    Level: advanced

3495:    Note:
3496:    If the timestepper cannot evaluate the error in a particular step
3497:    (eg. in the first step or restart steps after event handling),
3498:    this routine returns wlte=-1.0 .

3500: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3501: @*/
3502: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3503: {
3504:   PetscFunctionBegin;
3511:   PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3512:   PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3513:   PetscFunctionReturn(PETSC_SUCCESS);
3514: }

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

3519:    Collective

3521:    Input Parameters:
3522: +  ts - time stepping context
3523: .  order - desired order of accuracy
3524: -  done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)

3526:    Output Parameter:
3527: .  U - state at the end of the current step

3529:    Level: advanced

3531:    Notes:
3532:    This function cannot be called until all stages have been evaluated.

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

3536: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3537: @*/
3538: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3539: {
3540:   PetscFunctionBegin;
3544:   PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3545:   PetscFunctionReturn(PETSC_SUCCESS);
3546: }

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

3551:   Not collective

3553:   Input Parameter:
3554: . ts - time stepping context

3556:   Output Parameter:
3557: . initConditions - The function which computes an initial condition

3559:   Calling sequence of `initCondition`:
3560: .vb
3561:   PetscErrorCode initCondition(TS ts, Vec u)
3562: .ve
3563: + ts - The timestepping context
3564: - u  - The input vector in which the initial condition is stored

3566:    Level: advanced

3568: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3569: @*/
3570: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS, Vec))
3571: {
3572:   PetscFunctionBegin;
3575:   *initCondition = ts->ops->initcondition;
3576:   PetscFunctionReturn(PETSC_SUCCESS);
3577: }

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

3582:   Logically collective

3584:   Input Parameters:
3585: + ts  - time stepping context
3586: - initCondition - The function which computes an initial condition

3588:   Calling sequence of `initCondition`:
3589: $ PetscErrorCode initCondition(TS ts, Vec u)
3590: + ts - The timestepping context
3591: - u  - The input vector in which the initial condition is to be stored

3593:   Level: advanced

3595: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3596: @*/
3597: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS, Vec))
3598: {
3599:   PetscFunctionBegin;
3602:   ts->ops->initcondition = initCondition;
3603:   PetscFunctionReturn(PETSC_SUCCESS);
3604: }

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

3609:   Collective

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

3615:   Level: advanced

3617: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3618: @*/
3619: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3620: {
3621:   PetscFunctionBegin;
3624:   PetscTryTypeMethod(ts, initcondition, u);
3625:   PetscFunctionReturn(PETSC_SUCCESS);
3626: }

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

3631:   Not collective

3633:   Input Parameter:
3634: . ts - time stepping context

3636:   Output Parameter:
3637: . exactError - The function which computes the solution error

3639:   Calling sequence of `exactError`:
3640: $ PetscErrorCode exactError(TS ts, Vec u, Vec e)
3641: + ts - The timestepping context
3642: . u  - The approximate solution vector
3643: - e  - The vector in which the error is stored

3645:   Level: advanced

3647: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3648: @*/
3649: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS, Vec, Vec))
3650: {
3651:   PetscFunctionBegin;
3654:   *exactError = ts->ops->exacterror;
3655:   PetscFunctionReturn(PETSC_SUCCESS);
3656: }

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

3661:   Logically collective

3663:   Input Parameters:
3664: + ts - time stepping context
3665: - exactError - The function which computes the solution error

3667:   Calling sequence of `exactError`:
3668: $ PetscErrorCode exactError(TS ts, Vec u)
3669: + ts - The timestepping context
3670: . u  - The approximate solution vector
3671: - e  - The  vector in which the error is stored

3673:   Level: advanced

3675: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3676: @*/
3677: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS, Vec, Vec))
3678: {
3679:   PetscFunctionBegin;
3682:   ts->ops->exacterror = exactError;
3683:   PetscFunctionReturn(PETSC_SUCCESS);
3684: }

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

3689:   Collective

3691:   Input Parameters:
3692: + ts - time stepping context
3693: . u  - The approximate solution
3694: - e  - The `Vec` used to store the error

3696:   Level: advanced

3698: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3699: @*/
3700: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3701: {
3702:   PetscFunctionBegin;
3706:   PetscTryTypeMethod(ts, exacterror, u, e);
3707:   PetscFunctionReturn(PETSC_SUCCESS);
3708: }

3710: /*@
3711:    TSSolve - Steps the requested number of timesteps.

3713:    Collective

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

3720:    Level: beginner

3722:    Notes:
3723:    The final time returned by this function may be different from the time of the internally
3724:    held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
3725:    stepped over the final time.

3727: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3728: @*/
3729: PetscErrorCode TSSolve(TS ts, Vec u)
3730: {
3731:   Vec solution;

3733:   PetscFunctionBegin;

3737:   PetscCall(TSSetExactFinalTimeDefault(ts));
3738:   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 */
3739:     if (!ts->vec_sol || u == ts->vec_sol) {
3740:       PetscCall(VecDuplicate(u, &solution));
3741:       PetscCall(TSSetSolution(ts, solution));
3742:       PetscCall(VecDestroy(&solution)); /* grant ownership */
3743:     }
3744:     PetscCall(VecCopy(u, ts->vec_sol));
3745:     PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
3746:   } else if (u) PetscCall(TSSetSolution(ts, u));
3747:   PetscCall(TSSetUp(ts));
3748:   PetscCall(TSTrajectorySetUp(ts->trajectory, ts));

3750:   PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_MAX_INT, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
3751:   PetscCheck(ts->exact_final_time != TS_EXACTFINALTIME_UNSPECIFIED, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetExactFinalTime() or use -ts_exact_final_time <stepover,interpolate,matchstep> before calling TSSolve()");
3752:   PetscCheck(ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP || ts->adapt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Since TS is not adaptive you cannot use TS_EXACTFINALTIME_MATCHSTEP, suggest TS_EXACTFINALTIME_INTERPOLATE");
3753:   PetscCheck(!(ts->tspan && ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP), PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "You must use TS_EXACTFINALTIME_MATCHSTEP when using time span");

3755:   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 */
3756:     PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]));
3757:     ts->tspan->spanctr = 1;
3758:   }

3760:   if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));

3762:   /* reset number of steps only when the step is not restarted. ARKIMEX
3763:      restarts the step after an event. Resetting these counters in such case causes
3764:      TSTrajectory to incorrectly save the output files
3765:   */
3766:   /* reset time step and iteration counters */
3767:   if (!ts->steps) {
3768:     ts->ksp_its           = 0;
3769:     ts->snes_its          = 0;
3770:     ts->num_snes_failures = 0;
3771:     ts->reject            = 0;
3772:     ts->steprestart       = PETSC_TRUE;
3773:     ts->steprollback      = PETSC_FALSE;
3774:     ts->rhsjacobian.time  = PETSC_MIN_REAL;
3775:   }

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

3782:     if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
3783:     else maxdt = ts->max_time - ts->ptime;
3784:     ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
3785:   }
3786:   ts->reason = TS_CONVERGED_ITERATING;

3788:   {
3789:     PetscViewer       viewer;
3790:     PetscViewerFormat format;
3791:     PetscBool         flg;
3792:     static PetscBool  incall = PETSC_FALSE;

3794:     if (!incall) {
3795:       /* Estimate the convergence rate of the time discretization */
3796:       PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
3797:       if (flg) {
3798:         PetscConvEst conv;
3799:         DM           dm;
3800:         PetscReal   *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
3801:         PetscInt     Nf;
3802:         PetscBool    checkTemporal = PETSC_TRUE;

3804:         incall = PETSC_TRUE;
3805:         PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
3806:         PetscCall(TSGetDM(ts, &dm));
3807:         PetscCall(DMGetNumFields(dm, &Nf));
3808:         PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
3809:         PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
3810:         PetscCall(PetscConvEstUseTS(conv, checkTemporal));
3811:         PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
3812:         PetscCall(PetscConvEstSetFromOptions(conv));
3813:         PetscCall(PetscConvEstSetUp(conv));
3814:         PetscCall(PetscConvEstGetConvRate(conv, alpha));
3815:         PetscCall(PetscViewerPushFormat(viewer, format));
3816:         PetscCall(PetscConvEstRateView(conv, alpha, viewer));
3817:         PetscCall(PetscViewerPopFormat(viewer));
3818:         PetscCall(PetscViewerDestroy(&viewer));
3819:         PetscCall(PetscConvEstDestroy(&conv));
3820:         PetscCall(PetscFree(alpha));
3821:         incall = PETSC_FALSE;
3822:       }
3823:     }
3824:   }

3826:   PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));

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

3837:     if (!ts->steps) {
3838:       PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
3839:       PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
3840:     }

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

3868:     if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
3869:       PetscCall(TSInterpolate(ts, ts->max_time, u));
3870:       ts->solvetime = ts->max_time;
3871:       solution      = u;
3872:       PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
3873:     } else {
3874:       if (u) PetscCall(VecCopy(ts->vec_sol, u));
3875:       ts->solvetime = ts->ptime;
3876:       solution      = ts->vec_sol;
3877:     }
3878:   }

3880:   PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
3881:   PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
3882:   PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
3883:   if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
3884:   PetscFunctionReturn(PETSC_SUCCESS);
3885: }

3887: /*@
3888:    TSGetTime - Gets the time of the most recently completed step.

3890:    Not Collective

3892:    Input Parameter:
3893: .  ts - the `TS` context obtained from `TSCreate()`

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

3898:    Level: beginner

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

3904: .seealso: [](ch_ts), TS`, ``TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
3905: @*/
3906: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
3907: {
3908:   PetscFunctionBegin;
3911:   *t = ts->ptime;
3912:   PetscFunctionReturn(PETSC_SUCCESS);
3913: }

3915: /*@
3916:    TSGetPrevTime - Gets the starting time of the previously completed step.

3918:    Not Collective

3920:    Input Parameter:
3921: .  ts - the `TS` context obtained from `TSCreate()`

3923:    Output Parameter:
3924: .  t  - the previous time

3926:    Level: beginner

3928: .seealso: [](ch_ts), TS`, ``TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
3929: @*/
3930: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
3931: {
3932:   PetscFunctionBegin;
3935:   *t = ts->ptime_prev;
3936:   PetscFunctionReturn(PETSC_SUCCESS);
3937: }

3939: /*@
3940:    TSSetTime - Allows one to reset the time.

3942:    Logically Collective

3944:    Input Parameters:
3945: +  ts - the `TS` context obtained from `TSCreate()`
3946: -  time - the time

3948:    Level: intermediate

3950: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
3951: @*/
3952: PetscErrorCode TSSetTime(TS ts, PetscReal t)
3953: {
3954:   PetscFunctionBegin;
3957:   ts->ptime = t;
3958:   PetscFunctionReturn(PETSC_SUCCESS);
3959: }

3961: /*@C
3962:    TSSetOptionsPrefix - Sets the prefix used for searching for all
3963:    TS options in the database.

3965:    Logically Collective

3967:    Input Parameters:
3968: +  ts     - The `TS` context
3969: -  prefix - The prefix to prepend to all option names

3971:    Level: advanced

3973:    Note:
3974:    A hyphen (-) must NOT be given at the beginning of the prefix name.
3975:    The first character of all runtime options is AUTOMATICALLY the
3976:    hyphen.

3978: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
3979: @*/
3980: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
3981: {
3982:   SNES snes;

3984:   PetscFunctionBegin;
3986:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
3987:   PetscCall(TSGetSNES(ts, &snes));
3988:   PetscCall(SNESSetOptionsPrefix(snes, prefix));
3989:   PetscFunctionReturn(PETSC_SUCCESS);
3990: }

3992: /*@C
3993:    TSAppendOptionsPrefix - Appends to the prefix used for searching for all
3994:    TS options in the database.

3996:    Logically Collective

3998:    Input Parameters:
3999: +  ts     - The `TS` context
4000: -  prefix - The prefix to prepend to all option names

4002:    Level: advanced

4004:    Note:
4005:    A hyphen (-) must NOT be given at the beginning of the prefix name.
4006:    The first character of all runtime options is AUTOMATICALLY the
4007:    hyphen.

4009: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4010: @*/
4011: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4012: {
4013:   SNES snes;

4015:   PetscFunctionBegin;
4017:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4018:   PetscCall(TSGetSNES(ts, &snes));
4019:   PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4020:   PetscFunctionReturn(PETSC_SUCCESS);
4021: }

4023: /*@C
4024:    TSGetOptionsPrefix - Sets the prefix used for searching for all
4025:    `TS` options in the database.

4027:    Not Collective

4029:    Input Parameter:
4030: .  ts - The `TS` context

4032:    Output Parameter:
4033: .  prefix - A pointer to the prefix string used

4035:    Level: intermediate

4037:    Fortran Note:
4038:    The user should pass in a string 'prefix' of
4039:    sufficient length to hold the prefix.

4041: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4042: @*/
4043: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4044: {
4045:   PetscFunctionBegin;
4048:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4049:   PetscFunctionReturn(PETSC_SUCCESS);
4050: }

4052: /*@C
4053:    TSGetRHSJacobian - Returns the Jacobian J at the present timestep.

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

4057:    Input Parameter:
4058: .  ts  - The `TS` context obtained from `TSCreate()`

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

4066:    Level: intermediate

4068:    Note:
4069:     You can pass in `NULL` for any return argument you do not need.

4071: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`

4073: @*/
4074: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobian *func, void **ctx)
4075: {
4076:   DM dm;

4078:   PetscFunctionBegin;
4079:   if (Amat || Pmat) {
4080:     SNES snes;
4081:     PetscCall(TSGetSNES(ts, &snes));
4082:     PetscCall(SNESSetUpMatrices(snes));
4083:     PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4084:   }
4085:   PetscCall(TSGetDM(ts, &dm));
4086:   PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4087:   PetscFunctionReturn(PETSC_SUCCESS);
4088: }

4090: /*@C
4091:    TSGetIJacobian - Returns the implicit Jacobian at the present timestep.

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

4095:    Input Parameter:
4096: .  ts  - The `TS` context obtained from `TSCreate()`

4098:    Output Parameters:
4099: +  Amat  - The (approximate) Jacobian of F(t,U,U_t)
4100: .  Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4101: .  f   - The function to compute the matrices
4102: - ctx - User-defined context for Jacobian evaluation routine

4104:    Level: advanced

4106:    Note:
4107:     You can pass in `NULL` for any return argument you do not need.

4109: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4110: @*/
4111: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobian *f, void **ctx)
4112: {
4113:   DM dm;

4115:   PetscFunctionBegin;
4116:   if (Amat || Pmat) {
4117:     SNES snes;
4118:     PetscCall(TSGetSNES(ts, &snes));
4119:     PetscCall(SNESSetUpMatrices(snes));
4120:     PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4121:   }
4122:   PetscCall(TSGetDM(ts, &dm));
4123:   PetscCall(DMTSGetIJacobian(dm, f, ctx));
4124:   PetscFunctionReturn(PETSC_SUCCESS);
4125: }

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

4131:    Logically Collective

4133:    Input Parameters:
4134: +  ts - the `TS` integrator object
4135: -  dm - the dm, cannot be `NULL`

4137:    Level: intermediate

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

4144: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4145: @*/
4146: PetscErrorCode TSSetDM(TS ts, DM dm)
4147: {
4148:   SNES snes;
4149:   DMTS tsdm;

4151:   PetscFunctionBegin;
4154:   PetscCall(PetscObjectReference((PetscObject)dm));
4155:   if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4156:     if (ts->dm->dmts && !dm->dmts) {
4157:       PetscCall(DMCopyDMTS(ts->dm, dm));
4158:       PetscCall(DMGetDMTS(ts->dm, &tsdm));
4159:       /* Grant write privileges to the replacement DM */
4160:       if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4161:     }
4162:     PetscCall(DMDestroy(&ts->dm));
4163:   }
4164:   ts->dm = dm;

4166:   PetscCall(TSGetSNES(ts, &snes));
4167:   PetscCall(SNESSetDM(snes, dm));
4168:   PetscFunctionReturn(PETSC_SUCCESS);
4169: }

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

4174:    Not Collective

4176:    Input Parameter:
4177: . ts - the `TS`

4179:    Output Parameter:
4180: .  dm - the `DM`

4182:    Level: intermediate

4184: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4185: @*/
4186: PetscErrorCode TSGetDM(TS ts, DM *dm)
4187: {
4188:   PetscFunctionBegin;
4190:   if (!ts->dm) {
4191:     PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4192:     if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4193:   }
4194:   *dm = ts->dm;
4195:   PetscFunctionReturn(PETSC_SUCCESS);
4196: }

4198: /*@
4199:    SNESTSFormFunction - Function to evaluate nonlinear residual

4201:    Logically Collective

4203:    Input Parameters:
4204: + snes - nonlinear solver
4205: . U - the current state at which to evaluate the residual
4206: - ctx - user context, must be a TS

4208:    Output Parameter:
4209: . F - the nonlinear residual

4211:    Level: advanced

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

4217: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4218: @*/
4219: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4220: {
4221:   TS ts = (TS)ctx;

4223:   PetscFunctionBegin;
4228:   PetscCall((ts->ops->snesfunction)(snes, U, F, ts));
4229:   PetscFunctionReturn(PETSC_SUCCESS);
4230: }

4232: /*@
4233:    SNESTSFormJacobian - Function to evaluate the Jacobian

4235:    Collective

4237:    Input Parameters:
4238: + snes - nonlinear solver
4239: . U - the current state at which to evaluate the residual
4240: - ctx - user context, must be a `TS`

4242:    Output Parameters:
4243: + A - the Jacobian
4244: - B - the preconditioning matrix (may be the same as A)

4246:    Level: developer

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

4251: .seealso: [](ch_ts), `SNESSetJacobian()`
4252: @*/
4253: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4254: {
4255:   TS ts = (TS)ctx;

4257:   PetscFunctionBegin;
4265:   PetscCall((ts->ops->snesjacobian)(snes, U, A, B, ts));
4266:   PetscFunctionReturn(PETSC_SUCCESS);
4267: }

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

4272:    Collective

4274:    Input Parameters:
4275: +  ts - time stepping context
4276: .  t - time at which to evaluate
4277: .  U - state at which to evaluate
4278: -  ctx - context

4280:    Output Parameter:
4281: .  F - right hand side

4283:    Level: intermediate

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

4289: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4290: @*/
4291: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4292: {
4293:   Mat Arhs, Brhs;

4295:   PetscFunctionBegin;
4296:   PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4297:   /* undo the damage caused by shifting */
4298:   PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4299:   PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4300:   PetscCall(MatMult(Arhs, U, F));
4301:   PetscFunctionReturn(PETSC_SUCCESS);
4302: }

4304: /*@C
4305:    TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.

4307:    Collective

4309:    Input Parameters:
4310: +  ts - time stepping context
4311: .  t - time at which to evaluate
4312: .  U - state at which to evaluate
4313: -  ctx - context

4315:    Output Parameters:
4316: +  A - pointer to operator
4317: -  B - pointer to preconditioning matrix

4319:    Level: intermediate

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

4324: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4325: @*/
4326: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4327: {
4328:   PetscFunctionBegin;
4329:   PetscFunctionReturn(PETSC_SUCCESS);
4330: }

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

4335:    Collective

4337:    Input Parameters:
4338: +  ts - time stepping context
4339: .  t - time at which to evaluate
4340: .  U - state at which to evaluate
4341: .  Udot - time derivative of state vector
4342: -  ctx - context

4344:    Output Parameter:
4345: .  F - left hand side

4347:    Level: intermediate

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

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

4357: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4358: @*/
4359: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4360: {
4361:   Mat A, B;

4363:   PetscFunctionBegin;
4364:   PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4365:   PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4366:   PetscCall(MatMult(A, Udot, F));
4367:   PetscFunctionReturn(PETSC_SUCCESS);
4368: }

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

4373:    Collective

4375:    Input Parameters:
4376: +  ts - time stepping context
4377: .  t - time at which to evaluate
4378: .  U - state at which to evaluate
4379: .  Udot - time derivative of state vector
4380: .  shift - shift to apply
4381: -  ctx - context

4383:    Output Parameters:
4384: +  A - pointer to operator
4385: -  B - pointer to preconditioning matrix

4387:    Level: advanced

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

4392:    It is only appropriate for problems of the form

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

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

4400: $    shift*M + J

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

4405: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4406: @*/
4407: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4408: {
4409:   PetscFunctionBegin;
4410:   PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4411:   ts->ijacobian.shift = shift;
4412:   PetscFunctionReturn(PETSC_SUCCESS);
4413: }

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

4418:    Not Collective

4420:    Input Parameter:
4421: .  ts - the `TS` context

4423:    Output Parameter:
4424: .  equation_type - see `TSEquationType`

4426:    Level: beginner

4428: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4429: @*/
4430: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4431: {
4432:   PetscFunctionBegin;
4435:   *equation_type = ts->equation_type;
4436:   PetscFunctionReturn(PETSC_SUCCESS);
4437: }

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

4442:    Not Collective

4444:    Input Parameters:
4445: +  ts - the `TS` context
4446: -  equation_type - see `TSEquationType`

4448:    Level: advanced

4450: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4451: @*/
4452: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4453: {
4454:   PetscFunctionBegin;
4456:   ts->equation_type = equation_type;
4457:   PetscFunctionReturn(PETSC_SUCCESS);
4458: }

4460: /*@
4461:    TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.

4463:    Not Collective

4465:    Input Parameter:
4466: .  ts - the `TS` context

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

4472:    Level: beginner

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

4477: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4478: @*/
4479: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4480: {
4481:   PetscFunctionBegin;
4484:   *reason = ts->reason;
4485:   PetscFunctionReturn(PETSC_SUCCESS);
4486: }

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

4491:    Logically Collective; reason must contain common value

4493:    Input Parameters:
4494: +  ts - the `TS` context
4495: -  reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4496:             manual pages for the individual convergence tests for complete lists

4498:    Level: advanced

4500:    Note:
4501:    Can only be called while `TSSolve()` is active.

4503: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4504: @*/
4505: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4506: {
4507:   PetscFunctionBegin;
4509:   ts->reason = reason;
4510:   PetscFunctionReturn(PETSC_SUCCESS);
4511: }

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

4516:    Not Collective

4518:    Input Parameter:
4519: .  ts - the `TS` context

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

4524:    Level: beginner

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

4529: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4530: @*/
4531: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4532: {
4533:   PetscFunctionBegin;
4536:   *ftime = ts->solvetime;
4537:   PetscFunctionReturn(PETSC_SUCCESS);
4538: }

4540: /*@
4541:    TSGetSNESIterations - Gets the total number of nonlinear iterations
4542:    used by the time integrator.

4544:    Not Collective

4546:    Input Parameter:
4547: .  ts - `TS` context

4549:    Output Parameter:
4550: .  nits - number of nonlinear iterations

4552:    Level: intermediate

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

4557: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4558: @*/
4559: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4560: {
4561:   PetscFunctionBegin;
4564:   *nits = ts->snes_its;
4565:   PetscFunctionReturn(PETSC_SUCCESS);
4566: }

4568: /*@
4569:    TSGetKSPIterations - Gets the total number of linear iterations
4570:    used by the time integrator.

4572:    Not Collective

4574:    Input Parameter:
4575: .  ts - `TS` context

4577:    Output Parameter:
4578: .  lits - number of linear iterations

4580:    Level: intermediate

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

4585: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()`
4586: @*/
4587: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4588: {
4589:   PetscFunctionBegin;
4592:   *lits = ts->ksp_its;
4593:   PetscFunctionReturn(PETSC_SUCCESS);
4594: }

4596: /*@
4597:    TSGetStepRejections - Gets the total number of rejected steps.

4599:    Not Collective

4601:    Input Parameter:
4602: .  ts - `TS` context

4604:    Output Parameter:
4605: .  rejects - number of steps rejected

4607:    Level: intermediate

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

4612: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4613: @*/
4614: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4615: {
4616:   PetscFunctionBegin;
4619:   *rejects = ts->reject;
4620:   PetscFunctionReturn(PETSC_SUCCESS);
4621: }

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

4626:    Not Collective

4628:    Input Parameter:
4629: .  ts - `TS` context

4631:    Output Parameter:
4632: .  fails - number of failed nonlinear solves

4634:    Level: intermediate

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

4639: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4640: @*/
4641: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4642: {
4643:   PetscFunctionBegin;
4646:   *fails = ts->num_snes_failures;
4647:   PetscFunctionReturn(PETSC_SUCCESS);
4648: }

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

4653:    Not Collective

4655:    Input Parameters:
4656: +  ts - `TS` context
4657: -  rejects - maximum number of rejected steps, pass -1 for unlimited

4659:    Options Database Key:
4660: .  -ts_max_reject - Maximum number of step rejections before a step fails

4662:    Level: intermediate

4664: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4665: @*/
4666: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4667: {
4668:   PetscFunctionBegin;
4670:   ts->max_reject = rejects;
4671:   PetscFunctionReturn(PETSC_SUCCESS);
4672: }

4674: /*@
4675:    TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves

4677:    Not Collective

4679:    Input Parameters:
4680: +  ts - `TS` context
4681: -  fails - maximum number of failed nonlinear solves, pass -1 for unlimited

4683:    Options Database Key:
4684: .  -ts_max_snes_failures - Maximum number of nonlinear solve failures

4686:    Level: intermediate

4688: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4689: @*/
4690: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4691: {
4692:   PetscFunctionBegin;
4694:   ts->max_snes_failures = fails;
4695:   PetscFunctionReturn(PETSC_SUCCESS);
4696: }

4698: /*@
4699:    TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`

4701:    Not Collective

4703:    Input Parameters:
4704: +  ts - `TS` context
4705: -  err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure

4707:    Options Database Key:
4708: .  -ts_error_if_step_fails - Error if no step succeeds

4710:    Level: intermediate

4712: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4713: @*/
4714: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
4715: {
4716:   PetscFunctionBegin;
4718:   ts->errorifstepfailed = err;
4719:   PetscFunctionReturn(PETSC_SUCCESS);
4720: }

4722: /*@
4723:    TSGetAdapt - Get the adaptive controller context for the current method

4725:    Collective on `ts` if controller has not been created yet

4727:    Input Parameter:
4728: .  ts - time stepping context

4730:    Output Parameter:
4731: .  adapt - adaptive controller

4733:    Level: intermediate

4735: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
4736: @*/
4737: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
4738: {
4739:   PetscFunctionBegin;
4742:   if (!ts->adapt) {
4743:     PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
4744:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
4745:   }
4746:   *adapt = ts->adapt;
4747:   PetscFunctionReturn(PETSC_SUCCESS);
4748: }

4750: /*@
4751:    TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller

4753:    Logically Collective

4755:    Input Parameters:
4756: +  ts - time integration context
4757: .  atol - scalar absolute tolerances, `PETSC_DECIDE` to leave current value
4758: .  vatol - vector of absolute tolerances or `NULL`, used in preference to atol if present
4759: .  rtol - scalar relative tolerances, `PETSC_DECIDE` to leave current value
4760: -  vrtol - vector of relative tolerances or `NULL`, used in preference to atol if present

4762:    Options Database Keys:
4763: +  -ts_rtol <rtol> - relative tolerance for local truncation error
4764: -  -ts_atol <atol> - Absolute tolerance for local truncation error

4766:    Level: beginner

4768:    Notes:
4769:    With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
4770:    (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
4771:    computed only for the differential or the algebraic part then this can be done using the vector of
4772:    tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
4773:    differential part and infinity for the algebraic part, the LTE calculation will include only the
4774:    differential variables.

4776: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
4777: @*/
4778: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
4779: {
4780:   PetscFunctionBegin;
4781:   if (atol != (PetscReal)PETSC_DECIDE && atol != (PetscReal)PETSC_DEFAULT) ts->atol = atol;
4782:   if (vatol) {
4783:     PetscCall(PetscObjectReference((PetscObject)vatol));
4784:     PetscCall(VecDestroy(&ts->vatol));
4785:     ts->vatol = vatol;
4786:   }
4787:   if (rtol != (PetscReal)PETSC_DECIDE && rtol != (PetscReal)PETSC_DEFAULT) ts->rtol = rtol;
4788:   if (vrtol) {
4789:     PetscCall(PetscObjectReference((PetscObject)vrtol));
4790:     PetscCall(VecDestroy(&ts->vrtol));
4791:     ts->vrtol = vrtol;
4792:   }
4793:   PetscFunctionReturn(PETSC_SUCCESS);
4794: }

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

4799:    Logically Collective

4801:    Input Parameter:
4802: .  ts - time integration context

4804:    Output Parameters:
4805: +  atol - scalar absolute tolerances, `NULL` to ignore
4806: .  vatol - vector of absolute tolerances, `NULL` to ignore
4807: .  rtol - scalar relative tolerances, `NULL` to ignore
4808: -  vrtol - vector of relative tolerances, `NULL` to ignore

4810:    Level: beginner

4812: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
4813: @*/
4814: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
4815: {
4816:   PetscFunctionBegin;
4817:   if (atol) *atol = ts->atol;
4818:   if (vatol) *vatol = ts->vatol;
4819:   if (rtol) *rtol = ts->rtol;
4820:   if (vrtol) *vrtol = ts->vrtol;
4821:   PetscFunctionReturn(PETSC_SUCCESS);
4822: }

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

4827:    Collective

4829:    Input Parameters:
4830: +  ts - time stepping context
4831: .  U - state vector, usually ts->vec_sol
4832: -  Y - state vector to be compared to U

4834:    Output Parameters:
4835: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
4836: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
4837: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

4839:    Level: developer

4841: .seealso: [](ch_ts), `TS`, `TSErrorWeightedNorm()`, `TSErrorWeightedNormInfinity()`
4842: @*/
4843: PetscErrorCode TSErrorWeightedNorm2(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4844: {
4845:   PetscInt           i, n, N, rstart;
4846:   PetscInt           n_loc, na_loc, nr_loc;
4847:   PetscReal          n_glb, na_glb, nr_glb;
4848:   const PetscScalar *u, *y;
4849:   PetscReal          sum, suma, sumr, gsum, gsuma, gsumr, diff;
4850:   PetscReal          tol, tola, tolr;
4851:   PetscReal          err_loc[6], err_glb[6];

4853:   PetscFunctionBegin;
4859:   PetscCheckSameComm(U, 2, Y, 3);
4863:   PetscCheck(U != Y, PetscObjectComm((PetscObject)U), PETSC_ERR_ARG_IDN, "U and Y cannot be the same vector");

4865:   PetscCall(VecGetSize(U, &N));
4866:   PetscCall(VecGetLocalSize(U, &n));
4867:   PetscCall(VecGetOwnershipRange(U, &rstart, NULL));
4868:   PetscCall(VecGetArrayRead(U, &u));
4869:   PetscCall(VecGetArrayRead(Y, &y));
4870:   sum    = 0.;
4871:   n_loc  = 0;
4872:   suma   = 0.;
4873:   na_loc = 0;
4874:   sumr   = 0.;
4875:   nr_loc = 0;
4876:   if (ts->vatol && ts->vrtol) {
4877:     const PetscScalar *atol, *rtol;
4878:     PetscCall(VecGetArrayRead(ts->vatol, &atol));
4879:     PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
4880:     for (i = 0; i < n; i++) {
4881:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4882:       diff = PetscAbsScalar(y[i] - u[i]);
4883:       tola = PetscRealPart(atol[i]);
4884:       if (tola > 0.) {
4885:         suma += PetscSqr(diff / tola);
4886:         na_loc++;
4887:       }
4888:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4889:       if (tolr > 0.) {
4890:         sumr += PetscSqr(diff / tolr);
4891:         nr_loc++;
4892:       }
4893:       tol = tola + tolr;
4894:       if (tol > 0.) {
4895:         sum += PetscSqr(diff / tol);
4896:         n_loc++;
4897:       }
4898:     }
4899:     PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
4900:     PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
4901:   } else if (ts->vatol) { /* vector atol, scalar rtol */
4902:     const PetscScalar *atol;
4903:     PetscCall(VecGetArrayRead(ts->vatol, &atol));
4904:     for (i = 0; i < n; i++) {
4905:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4906:       diff = PetscAbsScalar(y[i] - u[i]);
4907:       tola = PetscRealPart(atol[i]);
4908:       if (tola > 0.) {
4909:         suma += PetscSqr(diff / tola);
4910:         na_loc++;
4911:       }
4912:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4913:       if (tolr > 0.) {
4914:         sumr += PetscSqr(diff / tolr);
4915:         nr_loc++;
4916:       }
4917:       tol = tola + tolr;
4918:       if (tol > 0.) {
4919:         sum += PetscSqr(diff / tol);
4920:         n_loc++;
4921:       }
4922:     }
4923:     PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
4924:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
4925:     const PetscScalar *rtol;
4926:     PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
4927:     for (i = 0; i < n; i++) {
4928:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4929:       diff = PetscAbsScalar(y[i] - u[i]);
4930:       tola = ts->atol;
4931:       if (tola > 0.) {
4932:         suma += PetscSqr(diff / tola);
4933:         na_loc++;
4934:       }
4935:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4936:       if (tolr > 0.) {
4937:         sumr += PetscSqr(diff / tolr);
4938:         nr_loc++;
4939:       }
4940:       tol = tola + tolr;
4941:       if (tol > 0.) {
4942:         sum += PetscSqr(diff / tol);
4943:         n_loc++;
4944:       }
4945:     }
4946:     PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
4947:   } else { /* scalar atol, scalar rtol */
4948:     for (i = 0; i < n; i++) {
4949:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4950:       diff = PetscAbsScalar(y[i] - u[i]);
4951:       tola = ts->atol;
4952:       if (tola > 0.) {
4953:         suma += PetscSqr(diff / tola);
4954:         na_loc++;
4955:       }
4956:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4957:       if (tolr > 0.) {
4958:         sumr += PetscSqr(diff / tolr);
4959:         nr_loc++;
4960:       }
4961:       tol = tola + tolr;
4962:       if (tol > 0.) {
4963:         sum += PetscSqr(diff / tol);
4964:         n_loc++;
4965:       }
4966:     }
4967:   }
4968:   PetscCall(VecRestoreArrayRead(U, &u));
4969:   PetscCall(VecRestoreArrayRead(Y, &y));

4971:   err_loc[0] = sum;
4972:   err_loc[1] = suma;
4973:   err_loc[2] = sumr;
4974:   err_loc[3] = (PetscReal)n_loc;
4975:   err_loc[4] = (PetscReal)na_loc;
4976:   err_loc[5] = (PetscReal)nr_loc;

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

4980:   gsum   = err_glb[0];
4981:   gsuma  = err_glb[1];
4982:   gsumr  = err_glb[2];
4983:   n_glb  = err_glb[3];
4984:   na_glb = err_glb[4];
4985:   nr_glb = err_glb[5];

4987:   *norm = 0.;
4988:   if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
4989:   *norma = 0.;
4990:   if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
4991:   *normr = 0.;
4992:   if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);

4994:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
4995:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
4996:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
4997:   PetscFunctionReturn(PETSC_SUCCESS);
4998: }

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

5003:    Collective

5005:    Input Parameters:
5006: +  ts - time stepping context
5007: .  U - state vector, usually ts->vec_sol
5008: -  Y - state vector to be compared to U

5010:    Output Parameters:
5011: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5012: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5013: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5015:    Level: developer

5017: .seealso: [](ch_ts), `TS`, `TSErrorWeightedNorm()`, `TSErrorWeightedNorm2()`
5018: @*/
5019: PetscErrorCode TSErrorWeightedNormInfinity(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5020: {
5021:   PetscInt           i, n, N, rstart;
5022:   const PetscScalar *u, *y;
5023:   PetscReal          max, gmax, maxa, gmaxa, maxr, gmaxr;
5024:   PetscReal          tol, tola, tolr, diff;
5025:   PetscReal          err_loc[3], err_glb[3];

5027:   PetscFunctionBegin;
5033:   PetscCheckSameComm(U, 2, Y, 3);
5037:   PetscCheck(U != Y, PetscObjectComm((PetscObject)U), PETSC_ERR_ARG_IDN, "U and Y cannot be the same vector");

5039:   PetscCall(VecGetSize(U, &N));
5040:   PetscCall(VecGetLocalSize(U, &n));
5041:   PetscCall(VecGetOwnershipRange(U, &rstart, NULL));
5042:   PetscCall(VecGetArrayRead(U, &u));
5043:   PetscCall(VecGetArrayRead(Y, &y));

5045:   max  = 0.;
5046:   maxa = 0.;
5047:   maxr = 0.;

5049:   if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
5050:     const PetscScalar *atol, *rtol;
5051:     PetscCall(VecGetArrayRead(ts->vatol, &atol));
5052:     PetscCall(VecGetArrayRead(ts->vrtol, &rtol));

5054:     for (i = 0; i < n; i++) {
5055:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5056:       diff = PetscAbsScalar(y[i] - u[i]);
5057:       tola = PetscRealPart(atol[i]);
5058:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5059:       tol  = tola + tolr;
5060:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
5061:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
5062:       if (tol > 0.) max = PetscMax(max, diff / tol);
5063:     }
5064:     PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5065:     PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5066:   } else if (ts->vatol) { /* vector atol, scalar rtol */
5067:     const PetscScalar *atol;
5068:     PetscCall(VecGetArrayRead(ts->vatol, &atol));
5069:     for (i = 0; i < n; i++) {
5070:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5071:       diff = PetscAbsScalar(y[i] - u[i]);
5072:       tola = PetscRealPart(atol[i]);
5073:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5074:       tol  = tola + tolr;
5075:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
5076:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
5077:       if (tol > 0.) max = PetscMax(max, diff / tol);
5078:     }
5079:     PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5080:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
5081:     const PetscScalar *rtol;
5082:     PetscCall(VecGetArrayRead(ts->vrtol, &rtol));

5084:     for (i = 0; i < n; i++) {
5085:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5086:       diff = PetscAbsScalar(y[i] - u[i]);
5087:       tola = ts->atol;
5088:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5089:       tol  = tola + tolr;
5090:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
5091:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
5092:       if (tol > 0.) max = PetscMax(max, diff / tol);
5093:     }
5094:     PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5095:   } else { /* scalar atol, scalar rtol */

5097:     for (i = 0; i < n; i++) {
5098:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5099:       diff = PetscAbsScalar(y[i] - u[i]);
5100:       tola = ts->atol;
5101:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5102:       tol  = tola + tolr;
5103:       if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
5104:       if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
5105:       if (tol > 0.) max = PetscMax(max, diff / tol);
5106:     }
5107:   }
5108:   PetscCall(VecRestoreArrayRead(U, &u));
5109:   PetscCall(VecRestoreArrayRead(Y, &y));
5110:   err_loc[0] = max;
5111:   err_loc[1] = maxa;
5112:   err_loc[2] = maxr;
5113:   PetscCall(MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
5114:   gmax  = err_glb[0];
5115:   gmaxa = err_glb[1];
5116:   gmaxr = err_glb[2];

5118:   *norm  = gmax;
5119:   *norma = gmaxa;
5120:   *normr = gmaxr;
5121:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5122:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5123:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5124:   PetscFunctionReturn(PETSC_SUCCESS);
5125: }

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

5130:    Collective

5132:    Input Parameters:
5133: +  ts - time stepping context
5134: .  U - state vector, usually ts->vec_sol
5135: .  Y - state vector to be compared to U
5136: -  wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

5138:    Output Parameters:
5139: +  norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5140: .  norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5141: -  normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

5143:    Options Database Key:
5144: .  -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5146:    Level: developer

5148: .seealso: [](ch_ts), `TS`, `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`, `TSErrorWeightedENorm`
5149: @*/
5150: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5151: {
5152:   PetscFunctionBegin;
5153:   if (wnormtype == NORM_2) PetscCall(TSErrorWeightedNorm2(ts, U, Y, norm, norma, normr));
5154:   else if (wnormtype == NORM_INFINITY) PetscCall(TSErrorWeightedNormInfinity(ts, U, Y, norm, norma, normr));
5155:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5156:   PetscFunctionReturn(PETSC_SUCCESS);
5157: }

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

5162:    Collective

5164:    Input Parameters:
5165: +  ts - time stepping context
5166: .  E - error vector
5167: .  U - state vector, usually ts->vec_sol
5168: -  Y - state vector, previous time step

5170:    Output Parameters:
5171: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5172: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5173: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5175:    Level: developer

5177: .seealso: [](ch_ts), `TS`, `TSErrorWeightedENorm()`, `TSErrorWeightedENormInfinity()`
5178: @*/
5179: PetscErrorCode TSErrorWeightedENorm2(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5180: {
5181:   PetscInt           i, n, N, rstart;
5182:   PetscInt           n_loc, na_loc, nr_loc;
5183:   PetscReal          n_glb, na_glb, nr_glb;
5184:   const PetscScalar *e, *u, *y;
5185:   PetscReal          err, sum, suma, sumr, gsum, gsuma, gsumr;
5186:   PetscReal          tol, tola, tolr;
5187:   PetscReal          err_loc[6], err_glb[6];

5189:   PetscFunctionBegin;
5197:   PetscCheckSameComm(E, 2, U, 3);
5198:   PetscCheckSameComm(U, 3, Y, 4);

5203:   PetscCall(VecGetSize(E, &N));
5204:   PetscCall(VecGetLocalSize(E, &n));
5205:   PetscCall(VecGetOwnershipRange(E, &rstart, NULL));
5206:   PetscCall(VecGetArrayRead(E, &e));
5207:   PetscCall(VecGetArrayRead(U, &u));
5208:   PetscCall(VecGetArrayRead(Y, &y));
5209:   sum    = 0.;
5210:   n_loc  = 0;
5211:   suma   = 0.;
5212:   na_loc = 0;
5213:   sumr   = 0.;
5214:   nr_loc = 0;
5215:   if (ts->vatol && ts->vrtol) {
5216:     const PetscScalar *atol, *rtol;
5217:     PetscCall(VecGetArrayRead(ts->vatol, &atol));
5218:     PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
5219:     for (i = 0; i < n; i++) {
5220:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5221:       err  = PetscAbsScalar(e[i]);
5222:       tola = PetscRealPart(atol[i]);
5223:       if (tola > 0.) {
5224:         suma += PetscSqr(err / tola);
5225:         na_loc++;
5226:       }
5227:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5228:       if (tolr > 0.) {
5229:         sumr += PetscSqr(err / tolr);
5230:         nr_loc++;
5231:       }
5232:       tol = tola + tolr;
5233:       if (tol > 0.) {
5234:         sum += PetscSqr(err / tol);
5235:         n_loc++;
5236:       }
5237:     }
5238:     PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5239:     PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5240:   } else if (ts->vatol) { /* vector atol, scalar rtol */
5241:     const PetscScalar *atol;
5242:     PetscCall(VecGetArrayRead(ts->vatol, &atol));
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:       if (tola > 0.) {
5248:         suma += PetscSqr(err / tola);
5249:         na_loc++;
5250:       }
5251:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5252:       if (tolr > 0.) {
5253:         sumr += PetscSqr(err / tolr);
5254:         nr_loc++;
5255:       }
5256:       tol = tola + tolr;
5257:       if (tol > 0.) {
5258:         sum += PetscSqr(err / tol);
5259:         n_loc++;
5260:       }
5261:     }
5262:     PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5263:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
5264:     const PetscScalar *rtol;
5265:     PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
5266:     for (i = 0; i < n; i++) {
5267:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5268:       err  = PetscAbsScalar(e[i]);
5269:       tola = ts->atol;
5270:       if (tola > 0.) {
5271:         suma += PetscSqr(err / tola);
5272:         na_loc++;
5273:       }
5274:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5275:       if (tolr > 0.) {
5276:         sumr += PetscSqr(err / tolr);
5277:         nr_loc++;
5278:       }
5279:       tol = tola + tolr;
5280:       if (tol > 0.) {
5281:         sum += PetscSqr(err / tol);
5282:         n_loc++;
5283:       }
5284:     }
5285:     PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5286:   } else { /* scalar atol, scalar rtol */
5287:     for (i = 0; i < n; i++) {
5288:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5289:       err  = PetscAbsScalar(e[i]);
5290:       tola = ts->atol;
5291:       if (tola > 0.) {
5292:         suma += PetscSqr(err / tola);
5293:         na_loc++;
5294:       }
5295:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5296:       if (tolr > 0.) {
5297:         sumr += PetscSqr(err / tolr);
5298:         nr_loc++;
5299:       }
5300:       tol = tola + tolr;
5301:       if (tol > 0.) {
5302:         sum += PetscSqr(err / tol);
5303:         n_loc++;
5304:       }
5305:     }
5306:   }
5307:   PetscCall(VecRestoreArrayRead(E, &e));
5308:   PetscCall(VecRestoreArrayRead(U, &u));
5309:   PetscCall(VecRestoreArrayRead(Y, &y));

5311:   err_loc[0] = sum;
5312:   err_loc[1] = suma;
5313:   err_loc[2] = sumr;
5314:   err_loc[3] = (PetscReal)n_loc;
5315:   err_loc[4] = (PetscReal)na_loc;
5316:   err_loc[5] = (PetscReal)nr_loc;

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

5320:   gsum   = err_glb[0];
5321:   gsuma  = err_glb[1];
5322:   gsumr  = err_glb[2];
5323:   n_glb  = err_glb[3];
5324:   na_glb = err_glb[4];
5325:   nr_glb = err_glb[5];

5327:   *norm = 0.;
5328:   if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
5329:   *norma = 0.;
5330:   if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
5331:   *normr = 0.;
5332:   if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);

5334:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5335:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5336:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5337:   PetscFunctionReturn(PETSC_SUCCESS);
5338: }

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

5343:    Collective

5345:    Input Parameters:
5346: +  ts - time stepping context
5347: .  E - error vector
5348: .  U - state vector, usually ts->vec_sol
5349: -  Y - state vector, previous time step

5351:    Output Parameters:
5352: +  norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5353: .  norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5354: -  normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances

5356:    Level: developer

5358: .seealso: [](ch_ts), `TS`, `TSErrorWeightedENorm()`, `TSErrorWeightedENorm2()`
5359: @*/
5360: PetscErrorCode TSErrorWeightedENormInfinity(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5361: {
5362:   PetscInt           i, n, N, rstart;
5363:   const PetscScalar *e, *u, *y;
5364:   PetscReal          err, max, gmax, maxa, gmaxa, maxr, gmaxr;
5365:   PetscReal          tol, tola, tolr;
5366:   PetscReal          err_loc[3], err_glb[3];

5368:   PetscFunctionBegin;
5376:   PetscCheckSameComm(E, 2, U, 3);
5377:   PetscCheckSameComm(U, 3, Y, 4);

5382:   PetscCall(VecGetSize(E, &N));
5383:   PetscCall(VecGetLocalSize(E, &n));
5384:   PetscCall(VecGetOwnershipRange(E, &rstart, NULL));
5385:   PetscCall(VecGetArrayRead(E, &e));
5386:   PetscCall(VecGetArrayRead(U, &u));
5387:   PetscCall(VecGetArrayRead(Y, &y));

5389:   max  = 0.;
5390:   maxa = 0.;
5391:   maxr = 0.;

5393:   if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
5394:     const PetscScalar *atol, *rtol;
5395:     PetscCall(VecGetArrayRead(ts->vatol, &atol));
5396:     PetscCall(VecGetArrayRead(ts->vrtol, &rtol));

5398:     for (i = 0; i < n; i++) {
5399:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5400:       err  = PetscAbsScalar(e[i]);
5401:       tola = PetscRealPart(atol[i]);
5402:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5403:       tol  = tola + tolr;
5404:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5405:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5406:       if (tol > 0.) max = PetscMax(max, err / tol);
5407:     }
5408:     PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5409:     PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5410:   } else if (ts->vatol) { /* vector atol, scalar rtol */
5411:     const PetscScalar *atol;
5412:     PetscCall(VecGetArrayRead(ts->vatol, &atol));
5413:     for (i = 0; i < n; i++) {
5414:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5415:       err  = PetscAbsScalar(e[i]);
5416:       tola = PetscRealPart(atol[i]);
5417:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5418:       tol  = tola + tolr;
5419:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5420:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5421:       if (tol > 0.) max = PetscMax(max, err / tol);
5422:     }
5423:     PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5424:   } else if (ts->vrtol) { /* scalar atol, vector rtol */
5425:     const PetscScalar *rtol;
5426:     PetscCall(VecGetArrayRead(ts->vrtol, &rtol));

5428:     for (i = 0; i < n; i++) {
5429:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5430:       err  = PetscAbsScalar(e[i]);
5431:       tola = ts->atol;
5432:       tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5433:       tol  = tola + tolr;
5434:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5435:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5436:       if (tol > 0.) max = PetscMax(max, err / tol);
5437:     }
5438:     PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5439:   } else { /* scalar atol, scalar rtol */

5441:     for (i = 0; i < n; i++) {
5442:       SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5443:       err  = PetscAbsScalar(e[i]);
5444:       tola = ts->atol;
5445:       tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5446:       tol  = tola + tolr;
5447:       if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5448:       if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5449:       if (tol > 0.) max = PetscMax(max, err / tol);
5450:     }
5451:   }
5452:   PetscCall(VecRestoreArrayRead(E, &e));
5453:   PetscCall(VecRestoreArrayRead(U, &u));
5454:   PetscCall(VecRestoreArrayRead(Y, &y));
5455:   err_loc[0] = max;
5456:   err_loc[1] = maxa;
5457:   err_loc[2] = maxr;
5458:   PetscCall(MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
5459:   gmax  = err_glb[0];
5460:   gmaxa = err_glb[1];
5461:   gmaxr = err_glb[2];

5463:   *norm  = gmax;
5464:   *norma = gmaxa;
5465:   *normr = gmaxr;
5466:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5467:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5468:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5469:   PetscFunctionReturn(PETSC_SUCCESS);
5470: }

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

5475:    Collective

5477:    Input Parameters:
5478: +  ts - time stepping context
5479: .  E - error vector
5480: .  U - state vector, usually ts->vec_sol
5481: .  Y - state vector, previous time step
5482: -  wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

5484:    Output Parameters:
5485: +  norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5486: .  norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5487: -  normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

5489:    Options Database Key:
5490: .  -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5492:    Level: developer

5494: .seealso: [](ch_ts), `TS`, `TSErrorWeightedENormInfinity()`, `TSErrorWeightedENorm2()`, `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`
5495: @*/
5496: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5497: {
5498:   PetscFunctionBegin;
5499:   if (wnormtype == NORM_2) PetscCall(TSErrorWeightedENorm2(ts, E, U, Y, norm, norma, normr));
5500:   else if (wnormtype == NORM_INFINITY) PetscCall(TSErrorWeightedENormInfinity(ts, E, U, Y, norm, norma, normr));
5501:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5502:   PetscFunctionReturn(PETSC_SUCCESS);
5503: }

5505: /*@
5506:    TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler

5508:    Logically Collective

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

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

5517:    Level: intermediate

5519: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5520: @*/
5521: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5522: {
5523:   PetscFunctionBegin;
5525:   ts->cfltime_local = cfltime;
5526:   ts->cfltime       = -1.;
5527:   PetscFunctionReturn(PETSC_SUCCESS);
5528: }

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

5533:    Collective

5535:    Input Parameter:
5536: .  ts - time stepping context

5538:    Output Parameter:
5539: .  cfltime - maximum stable time step for forward Euler

5541:    Level: advanced

5543: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5544: @*/
5545: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5546: {
5547:   PetscFunctionBegin;
5548:   if (ts->cfltime < 0) PetscCall(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5549:   *cfltime = ts->cfltime;
5550:   PetscFunctionReturn(PETSC_SUCCESS);
5551: }

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

5556:    Input Parameters:
5557: +  ts   - the `TS` context.
5558: .  xl   - lower bound.
5559: -  xu   - upper bound.

5561:    Level: advanced

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

5567: .seealso: [](ch_ts), `TS`
5568: @*/
5569: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5570: {
5571:   SNES snes;

5573:   PetscFunctionBegin;
5574:   PetscCall(TSGetSNES(ts, &snes));
5575:   PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5576:   PetscFunctionReturn(PETSC_SUCCESS);
5577: }

5579: /*@
5580:    TSComputeLinearStability - computes the linear stability function at a point

5582:    Collective

5584:    Input Parameters:
5585: +  ts - the `TS` context
5586: .  xr - real part of input argument
5587: -  xi - imaginary part of input argument

5589:    Output Parameters:
5590: +  yr - real part of function value
5591: -  yi - imaginary part of function value

5593:    Level: developer

5595: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5596: @*/
5597: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5598: {
5599:   PetscFunctionBegin;
5601:   PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5602:   PetscFunctionReturn(PETSC_SUCCESS);
5603: }

5605: /*@
5606:    TSRestartStep - Flags the solver to restart the next step

5608:    Collective

5610:    Input Parameter:
5611: .  ts - the `TS` context obtained from `TSCreate()`

5613:    Level: advanced

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

5623: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5624: @*/
5625: PetscErrorCode TSRestartStep(TS ts)
5626: {
5627:   PetscFunctionBegin;
5629:   ts->steprestart = PETSC_TRUE;
5630:   PetscFunctionReturn(PETSC_SUCCESS);
5631: }

5633: /*@
5634:    TSRollBack - Rolls back one time step

5636:    Collective

5638:    Input Parameter:
5639: .  ts - the `TS` context obtained from `TSCreate()`

5641:    Level: advanced

5643: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5644: @*/
5645: PetscErrorCode TSRollBack(TS ts)
5646: {
5647:   PetscFunctionBegin;
5649:   PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5650:   PetscUseTypeMethod(ts, rollback);
5651:   ts->time_step  = ts->ptime - ts->ptime_prev;
5652:   ts->ptime      = ts->ptime_prev;
5653:   ts->ptime_prev = ts->ptime_prev_rollback;
5654:   ts->steps--;
5655:   ts->steprollback = PETSC_TRUE;
5656:   PetscFunctionReturn(PETSC_SUCCESS);
5657: }

5659: /*@
5660:    TSGetStages - Get the number of stages and stage values

5662:    Input Parameter:
5663: .  ts - the `TS` context obtained from `TSCreate()`

5665:    Output Parameters:
5666: +  ns - the number of stages
5667: -  Y - the current stage vectors

5669:    Level: advanced

5671:    Note:
5672:    Both `ns` and `Y` can be `NULL`.

5674: .seealso: [](ch_ts), `TS`, `TSCreate()`
5675: @*/
5676: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5677: {
5678:   PetscFunctionBegin;
5682:   if (!ts->ops->getstages) {
5683:     if (ns) *ns = 0;
5684:     if (Y) *Y = NULL;
5685:   } else PetscUseTypeMethod(ts, getstages, ns, Y);
5686:   PetscFunctionReturn(PETSC_SUCCESS);
5687: }

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

5692:   Collective

5694:   Input Parameters:
5695: + ts - the `TS` context
5696: . t - current timestep
5697: . U - state vector
5698: . Udot - time derivative of state vector
5699: . shift - shift to apply, see note below
5700: - ctx - an optional user context

5702:   Output Parameters:
5703: + J - Jacobian matrix (not altered in this routine)
5704: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)

5706:   Level: intermediate

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

5711:   dF/dU + shift*dF/dUdot

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

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

5720: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5721: @*/
5722: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5723: {
5724:   SNES          snes;
5725:   MatFDColoring color;
5726:   PetscBool     hascolor, matcolor = PETSC_FALSE;

5728:   PetscFunctionBegin;
5729:   PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5730:   PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5731:   if (!color) {
5732:     DM         dm;
5733:     ISColoring iscoloring;

5735:     PetscCall(TSGetDM(ts, &dm));
5736:     PetscCall(DMHasColoring(dm, &hascolor));
5737:     if (hascolor && !matcolor) {
5738:       PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5739:       PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5740:       PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5741:       PetscCall(MatFDColoringSetFromOptions(color));
5742:       PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5743:       PetscCall(ISColoringDestroy(&iscoloring));
5744:     } else {
5745:       MatColoring mc;

5747:       PetscCall(MatColoringCreate(B, &mc));
5748:       PetscCall(MatColoringSetDistance(mc, 2));
5749:       PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5750:       PetscCall(MatColoringSetFromOptions(mc));
5751:       PetscCall(MatColoringApply(mc, &iscoloring));
5752:       PetscCall(MatColoringDestroy(&mc));
5753:       PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5754:       PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5755:       PetscCall(MatFDColoringSetFromOptions(color));
5756:       PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5757:       PetscCall(ISColoringDestroy(&iscoloring));
5758:     }
5759:     PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5760:     PetscCall(PetscObjectDereference((PetscObject)color));
5761:   }
5762:   PetscCall(TSGetSNES(ts, &snes));
5763:   PetscCall(MatFDColoringApply(B, color, U, snes));
5764:   if (J != B) {
5765:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5766:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5767:   }
5768:   PetscFunctionReturn(PETSC_SUCCESS);
5769: }

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

5774:     Input Parameters:
5775: +    ts - the `TS` context
5776: -    func - function called within `TSFunctionDomainError()`

5778:     Calling sequence of `func`:
5779: $   PetscErrorCode func(TS ts, PetscReal time, Vec state, PetscBool reject)
5780: +   ts - the TS context
5781: .   time - the current time (of the stage)
5782: .   state - the state to check if it is valid
5783: -   reject - (output parameter) `PETSC_FALSE` if the state is acceptable, `PETSC_TRUE` if not acceptable

5785:     Level: intermediate

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

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

5797: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5798: @*/

5800: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS, PetscReal, Vec, PetscBool *))
5801: {
5802:   PetscFunctionBegin;
5804:   ts->functiondomainerror = func;
5805:   PetscFunctionReturn(PETSC_SUCCESS);
5806: }

5808: /*@
5809:     TSFunctionDomainError - Checks if the current state is valid

5811:     Input Parameters:
5812: +    ts - the `TS` context
5813: .    stagetime - time of the simulation
5814: -    Y - state vector to check.

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

5819:     Level: developer

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

5825: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5826: @*/
5827: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5828: {
5829:   PetscFunctionBegin;
5831:   *accept = PETSC_TRUE;
5832:   if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5833:   PetscFunctionReturn(PETSC_SUCCESS);
5834: }

5836: /*@C
5837:   TSClone - This function clones a time step `TS` object.

5839:   Collective

5841:   Input Parameter:
5842: . tsin    - The input `TS`

5844:   Output Parameter:
5845: . tsout   - The output `TS` (cloned)

5847:   Level: developer

5849:   Notes:
5850:   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.
5851:   It will likely be replaced in the future with a mechanism of switching methods on the fly.

5853:   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
5854: .vb
5855:  SNES snes_dup = NULL;
5856:  TSGetSNES(ts,&snes_dup);
5857:  TSSetSNES(ts,snes_dup);
5858: .ve

5860: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5861: @*/
5862: PetscErrorCode TSClone(TS tsin, TS *tsout)
5863: {
5864:   TS     t;
5865:   SNES   snes_start;
5866:   DM     dm;
5867:   TSType type;

5869:   PetscFunctionBegin;
5871:   *tsout = NULL;

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

5875:   /* General TS description */
5876:   t->numbermonitors    = 0;
5877:   t->monitorFrequency  = 1;
5878:   t->setupcalled       = 0;
5879:   t->ksp_its           = 0;
5880:   t->snes_its          = 0;
5881:   t->nwork             = 0;
5882:   t->rhsjacobian.time  = PETSC_MIN_REAL;
5883:   t->rhsjacobian.scale = 1.;
5884:   t->ijacobian.shift   = 1.;

5886:   PetscCall(TSGetSNES(tsin, &snes_start));
5887:   PetscCall(TSSetSNES(t, snes_start));

5889:   PetscCall(TSGetDM(tsin, &dm));
5890:   PetscCall(TSSetDM(t, dm));

5892:   t->adapt = tsin->adapt;
5893:   PetscCall(PetscObjectReference((PetscObject)t->adapt));

5895:   t->trajectory = tsin->trajectory;
5896:   PetscCall(PetscObjectReference((PetscObject)t->trajectory));

5898:   t->event = tsin->event;
5899:   if (t->event) t->event->refct++;

5901:   t->problem_type      = tsin->problem_type;
5902:   t->ptime             = tsin->ptime;
5903:   t->ptime_prev        = tsin->ptime_prev;
5904:   t->time_step         = tsin->time_step;
5905:   t->max_time          = tsin->max_time;
5906:   t->steps             = tsin->steps;
5907:   t->max_steps         = tsin->max_steps;
5908:   t->equation_type     = tsin->equation_type;
5909:   t->atol              = tsin->atol;
5910:   t->rtol              = tsin->rtol;
5911:   t->max_snes_failures = tsin->max_snes_failures;
5912:   t->max_reject        = tsin->max_reject;
5913:   t->errorifstepfailed = tsin->errorifstepfailed;

5915:   PetscCall(TSGetType(tsin, &type));
5916:   PetscCall(TSSetType(t, type));

5918:   t->vec_sol = NULL;

5920:   t->cfltime          = tsin->cfltime;
5921:   t->cfltime_local    = tsin->cfltime_local;
5922:   t->exact_final_time = tsin->exact_final_time;

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

5926:   if (((PetscObject)tsin)->fortran_func_pointers) {
5927:     PetscInt i;
5928:     PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers));
5929:     for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5930:   }
5931:   *tsout = t;
5932:   PetscFunctionReturn(PETSC_SUCCESS);
5933: }

5935: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5936: {
5937:   TS ts = (TS)ctx;

5939:   PetscFunctionBegin;
5940:   PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5941:   PetscFunctionReturn(PETSC_SUCCESS);
5942: }

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

5947:    Logically Collective

5949:     Input Parameter:
5950: .    TS - the time stepping routine

5952:    Output Parameter:
5953: .   flg - `PETSC_TRUE` if the multiply is likely correct

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

5958:    Level: advanced

5960:    Note:
5961:     This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian

5963: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5964: @*/
5965: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5966: {
5967:   Mat           J, B;
5968:   TSRHSJacobian func;
5969:   void         *ctx;

5971:   PetscFunctionBegin;
5972:   PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5973:   PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5974:   PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5975:   PetscFunctionReturn(PETSC_SUCCESS);
5976: }

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

5981:    Logically Collective

5983:     Input Parameter:
5984: .    TS - the time stepping routine

5986:    Output Parameter:
5987: .   flg - `PETSC_TRUE` if the multiply is likely correct

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

5992:    Level: advanced

5994:    Notes:
5995:     This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian

5997: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5998: @*/
5999: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
6000: {
6001:   Mat           J, B;
6002:   void         *ctx;
6003:   TSRHSJacobian func;

6005:   PetscFunctionBegin;
6006:   PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
6007:   PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
6008:   PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
6009:   PetscFunctionReturn(PETSC_SUCCESS);
6010: }

6012: /*@
6013:   TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.

6015:   Logically Collective

6017:   Input Parameters:
6018: +  ts - timestepping context
6019: -  use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

6021:   Options Database Key:
6022: .   -ts_use_splitrhsfunction - <true,false>

6024:   Level: intermediate

6026:   Note:
6027:   This is only for multirate methods

6029: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
6030: @*/
6031: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
6032: {
6033:   PetscFunctionBegin;
6035:   ts->use_splitrhsfunction = use_splitrhsfunction;
6036:   PetscFunctionReturn(PETSC_SUCCESS);
6037: }

6039: /*@
6040:   TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.

6042:   Not Collective

6044:   Input Parameter:
6045: .  ts - timestepping context

6047:   Output Parameter:
6048: .  use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

6050:   Level: intermediate

6052: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
6053: @*/
6054: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
6055: {
6056:   PetscFunctionBegin;
6058:   *use_splitrhsfunction = ts->use_splitrhsfunction;
6059:   PetscFunctionReturn(PETSC_SUCCESS);
6060: }

6062: /*@
6063:     TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.

6065:    Logically  Collective

6067:    Input Parameters:
6068: +  ts - the time-stepper
6069: -  str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)

6071:    Level: intermediate

6073:    Note:
6074:      When the relationship between the nonzero structures is known and supplied the solution process can be much faster

6076: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
6077:  @*/
6078: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
6079: {
6080:   PetscFunctionBegin;
6082:   ts->axpy_pattern = str;
6083:   PetscFunctionReturn(PETSC_SUCCESS);
6084: }

6086: /*@
6087:   TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span

6089:   Collective

6091:   Input Parameters:
6092: + ts - the time-stepper
6093: . n - number of the time points (>=2)
6094: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.

6096:   Options Database Key:
6097: . -ts_time_span <t0,...tf> - Sets the time span

6099:   Level: intermediate

6101:   Notes:
6102:   The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
6103:   `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
6104:   The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may
6105:   pressure the memory system when using a large number of span points.

6107: .seealso: [](ch_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()`
6108:  @*/
6109: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
6110: {
6111:   PetscFunctionBegin;
6113:   PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
6114:   if (ts->tspan && n != ts->tspan->num_span_times) {
6115:     PetscCall(PetscFree(ts->tspan->span_times));
6116:     PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
6117:     PetscCall(PetscMalloc1(n, &ts->tspan->span_times));
6118:   }
6119:   if (!ts->tspan) {
6120:     TSTimeSpan tspan;
6121:     PetscCall(PetscNew(&tspan));
6122:     PetscCall(PetscMalloc1(n, &tspan->span_times));
6123:     tspan->reltol = 1e-6;
6124:     tspan->abstol = 10 * PETSC_MACHINE_EPSILON;
6125:     ts->tspan     = tspan;
6126:   }
6127:   ts->tspan->num_span_times = n;
6128:   PetscCall(PetscArraycpy(ts->tspan->span_times, span_times, n));
6129:   PetscCall(TSSetTime(ts, ts->tspan->span_times[0]));
6130:   PetscCall(TSSetMaxTime(ts, ts->tspan->span_times[n - 1]));
6131:   PetscFunctionReturn(PETSC_SUCCESS);
6132: }

6134: /*@C
6135:   TSGetTimeSpan - gets the time span set with `TSSetTimeSpan()`

6137:   Not Collective

6139:   Input Parameter:
6140: . ts - the time-stepper

6142:   Output Parameters:
6143: + n - number of the time points (>=2)
6144: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.

6146:   Level: beginner

6148:   Note:
6149:   The values obtained are valid until the `TS` object is destroyed.

6151:   Both `n` and `span_times` can be `NULL`.

6153: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()`
6154:  @*/
6155: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal **span_times)
6156: {
6157:   PetscFunctionBegin;
6161:   if (!ts->tspan) {
6162:     if (n) *n = 0;
6163:     if (span_times) *span_times = NULL;
6164:   } else {
6165:     if (n) *n = ts->tspan->num_span_times;
6166:     if (span_times) *span_times = ts->tspan->span_times;
6167:   }
6168:   PetscFunctionReturn(PETSC_SUCCESS);
6169: }

6171: /*@
6172:    TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.

6174:    Input Parameter:
6175: .  ts - the `TS` context obtained from `TSCreate()`

6177:    Output Parameters:
6178: +  nsol - the number of solutions
6179: -  Sols - the solution vectors

6181:    Level: intermediate

6183:    Notes:
6184:     Both `nsol` and `Sols` can be `NULL`.

6186:     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()`.
6187:     For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span.

6189: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`
6190: @*/
6191: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
6192: {
6193:   PetscFunctionBegin;
6197:   if (!ts->tspan) {
6198:     if (nsol) *nsol = 0;
6199:     if (Sols) *Sols = NULL;
6200:   } else {
6201:     if (nsol) *nsol = ts->tspan->spanctr;
6202:     if (Sols) *Sols = ts->tspan->vecs_sol;
6203:   }
6204:   PetscFunctionReturn(PETSC_SUCCESS);
6205: }

6207: /*@C
6208:   TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.

6210:   Collective

6212:   Input Parameters:
6213: + ts - the `TS` context
6214: . J - Jacobian matrix (not altered in this routine)
6215: - B - newly computed Jacobian matrix to use with preconditioner

6217:   Level: intermediate

6219:   Notes:
6220:   This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
6221:   many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
6222:   and multiple fields are involved.

6224:   Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
6225:   structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
6226:   usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
6227:   `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.

6229: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
6230: @*/
6231: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
6232: {
6233:   MatColoring   mc            = NULL;
6234:   ISColoring    iscoloring    = NULL;
6235:   MatFDColoring matfdcoloring = NULL;

6237:   PetscFunctionBegin;
6238:   /* Generate new coloring after eliminating zeros in the matrix */
6239:   PetscCall(MatEliminateZeros(B));
6240:   PetscCall(MatColoringCreate(B, &mc));
6241:   PetscCall(MatColoringSetDistance(mc, 2));
6242:   PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
6243:   PetscCall(MatColoringSetFromOptions(mc));
6244:   PetscCall(MatColoringApply(mc, &iscoloring));
6245:   PetscCall(MatColoringDestroy(&mc));
6246:   /* Replace the old coloring with the new one */
6247:   PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
6248:   PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
6249:   PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
6250:   PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
6251:   PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
6252:   PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
6253:   PetscCall(ISColoringDestroy(&iscoloring));
6254:   PetscFunctionReturn(PETSC_SUCCESS);
6255: }