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, init_time and max_time are set
 38: . -ts_eval_times <t0,...tn>                                          - time points where solutions are computed and stored for each indicated time
 39: . -ts_max_steps <steps>                                              - maximum time-step number to execute until (possibly with nonzero starting value)
 40: . -ts_run_steps <steps>                                              - maximum number of time steps for TSSolve to take on each call
 41: . -ts_init_time <time>                                               - initial time to start computation
 42: . -ts_final_time <time>                                              - final time to compute to (deprecated: use `-ts_max_time`)
 43: . -ts_dt <dt>                                                        - initial time step
 44: . -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
 45: . -ts_max_snes_failures <maxfailures>                                - Maximum number of nonlinear solve failures allowed
 46: . -ts_max_reject <maxrejects>                                        - Maximum number of step rejections before step fails
 47: . -ts_error_if_step_fails <true,false>                               - Error if no step succeeds
 48: . -ts_rtol <rtol>                                                    - relative tolerance for local truncation error
 49: . -ts_atol <atol>                                                    - Absolute tolerance for local truncation error
 50: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view               - test the Jacobian at each iteration against finite difference with RHS function
 51: . -ts_rhs_jacobian_test_mult_transpose                               - test the Jacobian at each iteration against finite difference with RHS function
 52: . -ts_adjoint_solve <yes,no>                                         - After solving the ODE/DAE solve the adjoint problem (requires `-ts_save_trajectory`)
 53: . -ts_fd_color                                                       - Use finite differences with coloring to compute IJacobian
 54: . -ts_monitor                                                        - print information at each timestep
 55: . -ts_monitor_cancel                                                 - Cancel all monitors
 56: . -ts_monitor_wall_clock_time                                        - Monitor wall-clock time, KSP iterations, and SNES iterations per step
 57: . -ts_monitor_lg_solution                                            - Monitor solution graphically
 58: . -ts_monitor_lg_error                                               - Monitor error graphically
 59: . -ts_monitor_error                                                  - Monitors norm of error
 60: . -ts_monitor_lg_timestep                                            - Monitor timestep size graphically
 61: . -ts_monitor_lg_timestep_log                                        - Monitor log timestep size graphically
 62: . -ts_monitor_lg_snes_iterations                                     - Monitor number nonlinear iterations for each timestep graphically
 63: . -ts_monitor_lg_ksp_iterations                                      - Monitor number nonlinear iterations for each timestep graphically
 64: . -ts_monitor_sp_eig                                                 - Monitor eigenvalues of linearized operator graphically
 65: . -ts_monitor_draw_solution                                          - Monitor solution graphically
 66: . -ts_monitor_draw_solution_phase  <xleft,yleft,xright,yright>       - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
 67: . -ts_monitor_draw_error                                             - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
 68: . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
 69: . -ts_monitor_solution_interval <interval>                           - output once every interval (default=1) time steps. Use -1 to only output at the end of the simulation
 70: . -ts_monitor_solution_skip_initial                                  - skip writing of initial condition
 71: . -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)
 72: . -ts_monitor_solution_vtk_interval <interval>                       - output once every interval (default=1) time steps. Use -1 to only output at the end of the simulation
 73: - -ts_monitor_envelope                                               - determine maximum and minimum value of each component of the solution over the solution time

 75:   Level: beginner

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

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

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

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

101:   PetscFunctionBegin;

104:   PetscCall(TSRegisterAll());
105:   PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));

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

114:   /* Handle generic TS options */
115:   PetscCall(PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL));
116:   PetscCall(PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL));
117:   PetscCall(PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", eval_times, &num_eval_times, &flg));
118:   if (flg) PetscCall(TSSetTimeSpan(ts, num_eval_times, eval_times));
119:   num_eval_times = PETSC_STATIC_ARRAY_LENGTH(eval_times);
120:   PetscCall(PetscOptionsRealArray("-ts_eval_times", "Evaluation time points", "TSSetEvaluationTimes", eval_times, &num_eval_times, &opt));
121:   PetscCheck(flg != opt || (!flg && !opt), PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "May not provide -ts_time_span and -ts_eval_times simultaneously");
122:   if (opt) PetscCall(TSSetEvaluationTimes(ts, num_eval_times, eval_times));
123:   PetscCall(PetscOptionsInt("-ts_max_steps", "Maximum time step number to execute to (possibly with non-zero starting value)", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL));
124:   PetscCall(PetscOptionsInt("-ts_run_steps", "Maximum number of time steps to take on each call to TSSolve()", "TSSetRunSteps", ts->run_steps, &ts->run_steps, NULL));
125:   PetscCall(PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL));
126:   PetscCall(PetscOptionsReal("-ts_dt", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg));
127:   if (flg) PetscCall(TSSetTimeStep(ts, time_step));
128:   PetscCall(PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg));
129:   if (flg) PetscCall(TSSetExactFinalTime(ts, eftopt));
130:   PetscCall(PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, &flg));
131:   if (flg) PetscCall(TSSetMaxSNESFailures(ts, ts->max_snes_failures));
132:   PetscCall(PetscOptionsInt("-ts_max_reject", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, &flg));
133:   if (flg) PetscCall(TSSetMaxStepRejections(ts, ts->max_reject));
134:   PetscCall(PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL));
135:   PetscCall(PetscOptionsBoundedReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL, 0));
136:   PetscCall(PetscOptionsBoundedReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL, 0));

138:   PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL));
139:   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));
140:   PetscCall(PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL));
141: #if defined(PETSC_HAVE_SAWS)
142:   {
143:     PetscBool set;
144:     flg = PETSC_FALSE;
145:     PetscCall(PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set));
146:     if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)ts, flg));
147:   }
148: #endif

150:   /* Monitor options */
151:   PetscCall(PetscOptionsDeprecated("-ts_monitor_frequency", "-ts_dmswarm_monitor_moments_interval", "3.24", "Retired in favor of monitor-specific intervals (ts_dmswarm_monitor_moments was the only monitor to use ts_monitor_frequency)"));
152:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL));
153:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_wall_clock_time", "Monitor wall-clock time, KSP iterations, and SNES iterations per step", "TSMonitorWallClockTime", TSMonitorWallClockTime, TSMonitorWallClockTimeSetUp));
154:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL));
155:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, TSMonitorSolutionSetup));
156:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL));
157:   PetscCall(PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg));
158:   if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)ts, monfilename));

160:   PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt));
161:   if (opt) {
162:     PetscInt  howoften = 1;
163:     DM        dm;
164:     PetscBool net;

166:     PetscCall(PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL));
167:     PetscCall(TSGetDM(ts, &dm));
168:     PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net));
169:     if (net) {
170:       TSMonitorLGCtxNetwork ctx;
171:       PetscCall(TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx));
172:       PetscCall(TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxNetworkDestroy));
173:       PetscCall(PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL));
174:     } else {
175:       TSMonitorLGCtx ctx;
176:       PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
177:       PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
178:     }
179:   }

181:   PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt));
182:   if (opt) {
183:     TSMonitorLGCtx ctx;
184:     PetscInt       howoften = 1;

186:     PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL));
187:     PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
188:     PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
189:   }
190:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL));

192:   PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt));
193:   if (opt) {
194:     TSMonitorLGCtx ctx;
195:     PetscInt       howoften = 1;

197:     PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
198:     PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
199:     PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
200:   }
201:   PetscCall(PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt));
202:   if (opt) {
203:     TSMonitorLGCtx ctx;
204:     PetscInt       howoften = 1;

206:     PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
207:     PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
208:     PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
209:     ctx->semilogy = PETSC_TRUE;
210:   }

212:   PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt));
213:   if (opt) {
214:     TSMonitorLGCtx ctx;
215:     PetscInt       howoften = 1;

217:     PetscCall(PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL));
218:     PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
219:     PetscCall(TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
220:   }
221:   PetscCall(PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt));
222:   if (opt) {
223:     TSMonitorLGCtx ctx;
224:     PetscInt       howoften = 1;

226:     PetscCall(PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL));
227:     PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
228:     PetscCall(TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
229:   }
230:   PetscCall(PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt));
231:   if (opt) {
232:     TSMonitorSPEigCtx ctx;
233:     PetscInt          howoften = 1;

235:     PetscCall(PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL));
236:     PetscCall(TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
237:     PetscCall(TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscCtxDestroyFn *)TSMonitorSPEigCtxDestroy));
238:   }
239:   PetscCall(PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase space from the DMSwarm", "TSMonitorSPSwarm", &opt));
240:   if (opt) {
241:     TSMonitorSPCtx ctx;
242:     PetscInt       howoften = 1, retain = 0;
243:     PetscBool      phase = PETSC_TRUE, create = PETSC_TRUE, multispecies = PETSC_FALSE;

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

265:     for (PetscInt i = 0; i < ts->numbermonitors; ++i)
266:       if (ts->monitor[i] == TSMonitorHGSwarmSolution) {
267:         create = PETSC_FALSE;
268:         break;
269:       }
270:     if (create) {
271:       DM       sw, dm;
272:       PetscInt Nc, Nb;

274:       PetscCall(TSGetDM(ts, &sw));
275:       PetscCall(DMSwarmGetCellDM(sw, &dm));
276:       PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc));
277:       Nb = PetscMin(20, PetscMax(10, Nc));
278:       PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL));
279:       PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL));
280:       PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL));
281:       PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL));
282:       PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx));
283:       PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscCtxDestroyFn *)TSMonitorHGCtxDestroy));
284:     }
285:   }
286:   opt = PETSC_FALSE;
287:   PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt));
288:   if (opt) {
289:     TSMonitorDrawCtx ctx;
290:     PetscInt         howoften = 1;

292:     PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL));
293:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
294:     PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
295:   }
296:   opt = PETSC_FALSE;
297:   PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt));
298:   if (opt) {
299:     TSMonitorDrawCtx ctx;
300:     PetscReal        bounds[4];
301:     PetscInt         n = 4;
302:     PetscDraw        draw;
303:     PetscDrawAxis    axis;

305:     PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL));
306:     PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field");
307:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx));
308:     PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw));
309:     PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis));
310:     PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]));
311:     PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2"));
312:     PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
313:   }
314:   opt = PETSC_FALSE;
315:   PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt));
316:   if (opt) {
317:     TSMonitorDrawCtx ctx;
318:     PetscInt         howoften = 1;

320:     PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL));
321:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
322:     PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
323:   }
324:   opt = PETSC_FALSE;
325:   PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt));
326:   if (opt) {
327:     TSMonitorDrawCtx ctx;
328:     PetscInt         howoften = 1;

330:     PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL));
331:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
332:     PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
333:   }

335:   opt = PETSC_FALSE;
336:   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));
337:   if (flg) {
338:     TSMonitorVTKCtx ctx;

340:     PetscCall(TSMonitorSolutionVTKCtxCreate(monfilename, &ctx));
341:     PetscCall(PetscOptionsInt("-ts_monitor_solution_vtk_interval", "Save every interval time step (-1 for last step only)", NULL, ctx->interval, &ctx->interval, NULL));
342:     PetscCall(TSMonitorSet(ts, (PetscErrorCode (*)(TS, PetscInt, PetscReal, Vec, void *))TSMonitorSolutionVTK, ctx, (PetscCtxDestroyFn *)TSMonitorSolutionVTKDestroy));
343:   }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

446:   Collective

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

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

454:   Level: advanced

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

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

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

472:   Collective

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

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

481:   Level: intermediate

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

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

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

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

502:   Collective

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

507:   Level: intermediate

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

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

525:   Collective

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

530:   Level: intermediate

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

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

546:   Collective

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

553:   Output Parameters:
554: + A - Jacobian matrix
555: - B - optional matrix used to compute the preconditioner, often the same as `A`

557:   Level: developer

559:   Note:
560:   Most users should not need to explicitly call this routine, as it
561:   is used internally within the ODE integrators.

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

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

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

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

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

609:   Collective

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

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

619:   Level: developer

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

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

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

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

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

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

658:   Collective

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

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

667:   Level: developer

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

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

688:   Collective

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

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

697:   Level: developer

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

707:   PetscFunctionBegin;
710:   PetscCall(TSGetDM(ts, &dm));
711:   PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));

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

717: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
718: {
719:   Mat            A, B;
720:   TSIJacobianFn *ijacobian;

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

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

769:   Collective

771:   Input Parameters:
772: + ts   - the `TS` context
773: . t    - current time
774: . U    - state vector
775: . Udot - time derivative of state vector
776: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSFunction should be kept separate

778:   Output Parameter:
779: . Y - right-hand side

781:   Level: developer

783:   Note:
784:   Most users should not need to explicitly call this routine, as it
785:   is used internally within the nonlinear solvers.

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

790: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
791: @*/
792: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
793: {
794:   TSIFunctionFn   *ifunction;
795:   TSRHSFunctionFn *rhsfunction;
796:   void            *ctx;
797:   DM               dm;

799:   PetscFunctionBegin;

805:   PetscCall(TSGetDM(ts, &dm));
806:   PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
807:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));

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

811:   PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, Udot, Y));
812:   if (ifunction) {
813:     PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
814:     ts->ifuncs++;
815:   }
816:   if (imex) {
817:     if (!ifunction) PetscCall(VecCopy(Udot, Y));
818:   } else if (rhsfunction) {
819:     if (ifunction) {
820:       Vec Frhs;

822:       PetscCall(DMGetGlobalVector(dm, &Frhs));
823:       PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
824:       PetscCall(VecAXPY(Y, -1, Frhs));
825:       PetscCall(DMRestoreGlobalVector(dm, &Frhs));
826:     } else {
827:       PetscCall(TSComputeRHSFunction(ts, t, U, Y));
828:       PetscCall(VecAYPX(Y, -1, Udot));
829:     }
830:   }
831:   PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, Udot, Y));
832:   PetscFunctionReturn(PETSC_SUCCESS);
833: }

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

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

841: */
842: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
843: {
844:   PetscFunctionBegin;
846:   PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
847:   PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");

849:   if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
850:   if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
851:   if (B && B == ts->Brhs && A != B) {
852:     if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
853:     if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
854:   }
855:   ts->rhsjacobian.shift = 0;
856:   ts->rhsjacobian.scale = 1.;
857:   PetscFunctionReturn(PETSC_SUCCESS);
858: }

860: /*@
861:   TSComputeIJacobian - Evaluates the Jacobian of the DAE

863:   Collective

865:   Input Parameters:
866: + ts    - the `TS` context
867: . t     - current timestep
868: . U     - state vector
869: . Udot  - time derivative of state vector
870: . shift - shift to apply, see note below
871: - imex  - flag indicates if the method is `TSARKIMEX` so that the RHSJacobian should be kept separate

873:   Output Parameters:
874: + A - Jacobian matrix
875: - B - matrix from which the preconditioner is constructed; often the same as `A`

877:   Level: developer

879:   Notes:
880:   If F(t,U,Udot)=0 is the DAE, the required Jacobian is
881: .vb
882:    dF/dU + shift*dF/dUdot
883: .ve
884:   Most users should not need to explicitly call this routine, as it
885:   is used internally within the nonlinear solvers.

887: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
888: @*/
889: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
890: {
891:   TSIJacobianFn   *ijacobian;
892:   TSRHSJacobianFn *rhsjacobian;
893:   DM               dm;
894:   void            *ctx;

896:   PetscFunctionBegin;

903:   PetscCall(TSGetDM(ts, &dm));
904:   PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
905:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));

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

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

