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: /* Logging support */
 11: PetscClassId  TS_CLASSID, DMTS_CLASSID;
 12: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;

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

 16: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt, TSAdaptType default_type)
 17: {
 18:   PetscFunctionBegin;
 20:   PetscAssertPointer(default_type, 2);
 21:   if (!((PetscObject)adapt)->type_name) PetscCall(TSAdaptSetType(adapt, default_type));
 22:   PetscFunctionReturn(PETSC_SUCCESS);
 23: }

 25: /*@
 26:   TSSetFromOptions - Sets various `TS` parameters from the options database

 28:   Collective

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

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

 70:   Level: beginner

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

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

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

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

 96:   PetscFunctionBegin;

 99:   PetscCall(TSRegisterAll());
100:   PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

408:   /* Handle specific TS options */
409:   PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);

411:   /* Handle TSAdapt options */
412:   PetscCall(TSGetAdapt(ts, &ts->adapt));
413:   PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
414:   PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject));

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

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

423:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
424:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject));
425:   PetscOptionsEnd();

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

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

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

442:   Collective

444:   Input Parameter:
445: . ts - the `TS` context obtained from `TSCreate()`

447:   Output Parameter:
448: . tr - the `TSTrajectory` object, if it exists

450:   Level: advanced

452:   Note:
453:   This routine should be called after all `TS` options have been set

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

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

468:   Collective

470:   Input Parameter:
471: . ts - the `TS` context obtained from `TSCreate()`

473:   Options Database Keys:
474: + -ts_save_trajectory      - saves the trajectory to a file
475: - -ts_trajectory_type type - set trajectory type

477:   Level: intermediate

479:   Notes:
480:   This routine should be called after all `TS` options have been set

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

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

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

498:   Collective

500:   Input Parameter:
501: . ts - the `TS` context obtained from `TSCreate()`

503:   Level: intermediate

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

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

521:   Collective

523:   Input Parameter:
524: . ts - the `TS` context obtained from `TSCreate()`

526:   Level: intermediate

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

538: /*@
539:   TSComputeRHSJacobian - Computes the Jacobian matrix that has been
540:   set with `TSSetRHSJacobian()`.

542:   Collective

544:   Input Parameters:
545: + ts - the `TS` context
546: . t  - current timestep
547: - U  - input vector

549:   Output Parameters:
550: + A - Jacobian matrix
551: - B - optional preconditioning matrix

553:   Level: developer

555:   Note:
556:   Most users should not need to explicitly call this routine, as it
557:   is used internally within the nonlinear solvers.

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

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

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

584:   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);
585:   if (rhsjacobianfunc) {
586:     PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
587:     PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
588:     ts->rhsjacs++;
589:     PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
590:   } else {
591:     PetscCall(MatZeroEntries(A));
592:     if (B && A != B) PetscCall(MatZeroEntries(B));
593:   }
594:   ts->rhsjacobian.time  = t;
595:   ts->rhsjacobian.shift = 0;
596:   ts->rhsjacobian.scale = 1.;
597:   PetscCall(PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid));
598:   PetscCall(PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate));
599:   PetscFunctionReturn(PETSC_SUCCESS);
600: }

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

605:   Collective

607:   Input Parameters:
608: + ts - the `TS` context
609: . t  - current time
610: - U  - state vector

612:   Output Parameter:
613: . y - right hand side

615:   Level: developer

617:   Note:
618:   Most users should not need to explicitly call this routine, as it
619:   is used internally within the nonlinear solvers.

621: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
622: @*/
623: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
624: {
625:   TSRHSFunctionFn *rhsfunction;
626:   TSIFunctionFn   *ifunction;
627:   void            *ctx;
628:   DM               dm;

630:   PetscFunctionBegin;
634:   PetscCall(TSGetDM(ts, &dm));
635:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx));
636:   PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));

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

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

651: /*@
652:   TSComputeSolutionFunction - Evaluates the solution function.

654:   Collective

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

660:   Output Parameter:
661: . U - the solution

663:   Level: developer

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

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

684:   Collective

686:   Input Parameters:
687: + ts - the `TS` context
688: - t  - current time

690:   Output Parameter:
691: . U - the function value

693:   Level: developer

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

703:   PetscFunctionBegin;
706:   PetscCall(TSGetDM(ts, &dm));
707:   PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));

709:   if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
710:   PetscFunctionReturn(PETSC_SUCCESS);
711: }

713: static PetscErrorCode TSGetRHSVec_Private(TS ts, Vec *Frhs)
714: {
715:   Vec F;

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

725: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
726: {
727:   Mat            A, B;
728:   TSIJacobianFn *ijacobian;

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

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

777:   Collective

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

786:   Output Parameter:
787: . Y - right hand side

789:   Level: developer

791:   Note:
792:   Most users should not need to explicitly call this routine, as it
793:   is used internally within the nonlinear solvers.

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

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

807:   PetscFunctionBegin;

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

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

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

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

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

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

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

866: /*@
867:   TSComputeIJacobian - Evaluates the Jacobian of the DAE

869:   Collective

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

879:   Output Parameters:
880: + A - Jacobian matrix
881: - B - matrix from which the preconditioner is constructed; often the same as `A`

883:   Level: developer

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

893: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
894: @*/
895: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
896: {
897:   TSIJacobianFn   *ijacobian;
898:   TSRHSJacobianFn *rhsjacobian;
899:   DM               dm;
900:   void            *ctx;

902:   PetscFunctionBegin;

909:   PetscCall(TSGetDM(ts, &dm));
910:   PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
911:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));

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

915:   PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
916:   if (ijacobian) {
917:     PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
918:     ts->ijacs++;
919:   }
920:   if (imex) {
921:     if (!ijacobian) { /* system was written as Udot = G(t,U) */
922:       PetscBool assembled;
923:       if (rhsjacobian) {
924:         Mat Arhs = NULL;
925:         PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
926:         if (A == Arhs) {
927:           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 */
928:           ts->rhsjacobian.time = PETSC_MIN_REAL;
929:         }
930:       }
931:       PetscCall(MatZeroEntries(A));
932:       PetscCall(MatAssembled(A, &assembled));
933:       if (!assembled) {
934:         PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
935:         PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
936:       }
937:       PetscCall(MatShift(A, shift));
938:       if (A != B) {
939:         PetscCall(MatZeroEntries(B));
940:         PetscCall(MatAssembled(B, &assembled));
941:         if (!assembled) {
942:           PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
943:           PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
944:         }
945:         PetscCall(MatShift(B, shift));
946:       }
947:     }
948:   } else {
949:     Mat Arhs = NULL, Brhs = NULL;

951:     /* RHSJacobian needs to be converted to part of IJacobian if exists */
952:     if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
953:     if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
954:       PetscObjectState Ustate;
955:       PetscObjectId    Uid;
956:       TSRHSFunctionFn *rhsfunction;

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

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

1006: /*@C
1007:   TSSetRHSFunction - Sets the routine for evaluating the function,
1008:   where U_t = G(t,u).

1010:   Logically Collective

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

1018:   Level: beginner

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

1023: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1024: @*/
1025: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, void *ctx)
1026: {
1027:   SNES snes;
1028:   Vec  ralloc = NULL;
1029:   DM   dm;

1031:   PetscFunctionBegin;

1035:   PetscCall(TSGetDM(ts, &dm));
1036:   PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1037:   PetscCall(TSGetSNES(ts, &snes));
1038:   if (!r && !ts->dm && ts->vec_sol) {
1039:     PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1040:     r = ralloc;
1041:   }
1042:   PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1043:   PetscCall(VecDestroy(&ralloc));
1044:   PetscFunctionReturn(PETSC_SUCCESS);
1045: }

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

1050:   Logically Collective

1052:   Input Parameters:
1053: + ts  - the `TS` context obtained from `TSCreate()`
1054: . f   - routine for evaluating the solution
1055: - ctx - [optional] user-defined context for private data for the
1056:           function evaluation routine (may be `NULL`)

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

1062:   Level: intermediate

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

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

1071: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1072: @*/
1073: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, void *ctx)
1074: {
1075:   DM dm;

1077:   PetscFunctionBegin;
1079:   PetscCall(TSGetDM(ts, &dm));
1080:   PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1081:   PetscFunctionReturn(PETSC_SUCCESS);
1082: }

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

1087:   Logically Collective

1089:   Input Parameters:
1090: + ts   - the `TS` context obtained from `TSCreate()`
1091: . func - routine for evaluating the forcing function
1092: - ctx  - [optional] user-defined context for private data for the function evaluation routine
1093:          (may be `NULL`)

1095:   Level: intermediate

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

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

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

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

1109: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
1110: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1111: @*/
1112: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, void *ctx)
1113: {
1114:   DM dm;

1116:   PetscFunctionBegin;
1118:   PetscCall(TSGetDM(ts, &dm));
1119:   PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1120:   PetscFunctionReturn(PETSC_SUCCESS);
1121: }

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

1127:   Logically Collective

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

1136:   Level: beginner

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

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

1144: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
1145: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
1146: @*/
1147: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, void *ctx)
1148: {
1149:   SNES           snes;
1150:   DM             dm;
1151:   TSIJacobianFn *ijacobian;

1153:   PetscFunctionBegin;
1157:   if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1158:   if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);

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

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

1181:   Logically Collective

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

1189:   Level: beginner

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

1194: .seealso: [](ch_ts), `TS`, `TSIFunctionFn`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`,
1195: `TSSetIJacobian()`
1196: @*/
1197: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunctionFn *f, void *ctx)
1198: {
1199:   SNES snes;
1200:   Vec  ralloc = NULL;
1201:   DM   dm;

1203:   PetscFunctionBegin;

1207:   PetscCall(TSGetDM(ts, &dm));
1208:   PetscCall(DMTSSetIFunction(dm, f, ctx));

1210:   PetscCall(TSGetSNES(ts, &snes));
1211:   if (!r && !ts->dm && ts->vec_sol) {
1212:     PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1213:     r = ralloc;
1214:   }
1215:   PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1216:   PetscCall(VecDestroy(&ralloc));
1217:   PetscFunctionReturn(PETSC_SUCCESS);
1218: }

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

1223:   Not Collective

1225:   Input Parameter:
1226: . ts - the `TS` context

1228:   Output Parameters:
1229: + r    - vector to hold residual (or `NULL`)
1230: . func - the function to compute residual (or `NULL`)
1231: - ctx  - the function context (or `NULL`)

1233:   Level: advanced

1235: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1236: @*/
1237: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunctionFn **func, void **ctx)
1238: {
1239:   SNES snes;
1240:   DM   dm;

1242:   PetscFunctionBegin;
1244:   PetscCall(TSGetSNES(ts, &snes));
1245:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1246:   PetscCall(TSGetDM(ts, &dm));
1247:   PetscCall(DMTSGetIFunction(dm, func, ctx));
1248:   PetscFunctionReturn(PETSC_SUCCESS);
1249: }

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

1254:   Not Collective

1256:   Input Parameter:
1257: . ts - the `TS` context

1259:   Output Parameters:
1260: + r    - vector to hold computed right hand side (or `NULL`)
1261: . func - the function to compute right hand side (or `NULL`)
1262: - ctx  - the function context (or `NULL`)

1264:   Level: advanced

1266: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1267: @*/
1268: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunctionFn **func, void **ctx)
1269: {
1270:   SNES snes;
1271:   DM   dm;

1273:   PetscFunctionBegin;
1275:   PetscCall(TSGetSNES(ts, &snes));
1276:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1277:   PetscCall(TSGetDM(ts, &dm));
1278:   PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1279:   PetscFunctionReturn(PETSC_SUCCESS);
1280: }

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

1286:   Logically Collective

1288:   Input Parameters:
1289: + ts   - the `TS` context obtained from `TSCreate()`
1290: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
1291: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1292: . f    - the Jacobian evaluation routine
1293: - ctx  - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)

1295:   Level: beginner

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

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

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

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

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

1315: .seealso: [](ch_ts), `TS`, `TSIJacobianFn`, `TSSetIFunction()`, `TSSetRHSJacobian()`,
1316: `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1317: @*/
1318: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobianFn *f, void *ctx)
1319: {
1320:   SNES snes;
1321:   DM   dm;

1323:   PetscFunctionBegin;
1327:   if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1328:   if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);

1330:   PetscCall(TSGetDM(ts, &dm));
1331:   PetscCall(DMTSSetIJacobian(dm, f, ctx));

1333:   PetscCall(TSGetSNES(ts, &snes));
1334:   PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1335:   PetscFunctionReturn(PETSC_SUCCESS);
1336: }

1338: /*@
1339:   TSRHSJacobianSetReuse - restore the RHS Jacobian before calling the user-provided `TSRHSJacobianFn` function again

1341:   Logically Collective

1343:   Input Parameters:
1344: + ts    - `TS` context obtained from `TSCreate()`
1345: - reuse - `PETSC_TRUE` if the RHS Jacobian

1347:   Level: intermediate

1349:   Notes:
1350:   Without this flag, `TS` will change the sign and shift the RHS Jacobian for a
1351:   finite-time-step implicit solve, in which case the user function will need to recompute the
1352:   entire Jacobian.  The `reuse `flag must be set if the evaluation function assumes that the
1353:   matrix entries have not been changed by the `TS`.

1355: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1356: @*/
1357: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1358: {
1359:   PetscFunctionBegin;
1360:   ts->rhsjacobian.reuse = reuse;
1361:   PetscFunctionReturn(PETSC_SUCCESS);
1362: }

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

1367:   Logically Collective

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

1375:   Level: beginner

1377: .seealso: [](ch_ts), `TS`, `TSI2FunctionFn`, `TSSetI2Jacobian()`, `TSSetIFunction()`,
1378: `TSCreate()`, `TSSetRHSFunction()`
1379: @*/
1380: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2FunctionFn *fun, void *ctx)
1381: {
1382:   DM dm;

1384:   PetscFunctionBegin;
1387:   PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1388:   PetscCall(TSGetDM(ts, &dm));
1389:   PetscCall(DMTSSetI2Function(dm, fun, ctx));
1390:   PetscFunctionReturn(PETSC_SUCCESS);
1391: }

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

1396:   Not Collective

1398:   Input Parameter:
1399: . ts - the `TS` context