945:     /* RHSJacobian needs to be converted to part of IJacobian if exists */
946:     if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
947:     if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
948:       PetscObjectState Ustate;
949:       PetscObjectId    Uid;
950:       TSRHSFunctionFn *rhsfunction;

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

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

1000: /*@C
1001:   TSSetRHSFunction - Sets the routine for evaluating the function,
1002:   where U_t = G(t,u).

1004:   Logically Collective

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

1012:   Level: beginner

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

1017: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1018: @*/
1019: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, void *ctx)
1020: {
1021:   SNES snes;
1022:   Vec  ralloc = NULL;
1023:   DM   dm;

1025:   PetscFunctionBegin;

1029:   PetscCall(TSGetDM(ts, &dm));
1030:   PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1031:   PetscCall(TSGetSNES(ts, &snes));
1032:   if (!r && !ts->dm && ts->vec_sol) {
1033:     PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1034:     r = ralloc;
1035:   }
1036:   PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1037:   PetscCall(VecDestroy(&ralloc));
1038:   PetscFunctionReturn(PETSC_SUCCESS);
1039: }

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

1044:   Logically Collective

1046:   Input Parameters:
1047: + ts  - the `TS` context obtained from `TSCreate()`
1048: . f   - routine for evaluating the solution
1049: - ctx - [optional] user-defined context for private data for the
1050:           function evaluation routine (may be `NULL`)

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

1056:   Level: intermediate

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

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

1065: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1066: @*/
1067: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, void *ctx)
1068: {
1069:   DM dm;

1071:   PetscFunctionBegin;
1073:   PetscCall(TSGetDM(ts, &dm));
1074:   PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1075:   PetscFunctionReturn(PETSC_SUCCESS);
1076: }

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

1081:   Logically Collective

1083:   Input Parameters:
1084: + ts   - the `TS` context obtained from `TSCreate()`
1085: . func - routine for evaluating the forcing function
1086: - ctx  - [optional] user-defined context for private data for the function evaluation routine
1087:          (may be `NULL`)

1089:   Level: intermediate

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

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

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

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

1103: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
1104: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1105: @*/
1106: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, void *ctx)
1107: {
1108:   DM dm;

1110:   PetscFunctionBegin;
1112:   PetscCall(TSGetDM(ts, &dm));
1113:   PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1114:   PetscFunctionReturn(PETSC_SUCCESS);
1115: }

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

1121:   Logically Collective

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

1130:   Level: beginner

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

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

1138: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
1139: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
1140: @*/
1141: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, void *ctx)
1142: {
1143:   SNES           snes;
1144:   DM             dm;
1145:   TSIJacobianFn *ijacobian;

1147:   PetscFunctionBegin;
1151:   if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1152:   if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);

1154:   PetscCall(TSGetDM(ts, &dm));
1155:   PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1156:   PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1157:   PetscCall(TSGetSNES(ts, &snes));
1158:   if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1159:   if (Amat) {
1160:     PetscCall(PetscObjectReference((PetscObject)Amat));
1161:     PetscCall(MatDestroy(&ts->Arhs));
1162:     ts->Arhs = Amat;
1163:   }
1164:   if (Pmat) {
1165:     PetscCall(PetscObjectReference((PetscObject)Pmat));
1166:     PetscCall(MatDestroy(&ts->Brhs));
1167:     ts->Brhs = Pmat;
1168:   }
1169:   PetscFunctionReturn(PETSC_SUCCESS);
1170: }

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

1175:   Logically Collective

1177:   Input Parameters:
1178: + ts  - the `TS` context obtained from `TSCreate()`
1179: . r   - vector to hold the residual (or `NULL` to have it created internally)
1180: . f   - the function evaluation routine
1181: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)

1183:   Level: beginner

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

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

1197:   PetscFunctionBegin;

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

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

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

1217:   Not Collective

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

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

1227:   Level: advanced

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

1236:   PetscFunctionBegin;
1238:   PetscCall(TSGetSNES(ts, &snes));
1239:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1240:   PetscCall(TSGetDM(ts, &dm));
1241:   PetscCall(DMTSGetIFunction(dm, func, ctx));
1242:   PetscFunctionReturn(PETSC_SUCCESS);
1243: }

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

1248:   Not Collective

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

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

1258:   Level: advanced

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

1267:   PetscFunctionBegin;
1269:   PetscCall(TSGetSNES(ts, &snes));
1270:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1271:   PetscCall(TSGetDM(ts, &dm));
1272:   PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1273:   PetscFunctionReturn(PETSC_SUCCESS);
1274: }

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

1280:   Logically Collective

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

1289:   Level: beginner

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

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

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

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

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

1309:   In case `TSSetRHSJacobian()` is also used in conjunction with a fully-implicit solver,
1310:   multilevel linear solvers, e.g. `PCMG`, will likely not work due to the way `TS` handles rhs matrices.

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

1320:   PetscFunctionBegin;
1324:   if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1325:   if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);

1327:   PetscCall(TSGetDM(ts, &dm));
1328:   PetscCall(DMTSSetIJacobian(dm, f, ctx));

1330:   PetscCall(TSGetSNES(ts, &snes));
1331:   PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1332:   PetscFunctionReturn(PETSC_SUCCESS);
1333: }

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

1338:   Logically Collective

1340:   Input Parameters:
1341: + ts    - `TS` context obtained from `TSCreate()`
1342: - reuse - `PETSC_TRUE` if the RHS Jacobian

1344:   Level: intermediate

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

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

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

1364:   Logically Collective

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

1372:   Level: beginner

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

1381:   PetscFunctionBegin;
1384:   PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1385:   PetscCall(TSGetDM(ts, &dm));
1386:   PetscCall(DMTSSetI2Function(dm, fun, ctx));
1387:   PetscFunctionReturn(PETSC_SUCCESS);
1388: }

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

1393:   Not Collective

1395:   Input Parameter:
1396: . ts - the `TS` context

1398:   Output Parameters:
1399: + r   - vector to hold residual (or `NULL`)
1400: . fun - the function to compute residual (or `NULL`)
1401: - ctx - the function context (or `NULL`)

1403:   Level: advanced

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

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

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

1425:   Logically Collective

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

1434:   Level: beginner

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

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

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

1450:   PetscFunctionBegin;
1454:   PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1455:   PetscCall(TSGetDM(ts, &dm));
1456:   PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1457:   PetscFunctionReturn(PETSC_SUCCESS);
1458: }

1460: /*@C
1461:   TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.

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

1465:   Input Parameter:
1466: . ts - The `TS` context obtained from `TSCreate()`

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

1474:   Level: advanced

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

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

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

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

1498:   Collective

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

1507:   Output Parameter:
1508: . F - the residual vector

1510:   Level: developer

1512:   Note:
1513:   Most users should not need to explicitly call this routine, as it
1514:   is used internally within the nonlinear solvers.

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

1525:   PetscFunctionBegin;

1532:   PetscCall(TSGetDM(ts, &dm));
1533:   PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1534:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));

1536:   if (!I2Function) {
1537:     PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1538:     PetscFunctionReturn(PETSC_SUCCESS);
1539:   }

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

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

1545:   if (rhsfunction) {
1546:     Vec Frhs;

1548:     PetscCall(DMGetGlobalVector(dm, &Frhs));
1549:     PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1550:     PetscCall(VecAXPY(F, -1, Frhs));
1551:     PetscCall(DMRestoreGlobalVector(dm, &Frhs));
1552:   }

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

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

1561:   Collective

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

1572:   Output Parameters:
1573: + J - Jacobian matrix
1574: - P - optional matrix used to construct the preconditioner

1576:   Level: developer

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

1581: $$
1582:   dF/dU + shiftV*dF/dV + shiftA*dF/dA
1583: $$

1585:   Most users should not need to explicitly call this routine, as it
1586:   is used internally within the ODE integrators.

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

1597:   PetscFunctionBegin;

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

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

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

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

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

1631:   Logically Collective

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

1638:   Level: advanced

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

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

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

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

1666:   Logically Collective

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

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

1675:   Level: developer

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

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

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

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

1704:   Logically Collective

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

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

1712:   Level: developer

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

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

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

1733:   Logically Collective

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

1740:   Level: beginner

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

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

1761:   Not Collective

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

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

1770:   Level: intermediate

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

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

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

1793:   Collective

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

1800:   Level: intermediate

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

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

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

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

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

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

1844:   Collective

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

1851:   Level: intermediate

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

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

1866:   Collective

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

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

1875:   Level: beginner

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

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

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

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

1901:   PetscFunctionBegin;
1903:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1905:   PetscCheckSameComm(ts, 1, viewer, 2);

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

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

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

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

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

2009:   PetscCall(PetscViewerASCIIPushTab(viewer));
2010:   PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &issundials));
2011:   PetscCall(PetscViewerASCIIPopTab(viewer));
2012:   PetscFunctionReturn(PETSC_SUCCESS);
2013: }

2015: /*@
2016:   TSSetApplicationContext - Sets an optional user-defined context for the timesteppers that may be accessed, for example inside the user provided
2017:   `TS` callbacks with `TSGetApplicationContext()`

2019:   Logically Collective

2021:   Input Parameters:
2022: + ts  - the `TS` context obtained from `TSCreate()`
2023: - ctx - user context

2025:   Level: intermediate

2027:   Fortran Note:
2028:   This only works when `ctx` is a Fortran derived type (it cannot be a `PetscObject`), we recommend writing a Fortran interface definition for this
2029:   function that tells the Fortran compiler the derived data type that is passed in as the `ctx` argument. See `TSGetApplicationContext()` for
2030:   an example.

2032: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2033: @*/
2034: PetscErrorCode TSSetApplicationContext(TS ts, PeCtx ctx)
2035: {
2036:   PetscFunctionBegin;
2038:   ts->ctx = ctx;
2039:   PetscFunctionReturn(PETSC_SUCCESS);
2040: }

2042: /*@
2043:   TSGetApplicationContext - Gets the user-defined context for the
2044:   timestepper that was set with `TSSetApplicationContext()`

2046:   Not Collective

2048:   Input Parameter:
2049: . ts - the `TS` context obtained from `TSCreate()`

2051:   Output Parameter:
2052: . ctx - a pointer to the user context

2054:   Level: intermediate

2056:   Fortran Notes:
2057:   This only works when the context is a Fortran derived type (it cannot be a `PetscObject`) and you **must** write a Fortran interface definition for this
2058:   function that tells the Fortran compiler the derived data type that is returned as the `ctx` argument. For example,
2059: .vb
2060:   Interface TSGetApplicationContext
2061:     Subroutine TSGetApplicationContext(ts,ctx,ierr)
2062:   #include <petsc/finclude/petscts.h>
2063:       use petscts
2064:       TS ts
2065:       type(tUsertype), pointer :: ctx
2066:       PetscErrorCode ierr
2067:     End Subroutine
2068:   End Interface TSGetApplicationContext
2069: .ve

2071:   The prototype for `ctx` must be
2072: .vb
2073:   type(tUsertype), pointer :: ctx
2074: .ve

2076: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2077: @*/
2078: PetscErrorCode TSGetApplicationContext(TS ts, void *ctx)
2079: {
2080:   PetscFunctionBegin;
2082:   *(void **)ctx = ts->ctx;
2083:   PetscFunctionReturn(PETSC_SUCCESS);
2084: }

2086: /*@
2087:   TSGetStepNumber - Gets the number of time steps completed.

2089:   Not Collective

2091:   Input Parameter:
2092: . ts - the `TS` context obtained from `TSCreate()`

2094:   Output Parameter:
2095: . steps - number of steps completed so far

2097:   Level: intermediate

2099: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2100: @*/
2101: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2102: {
2103:   PetscFunctionBegin;
2105:   PetscAssertPointer(steps, 2);
2106:   *steps = ts->steps;
2107:   PetscFunctionReturn(PETSC_SUCCESS);
2108: }

2110: /*@
2111:   TSSetStepNumber - Sets the number of steps completed.

2113:   Logically Collective

2115:   Input Parameters:
2116: + ts    - the `TS` context
2117: - steps - number of steps completed so far

2119:   Level: developer

2121:   Note:
2122:   For most uses of the `TS` solvers the user need not explicitly call
2123:   `TSSetStepNumber()`, as the step counter is appropriately updated in
2124:   `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2125:   reinitialize timestepping by setting the step counter to zero (and time
2126:   to the initial time) to solve a similar problem with different initial
2127:   conditions or parameters. Other possible use case is to continue
2128:   timestepping from a previously interrupted run in such a way that `TS`
2129:   monitors will be called with a initial nonzero step counter.

2131: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2132: @*/
2133: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2134: {
2135:   PetscFunctionBegin;
2138:   PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2139:   ts->steps = steps;
2140:   PetscFunctionReturn(PETSC_SUCCESS);
2141: }

2143: /*@
2144:   TSSetTimeStep - Allows one to reset the timestep at any time,
2145:   useful for simple pseudo-timestepping codes.

2147:   Logically Collective

2149:   Input Parameters:
2150: + ts        - the `TS` context obtained from `TSCreate()`
2151: - time_step - the size of the timestep

2153:   Level: intermediate

2155: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2156: @*/
2157: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2158: {
2159:   PetscFunctionBegin;
2162:   ts->time_step = time_step;
2163:   PetscFunctionReturn(PETSC_SUCCESS);
2164: }

2166: /*@
2167:   TSSetExactFinalTime - Determines whether to adapt the final time step to
2168:   match the exact final time, to interpolate the solution to the exact final time,
2169:   or to just return at the final time `TS` computed (which may be slightly larger
2170:   than the requested final time).

2172:   Logically Collective

2174:   Input Parameters:
2175: + ts     - the time-step context
2176: - eftopt - exact final time option
2177: .vb
2178:   TS_EXACTFINALTIME_STEPOVER    - Don't do anything if final time is exceeded, just use it
2179:   TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time if the final time is exceeded
2180:   TS_EXACTFINALTIME_MATCHSTEP   - Adapt final time step to ensure the computed final time exactly equals the requested final time
2181: .ve

2183:   Options Database Key:
2184: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step approach at runtime

2186:   Level: beginner

2188:   Note:
2189:   If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2190:   then the final time you selected.

2192: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2193: @*/
2194: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2195: {
2196:   PetscFunctionBegin;
2199:   ts->exact_final_time = eftopt;
2200:   PetscFunctionReturn(PETSC_SUCCESS);
2201: }

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

2206:   Not Collective

2208:   Input Parameter:
2209: . ts - the `TS` context

2211:   Output Parameter:
2212: . eftopt - exact final time option

2214:   Level: beginner

2216: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2217: @*/
2218: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2219: {
2220:   PetscFunctionBegin;
2222:   PetscAssertPointer(eftopt, 2);
2223:   *eftopt = ts->exact_final_time;
2224:   PetscFunctionReturn(PETSC_SUCCESS);
2225: }

2227: /*@
2228:   TSGetTimeStep - Gets the current timestep size.

2230:   Not Collective

2232:   Input Parameter:
2233: . ts - the `TS` context obtained from `TSCreate()`

2235:   Output Parameter:
2236: . dt - the current timestep size

2238:   Level: intermediate

2240: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2241: @*/
2242: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2243: {
2244:   PetscFunctionBegin;
2246:   PetscAssertPointer(dt, 2);
2247:   *dt = ts->time_step;
2248:   PetscFunctionReturn(PETSC_SUCCESS);
2249: }

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

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

2259:   Input Parameter:
2260: . ts - the `TS` context obtained from `TSCreate()`

2262:   Output Parameter:
2263: . v - the vector containing the solution

2265:   Level: intermediate

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

2271: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2272: @*/
2273: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2274: {
2275:   PetscFunctionBegin;
2277:   PetscAssertPointer(v, 2);
2278:   *v = ts->vec_sol;
2279:   PetscFunctionReturn(PETSC_SUCCESS);
2280: }

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

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

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

2299:   Level: advanced

2301: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2302: @*/
2303: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2304: {
2305:   PetscFunctionBegin;
2307:   if (!ts->ops->getsolutioncomponents) *n = 0;
2308:   else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2309:   PetscFunctionReturn(PETSC_SUCCESS);
2310: }

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

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

2318:   Input Parameters:
2319: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2320: - v  - the vector containing the auxiliary solution

2322:   Level: intermediate

2324: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2325: @*/
2326: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2327: {
2328:   PetscFunctionBegin;
2330:   if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2331:   else PetscCall(VecZeroEntries(*v));
2332:   PetscFunctionReturn(PETSC_SUCCESS);
2333: }

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

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

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

2346:   Level: intermediate

2348:   Note:
2349:   MUST call after `TSSetUp()`

2351: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2352: @*/
2353: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2354: {
2355:   PetscFunctionBegin;
2357:   if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2358:   else PetscCall(VecZeroEntries(*v));
2359:   PetscFunctionReturn(PETSC_SUCCESS);
2360: }

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

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

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

2373:   Level: intermediate

2375: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
2376: @*/
2377: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2378: {
2379:   PetscFunctionBegin;
2381:   PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2382:   PetscTryTypeMethod(ts, settimeerror, v);
2383:   PetscFunctionReturn(PETSC_SUCCESS);
2384: }

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

2390:   Not collective

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

2401:   Level: beginner

2403: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2404: @*/
2405: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2406: {
2407:   PetscFunctionBegin;
2409:   ts->problem_type = type;
2410:   if (type == TS_LINEAR) {
2411:     SNES snes;
2412:     PetscCall(TSGetSNES(ts, &snes));
2413:     PetscCall(SNESSetType(snes, SNESKSPONLY));
2414:   }
2415:   PetscFunctionReturn(PETSC_SUCCESS);
2416: }

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

2421:   Not collective

2423:   Input Parameter:
2424: . ts - The `TS`

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

2434:   Level: beginner

2436: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2437: @*/
2438: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2439: {
2440:   PetscFunctionBegin;
2442:   PetscAssertPointer(type, 2);
2443:   *type = ts->problem_type;
2444:   PetscFunctionReturn(PETSC_SUCCESS);
2445: }

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

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

2458:   PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2459:   if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2460:   else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2461:   PetscFunctionReturn(PETSC_SUCCESS);
2462: }

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

2467:   Collective

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

2472:   Level: advanced

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

2481: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2482: @*/
2483: PetscErrorCode TSSetUp(TS ts)
2484: {
2485:   DM dm;
2486:   PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2487:   PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2488:   TSIFunctionFn   *ifun;
2489:   TSIJacobianFn   *ijac;
2490:   TSI2JacobianFn  *i2jac;
2491:   TSRHSJacobianFn *rhsjac;

2493:   PetscFunctionBegin;
2495:   if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);

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

2502:   if (!ts->vec_sol) {
2503:     PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2504:     PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2505:   }

2507:   if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2508:     PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2509:     ts->Jacp = ts->Jacprhs;
2510:   }

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

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

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

2542:   PetscTryTypeMethod(ts, setup);

2544:   PetscCall(TSSetExactFinalTimeDefault(ts));

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

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

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

2565:   ts->setupcalled = PETSC_TRUE;
2566:   PetscFunctionReturn(PETSC_SUCCESS);
2567: }

2569: /*@
2570:   TSReset - Resets a `TS` context to the state it was in before `TSSetUp()` was called and removes any allocated `Vec` and `Mat` from its data structures

2572:   Collective

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

2577:   Level: developer

2579:   Notes:
2580:   Any options set on the `TS` object, including those set with `TSSetFromOptions()` remain.

2582:   See also `TSSetResize()` to change the size of the system being integrated (for example by adaptive mesh refinement) during the time integration.

2584: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSetResize()`
2585: @*/
2586: PetscErrorCode TSReset(TS ts)
2587: {
2588:   TS_RHSSplitLink ilink = ts->tsrhssplit, next;

2590:   PetscFunctionBegin;

2593:   PetscTryTypeMethod(ts, reset);
2594:   if (ts->snes) PetscCall(SNESReset(ts->snes));
2595:   if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));

2597:   PetscCall(MatDestroy(&ts->Arhs));
2598:   PetscCall(MatDestroy(&ts->Brhs));
2599:   PetscCall(VecDestroy(&ts->Frhs));
2600:   PetscCall(VecDestroy(&ts->vec_sol));
2601:   PetscCall(VecDestroy(&ts->vec_sol0));
2602:   PetscCall(VecDestroy(&ts->vec_dot));
2603:   PetscCall(VecDestroy(&ts->vatol));
2604:   PetscCall(VecDestroy(&ts->vrtol));
2605:   PetscCall(VecDestroyVecs(ts->nwork, &ts->work));

2607:   PetscCall(MatDestroy(&ts->Jacprhs));
2608:   PetscCall(MatDestroy(&ts->Jacp));
2609:   if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2610:   if (ts->quadraturets) {
2611:     PetscCall(TSReset(ts->quadraturets));
2612:     PetscCall(VecDestroy(&ts->vec_costintegrand));
2613:   }
2614:   while (ilink) {
2615:     next = ilink->next;
2616:     PetscCall(TSDestroy(&ilink->ts));
2617:     PetscCall(PetscFree(ilink->splitname));
2618:     PetscCall(ISDestroy(&ilink->is));
2619:     PetscCall(PetscFree(ilink));
2620:     ilink = next;
2621:   }
2622:   ts->tsrhssplit     = NULL;
2623:   ts->num_rhs_splits = 0;
2624:   if (ts->eval_times) {
2625:     PetscCall(PetscFree(ts->eval_times->time_points));
2626:     PetscCall(PetscFree(ts->eval_times->sol_times));
2627:     PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
2628:     PetscCall(PetscFree(ts->eval_times));
2629:   }
2630:   ts->rhsjacobian.time  = PETSC_MIN_REAL;
2631:   ts->rhsjacobian.scale = 1.0;
2632:   ts->ijacobian.shift   = 1.0;
2633:   ts->setupcalled       = PETSC_FALSE;
2634:   PetscFunctionReturn(PETSC_SUCCESS);
2635: }

2637: static PetscErrorCode TSResizeReset(TS);

2639: /*@
2640:   TSDestroy - Destroys the timestepper context that was created
2641:   with `TSCreate()`.

2643:   Collective

2645:   Input Parameter:
2646: . ts - the `TS` context obtained from `TSCreate()`

2648:   Level: beginner

2650: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2651: @*/
2652: PetscErrorCode TSDestroy(TS *ts)
2653: {
2654:   PetscFunctionBegin;
2655:   if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2657:   if (--((PetscObject)*ts)->refct > 0) {
2658:     *ts = NULL;
2659:     PetscFunctionReturn(PETSC_SUCCESS);
2660:   }

2662:   PetscCall(TSReset(*ts));
2663:   PetscCall(TSAdjointReset(*ts));
2664:   if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2665:   PetscCall(TSResizeReset(*ts));

2667:   /* if memory was published with SAWs then destroy it */
2668:   PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2669:   PetscTryTypeMethod(*ts, destroy);

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

2673:   PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2674:   PetscCall(TSEventDestroy(&(*ts)->event));

2676:   PetscCall(SNESDestroy(&(*ts)->snes));
2677:   PetscCall(SNESDestroy(&(*ts)->snesrhssplit));
2678:   PetscCall(DMDestroy(&(*ts)->dm));
2679:   PetscCall(TSMonitorCancel(*ts));
2680:   PetscCall(TSAdjointMonitorCancel(*ts));

2682:   PetscCall(TSDestroy(&(*ts)->quadraturets));
2683:   PetscCall(PetscHeaderDestroy(ts));
2684:   PetscFunctionReturn(PETSC_SUCCESS);
2685: }

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

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

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

2696:   Output Parameter:
2697: . snes - the nonlinear solver context

2699:   Level: beginner

2701:   Notes:
2702:   The user can then directly manipulate the `SNES` context to set various
2703:   options, etc.  Likewise, the user can then extract and manipulate the
2704:   `KSP`, and `PC` contexts as well.

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

2709: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2710: @*/
2711: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2712: {
2713:   PetscFunctionBegin;
2715:   PetscAssertPointer(snes, 2);
2716:   if (!ts->snes) {
2717:     PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2718:     PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2719:     PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2720:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2721:     if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2722:     if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2723:   }
2724:   *snes = ts->snes;
2725:   PetscFunctionReturn(PETSC_SUCCESS);
2726: }

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

2731:   Collective

2733:   Input Parameters:
2734: + ts   - the `TS` context obtained from `TSCreate()`
2735: - snes - the nonlinear solver context

2737:   Level: developer

2739:   Note:
2740:   Most users should have the `TS` created by calling `TSGetSNES()`

2742: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2743: @*/
2744: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2745: {
2746:   PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);

2748:   PetscFunctionBegin;
2751:   PetscCall(PetscObjectReference((PetscObject)snes));
2752:   PetscCall(SNESDestroy(&ts->snes));

2754:   ts->snes = snes;

2756:   PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2757:   PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2758:   if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2759:   PetscFunctionReturn(PETSC_SUCCESS);
2760: }

2762: /*@
2763:   TSGetKSP - Returns the `KSP` (linear solver) associated with
2764:   a `TS` (timestepper) context.

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

2768:   Input Parameter:
2769: . ts - the `TS` context obtained from `TSCreate()`

2771:   Output Parameter:
2772: . ksp - the nonlinear solver context

2774:   Level: beginner

2776:   Notes:
2777:   The user can then directly manipulate the `KSP` context to set various
2778:   options, etc.  Likewise, the user can then extract and manipulate the
2779:   `PC` context as well.

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

2784: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2785: @*/
2786: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2787: {
2788:   SNES snes;

2790:   PetscFunctionBegin;
2792:   PetscAssertPointer(ksp, 2);
2793:   PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2794:   PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
2795:   PetscCall(TSGetSNES(ts, &snes));
2796:   PetscCall(SNESGetKSP(snes, ksp));
2797:   PetscFunctionReturn(PETSC_SUCCESS);
2798: }

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

2802: /*@
2803:   TSSetMaxSteps - Sets the maximum number of steps to use.

2805:   Logically Collective

2807:   Input Parameters:
2808: + ts       - the `TS` context obtained from `TSCreate()`
2809: - maxsteps - maximum number of steps to use

2811:   Options Database Key:
2812: . -ts_max_steps <maxsteps> - Sets maxsteps

2814:   Level: intermediate

2816:   Note:
2817:   Use `PETSC_DETERMINE` to reset the maximum number of steps to the default from when the object's type was set

2819:   The default maximum number of steps is 5,000

2821:   Fortran Note:
2822:   Use `PETSC_DETERMINE_INTEGER`

2824: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2825: @*/
2826: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2827: {
2828:   PetscFunctionBegin;
2831:   if (maxsteps == PETSC_DETERMINE) {
2832:     ts->max_steps = ts->default_max_steps;
2833:   } else {
2834:     PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2835:     ts->max_steps = maxsteps;
2836:   }
2837:   PetscFunctionReturn(PETSC_SUCCESS);
2838: }

2840: /*@
2841:   TSGetMaxSteps - Gets the maximum number of steps to use.

2843:   Not Collective

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

2848:   Output Parameter:
2849: . maxsteps - maximum number of steps to use

2851:   Level: advanced

2853: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2854: @*/
2855: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2856: {
2857:   PetscFunctionBegin;
2859:   PetscAssertPointer(maxsteps, 2);
2860:   *maxsteps = ts->max_steps;
2861:   PetscFunctionReturn(PETSC_SUCCESS);
2862: }