1401:   Output Parameters:
1402: + r   - vector to hold residual (or `NULL`)
1403: . fun - the function to compute residual (or `NULL`)
1404: - ctx - the function context (or `NULL`)

1406:   Level: advanced

1408: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1409: @*/
1410: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2FunctionFn **fun, void **ctx)
1411: {
1412:   SNES snes;
1413:   DM   dm;

1415:   PetscFunctionBegin;
1417:   PetscCall(TSGetSNES(ts, &snes));
1418:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1419:   PetscCall(TSGetDM(ts, &dm));
1420:   PetscCall(DMTSGetI2Function(dm, fun, ctx));
1421:   PetscFunctionReturn(PETSC_SUCCESS);
1422: }

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

1428:   Logically Collective

1430:   Input Parameters:
1431: + ts  - the `TS` context obtained from `TSCreate()`
1432: . J   - matrix to hold the Jacobian values
1433: . P   - matrix for constructing the preconditioner (may be same as `J`)
1434: . jac - the Jacobian evaluation routine, see `TSI2JacobianFn` for the calling sequence
1435: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)

1437:   Level: beginner

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

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

1447: .seealso: [](ch_ts), `TS`, `TSI2JacobianFn`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1448: @*/
1449: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2JacobianFn *jac, void *ctx)
1450: {
1451:   DM dm;

1453:   PetscFunctionBegin;
1457:   PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1458:   PetscCall(TSGetDM(ts, &dm));
1459:   PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1460:   PetscFunctionReturn(PETSC_SUCCESS);
1461: }

1463: /*@C
1464:   TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.

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

1468:   Input Parameter:
1469: . ts - The `TS` context obtained from `TSCreate()`

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

1477:   Level: advanced

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

1482: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1483: @*/
1484: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2JacobianFn **jac, void **ctx)
1485: {
1486:   SNES snes;
1487:   DM   dm;

1489:   PetscFunctionBegin;
1490:   PetscCall(TSGetSNES(ts, &snes));
1491:   PetscCall(SNESSetUpMatrices(snes));
1492:   PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1493:   PetscCall(TSGetDM(ts, &dm));
1494:   PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1495:   PetscFunctionReturn(PETSC_SUCCESS);
1496: }

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

1501:   Collective

1503:   Input Parameters:
1504: + ts - the `TS` context
1505: . t  - current time
1506: . U  - state vector
1507: . V  - time derivative of state vector (U_t)
1508: - A  - second time derivative of state vector (U_tt)

1510:   Output Parameter:
1511: . F - the residual vector

1513:   Level: developer

1515:   Note:
1516:   Most users should not need to explicitly call this routine, as it
1517:   is used internally within the nonlinear solvers.

1519: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1520: @*/
1521: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1522: {
1523:   DM               dm;
1524:   TSI2FunctionFn  *I2Function;
1525:   void            *ctx;
1526:   TSRHSFunctionFn *rhsfunction;

1528:   PetscFunctionBegin;

1535:   PetscCall(TSGetDM(ts, &dm));
1536:   PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1537:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));

1539:   if (!I2Function) {
1540:     PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1541:     PetscFunctionReturn(PETSC_SUCCESS);
1542:   }

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

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

1548:   if (rhsfunction) {
1549:     Vec Frhs;
1550:     PetscCall(TSGetRHSVec_Private(ts, &Frhs));
1551:     PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1552:     PetscCall(VecAXPY(F, -1, Frhs));
1553:   }

1555:   PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, V, F));
1556:   PetscFunctionReturn(PETSC_SUCCESS);
1557: }

1559: /*@
1560:   TSComputeI2Jacobian - Evaluates the Jacobian of the DAE

1562:   Collective

1564:   Input Parameters:
1565: + ts     - the `TS` context
1566: . t      - current timestep
1567: . U      - state vector
1568: . V      - time derivative of state vector
1569: . A      - second time derivative of state vector
1570: . shiftV - shift to apply, see note below
1571: - shiftA - shift to apply, see note below

1573:   Output Parameters:
1574: + J - Jacobian matrix
1575: - P - optional preconditioning matrix

1577:   Level: developer

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

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

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

1587: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1588: @*/
1589: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1590: {
1591:   DM               dm;
1592:   TSI2JacobianFn  *I2Jacobian;
1593:   void            *ctx;
1594:   TSRHSJacobianFn *rhsjacobian;

1596:   PetscFunctionBegin;

1604:   PetscCall(TSGetDM(ts, &dm));
1605:   PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1606:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));

1608:   if (!I2Jacobian) {
1609:     PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1610:     PetscFunctionReturn(PETSC_SUCCESS);
1611:   }

1613:   PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, J, P));
1614:   PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1615:   if (rhsjacobian) {
1616:     Mat Jrhs, Prhs;
1617:     PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1618:     PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1619:     PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1620:     if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1621:   }

1623:   PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, J, P));
1624:   PetscFunctionReturn(PETSC_SUCCESS);
1625: }

1627: /*@C
1628:   TSSetTransientVariable - sets function to transform from state to transient variables

1630:   Logically Collective

1632:   Input Parameters:
1633: + ts   - time stepping context on which to change the transient variable
1634: . tvar - a function that transforms to transient variables, see `TSTransientVariableFn` for the calling sequence
1635: - ctx  - a context for tvar

1637:   Level: advanced

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

1649: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSTransientVariableFn`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1650: @*/
1651: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariableFn *tvar, void *ctx)
1652: {
1653:   DM dm;

1655:   PetscFunctionBegin;
1657:   PetscCall(TSGetDM(ts, &dm));
1658:   PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1659:   PetscFunctionReturn(PETSC_SUCCESS);
1660: }

1662: /*@
1663:   TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables

1665:   Logically Collective

1667:   Input Parameters:
1668: + ts - TS on which to compute
1669: - U  - state vector to be transformed to transient variables

1671:   Output Parameter:
1672: . C - transient (conservative) variable

1674:   Level: developer

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

1681: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1682: @*/
1683: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1684: {
1685:   DM   dm;
1686:   DMTS dmts;

1688:   PetscFunctionBegin;
1691:   PetscCall(TSGetDM(ts, &dm));
1692:   PetscCall(DMGetDMTS(dm, &dmts));
1693:   if (dmts->ops->transientvar) {
1695:     PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1696:   }
1697:   PetscFunctionReturn(PETSC_SUCCESS);
1698: }

1700: /*@
1701:   TSHasTransientVariable - determine whether transient variables have been set

1703:   Logically Collective

1705:   Input Parameter:
1706: . ts - `TS` on which to compute

1708:   Output Parameter:
1709: . has - `PETSC_TRUE` if transient variables have been set

1711:   Level: developer

1713: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1714: @*/
1715: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1716: {
1717:   DM   dm;
1718:   DMTS dmts;

1720:   PetscFunctionBegin;
1722:   PetscCall(TSGetDM(ts, &dm));
1723:   PetscCall(DMGetDMTS(dm, &dmts));
1724:   *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1725:   PetscFunctionReturn(PETSC_SUCCESS);
1726: }

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

1732:   Logically Collective

1734:   Input Parameters:
1735: + ts - the `TS` context obtained from `TSCreate()`
1736: . u  - the solution vector
1737: - v  - the time derivative vector

1739:   Level: beginner

1741: .seealso: [](ch_ts), `TS`
1742: @*/
1743: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1744: {
1745:   PetscFunctionBegin;
1749:   PetscCall(TSSetSolution(ts, u));
1750:   PetscCall(PetscObjectReference((PetscObject)v));
1751:   PetscCall(VecDestroy(&ts->vec_dot));
1752:   ts->vec_dot = v;
1753:   PetscFunctionReturn(PETSC_SUCCESS);
1754: }

1756: /*@
1757:   TS2GetSolution - Returns the solution and time derivative at the present timestep
1758:   for second order equations.

1760:   Not Collective

1762:   Input Parameter:
1763: . ts - the `TS` context obtained from `TSCreate()`

1765:   Output Parameters:
1766: + u - the vector containing the solution
1767: - v - the vector containing the time derivative

1769:   Level: intermediate

1771:   Notes:
1772:   It is valid to call this routine inside the function
1773:   that you are evaluating in order to move to the new timestep. This vector not
1774:   changed until the solution at the next timestep has been calculated.

1776: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1777: @*/
1778: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1779: {
1780:   PetscFunctionBegin;
1782:   if (u) PetscAssertPointer(u, 2);
1783:   if (v) PetscAssertPointer(v, 3);
1784:   if (u) *u = ts->vec_sol;
1785:   if (v) *v = ts->vec_dot;
1786:   PetscFunctionReturn(PETSC_SUCCESS);
1787: }

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

1792:   Collective