2864: /*@
2865:   TSSetRunSteps - Sets the maximum number of steps to take in each call to `TSSolve()`.

2867:   If the step count when `TSSolve()` is `start_step`, this will stop the simulation once `current_step - start_step >= run_steps`.
2868:   Comparatively, `TSSetMaxSteps()` will stop if `current_step >= max_steps`.
2869:   The simulation will stop when either condition is reached.

2871:   Logically Collective

2873:   Input Parameters:
2874: + ts       - the `TS` context obtained from `TSCreate()`
2875: - runsteps - maximum number of steps to take in each call to `TSSolve()`;

2877:   Options Database Key:
2878: . -ts_run_steps <runsteps> - Sets runsteps

2880:   Level: intermediate

2882:   Note:
2883:   The default is `PETSC_UNLIMITED`

2885: .seealso: [](ch_ts), `TS`, `TSGetRunSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`, `TSSetMaxSteps()`
2886: @*/
2887: PetscErrorCode TSSetRunSteps(TS ts, PetscInt runsteps)
2888: {
2889:   PetscFunctionBegin;
2892:   if (runsteps == PETSC_DETERMINE) {
2893:     ts->run_steps = PETSC_UNLIMITED;
2894:   } else {
2895:     PetscCheck(runsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Max number of steps to take in each call to TSSolve must be non-negative");
2896:     ts->run_steps = runsteps;
2897:   }
2898:   PetscFunctionReturn(PETSC_SUCCESS);
2899: }

2901: /*@
2902:   TSGetRunSteps - Gets the maximum number of steps to take in each call to `TSSolve()`.

2904:   Not Collective

2906:   Input Parameter:
2907: . ts - the `TS` context obtained from `TSCreate()`

2909:   Output Parameter:
2910: . runsteps - maximum number of steps to take in each call to `TSSolve`.

2912:   Level: advanced

2914: .seealso: [](ch_ts), `TS`, `TSSetRunSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`, `TSGetMaxSteps()`
2915: @*/
2916: PetscErrorCode TSGetRunSteps(TS ts, PetscInt *runsteps)
2917: {
2918:   PetscFunctionBegin;
2920:   PetscAssertPointer(runsteps, 2);
2921:   *runsteps = ts->run_steps;
2922:   PetscFunctionReturn(PETSC_SUCCESS);
2923: }

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

2928:   Logically Collective

2930:   Input Parameters:
2931: + ts      - the `TS` context obtained from `TSCreate()`
2932: - maxtime - final time to step to

2934:   Options Database Key:
2935: . -ts_max_time <maxtime> - Sets maxtime

2937:   Level: intermediate

2939:   Notes:
2940:   Use `PETSC_DETERMINE` to reset the maximum time to the default from when the object's type was set

2942:   The default maximum time is 5.0

2944:   Fortran Note:
2945:   Use `PETSC_DETERMINE_REAL`

2947: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2948: @*/
2949: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2950: {
2951:   PetscFunctionBegin;
2954:   if (maxtime == PETSC_DETERMINE) {
2955:     ts->max_time = ts->default_max_time;
2956:   } else {
2957:     ts->max_time = maxtime;
2958:   }
2959:   PetscFunctionReturn(PETSC_SUCCESS);
2960: }

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

2965:   Not Collective

2967:   Input Parameter:
2968: . ts - the `TS` context obtained from `TSCreate()`

2970:   Output Parameter:
2971: . maxtime - final time to step to

2973:   Level: advanced

2975: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2976: @*/
2977: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2978: {
2979:   PetscFunctionBegin;
2981:   PetscAssertPointer(maxtime, 2);
2982:   *maxtime = ts->max_time;
2983:   PetscFunctionReturn(PETSC_SUCCESS);
2984: }

2986: // PetscClangLinter pragma disable: -fdoc-*
2987: /*@
2988:   TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.

2990:   Level: deprecated

2992: @*/
2993: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2994: {
2995:   PetscFunctionBegin;
2997:   PetscCall(TSSetTime(ts, initial_time));
2998:   PetscCall(TSSetTimeStep(ts, time_step));
2999:   PetscFunctionReturn(PETSC_SUCCESS);
3000: }

3002: // PetscClangLinter pragma disable: -fdoc-*
3003: /*@
3004:   TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.

3006:   Level: deprecated

3008: @*/
3009: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
3010: {
3011:   PetscFunctionBegin;
3013:   if (maxsteps) {
3014:     PetscAssertPointer(maxsteps, 2);
3015:     *maxsteps = ts->max_steps;
3016:   }
3017:   if (maxtime) {
3018:     PetscAssertPointer(maxtime, 3);
3019:     *maxtime = ts->max_time;
3020:   }
3021:   PetscFunctionReturn(PETSC_SUCCESS);
3022: }

3024: // PetscClangLinter pragma disable: -fdoc-*
3025: /*@
3026:   TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.

3028:   Level: deprecated

3030: @*/
3031: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
3032: {
3033:   PetscFunctionBegin;
3034:   if (maxsteps != PETSC_CURRENT) PetscCall(TSSetMaxSteps(ts, maxsteps));
3035:   if (maxtime != (PetscReal)PETSC_CURRENT) PetscCall(TSSetMaxTime(ts, maxtime));
3036:   PetscFunctionReturn(PETSC_SUCCESS);
3037: }

3039: // PetscClangLinter pragma disable: -fdoc-*
3040: /*@
3041:   TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.

3043:   Level: deprecated

3045: @*/
3046: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
3047: {
3048:   return TSGetStepNumber(ts, steps);
3049: }

3051: // PetscClangLinter pragma disable: -fdoc-*
3052: /*@
3053:   TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.

3055:   Level: deprecated

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

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

3067:   Logically Collective

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

3073:   Level: beginner

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

3081:   PetscFunctionBegin;
3084:   PetscCall(PetscObjectReference((PetscObject)u));
3085:   PetscCall(VecDestroy(&ts->vec_sol));
3086:   ts->vec_sol = u;

3088:   PetscCall(TSGetDM(ts, &dm));
3089:   PetscCall(DMShellSetGlobalVector(dm, u));
3090:   PetscFunctionReturn(PETSC_SUCCESS);
3091: }

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

3097:   Logically Collective

3099:   Input Parameters:
3100: + ts   - The `TS` context obtained from `TSCreate()`
3101: - func - The function

3103:   Calling sequence of `func`:
3104: . ts - the `TS` context

3106:   Level: intermediate

3108: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3109: @*/
3110: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
3111: {
3112:   PetscFunctionBegin;
3114:   ts->prestep = func;
3115:   PetscFunctionReturn(PETSC_SUCCESS);
3116: }

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

3121:   Collective

3123:   Input Parameter:
3124: . ts - The `TS` context obtained from `TSCreate()`

3126:   Level: developer

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

3132: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3133: @*/
3134: PetscErrorCode TSPreStep(TS ts)
3135: {
3136:   PetscFunctionBegin;
3138:   if (ts->prestep) {
3139:     Vec              U;
3140:     PetscObjectId    idprev;
3141:     PetscBool        sameObject;
3142:     PetscObjectState sprev, spost;

3144:     PetscCall(TSGetSolution(ts, &U));
3145:     PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3146:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3147:     PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3148:     PetscCall(TSGetSolution(ts, &U));
3149:     PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3150:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3151:     if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3152:   }
3153:   PetscFunctionReturn(PETSC_SUCCESS);
3154: }

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

3160:   Logically Collective

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

3166:   Calling sequence of `func`:
3167: + ts        - the `TS` context
3168: - stagetime - the stage time

3170:   Level: intermediate

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

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

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

3191:   Logically Collective

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

3197:   Calling sequence of `func`:
3198: + ts         - the `TS` context
3199: . stagetime  - the stage time
3200: . stageindex - the stage index
3201: - Y          - Array of vectors (of size = total number of stages) with the stage solutions

3203:   Level: intermediate

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

3210: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3211: @*/
3212: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
3213: {
3214:   PetscFunctionBegin;
3216:   ts->poststage = func;
3217:   PetscFunctionReturn(PETSC_SUCCESS);
3218: }

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

3224:   Logically Collective

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

3230:   Calling sequence of `func`:
3231: . ts - the `TS` context

3233:   Level: intermediate

3235:   Note:
3236:   The function set by `TSSetPostEvaluate()` is called after the solution is evaluated, or after the step rollback.
3237:   Inside the `func` callback, the solution vector can be obtained with `TSGetSolution()`, and modified, if need be.
3238:   The time step can be obtained with `TSGetTimeStep()`, and the time at the start of the step - via `TSGetTime()`.
3239:   The potential changes to the solution vector introduced by event handling (`postevent()`) are not relevant for `TSSetPostEvaluate()`,
3240:   but are relevant for `TSSetPostStep()`, according to the function call scheme in `TSSolve()`, as shown below
3241: .vb
3242:   ...
3243:   Step()
3244:   PostEvaluate()
3245:   EventHandling()
3246:   step_rollback ? PostEvaluate() : PostStep()
3247:   ...
3248: .ve
3249:   where EventHandling() may result in one of the following three outcomes
3250: .vb
3251:   (1) | successful step | solution intact
3252:   (2) | successful step | solution modified by `postevent()`
3253:   (3) | step_rollback   | solution rolled back
3254: .ve

3256: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3257: @*/
3258: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
3259: {
3260:   PetscFunctionBegin;
3262:   ts->postevaluate = func;
3263:   PetscFunctionReturn(PETSC_SUCCESS);
3264: }

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

3269:   Collective

3271:   Input Parameters:
3272: + ts        - The `TS` context obtained from `TSCreate()`
3273: - stagetime - The absolute time of the current stage

3275:   Level: developer

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

3281: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3282: @*/
3283: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3284: {
3285:   PetscFunctionBegin;
3287:   if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3288:   PetscFunctionReturn(PETSC_SUCCESS);
3289: }

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

3294:   Collective

3296:   Input Parameters:
3297: + ts         - The `TS` context obtained from `TSCreate()`
3298: . stagetime  - The absolute time of the current stage
3299: . stageindex - Stage number
3300: - Y          - Array of vectors (of size = total number of stages) with the stage solutions

3302:   Level: developer

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

3308: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3309: @*/
3310: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3311: {
3312:   PetscFunctionBegin;
3314:   if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3315:   PetscFunctionReturn(PETSC_SUCCESS);
3316: }

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

3321:   Collective

3323:   Input Parameter:
3324: . ts - The `TS` context obtained from `TSCreate()`

3326:   Level: developer

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

3332: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3333: @*/
3334: PetscErrorCode TSPostEvaluate(TS ts)
3335: {
3336:   PetscFunctionBegin;
3338:   if (ts->postevaluate) {
3339:     Vec              U;
3340:     PetscObjectState sprev, spost;

3342:     PetscCall(TSGetSolution(ts, &U));
3343:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3344:     PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3345:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3346:     if (sprev != spost) PetscCall(TSRestartStep(ts));
3347:   }
3348:   PetscFunctionReturn(PETSC_SUCCESS);
3349: }

3351: /*@C
3352:   TSSetPostStep - Sets the general-purpose function
3353:   called once at the end of each successful time step.

3355:   Logically Collective

3357:   Input Parameters:
3358: + ts   - The `TS` context obtained from `TSCreate()`
3359: - func - The function

3361:   Calling sequence of `func`:
3362: . ts - the `TS` context

3364:   Level: intermediate

3366:   Note:
3367:   The function set by `TSSetPostStep()` is called after each successful step. If the event handler locates an event at the
3368:   given step, and `postevent()` modifies the solution vector, the solution vector obtained by `TSGetSolution()` inside `func` will
3369:   contain the changes. To get the solution without these changes, use `TSSetPostEvaluate()` to set the appropriate callback.
3370:   The scheme of the relevant function calls in `TSSolve()` is shown below
3371: .vb
3372:   ...
3373:   Step()
3374:   PostEvaluate()
3375:   EventHandling()
3376:   step_rollback ? PostEvaluate() : PostStep()
3377:   ...
3378: .ve
3379:   where EventHandling() may result in one of the following three outcomes
3380: .vb
3381:   (1) | successful step | solution intact
3382:   (2) | successful step | solution modified by `postevent()`
3383:   (3) | step_rollback   | solution rolled back
3384: .ve

3386: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3387: @*/
3388: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
3389: {
3390:   PetscFunctionBegin;
3392:   ts->poststep = func;
3393:   PetscFunctionReturn(PETSC_SUCCESS);
3394: }

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

3399:   Collective

3401:   Input Parameter:
3402: . ts - The `TS` context obtained from `TSCreate()`

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

3408:   Level: developer

3410: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPostStep()`
3411: @*/
3412: PetscErrorCode TSPostStep(TS ts)
3413: {
3414:   PetscFunctionBegin;
3416:   if (ts->poststep) {
3417:     Vec              U;
3418:     PetscObjectId    idprev;
3419:     PetscBool        sameObject;
3420:     PetscObjectState sprev, spost;

3422:     PetscCall(TSGetSolution(ts, &U));
3423:     PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3424:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3425:     PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3426:     PetscCall(TSGetSolution(ts, &U));
3427:     PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3428:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3429:     if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3430:   }
3431:   PetscFunctionReturn(PETSC_SUCCESS);
3432: }

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

3437:   Collective

3439:   Input Parameters:
3440: + ts - time stepping context
3441: - t  - time to interpolate to

3443:   Output Parameter:
3444: . U - state at given time

3446:   Level: intermediate

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

3451: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3452: @*/
3453: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3454: {
3455:   PetscFunctionBegin;
3458:   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);
3459:   PetscUseTypeMethod(ts, interpolate, t, U);
3460:   PetscFunctionReturn(PETSC_SUCCESS);
3461: }

3463: /*@
3464:   TSStep - Steps one time step

3466:   Collective

3468:   Input Parameter:
3469: . ts - the `TS` context obtained from `TSCreate()`

3471:   Level: developer

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

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

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

3482: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3483: @*/
3484: PetscErrorCode TSStep(TS ts)
3485: {
3486:   static PetscBool cite = PETSC_FALSE;
3487:   PetscReal        ptime;

3489:   PetscFunctionBegin;
3491:   PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3492:                                    "  title         = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3493:                                    "  author        = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3494:                                    "  journal       = {arXiv e-preprints},\n"
3495:                                    "  eprint        = {1806.01437},\n"
3496:                                    "  archivePrefix = {arXiv},\n"
3497:                                    "  year          = {2018}\n}\n",
3498:                                    &cite));
3499:   PetscCall(TSSetUp(ts));
3500:   PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3501:   if (ts->eval_times)
3502:     ts->eval_times->worktol = 0; /* In each step of TSSolve() 'eval_times->worktol' will be meaningfully defined (later) only once:
3503:                                                    in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */

3505:   PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->run_steps != PETSC_INT_MAX || ts->max_steps != PETSC_INT_MAX, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime(), TSSetMaxSteps(), or TSSetRunSteps() or use -ts_max_time <time>, -ts_max_steps <steps>, -ts_run_steps <steps>");
3506:   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()");
3507:   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");

3509:   if (!ts->vec_sol0) PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sol0));
3510:   PetscCall(VecCopy(ts->vec_sol, ts->vec_sol0));
3511:   ts->time_step0 = ts->time_step;

3513:   if (!ts->steps) ts->ptime_prev = ts->ptime;
3514:   ptime = ts->ptime;

3516:   ts->ptime_prev_rollback = ts->ptime_prev;
3517:   ts->reason              = TS_CONVERGED_ITERATING;

3519:   PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3520:   PetscUseTypeMethod(ts, step);
3521:   PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));

3523:   if (ts->reason >= 0) {
3524:     ts->ptime_prev = ptime;
3525:     ts->steps++;
3526:     ts->steprollback = PETSC_FALSE;
3527:     ts->steprestart  = PETSC_FALSE;
3528:     ts->stepresize   = PETSC_FALSE;
3529:   }

3531:   if (ts->reason < 0 && ts->errorifstepfailed) {
3532:     PetscCall(TSMonitorCancel(ts));
3533:     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 use unlimited to attempt recovery", TSConvergedReasons[ts->reason]);
3534:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3535:   }
3536:   PetscFunctionReturn(PETSC_SUCCESS);
3537: }

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

3543:   Collective

3545:   Input Parameters:
3546: + ts        - time stepping context
3547: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

3549:   Input/Output Parameter:
3550: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3551:            on output, the actual order of the error evaluation

3553:   Output Parameter:
3554: . wlte - the weighted local truncation error norm

3556:   Level: advanced

3558:   Note:
3559:   If the timestepper cannot evaluate the error in a particular step
3560:   (eg. in the first step or restart steps after event handling),
3561:   this routine returns wlte=-1.0 .

3563: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3564: @*/
3565: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3566: {
3567:   PetscFunctionBegin;
3571:   if (order) PetscAssertPointer(order, 3);
3573:   PetscAssertPointer(wlte, 4);
3574:   PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3575:   PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3576:   PetscFunctionReturn(PETSC_SUCCESS);
3577: }

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

3582:   Collective

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

3589:   Output Parameter:
3590: . U - state at the end of the current step

3592:   Level: advanced

3594:   Notes:
3595:   This function cannot be called until all stages have been evaluated.

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

3599: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3600: @*/
3601: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3602: {
3603:   PetscFunctionBegin;
3607:   PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3608:   PetscFunctionReturn(PETSC_SUCCESS);
3609: }

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

3614:   Not collective

3616:   Input Parameter:
3617: . ts - time stepping context

3619:   Output Parameter:
3620: . initCondition - The function which computes an initial condition

3622:   Calling sequence of `initCondition`:
3623: + ts - The timestepping context
3624: - u  - The input vector in which the initial condition is stored

3626:   Level: advanced

3628: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3629: @*/
3630: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3631: {
3632:   PetscFunctionBegin;
3634:   PetscAssertPointer(initCondition, 2);
3635:   *initCondition = ts->ops->initcondition;
3636:   PetscFunctionReturn(PETSC_SUCCESS);
3637: }

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

3642:   Logically collective

3644:   Input Parameters:
3645: + ts            - time stepping context
3646: - initCondition - The function which computes an initial condition

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

3652:   Level: advanced

3654: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3655: @*/
3656: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3657: {
3658:   PetscFunctionBegin;
3661:   ts->ops->initcondition = initCondition;
3662:   PetscFunctionReturn(PETSC_SUCCESS);
3663: }

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

3668:   Collective

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

3674:   Level: advanced

3676: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3677: @*/
3678: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3679: {
3680:   PetscFunctionBegin;
3683:   PetscTryTypeMethod(ts, initcondition, u);
3684:   PetscFunctionReturn(PETSC_SUCCESS);
3685: }

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

3690:   Not collective

3692:   Input Parameter:
3693: . ts - time stepping context

3695:   Output Parameter:
3696: . exactError - The function which computes the solution error

3698:   Calling sequence of `exactError`:
3699: + ts - The timestepping context
3700: . u  - The approximate solution vector
3701: - e  - The vector in which the error is stored

3703:   Level: advanced

3705: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3706: @*/
3707: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3708: {
3709:   PetscFunctionBegin;
3711:   PetscAssertPointer(exactError, 2);
3712:   *exactError = ts->ops->exacterror;
3713:   PetscFunctionReturn(PETSC_SUCCESS);
3714: }

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

3719:   Logically collective

3721:   Input Parameters:
3722: + ts         - time stepping context
3723: - exactError - The function which computes the solution error

3725:   Calling sequence of `exactError`:
3726: + ts - The timestepping context
3727: . u  - The approximate solution vector
3728: - e  - The  vector in which the error is stored

3730:   Level: advanced

3732: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3733: @*/
3734: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3735: {
3736:   PetscFunctionBegin;
3739:   ts->ops->exacterror = exactError;
3740:   PetscFunctionReturn(PETSC_SUCCESS);
3741: }

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

3746:   Collective

3748:   Input Parameters:
3749: + ts - time stepping context
3750: . u  - The approximate solution
3751: - e  - The `Vec` used to store the error

3753:   Level: advanced

3755: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3756: @*/
3757: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3758: {
3759:   PetscFunctionBegin;
3763:   PetscTryTypeMethod(ts, exacterror, u, e);
3764:   PetscFunctionReturn(PETSC_SUCCESS);
3765: }

3767: /*@C
3768:   TSSetResize - Sets the resize callbacks.

3770:   Logically Collective

3772:   Input Parameters:
3773: + ts       - The `TS` context obtained from `TSCreate()`
3774: . rollback - Whether a resize will restart the step
3775: . setup    - The setup function
3776: . transfer - The transfer function
3777: - ctx      - [optional] The user-defined context

3779:   Calling sequence of `setup`:
3780: + ts     - the `TS` context
3781: . step   - the current step
3782: . time   - the current time
3783: . state  - the current vector of state
3784: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3785: - ctx    - user defined context

3787:   Calling sequence of `transfer`:
3788: + ts      - the `TS` context
3789: . nv      - the number of vectors to be transferred
3790: . vecsin  - array of vectors to be transferred
3791: . vecsout - array of transferred vectors
3792: - ctx     - user defined context

3794:   Notes:
3795:   The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
3796:   depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
3797:   and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
3798:   that the size of the discrete problem has changed.
3799:   In both cases, the solver will collect the needed vectors that will be
3800:   transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
3801:   current solution vector, and other vectors needed by the specific solver used.
3802:   For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3803:   Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3804:   will be automatically reset if the sizes are changed and they must be specified again by the user
3805:   inside the `transfer` function.
3806:   The input and output arrays passed to `transfer` are allocated by PETSc.
3807:   Vectors in `vecsout` must be created by the user.
3808:   Ownership of vectors in `vecsout` is transferred to PETSc.

3810:   Level: advanced

3812: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3813: @*/
3814: PetscErrorCode TSSetResize(TS ts, PetscBool rollback, 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)
3815: {
3816:   PetscFunctionBegin;
3818:   ts->resizerollback = rollback;
3819:   ts->resizesetup    = setup;
3820:   ts->resizetransfer = transfer;
3821:   ts->resizectx      = ctx;
3822:   PetscFunctionReturn(PETSC_SUCCESS);
3823: }

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

3828:   Collective

3830:   Input Parameters:
3831: + ts   - The `TS` context obtained from `TSCreate()`
3832: - 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.

3834:   Level: developer

3836:   Note:
3837:   `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3838:    used within time stepping implementations,
3839:    so most users would not generally call this routine themselves.

3841: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3842: @*/
3843: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3844: {
3845:   PetscFunctionBegin;
3847:   PetscTryTypeMethod(ts, resizeregister, flg);
3848:   /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3849:   PetscFunctionReturn(PETSC_SUCCESS);
3850: }

3852: static PetscErrorCode TSResizeReset(TS ts)
3853: {
3854:   PetscFunctionBegin;
3856:   PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3857:   PetscFunctionReturn(PETSC_SUCCESS);
3858: }

3860: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3861: {
3862:   PetscFunctionBegin;
3865:   for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3866:   if (ts->resizetransfer) {
3867:     PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3868:     PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3869:   }
3870:   for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3871:   PetscFunctionReturn(PETSC_SUCCESS);
3872: }

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

3877:   Collective

3879:   Input Parameters:
3880: + ts   - The `TS` context obtained from `TSCreate()`
3881: . name - A string identifying the vector
3882: - vec  - The vector

3884:   Level: developer

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

3890: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3891: @*/
3892: PetscErrorCode TSResizeRegisterVec(TS ts, const char name[], Vec vec)
3893: {
3894:   PetscFunctionBegin;
3896:   PetscAssertPointer(name, 2);
3898:   PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3899:   PetscFunctionReturn(PETSC_SUCCESS);
3900: }

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

3905:   Collective

3907:   Input Parameters:
3908: + ts   - The `TS` context obtained from `TSCreate()`
3909: . name - A string identifying the vector
3910: - vec  - The vector

3912:   Level: developer

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

3918: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
3919: @*/
3920: PetscErrorCode TSResizeRetrieveVec(TS ts, const char name[], Vec *vec)
3921: {
3922:   PetscFunctionBegin;
3924:   PetscAssertPointer(name, 2);
3925:   PetscAssertPointer(vec, 3);
3926:   PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
3927:   PetscFunctionReturn(PETSC_SUCCESS);
3928: }

3930: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
3931: {
3932:   PetscInt        cnt;
3933:   PetscObjectList tmp;
3934:   Vec            *vecsin  = NULL;
3935:   const char    **namesin = NULL;

3937:   PetscFunctionBegin;
3938:   for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
3939:     if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
3940:   if (names) PetscCall(PetscMalloc1(cnt, &namesin));
3941:   if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
3942:   for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
3943:     if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
3944:       if (vecs) vecsin[cnt] = (Vec)tmp->obj;
3945:       if (names) namesin[cnt] = tmp->name;
3946:       cnt++;
3947:     }
3948:   }
3949:   if (nv) *nv = cnt;
3950:   if (names) *names = namesin;
3951:   if (vecs) *vecs = vecsin;
3952:   PetscFunctionReturn(PETSC_SUCCESS);
3953: }

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

3958:   Collective

3960:   Input Parameter:
3961: . ts - The `TS` context obtained from `TSCreate()`

3963:   Level: developer

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

3969: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3970: @*/
3971: PetscErrorCode TSResize(TS ts)
3972: {
3973:   PetscInt     nv      = 0;
3974:   const char **names   = NULL;
3975:   Vec         *vecsin  = NULL;
3976:   const char  *solname = "ts:vec_sol";

3978:   PetscFunctionBegin;
3980:   if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
3981:   if (ts->resizesetup) {
3982:     PetscCall(VecLockReadPush(ts->vec_sol));
3983:     PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &ts->stepresize, ts->resizectx));
3984:     PetscCall(VecLockReadPop(ts->vec_sol));
3985:     if (ts->stepresize) {
3986:       if (ts->resizerollback) {
3987:         PetscCall(TSRollBack(ts));
3988:         ts->time_step = ts->time_step0;
3989:       }
3990:       PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
3991:       PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
3992:     }
3993:   }

3995:   PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
3996:   if (nv) {
3997:     Vec *vecsout, vecsol;

3999:     /* Reset internal objects */
4000:     PetscCall(TSReset(ts));

4002:     /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
4003:     PetscCall(PetscCalloc1(nv, &vecsout));
4004:     PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
4005:     for (PetscInt i = 0; i < nv; i++) {
4006:       const char *name;
4007:       char       *oname;

4009:       PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
4010:       PetscCall(PetscStrallocpy(name, &oname));
4011:       PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
4012:       if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
4013:       PetscCall(PetscFree(oname));
4014:       PetscCall(VecDestroy(&vecsout[i]));
4015:     }
4016:     PetscCall(PetscFree(vecsout));
4017:     PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */

4019:     PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
4020:     if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
4021:     PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
4022:   }

4024:   PetscCall(PetscFree(names));
4025:   PetscCall(PetscFree(vecsin));
4026:   PetscCall(TSResizeReset(ts));
4027:   PetscFunctionReturn(PETSC_SUCCESS);
4028: }

4030: /*@
4031:   TSSolve - Steps the requested number of timesteps.

4033:   Collective

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

4040:   Level: beginner

4042:   Notes:
4043:   The final time returned by this function may be different from the time of the internally
4044:   held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
4045:   stepped over the final time.

4047: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
4048: @*/
4049: PetscErrorCode TSSolve(TS ts, Vec u)
4050: {
4051:   Vec solution;

4053:   PetscFunctionBegin;

4057:   PetscCall(TSSetExactFinalTimeDefault(ts));
4058:   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 */
4059:     if (!ts->vec_sol || u == ts->vec_sol) {
4060:       PetscCall(VecDuplicate(u, &solution));
4061:       PetscCall(TSSetSolution(ts, solution));
4062:       PetscCall(VecDestroy(&solution)); /* grant ownership */
4063:     }
4064:     PetscCall(VecCopy(u, ts->vec_sol));
4065:     PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
4066:   } else if (u) PetscCall(TSSetSolution(ts, u));
4067:   PetscCall(TSSetUp(ts));
4068:   PetscCall(TSTrajectorySetUp(ts->trajectory, ts));

4070:   PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->run_steps != PETSC_INT_MAX || ts->max_steps != PETSC_INT_MAX, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime(), TSSetMaxSteps(), or TSSetRunSteps() or use -ts_max_time <time>, -ts_max_steps <steps>, -ts_run_steps <steps>");
4071:   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()");
4072:   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");
4073:   PetscCheck(!(ts->eval_times && ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP), PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "You must use TS_EXACTFINALTIME_MATCHSTEP when using time span or evaluation times");

4075:   if (ts->eval_times) {
4076:     if (!ts->eval_times->sol_vecs) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
4077:     for (PetscInt i = 0; i < ts->eval_times->num_time_points; i++) {
4078:       PetscBool is_close = PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[i], ts->eval_times->reltol * ts->time_step + ts->eval_times->abstol, 0);
4079:       if (ts->ptime <= ts->eval_times->time_points[i] || is_close) {
4080:         ts->eval_times->time_point_idx = i;

4082:         PetscBool is_ptime_in_sol_times = PETSC_FALSE; // If current solution has already been saved, we should not save it again
4083:         if (ts->eval_times->sol_idx > 0) is_ptime_in_sol_times = PetscIsCloseAtTol(ts->ptime, ts->eval_times->sol_times[ts->eval_times->sol_idx - 1], ts->eval_times->reltol * ts->time_step + ts->eval_times->abstol, 0);
4084:         if (is_close && !is_ptime_in_sol_times) {
4085:           PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4086:           ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4087:           ts->eval_times->sol_idx++;
4088:           ts->eval_times->time_point_idx++;
4089:         }
4090:         break;
4091:       }
4092:     }
4093:   }

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

4097:   /* reset number of steps only when the step is not restarted. ARKIMEX
4098:      restarts the step after an event. Resetting these counters in such case causes
4099:      TSTrajectory to incorrectly save the output files
4100:   */
4101:   /* reset time step and iteration counters */
4102:   if (!ts->steps) {
4103:     ts->ksp_its           = 0;
4104:     ts->snes_its          = 0;
4105:     ts->num_snes_failures = 0;
4106:     ts->reject            = 0;
4107:     ts->steprestart       = PETSC_TRUE;
4108:     ts->steprollback      = PETSC_FALSE;
4109:     ts->stepresize        = PETSC_FALSE;
4110:     ts->rhsjacobian.time  = PETSC_MIN_REAL;
4111:   }

4113:   /* make sure initial time step does not overshoot final time or the next point in evaluation times */
4114:   if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
4115:     PetscReal maxdt;
4116:     PetscReal dt = ts->time_step;

4118:     if (ts->eval_times) maxdt = ts->eval_times->time_points[ts->eval_times->time_point_idx] - ts->ptime;
4119:     else maxdt = ts->max_time - ts->ptime;
4120:     ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
4121:   }
4122:   ts->reason = TS_CONVERGED_ITERATING;

4124:   {
4125:     PetscViewer       viewer;
4126:     PetscViewerFormat format;
4127:     PetscBool         flg;
4128:     static PetscBool  incall = PETSC_FALSE;

4130:     if (!incall) {
4131:       /* Estimate the convergence rate of the time discretization */
4132:       PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
4133:       if (flg) {
4134:         PetscConvEst conv;
4135:         DM           dm;
4136:         PetscReal   *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4137:         PetscInt     Nf;
4138:         PetscBool    checkTemporal = PETSC_TRUE;

4140:         incall = PETSC_TRUE;
4141:         PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
4142:         PetscCall(TSGetDM(ts, &dm));
4143:         PetscCall(DMGetNumFields(dm, &Nf));
4144:         PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
4145:         PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
4146:         PetscCall(PetscConvEstUseTS(conv, checkTemporal));
4147:         PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
4148:         PetscCall(PetscConvEstSetFromOptions(conv));
4149:         PetscCall(PetscConvEstSetUp(conv));
4150:         PetscCall(PetscConvEstGetConvRate(conv, alpha));
4151:         PetscCall(PetscViewerPushFormat(viewer, format));
4152:         PetscCall(PetscConvEstRateView(conv, alpha, viewer));
4153:         PetscCall(PetscViewerPopFormat(viewer));
4154:         PetscCall(PetscViewerDestroy(&viewer));
4155:         PetscCall(PetscConvEstDestroy(&conv));
4156:         PetscCall(PetscFree(alpha));
4157:         incall = PETSC_FALSE;
4158:       }
4159:     }
4160:   }

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

4164:   if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4165:     PetscUseTypeMethod(ts, solve);
4166:     if (u) PetscCall(VecCopy(ts->vec_sol, u));
4167:     ts->solvetime = ts->ptime;
4168:     solution      = ts->vec_sol;
4169:   } else { /* Step the requested number of timesteps. */
4170:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4171:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;

4173:     if (!ts->steps) {
4174:       PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4175:       PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4176:     }

4178:     ts->start_step = ts->steps; // records starting step
4179:     while (!ts->reason) {
4180:       PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4181:       if (!ts->steprollback || (ts->stepresize && ts->resizerollback)) PetscCall(TSPreStep(ts));
4182:       PetscCall(TSStep(ts));
4183:       if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4184:       if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4185:       if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4186:         if (ts->reason >= 0) ts->steps--;            /* Revert the step number changed by TSStep() */
4187:         PetscCall(TSForwardCostIntegral(ts));
4188:         if (ts->reason >= 0) ts->steps++;
4189:       }
4190:       if (ts->forward_solve) {            /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4191:         if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4192:         PetscCall(TSForwardStep(ts));
4193:         if (ts->reason >= 0) ts->steps++;
4194:       }
4195:       PetscCall(TSPostEvaluate(ts));
4196:       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. */
4197:       if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4198:       if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
4199:       /* check convergence */
4200:       if (!ts->reason) {
4201:         if ((ts->steps - ts->start_step) >= ts->run_steps) ts->reason = TS_CONVERGED_ITS;
4202:         else if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4203:         else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4204:       }
4205:       if (!ts->steprollback) {
4206:         PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4207:         PetscCall(TSPostStep(ts));
4208:         if (!ts->resizerollback) PetscCall(TSResize(ts));

4210:         if (ts->eval_times && ts->eval_times->time_point_idx < ts->eval_times->num_time_points && ts->reason >= 0) {
4211:           PetscCheck(ts->eval_times->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(eval_times->worktol > 0) in TSSolve()");
4212:           if (PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[ts->eval_times->time_point_idx], ts->eval_times->worktol, 0)) {
4213:             ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4214:             PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4215:             ts->eval_times->sol_idx++;
4216:             ts->eval_times->time_point_idx++;
4217:           }
4218:         }
4219:       }
4220:     }
4221:     PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));