1794:   Input Parameters:
1795: + ts     - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1796:            some related function before a call to `TSLoad()`.
1797: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`

1799:   Level: intermediate

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

1804: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1805: @*/
1806: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1807: {
1808:   PetscBool isbinary;
1809:   PetscInt  classid;
1810:   char      type[256];
1811:   DMTS      sdm;
1812:   DM        dm;

1814:   PetscFunctionBegin;
1817:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1818:   PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");

1820:   PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1821:   PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1822:   PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1823:   PetscCall(TSSetType(ts, type));
1824:   PetscTryTypeMethod(ts, load, viewer);
1825:   PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1826:   PetscCall(DMLoad(dm, viewer));
1827:   PetscCall(TSSetDM(ts, dm));
1828:   PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1829:   PetscCall(VecLoad(ts->vec_sol, viewer));
1830:   PetscCall(DMGetDMTS(ts->dm, &sdm));
1831:   PetscCall(DMTSLoad(sdm, viewer));
1832:   PetscFunctionReturn(PETSC_SUCCESS);
1833: }

1835: #include <petscdraw.h>
1836: #if defined(PETSC_HAVE_SAWS)
1837: #include <petscviewersaws.h>
1838: #endif

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

1843:   Collective

1845:   Input Parameters:
1846: + ts   - the `TS` context
1847: . obj  - Optional object that provides the prefix for the options database keys
1848: - name - command line option string to be passed by user

1850:   Level: intermediate

1852: .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1853: @*/
1854: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1855: {
1856:   PetscFunctionBegin;
1858:   PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1859:   PetscFunctionReturn(PETSC_SUCCESS);
1860: }

1862: /*@C
1863:   TSView - Prints the `TS` data structure.

1865:   Collective

1867:   Input Parameters:
1868: + ts     - the `TS` context obtained from `TSCreate()`
1869: - viewer - visualization context

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

1874:   Level: beginner

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

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

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

1889: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1890: @*/
1891: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1892: {
1893:   TSType    type;
1894:   PetscBool iascii, isstring, isundials, isbinary, isdraw;
1895:   DMTS      sdm;
1896: #if defined(PETSC_HAVE_SAWS)
1897:   PetscBool issaws;
1898: #endif

1900:   PetscFunctionBegin;
1902:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1904:   PetscCheckSameComm(ts, 1, viewer, 2);

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

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

1969:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1970:     PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1971:     PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
1972:     PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
1973:     PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
1974:     bottom = y - h;
1975:     PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1976:     PetscTryTypeMethod(ts, view, viewer);
1977:     if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1978:     if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
1979:     PetscCall(PetscDrawPopCurrentPoint(draw));
1980: #if defined(PETSC_HAVE_SAWS)
1981:   } else if (issaws) {
1982:     PetscMPIInt rank;
1983:     const char *name;

1985:     PetscCall(PetscObjectGetName((PetscObject)ts, &name));
1986:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1987:     if (!((PetscObject)ts)->amsmem && rank == 0) {
1988:       char dir[1024];

1990:       PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
1991:       PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
1992:       PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
1993:       PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
1994:       PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
1995:     }
1996:     PetscTryTypeMethod(ts, view, viewer);
1997: #endif
1998:   }
1999:   if (ts->snes && ts->usessnes) {
2000:     PetscCall(PetscViewerASCIIPushTab(viewer));
2001:     PetscCall(SNESView(ts->snes, viewer));
2002:     PetscCall(PetscViewerASCIIPopTab(viewer));
2003:   }
2004:   PetscCall(DMGetDMTS(ts->dm, &sdm));
2005:   PetscCall(DMTSView(sdm, viewer));

2007:   PetscCall(PetscViewerASCIIPushTab(viewer));
2008:   PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials));
2009:   PetscCall(PetscViewerASCIIPopTab(viewer));
2010:   PetscFunctionReturn(PETSC_SUCCESS);
2011: }

2013: /*@
2014:   TSSetApplicationContext - Sets an optional user-defined context for
2015:   the timesteppers.

2017:   Logically Collective

2019:   Input Parameters:
2020: + ts   - the `TS` context obtained from `TSCreate()`
2021: - usrP - user context

2023:   Level: intermediate

2025:   Fortran Notes:
2026:   You must write a Fortran interface definition for this
2027:   function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.

2029: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2030: @*/
2031: PetscErrorCode TSSetApplicationContext(TS ts, void *usrP)
2032: {
2033:   PetscFunctionBegin;
2035:   ts->user = usrP;
2036:   PetscFunctionReturn(PETSC_SUCCESS);
2037: }

2039: /*@
2040:   TSGetApplicationContext - Gets the user-defined context for the
2041:   timestepper that was set with `TSSetApplicationContext()`

2043:   Not Collective

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

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

2051:   Level: intermediate

2053:   Fortran Notes:
2054:   You must write a Fortran interface definition for this
2055:   function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.

2057: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2058: @*/
2059: PetscErrorCode TSGetApplicationContext(TS ts, void *usrP)
2060: {
2061:   PetscFunctionBegin;
2063:   *(void **)usrP = ts->user;
2064:   PetscFunctionReturn(PETSC_SUCCESS);
2065: }

2067: /*@
2068:   TSGetStepNumber - Gets the number of time steps completed.

2070:   Not Collective

2072:   Input Parameter:
2073: . ts - the `TS` context obtained from `TSCreate()`

2075:   Output Parameter:
2076: . steps - number of steps completed so far

2078:   Level: intermediate

2080: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2081: @*/
2082: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2083: {
2084:   PetscFunctionBegin;
2086:   PetscAssertPointer(steps, 2);
2087:   *steps = ts->steps;
2088:   PetscFunctionReturn(PETSC_SUCCESS);
2089: }

2091: /*@
2092:   TSSetStepNumber - Sets the number of steps completed.

2094:   Logically Collective

2096:   Input Parameters:
2097: + ts    - the `TS` context
2098: - steps - number of steps completed so far

2100:   Level: developer

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

2112: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2113: @*/
2114: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2115: {
2116:   PetscFunctionBegin;
2119:   PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2120:   ts->steps = steps;
2121:   PetscFunctionReturn(PETSC_SUCCESS);
2122: }

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

2128:   Logically Collective

2130:   Input Parameters:
2131: + ts        - the `TS` context obtained from `TSCreate()`
2132: - time_step - the size of the timestep

2134:   Level: intermediate

2136: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2137: @*/
2138: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2139: {
2140:   PetscFunctionBegin;
2143:   ts->time_step = time_step;
2144:   PetscFunctionReturn(PETSC_SUCCESS);
2145: }

2147: /*@
2148:   TSSetExactFinalTime - Determines whether to adapt the final time step to
2149:   match the exact final time, interpolate solution to the exact final time,
2150:   or just return at the final time `TS` computed.

2152:   Logically Collective

2154:   Input Parameters:
2155: + ts     - the time-step context
2156: - eftopt - exact final time option
2157: .vb
2158:   TS_EXACTFINALTIME_STEPOVER    - Don't do anything if final time is exceeded
2159:   TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2160:   TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2161: .ve

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

2166:   Level: beginner

2168:   Note:
2169:   If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2170:   then the final time you selected.

2172: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2173: @*/
2174: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2175: {
2176:   PetscFunctionBegin;
2179:   ts->exact_final_time = eftopt;
2180:   PetscFunctionReturn(PETSC_SUCCESS);
2181: }

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

2186:   Not Collective

2188:   Input Parameter:
2189: . ts - the `TS` context

2191:   Output Parameter:
2192: . eftopt - exact final time option

2194:   Level: beginner

2196: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2197: @*/
2198: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2199: {
2200:   PetscFunctionBegin;
2202:   PetscAssertPointer(eftopt, 2);
2203:   *eftopt = ts->exact_final_time;
2204:   PetscFunctionReturn(PETSC_SUCCESS);
2205: }

2207: /*@
2208:   TSGetTimeStep - Gets the current timestep size.

2210:   Not Collective

2212:   Input Parameter:
2213: . ts - the `TS` context obtained from `TSCreate()`

2215:   Output Parameter:
2216: . dt - the current timestep size

2218:   Level: intermediate

2220: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2221: @*/
2222: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2223: {
2224:   PetscFunctionBegin;
2226:   PetscAssertPointer(dt, 2);
2227:   *dt = ts->time_step;
2228:   PetscFunctionReturn(PETSC_SUCCESS);
2229: }

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

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

2239:   Input Parameter:
2240: . ts - the `TS` context obtained from `TSCreate()`

2242:   Output Parameter:
2243: . v - the vector containing the solution

2245:   Level: intermediate

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

2251: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2252: @*/
2253: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2254: {
2255:   PetscFunctionBegin;
2257:   PetscAssertPointer(v, 2);
2258:   *v = ts->vec_sol;
2259:   PetscFunctionReturn(PETSC_SUCCESS);
2260: }

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

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

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

2279:   Level: advanced

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

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

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

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

2302:   Level: intermediate

2304: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2305: @*/
2306: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2307: {
2308:   PetscFunctionBegin;
2310:   if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2311:   else PetscCall(VecZeroEntries(*v));
2312:   PetscFunctionReturn(PETSC_SUCCESS);
2313: }

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

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

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

2326:   Level: intermediate

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

2331: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2332: @*/
2333: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2334: {
2335:   PetscFunctionBegin;
2337:   if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2338:   else PetscCall(VecZeroEntries(*v));
2339:   PetscFunctionReturn(PETSC_SUCCESS);
2340: }

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

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

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

2353:   Level: intermediate

2355: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
2356: @*/
2357: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2358: {
2359:   PetscFunctionBegin;
2361:   PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2362:   PetscTryTypeMethod(ts, settimeerror, v);
2363:   PetscFunctionReturn(PETSC_SUCCESS);
2364: }

2366: /* ----- Routines to initialize and destroy a timestepper ---- */
2367: /*@
2368:   TSSetProblemType - Sets the type of problem to be solved.

2370:   Not collective

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

2381:   Level: beginner

2383: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2384: @*/
2385: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2386: {
2387:   PetscFunctionBegin;
2389:   ts->problem_type = type;
2390:   if (type == TS_LINEAR) {
2391:     SNES snes;
2392:     PetscCall(TSGetSNES(ts, &snes));
2393:     PetscCall(SNESSetType(snes, SNESKSPONLY));
2394:   }
2395:   PetscFunctionReturn(PETSC_SUCCESS);
2396: }

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

2401:   Not collective

2403:   Input Parameter:
2404: . ts - The `TS`

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

2414:   Level: beginner

2416: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2417: @*/
2418: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2419: {
2420:   PetscFunctionBegin;
2422:   PetscAssertPointer(type, 2);
2423:   *type = ts->problem_type;
2424:   PetscFunctionReturn(PETSC_SUCCESS);
2425: }

2427: /*
2428:     Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2429: */
2430: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2431: {
2432:   PetscBool isnone;

2434:   PetscFunctionBegin;
2435:   PetscCall(TSGetAdapt(ts, &ts->adapt));
2436:   PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));

2438:   PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2439:   if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2440:   else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2441:   PetscFunctionReturn(PETSC_SUCCESS);
2442: }

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

2447:   Collective

2449:   Input Parameter:
2450: . ts - the `TS` context obtained from `TSCreate()`

2452:   Level: advanced

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

2461: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2462: @*/
2463: PetscErrorCode TSSetUp(TS ts)
2464: {
2465:   DM dm;
2466:   PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2467:   PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2468:   TSIFunctionFn   *ifun;
2469:   TSIJacobianFn   *ijac;
2470:   TSI2JacobianFn  *i2jac;
2471:   TSRHSJacobianFn *rhsjac;

2473:   PetscFunctionBegin;
2475:   if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);

2477:   if (!((PetscObject)ts)->type_name) {
2478:     PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2479:     PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2480:   }

2482:   if (!ts->vec_sol) {
2483:     PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2484:     PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2485:   }

2487:   if (ts->tspan) {
2488:     if (!ts->tspan->vecs_sol) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2489:   }
2490:   if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2491:     PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2492:     ts->Jacp = ts->Jacprhs;
2493:   }

2495:   if (ts->quadraturets) {
2496:     PetscCall(TSSetUp(ts->quadraturets));
2497:     PetscCall(VecDestroy(&ts->vec_costintegrand));
2498:     PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2499:   }

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

2522:   PetscCall(TSGetAdapt(ts, &ts->adapt));
2523:   PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));

2525:   PetscTryTypeMethod(ts, setup);

2527:   PetscCall(TSSetExactFinalTimeDefault(ts));

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

2536:   /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2537:      Otherwise, the SNES will use coloring internally to form the Jacobian.
2538:    */
2539:   PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2540:   PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2541:   PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2542:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2543:   if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));

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

2548:   ts->setupcalled = PETSC_TRUE;
2549:   PetscFunctionReturn(PETSC_SUCCESS);
2550: }

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

2555:   Collective

2557:   Input Parameter:
2558: . ts - the `TS` context obtained from `TSCreate()`

2560:   Level: beginner

2562: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()`
2563: @*/
2564: PetscErrorCode TSReset(TS ts)
2565: {
2566:   TS_RHSSplitLink ilink = ts->tsrhssplit, next;

2568:   PetscFunctionBegin;

2571:   PetscTryTypeMethod(ts, reset);
2572:   if (ts->snes) PetscCall(SNESReset(ts->snes));
2573:   if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));

2575:   PetscCall(MatDestroy(&ts->Arhs));
2576:   PetscCall(MatDestroy(&ts->Brhs));
2577:   PetscCall(VecDestroy(&ts->Frhs));
2578:   PetscCall(VecDestroy(&ts->vec_sol));
2579:   PetscCall(VecDestroy(&ts->vec_dot));
2580:   PetscCall(VecDestroy(&ts->vatol));
2581:   PetscCall(VecDestroy(&ts->vrtol));
2582:   PetscCall(VecDestroyVecs(ts->nwork, &ts->work));

2584:   PetscCall(MatDestroy(&ts->Jacprhs));
2585:   PetscCall(MatDestroy(&ts->Jacp));
2586:   if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2587:   if (ts->quadraturets) {
2588:     PetscCall(TSReset(ts->quadraturets));
2589:     PetscCall(VecDestroy(&ts->vec_costintegrand));
2590:   }
2591:   while (ilink) {
2592:     next = ilink->next;
2593:     PetscCall(TSDestroy(&ilink->ts));
2594:     PetscCall(PetscFree(ilink->splitname));
2595:     PetscCall(ISDestroy(&ilink->is));
2596:     PetscCall(PetscFree(ilink));
2597:     ilink = next;
2598:   }
2599:   ts->tsrhssplit     = NULL;
2600:   ts->num_rhs_splits = 0;
2601:   if (ts->tspan) {
2602:     PetscCall(PetscFree(ts->tspan->span_times));
2603:     PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2604:     PetscCall(PetscFree(ts->tspan));
2605:   }
2606:   ts->setupcalled = PETSC_FALSE;
2607:   PetscFunctionReturn(PETSC_SUCCESS);
2608: }

2610: static PetscErrorCode TSResizeReset(TS);

2612: /*@C
2613:   TSDestroy - Destroys the timestepper context that was created
2614:   with `TSCreate()`.

2616:   Collective

2618:   Input Parameter:
2619: . ts - the `TS` context obtained from `TSCreate()`

2621:   Level: beginner

2623: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2624: @*/
2625: PetscErrorCode TSDestroy(TS *ts)
2626: {
2627:   PetscFunctionBegin;
2628:   if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2630:   if (--((PetscObject)*ts)->refct > 0) {
2631:     *ts = NULL;
2632:     PetscFunctionReturn(PETSC_SUCCESS);
2633:   }

2635:   PetscCall(TSReset(*ts));
2636:   PetscCall(TSAdjointReset(*ts));
2637:   if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2638:   PetscCall(TSResizeReset(*ts));

2640:   /* if memory was published with SAWs then destroy it */
2641:   PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2642:   PetscTryTypeMethod(*ts, destroy);

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

2646:   PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2647:   PetscCall(TSEventDestroy(&(*ts)->event));

2649:   PetscCall(SNESDestroy(&(*ts)->snes));
2650:   PetscCall(DMDestroy(&(*ts)->dm));
2651:   PetscCall(TSMonitorCancel(*ts));
2652:   PetscCall(TSAdjointMonitorCancel(*ts));

2654:   PetscCall(TSDestroy(&(*ts)->quadraturets));
2655:   PetscCall(PetscHeaderDestroy(ts));
2656:   PetscFunctionReturn(PETSC_SUCCESS);
2657: }

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

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

2665:   Input Parameter:
2666: . ts - the `TS` context obtained from `TSCreate()`

2668:   Output Parameter:
2669: . snes - the nonlinear solver context

2671:   Level: beginner

2673:   Notes:
2674:   The user can then directly manipulate the `SNES` context to set various
2675:   options, etc.  Likewise, the user can then extract and manipulate the
2676:   `KSP`, and `PC` contexts as well.

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

2681: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2682: @*/
2683: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2684: {
2685:   PetscFunctionBegin;
2687:   PetscAssertPointer(snes, 2);
2688:   if (!ts->snes) {
2689:     PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2690:     PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2691:     PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2692:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2693:     if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2694:     if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2695:   }
2696:   *snes = ts->snes;
2697:   PetscFunctionReturn(PETSC_SUCCESS);
2698: }

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

2703:   Collective

2705:   Input Parameters:
2706: + ts   - the `TS` context obtained from `TSCreate()`
2707: - snes - the nonlinear solver context

2709:   Level: developer

2711:   Note:
2712:   Most users should have the `TS` created by calling `TSGetSNES()`

2714: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2715: @*/
2716: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2717: {
2718:   PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);

2720:   PetscFunctionBegin;
2723:   PetscCall(PetscObjectReference((PetscObject)snes));
2724:   PetscCall(SNESDestroy(&ts->snes));

2726:   ts->snes = snes;

2728:   PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2729:   PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2730:   if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2731:   PetscFunctionReturn(PETSC_SUCCESS);
2732: }

2734: /*@
2735:   TSGetKSP - Returns the `KSP` (linear solver) associated with
2736:   a `TS` (timestepper) context.

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

2740:   Input Parameter:
2741: . ts - the `TS` context obtained from `TSCreate()`

2743:   Output Parameter:
2744: . ksp - the nonlinear solver context

2746:   Level: beginner

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

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

2756: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2757: @*/
2758: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2759: {
2760:   SNES snes;

2762:   PetscFunctionBegin;
2764:   PetscAssertPointer(ksp, 2);
2765:   PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2766:   PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
2767:   PetscCall(TSGetSNES(ts, &snes));
2768:   PetscCall(SNESGetKSP(snes, ksp));
2769:   PetscFunctionReturn(PETSC_SUCCESS);
2770: }

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

2774: /*@
2775:   TSSetMaxSteps - Sets the maximum number of steps to use.

2777:   Logically Collective

2779:   Input Parameters:
2780: + ts       - the `TS` context obtained from `TSCreate()`
2781: - maxsteps - maximum number of steps to use

2783:   Options Database Key:
2784: . -ts_max_steps <maxsteps> - Sets maxsteps

2786:   Level: intermediate

2788:   Note:
2789:   The default maximum number of steps is 5000

2791: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2792: @*/
2793: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2794: {
2795:   PetscFunctionBegin;
2798:   PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2799:   ts->max_steps = maxsteps;
2800:   PetscFunctionReturn(PETSC_SUCCESS);
2801: }

2803: /*@
2804:   TSGetMaxSteps - Gets the maximum number of steps to use.

2806:   Not Collective

2808:   Input Parameter:
2809: . ts - the `TS` context obtained from `TSCreate()`

2811:   Output Parameter:
2812: . maxsteps - maximum number of steps to use

2814:   Level: advanced

2816: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2817: @*/
2818: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2819: {
2820:   PetscFunctionBegin;
2822:   PetscAssertPointer(maxsteps, 2);
2823:   *maxsteps = ts->max_steps;
2824:   PetscFunctionReturn(PETSC_SUCCESS);
2825: }

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

2830:   Logically Collective

2832:   Input Parameters:
2833: + ts      - the `TS` context obtained from `TSCreate()`
2834: - maxtime - final time to step to

2836:   Options Database Key:
2837: . -ts_max_time <maxtime> - Sets maxtime

2839:   Level: intermediate

2841:   Notes:
2842:   The default maximum time is 5.0

2844: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2845: @*/
2846: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2847: {
2848:   PetscFunctionBegin;
2851:   ts->max_time = maxtime;
2852:   PetscFunctionReturn(PETSC_SUCCESS);
2853: }

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

2858:   Not Collective

2860:   Input Parameter:
2861: . ts - the `TS` context obtained from `TSCreate()`

2863:   Output Parameter:
2864: . maxtime - final time to step to

2866:   Level: advanced

2868: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2869: @*/
2870: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2871: {
2872:   PetscFunctionBegin;
2874:   PetscAssertPointer(maxtime, 2);
2875:   *maxtime = ts->max_time;
2876:   PetscFunctionReturn(PETSC_SUCCESS);
2877: }

2879: // PetscClangLinter pragma disable: -fdoc-*
2880: /*@
2881:   TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.

2883:   Level: deprecated

2885: @*/
2886: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2887: {
2888:   PetscFunctionBegin;
2890:   PetscCall(TSSetTime(ts, initial_time));
2891:   PetscCall(TSSetTimeStep(ts, time_step));
2892:   PetscFunctionReturn(PETSC_SUCCESS);
2893: }

2895: // PetscClangLinter pragma disable: -fdoc-*
2896: /*@
2897:   TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.

2899:   Level: deprecated

2901: @*/
2902: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2903: {
2904:   PetscFunctionBegin;
2906:   if (maxsteps) {
2907:     PetscAssertPointer(maxsteps, 2);
2908:     *maxsteps = ts->max_steps;
2909:   }
2910:   if (maxtime) {
2911:     PetscAssertPointer(maxtime, 3);
2912:     *maxtime = ts->max_time;
2913:   }
2914:   PetscFunctionReturn(PETSC_SUCCESS);
2915: }

2917: // PetscClangLinter pragma disable: -fdoc-*
2918: /*@
2919:   TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.

2921:   Level: deprecated

2923: @*/
2924: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2925: {
2926:   PetscFunctionBegin;
2930:   if (maxsteps >= 0) ts->max_steps = maxsteps;
2931:   if (maxtime != (PetscReal)PETSC_DEFAULT) ts->max_time = maxtime;
2932:   PetscFunctionReturn(PETSC_SUCCESS);
2933: }

2935: // PetscClangLinter pragma disable: -fdoc-*
2936: /*@
2937:   TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.

2939:   Level: deprecated

2941: @*/
2942: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
2943: {
2944:   return TSGetStepNumber(ts, steps);
2945: }

2947: // PetscClangLinter pragma disable: -fdoc-*
2948: /*@
2949:   TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.

2951:   Level: deprecated

2953: @*/
2954: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
2955: {
2956:   return TSGetStepNumber(ts, steps);
2957: }

2959: /*@
2960:   TSSetSolution - Sets the initial solution vector
2961:   for use by the `TS` routines.

2963:   Logically Collective

2965:   Input Parameters:
2966: + ts - the `TS` context obtained from `TSCreate()`
2967: - u  - the solution vector

2969:   Level: beginner

2971: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
2972: @*/
2973: PetscErrorCode TSSetSolution(TS ts, Vec u)
2974: {
2975:   DM dm;

2977:   PetscFunctionBegin;
2980:   PetscCall(PetscObjectReference((PetscObject)u));
2981:   PetscCall(VecDestroy(&ts->vec_sol));
2982:   ts->vec_sol = u;

2984:   PetscCall(TSGetDM(ts, &dm));
2985:   PetscCall(DMShellSetGlobalVector(dm, u));
2986:   PetscFunctionReturn(PETSC_SUCCESS);
2987: }

2989: /*@C
2990:   TSSetPreStep - Sets the general-purpose function
2991:   called once at the beginning of each time step.

2993:   Logically Collective

2995:   Input Parameters:
2996: + ts   - The `TS` context obtained from `TSCreate()`
2997: - func - The function

2999:   Calling sequence of `func`:
3000: . ts - the `TS` context

3002:   Level: intermediate

3004: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3005: @*/
3006: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
3007: {
3008:   PetscFunctionBegin;
3010:   ts->prestep = func;
3011:   PetscFunctionReturn(PETSC_SUCCESS);
3012: }

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

3017:   Collective

3019:   Input Parameter:
3020: . ts - The `TS` context obtained from `TSCreate()`

3022:   Level: developer

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

3028: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3029: @*/
3030: PetscErrorCode TSPreStep(TS ts)
3031: {
3032:   PetscFunctionBegin;
3034:   if (ts->prestep) {
3035:     Vec              U;
3036:     PetscObjectId    idprev;
3037:     PetscBool        sameObject;
3038:     PetscObjectState sprev, spost;

3040:     PetscCall(TSGetSolution(ts, &U));
3041:     PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3042:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3043:     PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3044:     PetscCall(TSGetSolution(ts, &U));
3045:     PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3046:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3047:     if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3048:   }
3049:   PetscFunctionReturn(PETSC_SUCCESS);
3050: }

3052: /*@C
3053:   TSSetPreStage - Sets the general-purpose function
3054:   called once at the beginning of each stage.

3056:   Logically Collective

3058:   Input Parameters:
3059: + ts   - The `TS` context obtained from `TSCreate()`
3060: - func - The function

3062:   Calling sequence of `func`:
3063: + ts        - the `TS` context
3064: - stagetime - the stage time

3066:   Level: intermediate

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

3073: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3074: @*/
3075: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime))
3076: {
3077:   PetscFunctionBegin;
3079:   ts->prestage = func;
3080:   PetscFunctionReturn(PETSC_SUCCESS);
3081: }

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

3087:   Logically Collective

3089:   Input Parameters:
3090: + ts   - The `TS` context obtained from `TSCreate()`
3091: - func - The function

3093:   Calling sequence of `func`:
3094: + ts         - the `TS` context
3095: . stagetime  - the stage time
3096: . stageindex - the stage index
3097: - Y          - Array of vectors (of size = total number of stages) with the stage solutions

3099:   Level: intermediate

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

3106: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3107: @*/
3108: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
3109: {
3110:   PetscFunctionBegin;
3112:   ts->poststage = func;
3113:   PetscFunctionReturn(PETSC_SUCCESS);
3114: }

3116: /*@C
3117:   TSSetPostEvaluate - Sets the general-purpose function
3118:   called once at the end of each step evaluation.

3120:   Logically Collective

3122:   Input Parameters:
3123: + ts   - The `TS` context obtained from `TSCreate()`
3124: - func - The function

3126:   Calling sequence of `func`:
3127: . ts - the `TS` context

3129:   Level: intermediate

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

3138: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3139: @*/
3140: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
3141: {
3142:   PetscFunctionBegin;
3144:   ts->postevaluate = func;
3145:   PetscFunctionReturn(PETSC_SUCCESS);
3146: }

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

3151:   Collective

3153:   Input Parameters:
3154: + ts        - The `TS` context obtained from `TSCreate()`
3155: - stagetime - The absolute time of the current stage

3157:   Level: developer

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

3163: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3164: @*/
3165: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3166: {
3167:   PetscFunctionBegin;
3169:   if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3170:   PetscFunctionReturn(PETSC_SUCCESS);
3171: }

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

3176:   Collective

3178:   Input Parameters:
3179: + ts         - The `TS` context obtained from `TSCreate()`
3180: . stagetime  - The absolute time of the current stage
3181: . stageindex - Stage number
3182: - Y          - Array of vectors (of size = total number of stages) with the stage solutions

3184:   Level: developer

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

3190: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3191: @*/
3192: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3193: {
3194:   PetscFunctionBegin;
3196:   if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3197:   PetscFunctionReturn(PETSC_SUCCESS);
3198: }

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

3203:   Collective

3205:   Input Parameter:
3206: . ts - The `TS` context obtained from `TSCreate()`

3208:   Level: developer

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

3214: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3215: @*/
3216: PetscErrorCode TSPostEvaluate(TS ts)
3217: {
3218:   PetscFunctionBegin;
3220:   if (ts->postevaluate) {
3221:     Vec              U;
3222:     PetscObjectState sprev, spost;

3224:     PetscCall(TSGetSolution(ts, &U));
3225:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3226:     PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3227:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3228:     if (sprev != spost) PetscCall(TSRestartStep(ts));
3229:   }
3230:   PetscFunctionReturn(PETSC_SUCCESS);
3231: }

3233: /*@C
3234:   TSSetPostStep - Sets the general-purpose function
3235:   called once at the end of each time step.

3237:   Logically Collective

3239:   Input Parameters:
3240: + ts   - The `TS` context obtained from `TSCreate()`
3241: - func - The function

3243:   Calling sequence of `func`:
3244: . ts - the `TS` context

3246:   Level: intermediate

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

3253: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3254: @*/
3255: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
3256: {
3257:   PetscFunctionBegin;
3259:   ts->poststep = func;
3260:   PetscFunctionReturn(PETSC_SUCCESS);
3261: }

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

3266:   Collective

3268:   Input Parameter:
3269: . ts - The `TS` context obtained from `TSCreate()`

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

3275:   Level: developer