4223:     if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4224:       if (!u) u = ts->vec_sol;
4225:       PetscCall(TSInterpolate(ts, ts->max_time, u));
4226:       ts->solvetime = ts->max_time;
4227:       solution      = u;
4228:       PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4229:     } else {
4230:       if (u) PetscCall(VecCopy(ts->vec_sol, u));
4231:       ts->solvetime = ts->ptime;
4232:       solution      = ts->vec_sol;
4233:     }
4234:   }

4236:   PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4237:   PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4238:   PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4239:   if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4240:   PetscFunctionReturn(PETSC_SUCCESS);
4241: }

4243: /*@
4244:   TSGetTime - Gets the time of the most recently completed step.

4246:   Not Collective

4248:   Input Parameter:
4249: . ts - the `TS` context obtained from `TSCreate()`

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

4254:   Level: beginner

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

4260: .seealso: [](ch_ts), `TS`, `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4261: @*/
4262: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4263: {
4264:   PetscFunctionBegin;
4266:   PetscAssertPointer(t, 2);
4267:   *t = ts->ptime;
4268:   PetscFunctionReturn(PETSC_SUCCESS);
4269: }

4271: /*@
4272:   TSGetPrevTime - Gets the starting time of the previously completed step.

4274:   Not Collective

4276:   Input Parameter:
4277: . ts - the `TS` context obtained from `TSCreate()`

4279:   Output Parameter:
4280: . t - the previous time

4282:   Level: beginner

4284: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4285: @*/
4286: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4287: {
4288:   PetscFunctionBegin;
4290:   PetscAssertPointer(t, 2);
4291:   *t = ts->ptime_prev;
4292:   PetscFunctionReturn(PETSC_SUCCESS);
4293: }

4295: /*@
4296:   TSSetTime - Allows one to reset the time.

4298:   Logically Collective

4300:   Input Parameters:
4301: + ts - the `TS` context obtained from `TSCreate()`
4302: - t  - the time

4304:   Level: intermediate

4306: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4307: @*/
4308: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4309: {
4310:   PetscFunctionBegin;
4313:   ts->ptime = t;
4314:   PetscFunctionReturn(PETSC_SUCCESS);
4315: }

4317: /*@
4318:   TSSetOptionsPrefix - Sets the prefix used for searching for all
4319:   TS options in the database.

4321:   Logically Collective

4323:   Input Parameters:
4324: + ts     - The `TS` context
4325: - prefix - The prefix to prepend to all option names

4327:   Level: advanced

4329:   Note:
4330:   A hyphen (-) must NOT be given at the beginning of the prefix name.
4331:   The first character of all runtime options is AUTOMATICALLY the
4332:   hyphen.

4334: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4335: @*/
4336: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4337: {
4338:   SNES snes;

4340:   PetscFunctionBegin;
4342:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4343:   PetscCall(TSGetSNES(ts, &snes));
4344:   PetscCall(SNESSetOptionsPrefix(snes, prefix));
4345:   PetscFunctionReturn(PETSC_SUCCESS);
4346: }

4348: /*@
4349:   TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4350:   TS options in the database.

4352:   Logically Collective

4354:   Input Parameters:
4355: + ts     - The `TS` context
4356: - prefix - The prefix to prepend to all option names

4358:   Level: advanced

4360:   Note:
4361:   A hyphen (-) must NOT be given at the beginning of the prefix name.
4362:   The first character of all runtime options is AUTOMATICALLY the
4363:   hyphen.

4365: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4366: @*/
4367: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4368: {
4369:   SNES snes;

4371:   PetscFunctionBegin;
4373:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4374:   PetscCall(TSGetSNES(ts, &snes));
4375:   PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4376:   PetscFunctionReturn(PETSC_SUCCESS);
4377: }

4379: /*@
4380:   TSGetOptionsPrefix - Sets the prefix used for searching for all
4381:   `TS` options in the database.

4383:   Not Collective

4385:   Input Parameter:
4386: . ts - The `TS` context

4388:   Output Parameter:
4389: . prefix - A pointer to the prefix string used

4391:   Level: intermediate

4393: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4394: @*/
4395: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4396: {
4397:   PetscFunctionBegin;
4399:   PetscAssertPointer(prefix, 2);
4400:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4401:   PetscFunctionReturn(PETSC_SUCCESS);
4402: }

4404: /*@C
4405:   TSGetRHSJacobian - Returns the Jacobian J at the present timestep.

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

4409:   Input Parameter:
4410: . ts - The `TS` context obtained from `TSCreate()`

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

4418:   Level: intermediate

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

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

4425: @*/
4426: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, void **ctx)
4427: {
4428:   DM dm;

4430:   PetscFunctionBegin;
4431:   if (Amat || Pmat) {
4432:     SNES snes;
4433:     PetscCall(TSGetSNES(ts, &snes));
4434:     PetscCall(SNESSetUpMatrices(snes));
4435:     PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4436:   }
4437:   PetscCall(TSGetDM(ts, &dm));
4438:   PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4439:   PetscFunctionReturn(PETSC_SUCCESS);
4440: }

4442: /*@C
4443:   TSGetIJacobian - Returns the implicit Jacobian at the present timestep.

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

4447:   Input Parameter:
4448: . ts - The `TS` context obtained from `TSCreate()`

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

4456:   Level: advanced

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

4461: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4462: @*/
4463: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, void **ctx)
4464: {
4465:   DM dm;

4467:   PetscFunctionBegin;
4468:   if (Amat || Pmat) {
4469:     SNES snes;
4470:     PetscCall(TSGetSNES(ts, &snes));
4471:     PetscCall(SNESSetUpMatrices(snes));
4472:     PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4473:   }
4474:   PetscCall(TSGetDM(ts, &dm));
4475:   PetscCall(DMTSGetIJacobian(dm, f, ctx));
4476:   PetscFunctionReturn(PETSC_SUCCESS);
4477: }

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

4483:   Logically Collective

4485:   Input Parameters:
4486: + ts - the `TS` integrator object
4487: - dm - the dm, cannot be `NULL`

4489:   Level: intermediate

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

4496: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4497: @*/
4498: PetscErrorCode TSSetDM(TS ts, DM dm)
4499: {
4500:   SNES snes;
4501:   DMTS tsdm;

4503:   PetscFunctionBegin;
4506:   PetscCall(PetscObjectReference((PetscObject)dm));
4507:   if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4508:     if (ts->dm->dmts && !dm->dmts) {
4509:       PetscCall(DMCopyDMTS(ts->dm, dm));
4510:       PetscCall(DMGetDMTS(ts->dm, &tsdm));
4511:       /* Grant write privileges to the replacement DM */
4512:       if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4513:     }
4514:     PetscCall(DMDestroy(&ts->dm));
4515:   }
4516:   ts->dm = dm;

4518:   PetscCall(TSGetSNES(ts, &snes));
4519:   PetscCall(SNESSetDM(snes, dm));
4520:   PetscFunctionReturn(PETSC_SUCCESS);
4521: }

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

4526:   Not Collective

4528:   Input Parameter:
4529: . ts - the `TS`

4531:   Output Parameter:
4532: . dm - the `DM`

4534:   Level: intermediate

4536: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4537: @*/
4538: PetscErrorCode TSGetDM(TS ts, DM *dm)
4539: {
4540:   PetscFunctionBegin;
4542:   if (!ts->dm) {
4543:     PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4544:     if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4545:   }
4546:   *dm = ts->dm;
4547:   PetscFunctionReturn(PETSC_SUCCESS);
4548: }

4550: /*@
4551:   SNESTSFormFunction - Function to evaluate nonlinear residual defined by an ODE solver algorithm implemented within `TS`

4553:   Logically Collective

4555:   Input Parameters:
4556: + snes - nonlinear solver
4557: . U    - the current state at which to evaluate the residual
4558: - ctx  - user context, must be a `TS`

4560:   Output Parameter:
4561: . F - the nonlinear residual

4563:   Level: developer

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

4569: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4570: @*/
4571: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4572: {
4573:   TS ts = (TS)ctx;

4575:   PetscFunctionBegin;
4580:   PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4581:   PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4582:   PetscFunctionReturn(PETSC_SUCCESS);
4583: }

4585: /*@
4586:   SNESTSFormJacobian - Function to evaluate the Jacobian defined by an ODE solver algorithm implemented within `TS`

4588:   Collective

4590:   Input Parameters:
4591: + snes - nonlinear solver
4592: . U    - the current state at which to evaluate the residual
4593: - ctx  - user context, must be a `TS`

4595:   Output Parameters:
4596: + A - the Jacobian
4597: - B - the matrix used to construct the preconditioner (often the same as `A`)

4599:   Level: developer

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

4604: .seealso: [](ch_ts), `SNESSetJacobian()`
4605: @*/
4606: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4607: {
4608:   TS ts = (TS)ctx;

4610:   PetscFunctionBegin;
4616:   PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4617:   PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4618:   PetscFunctionReturn(PETSC_SUCCESS);
4619: }

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

4624:   Collective

4626:   Input Parameters:
4627: + ts  - time stepping context
4628: . t   - time at which to evaluate
4629: . U   - state at which to evaluate
4630: - ctx - context

4632:   Output Parameter:
4633: . F - right-hand side

4635:   Level: intermediate

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

4641: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4642: @*/
4643: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4644: {
4645:   Mat Arhs, Brhs;

4647:   PetscFunctionBegin;
4648:   PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4649:   /* undo the damage caused by shifting */
4650:   PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4651:   PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4652:   PetscCall(MatMult(Arhs, U, F));
4653:   PetscFunctionReturn(PETSC_SUCCESS);
4654: }

4656: /*@C
4657:   TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.

4659:   Collective

4661:   Input Parameters:
4662: + ts  - time stepping context
4663: . t   - time at which to evaluate
4664: . U   - state at which to evaluate
4665: - ctx - context

4667:   Output Parameters:
4668: + A - Jacobian
4669: - B - matrix used to construct the preconditioner, often the same as `A`

4671:   Level: intermediate

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

4676: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4677: @*/
4678: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4679: {
4680:   PetscFunctionBegin;
4681:   PetscFunctionReturn(PETSC_SUCCESS);
4682: }

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

4687:   Collective

4689:   Input Parameters:
4690: + ts   - time stepping context
4691: . t    - time at which to evaluate
4692: . U    - state at which to evaluate
4693: . Udot - time derivative of state vector
4694: - ctx  - context

4696:   Output Parameter:
4697: . F - left hand side

4699:   Level: intermediate

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

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

4709: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4710: @*/
4711: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4712: {
4713:   Mat A, B;

4715:   PetscFunctionBegin;
4716:   PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4717:   PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4718:   PetscCall(MatMult(A, Udot, F));
4719:   PetscFunctionReturn(PETSC_SUCCESS);
4720: }

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

4725:   Collective

4727:   Input Parameters:
4728: + ts    - time stepping context
4729: . t     - time at which to evaluate
4730: . U     - state at which to evaluate
4731: . Udot  - time derivative of state vector
4732: . shift - shift to apply
4733: - ctx   - context

4735:   Output Parameters:
4736: + A - pointer to operator
4737: - B - pointer to matrix from which the preconditioner is built (often `A`)

4739:   Level: advanced

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

4744:   It is only appropriate for problems of the form

4746:   $$
4747:   M \dot{U} = F(U,t)
4748:   $$

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

4754:   $$
4755:   shift*M + J
4756:   $$

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

4761: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4762: @*/
4763: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4764: {
4765:   PetscFunctionBegin;
4766:   PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4767:   ts->ijacobian.shift = shift;
4768:   PetscFunctionReturn(PETSC_SUCCESS);
4769: }

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

4774:   Not Collective

4776:   Input Parameter:
4777: . ts - the `TS` context

4779:   Output Parameter:
4780: . equation_type - see `TSEquationType`

4782:   Level: beginner

4784: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4785: @*/
4786: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4787: {
4788:   PetscFunctionBegin;
4790:   PetscAssertPointer(equation_type, 2);
4791:   *equation_type = ts->equation_type;
4792:   PetscFunctionReturn(PETSC_SUCCESS);
4793: }

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

4798:   Not Collective

4800:   Input Parameters:
4801: + ts            - the `TS` context
4802: - equation_type - see `TSEquationType`

4804:   Level: advanced

4806: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4807: @*/
4808: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4809: {
4810:   PetscFunctionBegin;
4812:   ts->equation_type = equation_type;
4813:   PetscFunctionReturn(PETSC_SUCCESS);
4814: }

4816: /*@
4817:   TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.

4819:   Not Collective

4821:   Input Parameter:
4822: . ts - the `TS` context

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

4828:   Level: beginner

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

4833: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4834: @*/
4835: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4836: {
4837:   PetscFunctionBegin;
4839:   PetscAssertPointer(reason, 2);
4840:   *reason = ts->reason;
4841:   PetscFunctionReturn(PETSC_SUCCESS);
4842: }

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

4847:   Logically Collective; reason must contain common value

4849:   Input Parameters:
4850: + ts     - the `TS` context
4851: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4852:             manual pages for the individual convergence tests for complete lists

4854:   Level: advanced

4856:   Note:
4857:   Can only be called while `TSSolve()` is active.

4859: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4860: @*/
4861: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4862: {
4863:   PetscFunctionBegin;
4865:   ts->reason = reason;
4866:   PetscFunctionReturn(PETSC_SUCCESS);
4867: }

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

4872:   Not Collective

4874:   Input Parameter:
4875: . ts - the `TS` context

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

4880:   Level: beginner

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

4885: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4886: @*/
4887: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4888: {
4889:   PetscFunctionBegin;
4891:   PetscAssertPointer(ftime, 2);
4892:   *ftime = ts->solvetime;
4893:   PetscFunctionReturn(PETSC_SUCCESS);
4894: }

4896: /*@
4897:   TSGetSNESIterations - Gets the total number of nonlinear iterations
4898:   used by the time integrator.

4900:   Not Collective

4902:   Input Parameter:
4903: . ts - `TS` context

4905:   Output Parameter:
4906: . nits - number of nonlinear iterations

4908:   Level: intermediate

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

4913: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4914: @*/
4915: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4916: {
4917:   PetscFunctionBegin;
4919:   PetscAssertPointer(nits, 2);
4920:   *nits = ts->snes_its;
4921:   PetscFunctionReturn(PETSC_SUCCESS);
4922: }

4924: /*@
4925:   TSGetKSPIterations - Gets the total number of linear iterations
4926:   used by the time integrator.

4928:   Not Collective

4930:   Input Parameter:
4931: . ts - `TS` context

4933:   Output Parameter:
4934: . lits - number of linear iterations

4936:   Level: intermediate

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

4941: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`
4942: @*/
4943: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4944: {
4945:   PetscFunctionBegin;
4947:   PetscAssertPointer(lits, 2);
4948:   *lits = ts->ksp_its;
4949:   PetscFunctionReturn(PETSC_SUCCESS);
4950: }