3277: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()`
3278: @*/
3279: PetscErrorCode TSPostStep(TS ts)
3280: {
3281:   PetscFunctionBegin;
3283:   if (ts->poststep) {
3284:     Vec              U;
3285:     PetscObjectId    idprev;
3286:     PetscBool        sameObject;
3287:     PetscObjectState sprev, spost;

3289:     PetscCall(TSGetSolution(ts, &U));
3290:     PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3291:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3292:     PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3293:     PetscCall(TSGetSolution(ts, &U));
3294:     PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3295:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3296:     if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3297:   }
3298:   PetscFunctionReturn(PETSC_SUCCESS);
3299: }

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

3304:   Collective

3306:   Input Parameters:
3307: + ts - time stepping context
3308: - t  - time to interpolate to

3310:   Output Parameter:
3311: . U - state at given time

3313:   Level: intermediate

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

3318: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3319: @*/
3320: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3321: {
3322:   PetscFunctionBegin;
3325:   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);
3326:   PetscUseTypeMethod(ts, interpolate, t, U);
3327:   PetscFunctionReturn(PETSC_SUCCESS);
3328: }

3330: /*@
3331:   TSStep - Steps one time step

3333:   Collective

3335:   Input Parameter:
3336: . ts - the `TS` context obtained from `TSCreate()`

3338:   Level: developer

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

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

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

3349: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3350: @*/
3351: PetscErrorCode TSStep(TS ts)
3352: {
3353:   static PetscBool cite = PETSC_FALSE;
3354:   PetscReal        ptime;

3356:   PetscFunctionBegin;
3358:   PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3359:                                    "  title         = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3360:                                    "  author        = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3361:                                    "  journal       = {arXiv e-preprints},\n"
3362:                                    "  eprint        = {1806.01437},\n"
3363:                                    "  archivePrefix = {arXiv},\n"
3364:                                    "  year          = {2018}\n}\n",
3365:                                    &cite));
3366:   PetscCall(TSSetUp(ts));
3367:   PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3368:   if (ts->tspan)
3369:     ts->tspan->worktol = 0; /* In each step of TSSolve() 'tspan->worktol' will be meaningfully defined (later) only once:
3370:                                                    in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */

3372:   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>");
3373:   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()");
3374:   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");

3376:   if (!ts->steps) ts->ptime_prev = ts->ptime;
3377:   ptime                   = ts->ptime;
3378:   ts->ptime_prev_rollback = ts->ptime_prev;
3379:   ts->reason              = TS_CONVERGED_ITERATING;

3381:   PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3382:   PetscUseTypeMethod(ts, step);
3383:   PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));

3385:   if (ts->reason >= 0) {
3386:     ts->ptime_prev = ptime;
3387:     ts->steps++;
3388:     ts->steprollback = PETSC_FALSE;
3389:     ts->steprestart  = PETSC_FALSE;
3390:   }
3391:   if (!ts->reason) {
3392:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3393:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3394:   }

3396:   if (ts->reason < 0 && ts->errorifstepfailed) {
3397:     PetscCall(TSMonitorCancel(ts));
3398:     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]);
3399:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3400:   }
3401:   PetscFunctionReturn(PETSC_SUCCESS);
3402: }

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

3408:   Collective

3410:   Input Parameters:
3411: + ts        - time stepping context
3412: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

3414:   Input/Output Parameter:
3415: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3416:            on output, the actual order of the error evaluation

3418:   Output Parameter:
3419: . wlte - the weighted local truncation error norm

3421:   Level: advanced

3423:   Note:
3424:   If the timestepper cannot evaluate the error in a particular step
3425:   (eg. in the first step or restart steps after event handling),
3426:   this routine returns wlte=-1.0 .

3428: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3429: @*/
3430: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3431: {
3432:   PetscFunctionBegin;
3436:   if (order) PetscAssertPointer(order, 3);
3438:   PetscAssertPointer(wlte, 4);
3439:   PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3440:   PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3441:   PetscFunctionReturn(PETSC_SUCCESS);
3442: }

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

3447:   Collective

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

3454:   Output Parameter:
3455: . U - state at the end of the current step

3457:   Level: advanced

3459:   Notes:
3460:   This function cannot be called until all stages have been evaluated.

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

3464: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3465: @*/
3466: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3467: {
3468:   PetscFunctionBegin;
3472:   PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3473:   PetscFunctionReturn(PETSC_SUCCESS);
3474: }

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

3479:   Not collective

3481:   Input Parameter:
3482: . ts - time stepping context

3484:   Output Parameter:
3485: . initCondition - The function which computes an initial condition

3487:   Calling sequence of `initCondition`:
3488: + ts - The timestepping context
3489: - u  - The input vector in which the initial condition is stored

3491:   Level: advanced

3493: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3494: @*/
3495: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3496: {
3497:   PetscFunctionBegin;
3499:   PetscAssertPointer(initCondition, 2);
3500:   *initCondition = ts->ops->initcondition;
3501:   PetscFunctionReturn(PETSC_SUCCESS);
3502: }

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

3507:   Logically collective

3509:   Input Parameters:
3510: + ts            - time stepping context
3511: - initCondition - The function which computes an initial condition

3513:   Calling sequence of `initCondition`:
3514: + ts - The timestepping context
3515: - e  - The input vector in which the initial condition is to be stored

3517:   Level: advanced

3519: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3520: @*/
3521: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3522: {
3523:   PetscFunctionBegin;
3526:   ts->ops->initcondition = initCondition;
3527:   PetscFunctionReturn(PETSC_SUCCESS);
3528: }

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

3533:   Collective

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

3539:   Level: advanced

3541: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3542: @*/
3543: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3544: {
3545:   PetscFunctionBegin;
3548:   PetscTryTypeMethod(ts, initcondition, u);
3549:   PetscFunctionReturn(PETSC_SUCCESS);
3550: }

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

3555:   Not collective

3557:   Input Parameter:
3558: . ts - time stepping context

3560:   Output Parameter:
3561: . exactError - The function which computes the solution error

3563:   Calling sequence of `exactError`:
3564: + ts - The timestepping context
3565: . u  - The approximate solution vector
3566: - e  - The vector in which the error is stored

3568:   Level: advanced

3570: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3571: @*/
3572: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3573: {
3574:   PetscFunctionBegin;
3576:   PetscAssertPointer(exactError, 2);
3577:   *exactError = ts->ops->exacterror;
3578:   PetscFunctionReturn(PETSC_SUCCESS);
3579: }

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

3584:   Logically collective

3586:   Input Parameters:
3587: + ts         - time stepping context
3588: - exactError - The function which computes the solution error

3590:   Calling sequence of `exactError`:
3591: + ts - The timestepping context
3592: . u  - The approximate solution vector
3593: - e  - The  vector in which the error is stored

3595:   Level: advanced

3597: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3598: @*/
3599: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3600: {
3601:   PetscFunctionBegin;
3604:   ts->ops->exacterror = exactError;
3605:   PetscFunctionReturn(PETSC_SUCCESS);
3606: }

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

3611:   Collective

3613:   Input Parameters:
3614: + ts - time stepping context
3615: . u  - The approximate solution
3616: - e  - The `Vec` used to store the error

3618:   Level: advanced

3620: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3621: @*/
3622: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3623: {
3624:   PetscFunctionBegin;
3628:   PetscTryTypeMethod(ts, exacterror, u, e);
3629:   PetscFunctionReturn(PETSC_SUCCESS);
3630: }

3632: /*@C
3633:   TSSetResize - Sets the resize callbacks.

3635:   Logically Collective

3637:   Input Parameters:
3638: + ts       - The `TS` context obtained from `TSCreate()`
3639: . setup    - The setup function
3640: . transfer - The transfer function
3641: - ctx      - [optional] The user-defined context

3643:   Calling sequence of `setup`:
3644: + ts     - the `TS` context
3645: . step   - the current step
3646: . time   - the current time
3647: . state  - the current vector of state
3648: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3649: - ctx    - user defined context

3651:   Calling sequence of `transfer`:
3652: + ts      - the `TS` context
3653: . nv      - the number of vectors to be transferred
3654: . vecsin  - array of vectors to be transferred
3655: . vecsout - array of transferred vectors
3656: - ctx     - user defined context

3658:   Notes:
3659:   The `setup` function is called inside `TSSolve()` after `TSPostStep()` at the end of each time step
3660:   to determine if the problem size has changed.
3661:   If it is the case, the solver will collect the needed vectors that need to be
3662:   transferred from the old to the new sizes using `transfer`. These vectors will include the current
3663:   solution vector, and other vectors needed by the specific solver used.
3664:   For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3665:   Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3666:   will be automatically reset if the sizes are changed and they must be specified again by the user
3667:   inside the `transfer` function.
3668:   The input and output arrays passed to `transfer` are allocated by PETSc.
3669:   Vectors in `vecsout` must be created by the user.
3670:   Ownership of vectors in `vecsout` is transferred to PETSc.

3672:   Level: advanced

3674: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3675: @*/
3676: PetscErrorCode TSSetResize(TS ts, PetscErrorCode (*setup)(TS ts, PetscInt step, PetscReal time, Vec state, PetscBool *resize, void *ctx), PetscErrorCode (*transfer)(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx), void *ctx)
3677: {
3678:   PetscFunctionBegin;
3680:   ts->resizesetup    = setup;
3681:   ts->resizetransfer = transfer;
3682:   ts->resizectx      = ctx;
3683:   PetscFunctionReturn(PETSC_SUCCESS);
3684: }

3686: /*
3687:   TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.

3689:   Collective

3691:   Input Parameters:
3692: + ts   - The `TS` context obtained from `TSCreate()`
3693: - flg - If `PETSC_TRUE` each TS implementation (e.g. `TSBDF`) will register vectors to be transferred, if `PETSC_FALSE` vectors will be imported from transferred vectors.

3695:   Level: developer

3697:   Note:
3698:   `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3699:    used within time stepping implementations,
3700:    so most users would not generally call this routine themselves.

3702: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3703: @*/
3704: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3705: {
3706:   PetscFunctionBegin;
3708:   PetscTryTypeMethod(ts, resizeregister, flg);
3709:   /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3710:   PetscFunctionReturn(PETSC_SUCCESS);
3711: }

3713: static PetscErrorCode TSResizeReset(TS ts)
3714: {
3715:   PetscFunctionBegin;
3717:   PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3718:   PetscFunctionReturn(PETSC_SUCCESS);
3719: }

3721: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3722: {
3723:   PetscFunctionBegin;
3726:   for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3727:   if (ts->resizetransfer) {
3728:     PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3729:     PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3730:   }
3731:   for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3732:   PetscFunctionReturn(PETSC_SUCCESS);
3733: }

3735: /*@C
3736:   TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.

3738:   Collective

3740:   Input Parameters:
3741: + ts   - The `TS` context obtained from `TSCreate()`
3742: . name - A string identifying the vector
3743: - vec  - The vector

3745:   Level: developer

3747:   Note:
3748:   `TSResizeRegisterVec()` is typically used within time stepping implementations,
3749:   so most users would not generally call this routine themselves.

3751: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3752: @*/
3753: PetscErrorCode TSResizeRegisterVec(TS ts, const char *name, Vec vec)
3754: {
3755:   PetscFunctionBegin;
3757:   PetscAssertPointer(name, 2);
3759:   PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3760:   PetscFunctionReturn(PETSC_SUCCESS);
3761: }

3763: /*@C
3764:   TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.

3766:   Collective

3768:   Input Parameters:
3769: + ts   - The `TS` context obtained from `TSCreate()`
3770: . name - A string identifying the vector
3771: - vec  - The vector

3773:   Level: developer

3775:   Note:
3776:   `TSResizeRetrieveVec()` is typically used within time stepping implementations,
3777:   so most users would not generally call this routine themselves.

3779: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
3780: @*/
3781: PetscErrorCode TSResizeRetrieveVec(TS ts, const char *name, Vec *vec)
3782: {
3783:   PetscFunctionBegin;
3785:   PetscAssertPointer(name, 2);
3786:   PetscAssertPointer(vec, 3);
3787:   PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
3788:   PetscFunctionReturn(PETSC_SUCCESS);
3789: }

3791: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
3792: {
3793:   PetscInt        cnt;
3794:   PetscObjectList tmp;
3795:   Vec            *vecsin  = NULL;
3796:   const char    **namesin = NULL;

3798:   PetscFunctionBegin;
3799:   for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
3800:     if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
3801:   if (names) PetscCall(PetscMalloc1(cnt, &namesin));
3802:   if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
3803:   for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
3804:     if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
3805:       if (vecs) vecsin[cnt] = (Vec)tmp->obj;
3806:       if (names) namesin[cnt] = tmp->name;
3807:       cnt++;
3808:     }
3809:   }
3810:   if (nv) *nv = cnt;
3811:   if (names) *names = namesin;
3812:   if (vecs) *vecs = vecsin;
3813:   PetscFunctionReturn(PETSC_SUCCESS);
3814: }

3816: /*@
3817:   TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`

3819:   Collective

3821:   Input Parameter:
3822: . ts - The `TS` context obtained from `TSCreate()`

3824:   Level: developer

3826:   Note:
3827:   `TSResize()` is typically used within time stepping implementations,
3828:   so most users would not generally call this routine themselves.

3830: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3831: @*/
3832: PetscErrorCode TSResize(TS ts)
3833: {
3834:   PetscInt     nv      = 0;
3835:   const char **names   = NULL;
3836:   Vec         *vecsin  = NULL;
3837:   const char  *solname = "ts:vec_sol";

3839:   PetscFunctionBegin;
3841:   if (ts->resizesetup) {
3842:     PetscBool flg = PETSC_FALSE;

3844:     PetscCall(VecLockReadPush(ts->vec_sol));
3845:     PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &flg, ts->resizectx));
3846:     PetscCall(VecLockReadPop(ts->vec_sol));
3847:     if (flg) {
3848:       PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
3849:       PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
3850:     }
3851:   }

3853:   PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
3854:   if (nv) {
3855:     Vec *vecsout, vecsol;

3857:     /* Reset internal objects */
3858:     PetscCall(TSReset(ts));

3860:     /* Transfer needed vectors (users can call SetJacobian, SetDM here) */
3861:     PetscCall(PetscCalloc1(nv, &vecsout));
3862:     PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
3863:     for (PetscInt i = 0; i < nv; i++) {
3864:       PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
3865:       PetscCall(VecDestroy(&vecsout[i]));
3866:     }
3867:     PetscCall(PetscFree(vecsout));
3868:     PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */

3870:     PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
3871:     if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
3872:     PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
3873:   }

3875:   PetscCall(PetscFree(names));
3876:   PetscCall(PetscFree(vecsin));
3877:   PetscCall(TSResizeReset(ts));
3878:   PetscFunctionReturn(PETSC_SUCCESS);
3879: }

3881: /*@
3882:   TSSolve - Steps the requested number of timesteps.

3884:   Collective

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

3891:   Level: beginner

3893:   Notes:
3894:   The final time returned by this function may be different from the time of the internally
3895:   held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
3896:   stepped over the final time.

3898: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3899: @*/
3900: PetscErrorCode TSSolve(TS ts, Vec u)
3901: {
3902:   Vec solution;

3904:   PetscFunctionBegin;

3908:   PetscCall(TSSetExactFinalTimeDefault(ts));
3909:   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 */
3910:     if (!ts->vec_sol || u == ts->vec_sol) {
3911:       PetscCall(VecDuplicate(u, &solution));
3912:       PetscCall(TSSetSolution(ts, solution));
3913:       PetscCall(VecDestroy(&solution)); /* grant ownership */
3914:     }
3915:     PetscCall(VecCopy(u, ts->vec_sol));
3916:     PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
3917:   } else if (u) PetscCall(TSSetSolution(ts, u));
3918:   PetscCall(TSSetUp(ts));
3919:   PetscCall(TSTrajectorySetUp(ts->trajectory, ts));

3921:   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>");
3922:   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()");
3923:   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");
3924:   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");

3926:   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 */
3927:     PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]));
3928:     ts->tspan->spanctr = 1;
3929:   }

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

3933:   /* reset number of steps only when the step is not restarted. ARKIMEX
3934:      restarts the step after an event. Resetting these counters in such case causes
3935:      TSTrajectory to incorrectly save the output files
3936:   */
3937:   /* reset time step and iteration counters */
3938:   if (!ts->steps) {
3939:     ts->ksp_its           = 0;
3940:     ts->snes_its          = 0;
3941:     ts->num_snes_failures = 0;
3942:     ts->reject            = 0;
3943:     ts->steprestart       = PETSC_TRUE;
3944:     ts->steprollback      = PETSC_FALSE;
3945:     ts->rhsjacobian.time  = PETSC_MIN_REAL;
3946:   }

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

3953:     if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
3954:     else maxdt = ts->max_time - ts->ptime;
3955:     ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
3956:   }
3957:   ts->reason = TS_CONVERGED_ITERATING;

3959:   {
3960:     PetscViewer       viewer;
3961:     PetscViewerFormat format;
3962:     PetscBool         flg;
3963:     static PetscBool  incall = PETSC_FALSE;

3965:     if (!incall) {
3966:       /* Estimate the convergence rate of the time discretization */
3967:       PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
3968:       if (flg) {
3969:         PetscConvEst conv;
3970:         DM           dm;
3971:         PetscReal   *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
3972:         PetscInt     Nf;
3973:         PetscBool    checkTemporal = PETSC_TRUE;

3975:         incall = PETSC_TRUE;
3976:         PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
3977:         PetscCall(TSGetDM(ts, &dm));
3978:         PetscCall(DMGetNumFields(dm, &Nf));
3979:         PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
3980:         PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
3981:         PetscCall(PetscConvEstUseTS(conv, checkTemporal));
3982:         PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
3983:         PetscCall(PetscConvEstSetFromOptions(conv));
3984:         PetscCall(PetscConvEstSetUp(conv));
3985:         PetscCall(PetscConvEstGetConvRate(conv, alpha));
3986:         PetscCall(PetscViewerPushFormat(viewer, format));
3987:         PetscCall(PetscConvEstRateView(conv, alpha, viewer));
3988:         PetscCall(PetscViewerPopFormat(viewer));
3989:         PetscCall(PetscOptionsRestoreViewer(&viewer));
3990:         PetscCall(PetscConvEstDestroy(&conv));
3991:         PetscCall(PetscFree(alpha));
3992:         incall = PETSC_FALSE;
3993:       }
3994:     }
3995:   }

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