4952: /*@
4953:   TSGetStepRejections - Gets the total number of rejected steps.

4955:   Not Collective

4957:   Input Parameter:
4958: . ts - `TS` context

4960:   Output Parameter:
4961: . rejects - number of steps rejected

4963:   Level: intermediate

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

4968: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4969: @*/
4970: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4971: {
4972:   PetscFunctionBegin;
4974:   PetscAssertPointer(rejects, 2);
4975:   *rejects = ts->reject;
4976:   PetscFunctionReturn(PETSC_SUCCESS);
4977: }

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

4982:   Not Collective

4984:   Input Parameter:
4985: . ts - `TS` context

4987:   Output Parameter:
4988: . fails - number of failed nonlinear solves

4990:   Level: intermediate

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

4995: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4996: @*/
4997: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4998: {
4999:   PetscFunctionBegin;
5001:   PetscAssertPointer(fails, 2);
5002:   *fails = ts->num_snes_failures;
5003:   PetscFunctionReturn(PETSC_SUCCESS);
5004: }

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

5009:   Not Collective

5011:   Input Parameters:
5012: + ts      - `TS` context
5013: - rejects - maximum number of rejected steps, pass `PETSC_UNLIMITED` for unlimited

5015:   Options Database Key:
5016: . -ts_max_reject - Maximum number of step rejections before a step fails

5018:   Level: intermediate

5020:   Developer Note:
5021:   The options database name is incorrect.

5023: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
5024: @*/
5025: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
5026: {
5027:   PetscFunctionBegin;
5029:   if (rejects == PETSC_UNLIMITED || rejects == -1) {
5030:     ts->max_reject = PETSC_UNLIMITED;
5031:   } else {
5032:     PetscCheck(rejects >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of rejections");
5033:     ts->max_reject = rejects;
5034:   }
5035:   PetscFunctionReturn(PETSC_SUCCESS);
5036: }

5038: /*@
5039:   TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves

5041:   Not Collective

5043:   Input Parameters:
5044: + ts    - `TS` context
5045: - fails - maximum number of failed nonlinear solves, pass `PETSC_UNLIMITED` to allow any number of failures.

5047:   Options Database Key:
5048: . -ts_max_snes_failures - Maximum number of nonlinear solve failures

5050:   Level: intermediate

5052: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
5053: @*/
5054: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
5055: {
5056:   PetscFunctionBegin;
5058:   if (fails == PETSC_UNLIMITED || fails == -1) {
5059:     ts->max_snes_failures = PETSC_UNLIMITED;
5060:   } else {
5061:     PetscCheck(fails >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures");
5062:     ts->max_snes_failures = fails;
5063:   }
5064:   PetscFunctionReturn(PETSC_SUCCESS);
5065: }

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

5070:   Not Collective

5072:   Input Parameters:
5073: + ts  - `TS` context
5074: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure

5076:   Options Database Key:
5077: . -ts_error_if_step_fails - Error if no step succeeds

5079:   Level: intermediate

5081: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
5082: @*/
5083: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
5084: {
5085:   PetscFunctionBegin;
5087:   ts->errorifstepfailed = err;
5088:   PetscFunctionReturn(PETSC_SUCCESS);
5089: }

5091: /*@
5092:   TSGetAdapt - Get the adaptive controller context for the current method

5094:   Collective if controller has not yet been created

5096:   Input Parameter:
5097: . ts - time stepping context

5099:   Output Parameter:
5100: . adapt - adaptive controller

5102:   Level: intermediate

5104: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
5105: @*/
5106: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
5107: {
5108:   PetscFunctionBegin;
5110:   PetscAssertPointer(adapt, 2);
5111:   if (!ts->adapt) {
5112:     PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
5113:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
5114:   }
5115:   *adapt = ts->adapt;
5116:   PetscFunctionReturn(PETSC_SUCCESS);
5117: }

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

5122:   Logically Collective

5124:   Input Parameters:
5125: + ts    - time integration context
5126: . atol  - scalar absolute tolerances
5127: . vatol - vector of absolute tolerances or `NULL`, used in preference to `atol` if present
5128: . rtol  - scalar relative tolerances
5129: - vrtol - vector of relative tolerances or `NULL`, used in preference to `rtol` if present

5131:   Options Database Keys:
5132: + -ts_rtol <rtol> - relative tolerance for local truncation error
5133: - -ts_atol <atol> - Absolute tolerance for local truncation error

5135:   Level: beginner

5137:   Notes:
5138:   `PETSC_CURRENT` or `PETSC_DETERMINE` may be used for `atol` or `rtol` to indicate the current value
5139:   or the default value from when the object's type was set.

5141:   With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5142:   (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5143:   computed only for the differential or the algebraic part then this can be done using the vector of
5144:   tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
5145:   differential part and infinity for the algebraic part, the LTE calculation will include only the
5146:   differential variables.

5148:   Fortran Note:
5149:   Use `PETSC_CURRENT_INTEGER` or `PETSC_DETERMINE_INTEGER`.

5151: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
5152: @*/
5153: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
5154: {
5155:   PetscFunctionBegin;
5156:   if (atol == (PetscReal)PETSC_DETERMINE) {
5157:     ts->atol = ts->default_atol;
5158:   } else if (atol != (PetscReal)PETSC_CURRENT) {
5159:     PetscCheck(atol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)atol);
5160:     ts->atol = atol;
5161:   }

5163:   if (vatol) {
5164:     PetscCall(PetscObjectReference((PetscObject)vatol));
5165:     PetscCall(VecDestroy(&ts->vatol));
5166:     ts->vatol = vatol;
5167:   }

5169:   if (rtol == (PetscReal)PETSC_DETERMINE) {
5170:     ts->rtol = ts->default_rtol;
5171:   } else if (rtol != (PetscReal)PETSC_CURRENT) {
5172:     PetscCheck(rtol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative", (double)rtol);
5173:     ts->rtol = rtol;
5174:   }

5176:   if (vrtol) {
5177:     PetscCall(PetscObjectReference((PetscObject)vrtol));
5178:     PetscCall(VecDestroy(&ts->vrtol));
5179:     ts->vrtol = vrtol;
5180:   }
5181:   PetscFunctionReturn(PETSC_SUCCESS);
5182: }

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

5187:   Logically Collective

5189:   Input Parameter:
5190: . ts - time integration context

5192:   Output Parameters:
5193: + atol  - scalar absolute tolerances, `NULL` to ignore
5194: . vatol - vector of absolute tolerances, `NULL` to ignore
5195: . rtol  - scalar relative tolerances, `NULL` to ignore
5196: - vrtol - vector of relative tolerances, `NULL` to ignore

5198:   Level: beginner

5200: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
5201: @*/
5202: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
5203: {
5204:   PetscFunctionBegin;
5205:   if (atol) *atol = ts->atol;
5206:   if (vatol) *vatol = ts->vatol;
5207:   if (rtol) *rtol = ts->rtol;
5208:   if (vrtol) *vrtol = ts->vrtol;
5209:   PetscFunctionReturn(PETSC_SUCCESS);
5210: }

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

5215:   Collective

5217:   Input Parameters:
5218: + ts        - time stepping context
5219: . U         - state vector, usually ts->vec_sol
5220: . Y         - state vector to be compared to U
5221: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

5223:   Output Parameters:
5224: + norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5225: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5226: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

5228:   Options Database Key:
5229: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5231:   Level: developer

5233: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5234: @*/
5235: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5236: {
5237:   PetscInt norma_loc, norm_loc, normr_loc;

5239:   PetscFunctionBegin;
5241:   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));
5242:   if (wnormtype == NORM_2) {
5243:     if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5244:     if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5245:     if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5246:   }
5247:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5248:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5249:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5250:   PetscFunctionReturn(PETSC_SUCCESS);
5251: }

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

5256:   Collective

5258:   Input Parameters:
5259: + ts        - time stepping context
5260: . E         - error vector
5261: . U         - state vector, usually ts->vec_sol
5262: . Y         - state vector, previous time step
5263: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

5265:   Output Parameters:
5266: + norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5267: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5268: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

5270:   Options Database Key:
5271: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

5273:   Level: developer

5275: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5276: @*/
5277: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5278: {
5279:   PetscInt norma_loc, norm_loc, normr_loc;

5281:   PetscFunctionBegin;
5283:   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));
5284:   if (wnormtype == NORM_2) {
5285:     if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5286:     if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5287:     if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5288:   }
5289:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5290:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5291:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5292:   PetscFunctionReturn(PETSC_SUCCESS);
5293: }

5295: /*@
5296:   TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler

5298:   Logically Collective

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

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

5307:   Level: intermediate

5309: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5310: @*/
5311: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5312: {
5313:   PetscFunctionBegin;
5315:   ts->cfltime_local = cfltime;
5316:   ts->cfltime       = -1.;
5317:   PetscFunctionReturn(PETSC_SUCCESS);
5318: }

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

5323:   Collective

5325:   Input Parameter:
5326: . ts - time stepping context

5328:   Output Parameter:
5329: . cfltime - maximum stable time step for forward Euler

5331:   Level: advanced

5333: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5334: @*/
5335: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5336: {
5337:   PetscFunctionBegin;
5338:   if (ts->cfltime < 0) PetscCallMPI(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5339:   *cfltime = ts->cfltime;
5340:   PetscFunctionReturn(PETSC_SUCCESS);
5341: }

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

5346:   Input Parameters:
5347: + ts - the `TS` context.
5348: . xl - lower bound.
5349: - xu - upper bound.

5351:   Level: advanced

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

5357: .seealso: [](ch_ts), `TS`
5358: @*/
5359: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5360: {
5361:   SNES snes;

5363:   PetscFunctionBegin;
5364:   PetscCall(TSGetSNES(ts, &snes));
5365:   PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5366:   PetscFunctionReturn(PETSC_SUCCESS);
5367: }

5369: /*@
5370:   TSComputeLinearStability - computes the linear stability function at a point

5372:   Collective

5374:   Input Parameters:
5375: + ts - the `TS` context
5376: . xr - real part of input argument
5377: - xi - imaginary part of input argument

5379:   Output Parameters:
5380: + yr - real part of function value
5381: - yi - imaginary part of function value

5383:   Level: developer

5385: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5386: @*/
5387: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5388: {
5389:   PetscFunctionBegin;
5391:   PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5392:   PetscFunctionReturn(PETSC_SUCCESS);
5393: }

5395: /*@
5396:   TSRestartStep - Flags the solver to restart the next step

5398:   Collective

5400:   Input Parameter:
5401: . ts - the `TS` context obtained from `TSCreate()`

5403:   Level: advanced

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

5413: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5414: @*/
5415: PetscErrorCode TSRestartStep(TS ts)
5416: {
5417:   PetscFunctionBegin;
5419:   ts->steprestart = PETSC_TRUE;
5420:   PetscFunctionReturn(PETSC_SUCCESS);
5421: }

5423: /*@
5424:   TSRollBack - Rolls back one time step

5426:   Collective

5428:   Input Parameter:
5429: . ts - the `TS` context obtained from `TSCreate()`

5431:   Level: advanced

5433: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5434: @*/
5435: PetscErrorCode TSRollBack(TS ts)
5436: {
5437:   PetscFunctionBegin;
5439:   PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5440:   PetscTryTypeMethod(ts, rollback);
5441:   PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
5442:   ts->time_step  = ts->ptime - ts->ptime_prev;
5443:   ts->ptime      = ts->ptime_prev;
5444:   ts->ptime_prev = ts->ptime_prev_rollback;
5445:   ts->steps--;
5446:   ts->steprollback = PETSC_TRUE;
5447:   PetscFunctionReturn(PETSC_SUCCESS);
5448: }

5450: /*@
5451:   TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step

5453:   Not collective

5455:   Input Parameter:
5456: . ts - the `TS` context obtained from `TSCreate()`

5458:   Output Parameter:
5459: . flg - the rollback flag

5461:   Level: advanced

5463: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
5464: @*/
5465: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
5466: {
5467:   PetscFunctionBegin;
5469:   PetscAssertPointer(flg, 2);
5470:   *flg = ts->steprollback;
5471:   PetscFunctionReturn(PETSC_SUCCESS);
5472: }

5474: /*@
5475:   TSGetStepResize - Get the internal flag indicating if the current step is after a resize.

5477:   Not collective

5479:   Input Parameter:
5480: . ts - the `TS` context obtained from `TSCreate()`

5482:   Output Parameter:
5483: . flg - the resize flag

5485:   Level: advanced

5487: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetResize()`
5488: @*/
5489: PetscErrorCode TSGetStepResize(TS ts, PetscBool *flg)
5490: {
5491:   PetscFunctionBegin;
5493:   PetscAssertPointer(flg, 2);
5494:   *flg = ts->stepresize;
5495:   PetscFunctionReturn(PETSC_SUCCESS);
5496: }

5498: /*@
5499:   TSGetStages - Get the number of stages and stage values

5501:   Input Parameter:
5502: . ts - the `TS` context obtained from `TSCreate()`

5504:   Output Parameters:
5505: + ns - the number of stages
5506: - Y  - the current stage vectors

5508:   Level: advanced

5510:   Note:
5511:   Both `ns` and `Y` can be `NULL`.

5513: .seealso: [](ch_ts), `TS`, `TSCreate()`
5514: @*/
5515: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5516: {
5517:   PetscFunctionBegin;
5519:   if (ns) PetscAssertPointer(ns, 2);
5520:   if (Y) PetscAssertPointer(Y, 3);
5521:   if (!ts->ops->getstages) {
5522:     if (ns) *ns = 0;
5523:     if (Y) *Y = NULL;
5524:   } else PetscUseTypeMethod(ts, getstages, ns, Y);
5525:   PetscFunctionReturn(PETSC_SUCCESS);
5526: }

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

5531:   Collective

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

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

5545:   Level: intermediate

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

5550:   dF/dU + shift*dF/dUdot

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

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

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

5567:   PetscFunctionBegin;
5568:   PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5569:   PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5570:   if (!color) {
5571:     DM         dm;
5572:     ISColoring iscoloring;

5574:     PetscCall(TSGetDM(ts, &dm));
5575:     PetscCall(DMHasColoring(dm, &hascolor));
5576:     if (hascolor && !matcolor) {
5577:       PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5578:       PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5579:       PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5580:       PetscCall(MatFDColoringSetFromOptions(color));
5581:       PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5582:       PetscCall(ISColoringDestroy(&iscoloring));
5583:     } else {
5584:       MatColoring mc;

5586:       PetscCall(MatColoringCreate(B, &mc));
5587:       PetscCall(MatColoringSetDistance(mc, 2));
5588:       PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5589:       PetscCall(MatColoringSetFromOptions(mc));
5590:       PetscCall(MatColoringApply(mc, &iscoloring));
5591:       PetscCall(MatColoringDestroy(&mc));
5592:       PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5593:       PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5594:       PetscCall(MatFDColoringSetFromOptions(color));
5595:       PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5596:       PetscCall(ISColoringDestroy(&iscoloring));
5597:     }
5598:     PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5599:     PetscCall(PetscObjectDereference((PetscObject)color));
5600:   }
5601:   PetscCall(TSGetSNES(ts, &snes));
5602:   PetscCall(MatFDColoringApply(B, color, U, snes));
5603:   if (J != B) {
5604:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5605:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5606:   }
5607:   PetscFunctionReturn(PETSC_SUCCESS);
5608: }

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

5613:   Input Parameters:
5614: + ts   - the `TS` context
5615: - func - function called within `TSFunctionDomainError()`

5617:   Calling sequence of `func`:
5618: + ts     - the `TS` context
5619: . time   - the current time (of the stage)
5620: . state  - the state to check if it is valid
5621: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable

5623:   Level: intermediate

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

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

5635: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5636: @*/
5637: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5638: {
5639:   PetscFunctionBegin;
5641:   ts->functiondomainerror = func;
5642:   PetscFunctionReturn(PETSC_SUCCESS);
5643: }

5645: /*@
5646:   TSFunctionDomainError - Checks if the current state is valid

5648:   Input Parameters:
5649: + ts        - the `TS` context
5650: . stagetime - time of the simulation
5651: - Y         - state vector to check.

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

5656:   Level: developer

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

5662: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5663: @*/
5664: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5665: {
5666:   PetscFunctionBegin;
5668:   *accept = PETSC_TRUE;
5669:   if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5670:   PetscFunctionReturn(PETSC_SUCCESS);
5671: }

5673: /*@
5674:   TSClone - This function clones a time step `TS` object.

5676:   Collective

5678:   Input Parameter:
5679: . tsin - The input `TS`

5681:   Output Parameter:
5682: . tsout - The output `TS` (cloned)

5684:   Level: developer

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

5690:   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
5691: .vb
5692:  SNES snes_dup = NULL;
5693:  TSGetSNES(ts,&snes_dup);
5694:  TSSetSNES(ts,snes_dup);
5695: .ve

5697: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5698: @*/
5699: PetscErrorCode TSClone(TS tsin, TS *tsout)
5700: {
5701:   TS     t;
5702:   SNES   snes_start;
5703:   DM     dm;
5704:   TSType type;

5706:   PetscFunctionBegin;
5707:   PetscAssertPointer(tsin, 1);
5708:   *tsout = NULL;

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

5712:   /* General TS description */
5713:   t->numbermonitors    = 0;
5714:   t->setupcalled       = PETSC_FALSE;
5715:   t->ksp_its           = 0;
5716:   t->snes_its          = 0;
5717:   t->nwork             = 0;
5718:   t->rhsjacobian.time  = PETSC_MIN_REAL;
5719:   t->rhsjacobian.scale = 1.;
5720:   t->ijacobian.shift   = 1.;

5722:   PetscCall(TSGetSNES(tsin, &snes_start));
5723:   PetscCall(TSSetSNES(t, snes_start));

5725:   PetscCall(TSGetDM(tsin, &dm));
5726:   PetscCall(TSSetDM(t, dm));

5728:   t->adapt = tsin->adapt;
5729:   PetscCall(PetscObjectReference((PetscObject)t->adapt));

5731:   t->trajectory = tsin->trajectory;
5732:   PetscCall(PetscObjectReference((PetscObject)t->trajectory));

5734:   t->event = tsin->event;
5735:   if (t->event) t->event->refct++;

5737:   t->problem_type      = tsin->problem_type;
5738:   t->ptime             = tsin->ptime;
5739:   t->ptime_prev        = tsin->ptime_prev;
5740:   t->time_step         = tsin->time_step;
5741:   t->max_time          = tsin->max_time;
5742:   t->steps             = tsin->steps;
5743:   t->max_steps         = tsin->max_steps;
5744:   t->equation_type     = tsin->equation_type;
5745:   t->atol              = tsin->atol;
5746:   t->rtol              = tsin->rtol;
5747:   t->max_snes_failures = tsin->max_snes_failures;
5748:   t->max_reject        = tsin->max_reject;
5749:   t->errorifstepfailed = tsin->errorifstepfailed;

5751:   PetscCall(TSGetType(tsin, &type));
5752:   PetscCall(TSSetType(t, type));

5754:   t->vec_sol = NULL;

5756:   t->cfltime          = tsin->cfltime;
5757:   t->cfltime_local    = tsin->cfltime_local;
5758:   t->exact_final_time = tsin->exact_final_time;

5760:   t->ops[0] = tsin->ops[0];

5762:   if (((PetscObject)tsin)->fortran_func_pointers) {
5763:     PetscInt i;
5764:     PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers));
5765:     for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5766:   }
5767:   *tsout = t;
5768:   PetscFunctionReturn(PETSC_SUCCESS);
5769: }

5771: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5772: {
5773:   TS ts = (TS)ctx;

5775:   PetscFunctionBegin;
5776:   PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5777:   PetscFunctionReturn(PETSC_SUCCESS);
5778: }

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

5783:   Logically Collective

5785:   Input Parameter:
5786: . ts - the time stepping routine

5788:   Output Parameter:
5789: . flg - `PETSC_TRUE` if the multiply is likely correct

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

5794:   Level: advanced

5796:   Note:
5797:   This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian

5799: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5800: @*/
5801: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5802: {
5803:   Mat              J, B;
5804:   TSRHSJacobianFn *func;
5805:   void            *ctx;

5807:   PetscFunctionBegin;
5808:   PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5809:   PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5810:   PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5811:   PetscFunctionReturn(PETSC_SUCCESS);
5812: }

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

5817:   Logically Collective

5819:   Input Parameter:
5820: . ts - the time stepping routine

5822:   Output Parameter:
5823: . flg - `PETSC_TRUE` if the multiply is likely correct

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

5828:   Level: advanced

5830:   Notes:
5831:   This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian

5833: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5834: @*/
5835: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5836: {
5837:   Mat              J, B;
5838:   void            *ctx;
5839:   TSRHSJacobianFn *func;

5841:   PetscFunctionBegin;
5842:   PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5843:   PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5844:   PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5845:   PetscFunctionReturn(PETSC_SUCCESS);
5846: }

5848: /*@
5849:   TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.

5851:   Logically Collective

5853:   Input Parameters:
5854: + ts                   - timestepping context
5855: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

5857:   Options Database Key:
5858: . -ts_use_splitrhsfunction - <true,false>

5860:   Level: intermediate

5862:   Note:
5863:   This is only for multirate methods

5865: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5866: @*/
5867: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5868: {
5869:   PetscFunctionBegin;
5871:   ts->use_splitrhsfunction = use_splitrhsfunction;
5872:   PetscFunctionReturn(PETSC_SUCCESS);
5873: }

5875: /*@
5876:   TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.

5878:   Not Collective

5880:   Input Parameter:
5881: . ts - timestepping context

5883:   Output Parameter:
5884: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

5886:   Level: intermediate

5888: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5889: @*/
5890: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5891: {
5892:   PetscFunctionBegin;
5894:   *use_splitrhsfunction = ts->use_splitrhsfunction;
5895:   PetscFunctionReturn(PETSC_SUCCESS);
5896: }

5898: /*@
5899:   TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.

5901:   Logically  Collective

5903:   Input Parameters:
5904: + ts  - the time-stepper
5905: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)

5907:   Level: intermediate

5909:   Note:
5910:   When the relationship between the nonzero structures is known and supplied the solution process can be much faster

5912: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
5913:  @*/
5914: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5915: {
5916:   PetscFunctionBegin;
5918:   ts->axpy_pattern = str;
5919:   PetscFunctionReturn(PETSC_SUCCESS);
5920: }

5922: /*@
5923:   TSSetEvaluationTimes - sets the evaluation points. The solution will be computed and stored for each time requested

5925:   Collective

5927:   Input Parameters:
5928: + ts          - the time-stepper
5929: . n           - number of the time points
5930: - time_points - array of the time points, must be increasing

5932:   Options Database Key:
5933: . -ts_eval_times <t0,...tn> - Sets the evaluation times

5935:   Level: intermediate

5937:   Notes:
5938:   The elements in `time_points` must be all increasing. They correspond to the intermediate points to be saved.

5940:   `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.

5942:   The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using evaluation times may
5943:   pressure the memory system when using a large number of time points.

5945: .seealso: [](ch_ts), `TS`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`, `TSSetTimeSpan()`
5946:  @*/
5947: PetscErrorCode TSSetEvaluationTimes(TS ts, PetscInt n, PetscReal *time_points)
5948: {
5949:   PetscBool is_sorted;

5951:   PetscFunctionBegin;
5953:   if (ts->eval_times) { // Reset eval_times
5954:     ts->eval_times->sol_idx        = 0;
5955:     ts->eval_times->time_point_idx = 0;
5956:     if (n != ts->eval_times->num_time_points) {
5957:       PetscCall(PetscFree(ts->eval_times->time_points));
5958:       PetscCall(PetscFree(ts->eval_times->sol_times));
5959:       PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
5960:     } else {
5961:       PetscCall(PetscArrayzero(ts->eval_times->sol_times, n));
5962:       for (PetscInt i = 0; i < n; i++) PetscCall(VecZeroEntries(ts->eval_times->sol_vecs[i]));
5963:     }
5964:   } else { // Create/initialize eval_times
5965:     TSEvaluationTimes eval_times;
5966:     PetscCall(PetscNew(&eval_times));
5967:     PetscCall(PetscMalloc1(n, &eval_times->time_points));
5968:     PetscCall(PetscMalloc1(n, &eval_times->sol_times));
5969:     eval_times->reltol  = 1e-6;
5970:     eval_times->abstol  = 10 * PETSC_MACHINE_EPSILON;
5971:     eval_times->worktol = 0;
5972:     ts->eval_times      = eval_times;
5973:   }
5974:   ts->eval_times->num_time_points = n;
5975:   PetscCall(PetscSortedReal(n, time_points, &is_sorted));
5976:   PetscCheck(is_sorted, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "time_points array must be sorted");
5977:   PetscCall(PetscArraycpy(ts->eval_times->time_points, time_points, n));
5978:   // Note: ts->vec_sol not guaranteed to exist, so ts->eval_times->sol_vecs allocated at TSSolve time
5979:   PetscFunctionReturn(PETSC_SUCCESS);
5980: }

5982: /*@C
5983:   TSGetEvaluationTimes - gets the evaluation times set with `TSSetEvaluationTimes()`

5985:   Not Collective

5987:   Input Parameter:
5988: . ts - the time-stepper

5990:   Output Parameters:
5991: + n           - number of the time points
5992: - time_points - array of the time points

5994:   Level: beginner

5996:   Note:
5997:   The values obtained are valid until the `TS` object is destroyed.

5999:   Both `n` and `time_points` can be `NULL`.

6001:   Also used to see time points set by `TSSetTimeSpan()`.

6003: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6004:  @*/
6005: PetscErrorCode TSGetEvaluationTimes(TS ts, PetscInt *n, const PetscReal *time_points[])
6006: {
6007:   PetscFunctionBegin;
6009:   if (n) PetscAssertPointer(n, 2);
6010:   if (time_points) PetscAssertPointer(time_points, 3);
6011:   if (!ts->eval_times) {
6012:     if (n) *n = 0;
6013:     if (time_points) *time_points = NULL;
6014:   } else {
6015:     if (n) *n = ts->eval_times->num_time_points;
6016:     if (time_points) *time_points = ts->eval_times->time_points;
6017:   }
6018:   PetscFunctionReturn(PETSC_SUCCESS);
6019: }

6021: /*@C
6022:   TSGetEvaluationSolutions - Get the number of solutions and the solutions at the evaluation time points specified

6024:   Input Parameter:
6025: . ts - the `TS` context obtained from `TSCreate()`

6027:   Output Parameters:
6028: + nsol      - the number of solutions
6029: . sol_times - array of solution times corresponding to the solution vectors. See note below
6030: - Sols      - the solution vectors

6032:   Level: intermediate

6034:   Notes:
6035:   Both `nsol` and `Sols` can be `NULL`.

6037:   Some time points in the evaluation points may be skipped by `TS` so that `nsol` is less than the number of points specified by `TSSetEvaluationTimes()`.
6038:   For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain evaluation times.

6040:   Also used to see view solutions requested by `TSSetTimeSpan()`.

6042: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`
6043: @*/
6044: PetscErrorCode TSGetEvaluationSolutions(TS ts, PetscInt *nsol, const PetscReal *sol_times[], Vec **Sols)
6045: {
6046:   PetscFunctionBegin;
6048:   if (nsol) PetscAssertPointer(nsol, 2);
6049:   if (sol_times) PetscAssertPointer(sol_times, 3);
6050:   if (Sols) PetscAssertPointer(Sols, 4);
6051:   if (!ts->eval_times) {
6052:     if (nsol) *nsol = 0;
6053:     if (sol_times) *sol_times = NULL;
6054:     if (Sols) *Sols = NULL;
6055:   } else {
6056:     if (nsol) *nsol = ts->eval_times->sol_idx;
6057:     if (sol_times) *sol_times = ts->eval_times->sol_times;
6058:     if (Sols) *Sols = ts->eval_times->sol_vecs;
6059:   }
6060:   PetscFunctionReturn(PETSC_SUCCESS);
6061: }

6063: /*@
6064:   TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span

6066:   Collective

6068:   Input Parameters:
6069: + ts         - the time-stepper
6070: . n          - number of the time points (>=2)
6071: - span_times - array of the time points, must be increasing. The first element and the last element are the initial time and the final time respectively.

6073:   Options Database Key:
6074: . -ts_time_span <t0,...tf> - Sets the time span

6076:   Level: intermediate

6078:   Notes:
6079:   This function is identical to `TSSetEvaluationTimes()`, except that it also sets the initial time and final time for the `ts` to the first and last `span_times` entries.

6081:   The elements in `span_times` must be all increasing. They correspond to the intermediate points to be saved.

6083:   `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.

6085:   The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using time span may
6086:   pressure the memory system when using a large number of span points.

6088: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6089:  @*/
6090: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
6091: {
6092:   PetscFunctionBegin;
6094:   PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
6095:   PetscCall(TSSetEvaluationTimes(ts, n, span_times));
6096:   PetscCall(TSSetTime(ts, span_times[0]));
6097:   PetscCall(TSSetMaxTime(ts, span_times[n - 1]));
6098:   PetscFunctionReturn(PETSC_SUCCESS);
6099: }

6101: /*@
6102:   TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.

6104:   Collective

6106:   Input Parameters:
6107: + ts - the `TS` context
6108: . J  - Jacobian matrix (not altered in this routine)
6109: - B  - newly computed Jacobian matrix to use with preconditioner

6111:   Level: intermediate

6113:   Notes:
6114:   This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
6115:   many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
6116:   and multiple fields are involved.

6118:   Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
6119:   structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
6120:   usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
6121:   `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.

6123: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
6124: @*/
6125: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
6126: {
6127:   MatColoring   mc            = NULL;
6128:   ISColoring    iscoloring    = NULL;
6129:   MatFDColoring matfdcoloring = NULL;

6131:   PetscFunctionBegin;
6132:   /* Generate new coloring after eliminating zeros in the matrix */
6133:   PetscCall(MatEliminateZeros(B, PETSC_TRUE));
6134:   PetscCall(MatColoringCreate(B, &mc));
6135:   PetscCall(MatColoringSetDistance(mc, 2));
6136:   PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
6137:   PetscCall(MatColoringSetFromOptions(mc));
6138:   PetscCall(MatColoringApply(mc, &iscoloring));
6139:   PetscCall(MatColoringDestroy(&mc));
6140:   /* Replace the old coloring with the new one */
6141:   PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
6142:   PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
6143:   PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
6144:   PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
6145:   PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
6146:   PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
6147:   PetscCall(ISColoringDestroy(&iscoloring));
6148:   PetscFunctionReturn(PETSC_SUCCESS);
6149: }