3999:   if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4000:     PetscUseTypeMethod(ts, solve);
4001:     if (u) PetscCall(VecCopy(ts->vec_sol, u));
4002:     ts->solvetime = ts->ptime;
4003:     solution      = ts->vec_sol;
4004:   } else { /* Step the requested number of timesteps. */
4005:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4006:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;

4008:     if (!ts->steps) {
4009:       PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4010:       PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4011:     }

4013:     while (!ts->reason) {
4014:       PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4015:       if (!ts->steprollback) PetscCall(TSPreStep(ts));
4016:       PetscCall(TSStep(ts));
4017:       if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4018:       if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4019:       if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4020:         if (ts->reason >= 0) ts->steps--;            /* Revert the step number changed by TSStep() */
4021:         PetscCall(TSForwardCostIntegral(ts));
4022:         if (ts->reason >= 0) ts->steps++;
4023:       }
4024:       if (ts->forward_solve) {            /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4025:         if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4026:         PetscCall(TSForwardStep(ts));
4027:         if (ts->reason >= 0) ts->steps++;
4028:       }
4029:       PetscCall(TSPostEvaluate(ts));
4030:       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. */
4031:       if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4032:       if (!ts->steprollback) {
4033:         PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4034:         PetscCall(TSPostStep(ts));
4035:         PetscCall(TSResize(ts));

4037:         if (ts->tspan && ts->tspan->spanctr < ts->tspan->num_span_times) {
4038:           PetscCheck(ts->tspan->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(tspan->worktol > 0) in TSSolve()");
4039:           if (PetscIsCloseAtTol(ts->ptime, ts->tspan->span_times[ts->tspan->spanctr], ts->tspan->worktol, 0)) PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[ts->tspan->spanctr++]));
4040:         }
4041:       }
4042:     }
4043:     PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));

4045:     if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4046:       if (!u) u = ts->vec_sol;
4047:       PetscCall(TSInterpolate(ts, ts->max_time, u));
4048:       ts->solvetime = ts->max_time;
4049:       solution      = u;
4050:       PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4051:     } else {
4052:       if (u) PetscCall(VecCopy(ts->vec_sol, u));
4053:       ts->solvetime = ts->ptime;
4054:       solution      = ts->vec_sol;
4055:     }
4056:   }

4058:   PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4059:   PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4060:   PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4061:   if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4062:   PetscFunctionReturn(PETSC_SUCCESS);
4063: }

4065: /*@
4066:   TSGetTime - Gets the time of the most recently completed step.

4068:   Not Collective

4070:   Input Parameter:
4071: . ts - the `TS` context obtained from `TSCreate()`

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

4076:   Level: beginner

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

4082: .seealso: [](ch_ts), `TS`, ``TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4083: @*/
4084: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4085: {
4086:   PetscFunctionBegin;
4088:   PetscAssertPointer(t, 2);
4089:   *t = ts->ptime;
4090:   PetscFunctionReturn(PETSC_SUCCESS);
4091: }

4093: /*@
4094:   TSGetPrevTime - Gets the starting time of the previously completed step.

4096:   Not Collective

4098:   Input Parameter:
4099: . ts - the `TS` context obtained from `TSCreate()`

4101:   Output Parameter:
4102: . t - the previous time

4104:   Level: beginner

4106: .seealso: [](ch_ts), `TS`, ``TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4107: @*/
4108: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4109: {
4110:   PetscFunctionBegin;
4112:   PetscAssertPointer(t, 2);
4113:   *t = ts->ptime_prev;
4114:   PetscFunctionReturn(PETSC_SUCCESS);
4115: }

4117: /*@
4118:   TSSetTime - Allows one to reset the time.

4120:   Logically Collective

4122:   Input Parameters:
4123: + ts - the `TS` context obtained from `TSCreate()`
4124: - t  - the time

4126:   Level: intermediate

4128: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4129: @*/
4130: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4131: {
4132:   PetscFunctionBegin;
4135:   ts->ptime = t;
4136:   PetscFunctionReturn(PETSC_SUCCESS);
4137: }

4139: /*@C
4140:   TSSetOptionsPrefix - Sets the prefix used for searching for all
4141:   TS options in the database.

4143:   Logically Collective

4145:   Input Parameters:
4146: + ts     - The `TS` context
4147: - prefix - The prefix to prepend to all option names

4149:   Level: advanced

4151:   Note:
4152:   A hyphen (-) must NOT be given at the beginning of the prefix name.
4153:   The first character of all runtime options is AUTOMATICALLY the
4154:   hyphen.

4156: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4157: @*/
4158: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4159: {
4160:   SNES snes;

4162:   PetscFunctionBegin;
4164:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4165:   PetscCall(TSGetSNES(ts, &snes));
4166:   PetscCall(SNESSetOptionsPrefix(snes, prefix));
4167:   PetscFunctionReturn(PETSC_SUCCESS);
4168: }

4170: /*@C
4171:   TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4172:   TS options in the database.

4174:   Logically Collective

4176:   Input Parameters:
4177: + ts     - The `TS` context
4178: - prefix - The prefix to prepend to all option names

4180:   Level: advanced

4182:   Note:
4183:   A hyphen (-) must NOT be given at the beginning of the prefix name.
4184:   The first character of all runtime options is AUTOMATICALLY the
4185:   hyphen.

4187: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4188: @*/
4189: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4190: {
4191:   SNES snes;

4193:   PetscFunctionBegin;
4195:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4196:   PetscCall(TSGetSNES(ts, &snes));
4197:   PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4198:   PetscFunctionReturn(PETSC_SUCCESS);
4199: }

4201: /*@C
4202:   TSGetOptionsPrefix - Sets the prefix used for searching for all
4203:   `TS` options in the database.

4205:   Not Collective

4207:   Input Parameter:
4208: . ts - The `TS` context

4210:   Output Parameter:
4211: . prefix - A pointer to the prefix string used

4213:   Level: intermediate

4215:   Fortran Notes:
4216:   The user should pass in a string 'prefix' of
4217:   sufficient length to hold the prefix.

4219: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4220: @*/
4221: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4222: {
4223:   PetscFunctionBegin;
4225:   PetscAssertPointer(prefix, 2);
4226:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4227:   PetscFunctionReturn(PETSC_SUCCESS);
4228: }

4230: /*@C
4231:   TSGetRHSJacobian - Returns the Jacobian J at the present timestep.

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

4235:   Input Parameter:
4236: . ts - The `TS` context obtained from `TSCreate()`

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

4244:   Level: intermediate

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

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

4251: @*/
4252: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, void **ctx)
4253: {
4254:   DM dm;

4256:   PetscFunctionBegin;
4257:   if (Amat || Pmat) {
4258:     SNES snes;
4259:     PetscCall(TSGetSNES(ts, &snes));
4260:     PetscCall(SNESSetUpMatrices(snes));
4261:     PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4262:   }
4263:   PetscCall(TSGetDM(ts, &dm));
4264:   PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4265:   PetscFunctionReturn(PETSC_SUCCESS);
4266: }

4268: /*@C
4269:   TSGetIJacobian - Returns the implicit Jacobian at the present timestep.

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

4273:   Input Parameter:
4274: . ts - The `TS` context obtained from `TSCreate()`

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

4282:   Level: advanced

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

4287: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4288: @*/
4289: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, void **ctx)
4290: {
4291:   DM dm;

4293:   PetscFunctionBegin;
4294:   if (Amat || Pmat) {
4295:     SNES snes;
4296:     PetscCall(TSGetSNES(ts, &snes));
4297:     PetscCall(SNESSetUpMatrices(snes));
4298:     PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4299:   }
4300:   PetscCall(TSGetDM(ts, &dm));
4301:   PetscCall(DMTSGetIJacobian(dm, f, ctx));
4302:   PetscFunctionReturn(PETSC_SUCCESS);
4303: }

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

4309:   Logically Collective

4311:   Input Parameters:
4312: + ts - the `TS` integrator object
4313: - dm - the dm, cannot be `NULL`

4315:   Level: intermediate

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

4322: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4323: @*/
4324: PetscErrorCode TSSetDM(TS ts, DM dm)
4325: {
4326:   SNES snes;
4327:   DMTS tsdm;

4329:   PetscFunctionBegin;
4332:   PetscCall(PetscObjectReference((PetscObject)dm));
4333:   if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4334:     if (ts->dm->dmts && !dm->dmts) {
4335:       PetscCall(DMCopyDMTS(ts->dm, dm));
4336:       PetscCall(DMGetDMTS(ts->dm, &tsdm));
4337:       /* Grant write privileges to the replacement DM */
4338:       if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4339:     }
4340:     PetscCall(DMDestroy(&ts->dm));
4341:   }
4342:   ts->dm = dm;

4344:   PetscCall(TSGetSNES(ts, &snes));
4345:   PetscCall(SNESSetDM(snes, dm));
4346:   PetscFunctionReturn(PETSC_SUCCESS);
4347: }

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

4352:   Not Collective

4354:   Input Parameter:
4355: . ts - the `TS`

4357:   Output Parameter:
4358: . dm - the `DM`

4360:   Level: intermediate

4362: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4363: @*/
4364: PetscErrorCode TSGetDM(TS ts, DM *dm)
4365: {
4366:   PetscFunctionBegin;
4368:   if (!ts->dm) {
4369:     PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4370:     if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4371:   }
4372:   *dm = ts->dm;
4373:   PetscFunctionReturn(PETSC_SUCCESS);
4374: }

4376: /*@
4377:   SNESTSFormFunction - Function to evaluate nonlinear residual

4379:   Logically Collective

4381:   Input Parameters:
4382: + snes - nonlinear solver
4383: . U    - the current state at which to evaluate the residual
4384: - ctx  - user context, must be a TS

4386:   Output Parameter:
4387: . F - the nonlinear residual

4389:   Level: advanced

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

4395: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4396: @*/
4397: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4398: {
4399:   TS ts = (TS)ctx;

4401:   PetscFunctionBegin;
4406:   PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4407:   PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4408:   PetscFunctionReturn(PETSC_SUCCESS);
4409: }

4411: /*@
4412:   SNESTSFormJacobian - Function to evaluate the Jacobian

4414:   Collective

4416:   Input Parameters:
4417: + snes - nonlinear solver
4418: . U    - the current state at which to evaluate the residual
4419: - ctx  - user context, must be a `TS`

4421:   Output Parameters:
4422: + A - the Jacobian
4423: - B - the preconditioning matrix (may be the same as A)

4425:   Level: developer

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

4430: .seealso: [](ch_ts), `SNESSetJacobian()`
4431: @*/
4432: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4433: {
4434:   TS ts = (TS)ctx;

4436:   PetscFunctionBegin;
4442:   PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4443:   PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4444:   PetscFunctionReturn(PETSC_SUCCESS);
4445: }

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

4450:   Collective

4452:   Input Parameters:
4453: + ts  - time stepping context
4454: . t   - time at which to evaluate
4455: . U   - state at which to evaluate
4456: - ctx - context

4458:   Output Parameter:
4459: . F - right hand side

4461:   Level: intermediate

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

4467: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4468: @*/
4469: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4470: {
4471:   Mat Arhs, Brhs;

4473:   PetscFunctionBegin;
4474:   PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4475:   /* undo the damage caused by shifting */
4476:   PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4477:   PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4478:   PetscCall(MatMult(Arhs, U, F));
4479:   PetscFunctionReturn(PETSC_SUCCESS);
4480: }

4482: /*@C
4483:   TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.

4485:   Collective

4487:   Input Parameters:
4488: + ts  - time stepping context
4489: . t   - time at which to evaluate
4490: . U   - state at which to evaluate
4491: - ctx - context

4493:   Output Parameters:
4494: + A - pointer to operator
4495: - B - pointer to preconditioning matrix

4497:   Level: intermediate

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

4502: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4503: @*/
4504: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4505: {
4506:   PetscFunctionBegin;
4507:   PetscFunctionReturn(PETSC_SUCCESS);
4508: }

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

4513:   Collective

4515:   Input Parameters:
4516: + ts   - time stepping context
4517: . t    - time at which to evaluate
4518: . U    - state at which to evaluate
4519: . Udot - time derivative of state vector
4520: - ctx  - context

4522:   Output Parameter:
4523: . F - left hand side

4525:   Level: intermediate

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

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

4535: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4536: @*/
4537: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4538: {
4539:   Mat A, B;

4541:   PetscFunctionBegin;
4542:   PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4543:   PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4544:   PetscCall(MatMult(A, Udot, F));
4545:   PetscFunctionReturn(PETSC_SUCCESS);
4546: }

4548: /*@C
4549:   TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE

4551:   Collective

4553:   Input Parameters:
4554: + ts    - time stepping context
4555: . t     - time at which to evaluate
4556: . U     - state at which to evaluate
4557: . Udot  - time derivative of state vector
4558: . shift - shift to apply
4559: - ctx   - context

4561:   Output Parameters:
4562: + A - pointer to operator
4563: - B - pointer to matrix from which the preconditioner is built (often `A`)

4565:   Level: advanced

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

4570:   It is only appropriate for problems of the form

4572:   $$
4573:   M \dot{U} = F(U,t)
4574:   $$

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

4580:   $$
4581:   shift*M + J
4582:   $$

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

4587: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4588: @*/
4589: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4590: {
4591:   PetscFunctionBegin;
4592:   PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4593:   ts->ijacobian.shift = shift;
4594:   PetscFunctionReturn(PETSC_SUCCESS);
4595: }

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

4600:   Not Collective

4602:   Input Parameter:
4603: . ts - the `TS` context

4605:   Output Parameter:
4606: . equation_type - see `TSEquationType`

4608:   Level: beginner

4610: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4611: @*/
4612: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4613: {
4614:   PetscFunctionBegin;
4616:   PetscAssertPointer(equation_type, 2);
4617:   *equation_type = ts->equation_type;
4618:   PetscFunctionReturn(PETSC_SUCCESS);
4619: }

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

4624:   Not Collective

4626:   Input Parameters:
4627: + ts            - the `TS` context
4628: - equation_type - see `TSEquationType`

4630:   Level: advanced

4632: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4633: @*/
4634: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4635: {
4636:   PetscFunctionBegin;
4638:   ts->equation_type = equation_type;
4639:   PetscFunctionReturn(PETSC_SUCCESS);
4640: }

4642: /*@
4643:   TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.

4645:   Not Collective

4647:   Input Parameter:
4648: . ts - the `TS` context

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

4654:   Level: beginner

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

4659: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4660: @*/
4661: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4662: {
4663:   PetscFunctionBegin;
4665:   PetscAssertPointer(reason, 2);
4666:   *reason = ts->reason;
4667:   PetscFunctionReturn(PETSC_SUCCESS);
4668: }

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

4673:   Logically Collective; reason must contain common value

4675:   Input Parameters:
4676: + ts     - the `TS` context
4677: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4678:             manual pages for the individual convergence tests for complete lists

4680:   Level: advanced

4682:   Note:
4683:   Can only be called while `TSSolve()` is active.

4685: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4686: @*/
4687: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4688: {
4689:   PetscFunctionBegin;
4691:   ts->reason = reason;
4692:   PetscFunctionReturn(PETSC_SUCCESS);
4693: }

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

4698:   Not Collective

4700:   Input Parameter:
4701: . ts - the `TS` context

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

4706:   Level: beginner

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

4711: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4712: @*/
4713: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4714: {
4715:   PetscFunctionBegin;
4717:   PetscAssertPointer(ftime, 2);
4718:   *ftime = ts->solvetime;
4719:   PetscFunctionReturn(PETSC_SUCCESS);
4720: }

4722: /*@
4723:   TSGetSNESIterations - Gets the total number of nonlinear iterations
4724:   used by the time integrator.

4726:   Not Collective

4728:   Input Parameter:
4729: . ts - `TS` context

4731:   Output Parameter:
4732: . nits - number of nonlinear iterations

4734:   Level: intermediate

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

4739: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4740: @*/
4741: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4742: {
4743:   PetscFunctionBegin;
4745:   PetscAssertPointer(nits, 2);
4746:   *nits = ts->snes_its;
4747:   PetscFunctionReturn(PETSC_SUCCESS);
4748: }

4750: /*@
4751:   TSGetKSPIterations - Gets the total number of linear iterations
4752:   used by the time integrator.

4754:   Not Collective

4756:   Input Parameter:
4757: . ts - `TS` context

4759:   Output Parameter:
4760: . lits - number of linear iterations

4762:   Level: intermediate

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

4767: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()`
4768: @*/
4769: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4770: {
4771:   PetscFunctionBegin;
4773:   PetscAssertPointer(lits, 2);
4774:   *lits = ts->ksp_its;
4775:   PetscFunctionReturn(PETSC_SUCCESS);
4776: }

4778: /*@
4779:   TSGetStepRejections - Gets the total number of rejected steps.

4781:   Not Collective

4783:   Input Parameter:
4784: . ts - `TS` context

4786:   Output Parameter:
4787: . rejects - number of steps rejected

4789:   Level: intermediate

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

4794: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4795: @*/
4796: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4797: {
4798:   PetscFunctionBegin;
4800:   PetscAssertPointer(rejects, 2);
4801:   *rejects = ts->reject;
4802:   PetscFunctionReturn(PETSC_SUCCESS);
4803: }

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

4808:   Not Collective

4810:   Input Parameter:
4811: . ts - `TS` context

4813:   Output Parameter:
4814: . fails - number of failed nonlinear solves

4816:   Level: intermediate

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

4821: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4822: @*/
4823: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4824: {
4825:   PetscFunctionBegin;
4827:   PetscAssertPointer(fails, 2);
4828:   *fails = ts->num_snes_failures;
4829:   PetscFunctionReturn(PETSC_SUCCESS);
4830: }

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

4835:   Not Collective

4837:   Input Parameters:
4838: + ts      - `TS` context
4839: - rejects - maximum number of rejected steps, pass -1 for unlimited

4841:   Options Database Key:
4842: . -ts_max_reject - Maximum number of step rejections before a step fails

4844:   Level: intermediate

4846: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4847: @*/
4848: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4849: {
4850:   PetscFunctionBegin;
4852:   ts->max_reject = rejects;
4853:   PetscFunctionReturn(PETSC_SUCCESS);
4854: }

4856: /*@
4857:   TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves

4859:   Not Collective

4861:   Input Parameters:
4862: + ts    - `TS` context
4863: - fails - maximum number of failed nonlinear solves, pass -1 for unlimited

4865:   Options Database Key:
4866: . -ts_max_snes_failures - Maximum number of nonlinear solve failures

4868:   Level: intermediate

4870: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4871: @*/
4872: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4873: {
4874:   PetscFunctionBegin;
4876:   ts->max_snes_failures = fails;
4877:   PetscFunctionReturn(PETSC_SUCCESS);
4878: }

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

4883:   Not Collective

4885:   Input Parameters:
4886: + ts  - `TS` context
4887: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure

4889:   Options Database Key:
4890: . -ts_error_if_step_fails - Error if no step succeeds

4892:   Level: intermediate

4894: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
4895: @*/
4896: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
4897: {
4898:   PetscFunctionBegin;
4900:   ts->errorifstepfailed = err;
4901:   PetscFunctionReturn(PETSC_SUCCESS);
4902: }

4904: /*@
4905:   TSGetAdapt - Get the adaptive controller context for the current method

4907:   Collective if controller has not yet been created

4909:   Input Parameter:
4910: . ts - time stepping context

4912:   Output Parameter:
4913: . adapt - adaptive controller

4915:   Level: intermediate

4917: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
4918: @*/
4919: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
4920: {
4921:   PetscFunctionBegin;
4923:   PetscAssertPointer(adapt, 2);
4924:   if (!ts->adapt) {
4925:     PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
4926:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
4927:   }
4928:   *adapt = ts->adapt;
4929:   PetscFunctionReturn(PETSC_SUCCESS);
4930: }

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

4935:   Logically Collective

4937:   Input Parameters:
4938: + ts    - time integration context
4939: . atol  - scalar absolute tolerances, `PETSC_DECIDE` to leave current value
4940: . vatol - vector of absolute tolerances or `NULL`, used in preference to atol if present
4941: . rtol  - scalar relative tolerances, `PETSC_DECIDE` to leave current value
4942: - vrtol - vector of relative tolerances or `NULL`, used in preference to atol if present

4944:   Options Database Keys:
4945: + -ts_rtol <rtol> - relative tolerance for local truncation error
4946: - -ts_atol <atol> - Absolute tolerance for local truncation error

4948:   Level: beginner

4950:   Notes:
4951:   With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
4952:   (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
4953:   computed only for the differential or the algebraic part then this can be done using the vector of
4954:   tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
4955:   differential part and infinity for the algebraic part, the LTE calculation will include only the
4956:   differential variables.

4958: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
4959: @*/
4960: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
4961: {
4962:   PetscFunctionBegin;
4963:   if (atol != (PetscReal)PETSC_DECIDE && atol != (PetscReal)PETSC_DEFAULT) ts->atol = atol;
4964:   if (vatol) {
4965:     PetscCall(PetscObjectReference((PetscObject)vatol));
4966:     PetscCall(VecDestroy(&ts->vatol));
4967:     ts->vatol = vatol;
4968:   }
4969:   if (rtol != (PetscReal)PETSC_DECIDE && rtol != (PetscReal)PETSC_DEFAULT) ts->rtol = rtol;
4970:   if (vrtol) {
4971:     PetscCall(PetscObjectReference((PetscObject)vrtol));
4972:     PetscCall(VecDestroy(&ts->vrtol));
4973:     ts->vrtol = vrtol;
4974:   }
4975:   PetscFunctionReturn(PETSC_SUCCESS);
4976: }

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

4981:   Logically Collective

4983:   Input Parameter:
4984: . ts - time integration context

4986:   Output Parameters:
4987: + atol  - scalar absolute tolerances, `NULL` to ignore
4988: . vatol - vector of absolute tolerances, `NULL` to ignore
4989: . rtol  - scalar relative tolerances, `NULL` to ignore
4990: - vrtol - vector of relative tolerances, `NULL` to ignore

4992:   Level: beginner

4994: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
4995: @*/
4996: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
4997: {
4998:   PetscFunctionBegin;
4999:   if (atol) *atol = ts->atol;
5000:   if (vatol) *vatol = ts->vatol;
5001:   if (rtol) *rtol = ts->rtol;
5002:   if (vrtol) *vrtol = ts->vrtol;
5003:   PetscFunctionReturn(PETSC_SUCCESS);
5004: }

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

5009:   Collective

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

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

5022:   Options Database Key:
5023: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5025:   Level: developer

5027: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5028: @*/
5029: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5030: {
5031:   PetscInt norma_loc, norm_loc, normr_loc;

5033:   PetscFunctionBegin;
5035:   PetscCall(VecErrorWeightedNorms(U, Y, NULL, wnormtype, ts->atol, ts->vatol, ts->rtol, ts->vrtol, ts->adapt->ignore_max, norm, &norm_loc, norma, &norma_loc, normr, &normr_loc));
5036:   if (wnormtype == NORM_2) {
5037:     if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5038:     if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5039:     if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5040:   }
5041:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5042:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5043:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5044:   PetscFunctionReturn(PETSC_SUCCESS);
5045: }

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

5050:   Collective

5052:   Input Parameters:
5053: + ts        - time stepping context
5054: . E         - error vector
5055: . U         - state vector, usually ts->vec_sol
5056: . Y         - state vector, previous time step
5057: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

5059:   Output Parameters:
5060: + norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5061: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5062: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

5064:   Options Database Key:
5065: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5067:   Level: developer

5069: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5070: @*/
5071: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5072: {
5073:   PetscInt norma_loc, norm_loc, normr_loc;

5075:   PetscFunctionBegin;
5077:   PetscCall(VecErrorWeightedNorms(U, Y, E, wnormtype, ts->atol, ts->vatol, ts->rtol, ts->vrtol, ts->adapt->ignore_max, norm, &norm_loc, norma, &norma_loc, normr, &normr_loc));
5078:   if (wnormtype == NORM_2) {
5079:     if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5080:     if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5081:     if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5082:   }
5083:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5084:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5085:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5086:   PetscFunctionReturn(PETSC_SUCCESS);
5087: }

5089: /*@
5090:   TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler

5092:   Logically Collective

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

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

5101:   Level: intermediate

5103: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5104: @*/
5105: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5106: {
5107:   PetscFunctionBegin;
5109:   ts->cfltime_local = cfltime;
5110:   ts->cfltime       = -1.;
5111:   PetscFunctionReturn(PETSC_SUCCESS);
5112: }

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

5117:   Collective

5119:   Input Parameter:
5120: . ts - time stepping context

5122:   Output Parameter:
5123: . cfltime - maximum stable time step for forward Euler

5125:   Level: advanced

5127: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5128: @*/
5129: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5130: {
5131:   PetscFunctionBegin;
5132:   if (ts->cfltime < 0) PetscCall(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5133:   *cfltime = ts->cfltime;
5134:   PetscFunctionReturn(PETSC_SUCCESS);
5135: }

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

5140:   Input Parameters:
5141: + ts - the `TS` context.
5142: . xl - lower bound.
5143: - xu - upper bound.

5145:   Level: advanced

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

5151: .seealso: [](ch_ts), `TS`
5152: @*/
5153: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5154: {
5155:   SNES snes;

5157:   PetscFunctionBegin;
5158:   PetscCall(TSGetSNES(ts, &snes));
5159:   PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5160:   PetscFunctionReturn(PETSC_SUCCESS);
5161: }

5163: /*@
5164:   TSComputeLinearStability - computes the linear stability function at a point

5166:   Collective

5168:   Input Parameters:
5169: + ts - the `TS` context
5170: . xr - real part of input argument
5171: - xi - imaginary part of input argument

5173:   Output Parameters:
5174: + yr - real part of function value
5175: - yi - imaginary part of function value

5177:   Level: developer

5179: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5180: @*/
5181: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5182: {
5183:   PetscFunctionBegin;
5185:   PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5186:   PetscFunctionReturn(PETSC_SUCCESS);
5187: }

5189: /*@
5190:   TSRestartStep - Flags the solver to restart the next step

5192:   Collective

5194:   Input Parameter:
5195: . ts - the `TS` context obtained from `TSCreate()`

5197:   Level: advanced

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

5207: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5208: @*/
5209: PetscErrorCode TSRestartStep(TS ts)
5210: {
5211:   PetscFunctionBegin;
5213:   ts->steprestart = PETSC_TRUE;
5214:   PetscFunctionReturn(PETSC_SUCCESS);
5215: }

5217: /*@
5218:   TSRollBack - Rolls back one time step

5220:   Collective

5222:   Input Parameter:
5223: . ts - the `TS` context obtained from `TSCreate()`

5225:   Level: advanced

5227: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5228: @*/
5229: PetscErrorCode TSRollBack(TS ts)
5230: {
5231:   PetscFunctionBegin;
5233:   PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5234:   PetscUseTypeMethod(ts, rollback);
5235:   ts->time_step  = ts->ptime - ts->ptime_prev;
5236:   ts->ptime      = ts->ptime_prev;
5237:   ts->ptime_prev = ts->ptime_prev_rollback;
5238:   ts->steps--;
5239:   ts->steprollback = PETSC_TRUE;
5240:   PetscFunctionReturn(PETSC_SUCCESS);
5241: }

5243: /*@
5244:   TSGetStages - Get the number of stages and stage values

5246:   Input Parameter:
5247: . ts - the `TS` context obtained from `TSCreate()`

5249:   Output Parameters:
5250: + ns - the number of stages
5251: - Y  - the current stage vectors

5253:   Level: advanced

5255:   Note:
5256:   Both `ns` and `Y` can be `NULL`.

5258: .seealso: [](ch_ts), `TS`, `TSCreate()`
5259: @*/
5260: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5261: {
5262:   PetscFunctionBegin;
5264:   if (ns) PetscAssertPointer(ns, 2);
5265:   if (Y) PetscAssertPointer(Y, 3);
5266:   if (!ts->ops->getstages) {
5267:     if (ns) *ns = 0;
5268:     if (Y) *Y = NULL;
5269:   } else PetscUseTypeMethod(ts, getstages, ns, Y);
5270:   PetscFunctionReturn(PETSC_SUCCESS);
5271: }

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

5276:   Collective

5278:   Input Parameters:
5279: + ts    - the `TS` context
5280: . t     - current timestep
5281: . U     - state vector
5282: . Udot  - time derivative of state vector
5283: . shift - shift to apply, see note below
5284: - ctx   - an optional user context

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

5290:   Level: intermediate

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

5295:   dF/dU + shift*dF/dUdot

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

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

5304: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5305: @*/
5306: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5307: {
5308:   SNES          snes;
5309:   MatFDColoring color;
5310:   PetscBool     hascolor, matcolor = PETSC_FALSE;

5312:   PetscFunctionBegin;
5313:   PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5314:   PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5315:   if (!color) {
5316:     DM         dm;
5317:     ISColoring iscoloring;

5319:     PetscCall(TSGetDM(ts, &dm));
5320:     PetscCall(DMHasColoring(dm, &hascolor));
5321:     if (hascolor && !matcolor) {
5322:       PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5323:       PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5324:       PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5325:       PetscCall(MatFDColoringSetFromOptions(color));
5326:       PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5327:       PetscCall(ISColoringDestroy(&iscoloring));
5328:     } else {
5329:       MatColoring mc;

5331:       PetscCall(MatColoringCreate(B, &mc));
5332:       PetscCall(MatColoringSetDistance(mc, 2));
5333:       PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5334:       PetscCall(MatColoringSetFromOptions(mc));
5335:       PetscCall(MatColoringApply(mc, &iscoloring));
5336:       PetscCall(MatColoringDestroy(&mc));
5337:       PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5338:       PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5339:       PetscCall(MatFDColoringSetFromOptions(color));
5340:       PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5341:       PetscCall(ISColoringDestroy(&iscoloring));
5342:     }
5343:     PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5344:     PetscCall(PetscObjectDereference((PetscObject)color));
5345:   }
5346:   PetscCall(TSGetSNES(ts, &snes));
5347:   PetscCall(MatFDColoringApply(B, color, U, snes));
5348:   if (J != B) {
5349:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5350:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5351:   }
5352:   PetscFunctionReturn(PETSC_SUCCESS);
5353: }

5355: /*@C
5356:   TSSetFunctionDomainError - Set a function that tests if the current state vector is valid

5358:   Input Parameters:
5359: + ts   - the `TS` context
5360: - func - function called within `TSFunctionDomainError()`

5362:   Calling sequence of `func`:
5363: + ts     - the `TS` context
5364: . time   - the current time (of the stage)
5365: . state  - the state to check if it is valid
5366: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable

5368:   Level: intermediate

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

5376:   Developer Notes:
5377:   The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5378:   since one takes a function pointer and the other does not.

5380: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5381: @*/
5382: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5383: {
5384:   PetscFunctionBegin;
5386:   ts->functiondomainerror = func;
5387:   PetscFunctionReturn(PETSC_SUCCESS);
5388: }

5390: /*@
5391:   TSFunctionDomainError - Checks if the current state is valid

5393:   Input Parameters:
5394: + ts        - the `TS` context
5395: . stagetime - time of the simulation
5396: - Y         - state vector to check.

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

5401:   Level: developer

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

5407: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5408: @*/
5409: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5410: {
5411:   PetscFunctionBegin;
5413:   *accept = PETSC_TRUE;
5414:   if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5415:   PetscFunctionReturn(PETSC_SUCCESS);
5416: }

5418: /*@C
5419:   TSClone - This function clones a time step `TS` object.

5421:   Collective

5423:   Input Parameter:
5424: . tsin - The input `TS`

5426:   Output Parameter:
5427: . tsout - The output `TS` (cloned)

5429:   Level: developer

5431:   Notes:
5432:   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.
5433:   It will likely be replaced in the future with a mechanism of switching methods on the fly.

5435:   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
5436: .vb
5437:  SNES snes_dup = NULL;
5438:  TSGetSNES(ts,&snes_dup);
5439:  TSSetSNES(ts,snes_dup);
5440: .ve

5442: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5443: @*/
5444: PetscErrorCode TSClone(TS tsin, TS *tsout)
5445: {
5446:   TS     t;
5447:   SNES   snes_start;
5448:   DM     dm;
5449:   TSType type;

5451:   PetscFunctionBegin;
5452:   PetscAssertPointer(tsin, 1);
5453:   *tsout = NULL;

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

5457:   /* General TS description */
5458:   t->numbermonitors    = 0;
5459:   t->monitorFrequency  = 1;
5460:   t->setupcalled       = 0;
5461:   t->ksp_its           = 0;
5462:   t->snes_its          = 0;
5463:   t->nwork             = 0;
5464:   t->rhsjacobian.time  = PETSC_MIN_REAL;
5465:   t->rhsjacobian.scale = 1.;
5466:   t->ijacobian.shift   = 1.;

5468:   PetscCall(TSGetSNES(tsin, &snes_start));
5469:   PetscCall(TSSetSNES(t, snes_start));

5471:   PetscCall(TSGetDM(tsin, &dm));
5472:   PetscCall(TSSetDM(t, dm));

5474:   t->adapt = tsin->adapt;
5475:   PetscCall(PetscObjectReference((PetscObject)t->adapt));

5477:   t->trajectory = tsin->trajectory;
5478:   PetscCall(PetscObjectReference((PetscObject)t->trajectory));

5480:   t->event = tsin->event;
5481:   if (t->event) t->event->refct++;

5483:   t->problem_type      = tsin->problem_type;
5484:   t->ptime             = tsin->ptime;
5485:   t->ptime_prev        = tsin->ptime_prev;
5486:   t->time_step         = tsin->time_step;
5487:   t->max_time          = tsin->max_time;
5488:   t->steps             = tsin->steps;
5489:   t->max_steps         = tsin->max_steps;
5490:   t->equation_type     = tsin->equation_type;
5491:   t->atol              = tsin->atol;
5492:   t->rtol              = tsin->rtol;
5493:   t->max_snes_failures = tsin->max_snes_failures;
5494:   t->max_reject        = tsin->max_reject;
5495:   t->errorifstepfailed = tsin->errorifstepfailed;

5497:   PetscCall(TSGetType(tsin, &type));
5498:   PetscCall(TSSetType(t, type));

5500:   t->vec_sol = NULL;

5502:   t->cfltime          = tsin->cfltime;
5503:   t->cfltime_local    = tsin->cfltime_local;
5504:   t->exact_final_time = tsin->exact_final_time;

5506:   t->ops[0] = tsin->ops[0];

5508:   if (((PetscObject)tsin)->fortran_func_pointers) {
5509:     PetscInt i;
5510:     PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers));
5511:     for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5512:   }
5513:   *tsout = t;
5514:   PetscFunctionReturn(PETSC_SUCCESS);
5515: }

5517: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5518: {
5519:   TS ts = (TS)ctx;

5521:   PetscFunctionBegin;
5522:   PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5523:   PetscFunctionReturn(PETSC_SUCCESS);
5524: }

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

5529:   Logically Collective

5531:   Input Parameter:
5532: . ts - the time stepping routine

5534:   Output Parameter:
5535: . flg - `PETSC_TRUE` if the multiply is likely correct

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

5540:   Level: advanced

5542:   Note:
5543:   This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian

5545: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5546: @*/
5547: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5548: {
5549:   Mat              J, B;
5550:   TSRHSJacobianFn *func;
5551:   void            *ctx;

5553:   PetscFunctionBegin;
5554:   PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5555:   PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5556:   PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5557:   PetscFunctionReturn(PETSC_SUCCESS);
5558: }

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

5563:   Logically Collective

5565:   Input Parameter:
5566: . ts - the time stepping routine

5568:   Output Parameter:
5569: . flg - `PETSC_TRUE` if the multiply is likely correct

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

5574:   Level: advanced

5576:   Notes:
5577:   This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian

5579: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5580: @*/
5581: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5582: {
5583:   Mat              J, B;
5584:   void            *ctx;
5585:   TSRHSJacobianFn *func;

5587:   PetscFunctionBegin;
5588:   PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5589:   PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5590:   PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5591:   PetscFunctionReturn(PETSC_SUCCESS);
5592: }

5594: /*@
5595:   TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.

5597:   Logically Collective

5599:   Input Parameters:
5600: + ts                   - timestepping context
5601: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

5603:   Options Database Key:
5604: . -ts_use_splitrhsfunction - <true,false>

5606:   Level: intermediate

5608:   Note:
5609:   This is only for multirate methods

5611: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5612: @*/
5613: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5614: {
5615:   PetscFunctionBegin;
5617:   ts->use_splitrhsfunction = use_splitrhsfunction;
5618:   PetscFunctionReturn(PETSC_SUCCESS);
5619: }

5621: /*@
5622:   TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.

5624:   Not Collective

5626:   Input Parameter:
5627: . ts - timestepping context

5629:   Output Parameter:
5630: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

5632:   Level: intermediate

5634: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5635: @*/
5636: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5637: {
5638:   PetscFunctionBegin;
5640:   *use_splitrhsfunction = ts->use_splitrhsfunction;
5641:   PetscFunctionReturn(PETSC_SUCCESS);
5642: }

5644: /*@
5645:   TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.

5647:   Logically  Collective

5649:   Input Parameters:
5650: + ts  - the time-stepper
5651: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)

5653:   Level: intermediate

5655:   Note:
5656:   When the relationship between the nonzero structures is known and supplied the solution process can be much faster

5658: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
5659:  @*/
5660: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5661: {
5662:   PetscFunctionBegin;
5664:   ts->axpy_pattern = str;
5665:   PetscFunctionReturn(PETSC_SUCCESS);
5666: }

5668: /*@
5669:   TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span

5671:   Collective

5673:   Input Parameters:
5674: + ts         - the time-stepper
5675: . n          - number of the time points (>=2)
5676: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.

5678:   Options Database Key:
5679: . -ts_time_span <t0,...tf> - Sets the time span

5681:   Level: intermediate

5683:   Notes:
5684:   The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
5685:   `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
5686:   The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may
5687:   pressure the memory system when using a large number of span points.

5689: .seealso: [](ch_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()`
5690:  @*/
5691: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
5692: {
5693:   PetscFunctionBegin;
5695:   PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
5696:   if (ts->tspan && n != ts->tspan->num_span_times) {
5697:     PetscCall(PetscFree(ts->tspan->span_times));
5698:     PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
5699:     PetscCall(PetscMalloc1(n, &ts->tspan->span_times));
5700:   }
5701:   if (!ts->tspan) {
5702:     TSTimeSpan tspan;
5703:     PetscCall(PetscNew(&tspan));
5704:     PetscCall(PetscMalloc1(n, &tspan->span_times));
5705:     tspan->reltol  = 1e-6;
5706:     tspan->abstol  = 10 * PETSC_MACHINE_EPSILON;
5707:     tspan->worktol = 0;
5708:     ts->tspan      = tspan;
5709:   }
5710:   ts->tspan->num_span_times = n;
5711:   PetscCall(PetscArraycpy(ts->tspan->span_times, span_times, n));
5712:   PetscCall(TSSetTime(ts, ts->tspan->span_times[0]));
5713:   PetscCall(TSSetMaxTime(ts, ts->tspan->span_times[n - 1]));
5714:   PetscFunctionReturn(PETSC_SUCCESS);
5715: }

5717: /*@C
5718:   TSGetTimeSpan - gets the time span set with `TSSetTimeSpan()`

5720:   Not Collective

5722:   Input Parameter:
5723: . ts - the time-stepper

5725:   Output Parameters:
5726: + n          - number of the time points (>=2)
5727: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.

5729:   Level: beginner

5731:   Note:
5732:   The values obtained are valid until the `TS` object is destroyed.

5734:   Both `n` and `span_times` can be `NULL`.

5736: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()`
5737:  @*/
5738: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal **span_times)
5739: {
5740:   PetscFunctionBegin;
5742:   if (n) PetscAssertPointer(n, 2);
5743:   if (span_times) PetscAssertPointer(span_times, 3);
5744:   if (!ts->tspan) {
5745:     if (n) *n = 0;
5746:     if (span_times) *span_times = NULL;
5747:   } else {
5748:     if (n) *n = ts->tspan->num_span_times;
5749:     if (span_times) *span_times = ts->tspan->span_times;
5750:   }
5751:   PetscFunctionReturn(PETSC_SUCCESS);
5752: }

5754: /*@
5755:   TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.

5757:   Input Parameter:
5758: . ts - the `TS` context obtained from `TSCreate()`

5760:   Output Parameters:
5761: + nsol - the number of solutions
5762: - Sols - the solution vectors

5764:   Level: intermediate

5766:   Notes:
5767:   Both `nsol` and `Sols` can be `NULL`.

5769:   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()`.
5770:   For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span.

5772: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`
5773: @*/
5774: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
5775: {
5776:   PetscFunctionBegin;
5778:   if (nsol) PetscAssertPointer(nsol, 2);
5779:   if (Sols) PetscAssertPointer(Sols, 3);
5780:   if (!ts->tspan) {
5781:     if (nsol) *nsol = 0;
5782:     if (Sols) *Sols = NULL;
5783:   } else {
5784:     if (nsol) *nsol = ts->tspan->spanctr;
5785:     if (Sols) *Sols = ts->tspan->vecs_sol;
5786:   }
5787:   PetscFunctionReturn(PETSC_SUCCESS);
5788: }

5790: /*@C
5791:   TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.

5793:   Collective

5795:   Input Parameters:
5796: + ts - the `TS` context
5797: . J  - Jacobian matrix (not altered in this routine)
5798: - B  - newly computed Jacobian matrix to use with preconditioner

5800:   Level: intermediate

5802:   Notes:
5803:   This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
5804:   many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
5805:   and multiple fields are involved.

5807:   Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
5808:   structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
5809:   usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
5810:   `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.

5812: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5813: @*/
5814: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
5815: {
5816:   MatColoring   mc            = NULL;
5817:   ISColoring    iscoloring    = NULL;
5818:   MatFDColoring matfdcoloring = NULL;

5820:   PetscFunctionBegin;
5821:   /* Generate new coloring after eliminating zeros in the matrix */
5822:   PetscCall(MatEliminateZeros(B, PETSC_TRUE));
5823:   PetscCall(MatColoringCreate(B, &mc));
5824:   PetscCall(MatColoringSetDistance(mc, 2));
5825:   PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5826:   PetscCall(MatColoringSetFromOptions(mc));
5827:   PetscCall(MatColoringApply(mc, &iscoloring));
5828:   PetscCall(MatColoringDestroy(&mc));
5829:   /* Replace the old coloring with the new one */
5830:   PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
5831:   PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5832:   PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
5833:   PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
5834:   PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
5835:   PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
5836:   PetscCall(ISColoringDestroy(&iscoloring));
5837:   PetscFunctionReturn(PETSC_SUCCESS);
5838: }