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] = {0};
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: if (ts->usessnes && ts->snes) PetscCall(SNESMonitorCancel(ts->snes));
3534: 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]);
3535: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3536: }
3537: PetscFunctionReturn(PETSC_SUCCESS);
3538: }
3540: /*@
3541: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3542: at the end of a time step with a given order of accuracy.
3544: Collective
3546: Input Parameters:
3547: + ts - time stepping context
3548: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3550: Input/Output Parameter:
3551: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3552: on output, the actual order of the error evaluation
3554: Output Parameter:
3555: . wlte - the weighted local truncation error norm
3557: Level: advanced
3559: Note:
3560: If the timestepper cannot evaluate the error in a particular step
3561: (eg. in the first step or restart steps after event handling),
3562: this routine returns wlte=-1.0 .
3564: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3565: @*/
3566: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3567: {
3568: PetscFunctionBegin;
3572: if (order) PetscAssertPointer(order, 3);
3574: PetscAssertPointer(wlte, 4);
3575: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3576: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3577: PetscFunctionReturn(PETSC_SUCCESS);
3578: }
3580: /*@
3581: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3583: Collective
3585: Input Parameters:
3586: + ts - time stepping context
3587: . order - desired order of accuracy
3588: - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)
3590: Output Parameter:
3591: . U - state at the end of the current step
3593: Level: advanced
3595: Notes:
3596: This function cannot be called until all stages have been evaluated.
3598: 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.
3600: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3601: @*/
3602: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3603: {
3604: PetscFunctionBegin;
3608: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3609: PetscFunctionReturn(PETSC_SUCCESS);
3610: }
3612: /*@C
3613: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3615: Not collective
3617: Input Parameter:
3618: . ts - time stepping context
3620: Output Parameter:
3621: . initCondition - The function which computes an initial condition
3623: Calling sequence of `initCondition`:
3624: + ts - The timestepping context
3625: - u - The input vector in which the initial condition is stored
3627: Level: advanced
3629: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3630: @*/
3631: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3632: {
3633: PetscFunctionBegin;
3635: PetscAssertPointer(initCondition, 2);
3636: *initCondition = ts->ops->initcondition;
3637: PetscFunctionReturn(PETSC_SUCCESS);
3638: }
3640: /*@C
3641: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3643: Logically collective
3645: Input Parameters:
3646: + ts - time stepping context
3647: - initCondition - The function which computes an initial condition
3649: Calling sequence of `initCondition`:
3650: + ts - The timestepping context
3651: - e - The input vector in which the initial condition is to be stored
3653: Level: advanced
3655: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3656: @*/
3657: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3658: {
3659: PetscFunctionBegin;
3662: ts->ops->initcondition = initCondition;
3663: PetscFunctionReturn(PETSC_SUCCESS);
3664: }
3666: /*@
3667: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3669: Collective
3671: Input Parameters:
3672: + ts - time stepping context
3673: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3675: Level: advanced
3677: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3678: @*/
3679: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3680: {
3681: PetscFunctionBegin;
3684: PetscTryTypeMethod(ts, initcondition, u);
3685: PetscFunctionReturn(PETSC_SUCCESS);
3686: }
3688: /*@C
3689: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3691: Not collective
3693: Input Parameter:
3694: . ts - time stepping context
3696: Output Parameter:
3697: . exactError - The function which computes the solution error
3699: Calling sequence of `exactError`:
3700: + ts - The timestepping context
3701: . u - The approximate solution vector
3702: - e - The vector in which the error is stored
3704: Level: advanced
3706: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3707: @*/
3708: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3709: {
3710: PetscFunctionBegin;
3712: PetscAssertPointer(exactError, 2);
3713: *exactError = ts->ops->exacterror;
3714: PetscFunctionReturn(PETSC_SUCCESS);
3715: }
3717: /*@C
3718: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3720: Logically collective
3722: Input Parameters:
3723: + ts - time stepping context
3724: - exactError - The function which computes the solution error
3726: Calling sequence of `exactError`:
3727: + ts - The timestepping context
3728: . u - The approximate solution vector
3729: - e - The vector in which the error is stored
3731: Level: advanced
3733: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3734: @*/
3735: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3736: {
3737: PetscFunctionBegin;
3740: ts->ops->exacterror = exactError;
3741: PetscFunctionReturn(PETSC_SUCCESS);
3742: }
3744: /*@
3745: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3747: Collective
3749: Input Parameters:
3750: + ts - time stepping context
3751: . u - The approximate solution
3752: - e - The `Vec` used to store the error
3754: Level: advanced
3756: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3757: @*/
3758: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3759: {
3760: PetscFunctionBegin;
3764: PetscTryTypeMethod(ts, exacterror, u, e);
3765: PetscFunctionReturn(PETSC_SUCCESS);
3766: }
3768: /*@C
3769: TSSetResize - Sets the resize callbacks.
3771: Logically Collective
3773: Input Parameters:
3774: + ts - The `TS` context obtained from `TSCreate()`
3775: . rollback - Whether a resize will restart the step
3776: . setup - The setup function
3777: . transfer - The transfer function
3778: - ctx - [optional] The user-defined context
3780: Calling sequence of `setup`:
3781: + ts - the `TS` context
3782: . step - the current step
3783: . time - the current time
3784: . state - the current vector of state
3785: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3786: - ctx - user defined context
3788: Calling sequence of `transfer`:
3789: + ts - the `TS` context
3790: . nv - the number of vectors to be transferred
3791: . vecsin - array of vectors to be transferred
3792: . vecsout - array of transferred vectors
3793: - ctx - user defined context
3795: Notes:
3796: The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
3797: depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
3798: and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
3799: that the size of the discrete problem has changed.
3800: In both cases, the solver will collect the needed vectors that will be
3801: transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
3802: current solution vector, and other vectors needed by the specific solver used.
3803: For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3804: Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3805: will be automatically reset if the sizes are changed and they must be specified again by the user
3806: inside the `transfer` function.
3807: The input and output arrays passed to `transfer` are allocated by PETSc.
3808: Vectors in `vecsout` must be created by the user.
3809: Ownership of vectors in `vecsout` is transferred to PETSc.
3811: Level: advanced
3813: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3814: @*/
3815: 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)
3816: {
3817: PetscFunctionBegin;
3819: ts->resizerollback = rollback;
3820: ts->resizesetup = setup;
3821: ts->resizetransfer = transfer;
3822: ts->resizectx = ctx;
3823: PetscFunctionReturn(PETSC_SUCCESS);
3824: }
3826: /*
3827: TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.
3829: Collective
3831: Input Parameters:
3832: + ts - The `TS` context obtained from `TSCreate()`
3833: - 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.
3835: Level: developer
3837: Note:
3838: `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3839: used within time stepping implementations,
3840: so most users would not generally call this routine themselves.
3842: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3843: @*/
3844: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3845: {
3846: PetscFunctionBegin;
3848: PetscTryTypeMethod(ts, resizeregister, flg);
3849: /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3850: PetscFunctionReturn(PETSC_SUCCESS);
3851: }
3853: static PetscErrorCode TSResizeReset(TS ts)
3854: {
3855: PetscFunctionBegin;
3857: PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3858: PetscFunctionReturn(PETSC_SUCCESS);
3859: }
3861: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3862: {
3863: PetscFunctionBegin;
3866: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3867: if (ts->resizetransfer) {
3868: PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3869: PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3870: }
3871: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3872: PetscFunctionReturn(PETSC_SUCCESS);
3873: }
3875: /*@C
3876: TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.
3878: Collective
3880: Input Parameters:
3881: + ts - The `TS` context obtained from `TSCreate()`
3882: . name - A string identifying the vector
3883: - vec - The vector
3885: Level: developer
3887: Note:
3888: `TSResizeRegisterVec()` is typically used within time stepping implementations,
3889: so most users would not generally call this routine themselves.
3891: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3892: @*/
3893: PetscErrorCode TSResizeRegisterVec(TS ts, const char name[], Vec vec)
3894: {
3895: PetscFunctionBegin;
3897: PetscAssertPointer(name, 2);
3899: PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3900: PetscFunctionReturn(PETSC_SUCCESS);
3901: }
3903: /*@C
3904: TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.
3906: Collective
3908: Input Parameters:
3909: + ts - The `TS` context obtained from `TSCreate()`
3910: . name - A string identifying the vector
3911: - vec - The vector
3913: Level: developer
3915: Note:
3916: `TSResizeRetrieveVec()` is typically used within time stepping implementations,
3917: so most users would not generally call this routine themselves.
3919: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
3920: @*/
3921: PetscErrorCode TSResizeRetrieveVec(TS ts, const char name[], Vec *vec)
3922: {
3923: PetscFunctionBegin;
3925: PetscAssertPointer(name, 2);
3926: PetscAssertPointer(vec, 3);
3927: PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
3928: PetscFunctionReturn(PETSC_SUCCESS);
3929: }
3931: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
3932: {
3933: PetscInt cnt;
3934: PetscObjectList tmp;
3935: Vec *vecsin = NULL;
3936: const char **namesin = NULL;
3938: PetscFunctionBegin;
3939: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
3940: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
3941: if (names) PetscCall(PetscMalloc1(cnt, &namesin));
3942: if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
3943: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
3944: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
3945: if (vecs) vecsin[cnt] = (Vec)tmp->obj;
3946: if (names) namesin[cnt] = tmp->name;
3947: cnt++;
3948: }
3949: }
3950: if (nv) *nv = cnt;
3951: if (names) *names = namesin;
3952: if (vecs) *vecs = vecsin;
3953: PetscFunctionReturn(PETSC_SUCCESS);
3954: }
3956: /*@
3957: TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`
3959: Collective
3961: Input Parameter:
3962: . ts - The `TS` context obtained from `TSCreate()`
3964: Level: developer
3966: Note:
3967: `TSResize()` is typically used within time stepping implementations,
3968: so most users would not generally call this routine themselves.
3970: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3971: @*/
3972: PetscErrorCode TSResize(TS ts)
3973: {
3974: PetscInt nv = 0;
3975: const char **names = NULL;
3976: Vec *vecsin = NULL;
3977: const char *solname = "ts:vec_sol";
3979: PetscFunctionBegin;
3981: if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
3982: if (ts->resizesetup) {
3983: PetscCall(VecLockReadPush(ts->vec_sol));
3984: PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &ts->stepresize, ts->resizectx));
3985: PetscCall(VecLockReadPop(ts->vec_sol));
3986: if (ts->stepresize) {
3987: if (ts->resizerollback) {
3988: PetscCall(TSRollBack(ts));
3989: ts->time_step = ts->time_step0;
3990: }
3991: PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
3992: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
3993: }
3994: }
3996: PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
3997: if (nv) {
3998: Vec *vecsout, vecsol;
4000: /* Reset internal objects */
4001: PetscCall(TSReset(ts));
4003: /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
4004: PetscCall(PetscCalloc1(nv, &vecsout));
4005: PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
4006: for (PetscInt i = 0; i < nv; i++) {
4007: const char *name;
4008: char *oname;
4010: PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
4011: PetscCall(PetscStrallocpy(name, &oname));
4012: PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
4013: if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
4014: PetscCall(PetscFree(oname));
4015: PetscCall(VecDestroy(&vecsout[i]));
4016: }
4017: PetscCall(PetscFree(vecsout));
4018: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */
4020: PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
4021: if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
4022: PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
4023: }
4025: PetscCall(PetscFree(names));
4026: PetscCall(PetscFree(vecsin));
4027: PetscCall(TSResizeReset(ts));
4028: PetscFunctionReturn(PETSC_SUCCESS);
4029: }
4031: /*@
4032: TSSolve - Steps the requested number of timesteps.
4034: Collective
4036: Input Parameters:
4037: + ts - the `TS` context obtained from `TSCreate()`
4038: - u - the solution vector (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
4039: otherwise it must contain the initial conditions and will contain the solution at the final requested time
4041: Level: beginner
4043: Notes:
4044: The final time returned by this function may be different from the time of the internally
4045: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
4046: stepped over the final time.
4048: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
4049: @*/
4050: PetscErrorCode TSSolve(TS ts, Vec u)
4051: {
4052: Vec solution;
4054: PetscFunctionBegin;
4058: PetscCall(TSSetExactFinalTimeDefault(ts));
4059: 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 */
4060: if (!ts->vec_sol || u == ts->vec_sol) {
4061: PetscCall(VecDuplicate(u, &solution));
4062: PetscCall(TSSetSolution(ts, solution));
4063: PetscCall(VecDestroy(&solution)); /* grant ownership */
4064: }
4065: PetscCall(VecCopy(u, ts->vec_sol));
4066: PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
4067: } else if (u) PetscCall(TSSetSolution(ts, u));
4068: PetscCall(TSSetUp(ts));
4069: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
4071: 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>");
4072: 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()");
4073: 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");
4074: 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");
4076: if (ts->eval_times) {
4077: if (!ts->eval_times->sol_vecs) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
4078: for (PetscInt i = 0; i < ts->eval_times->num_time_points; i++) {
4079: PetscBool is_close = PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[i], ts->eval_times->reltol * ts->time_step + ts->eval_times->abstol, 0);
4080: if (ts->ptime <= ts->eval_times->time_points[i] || is_close) {
4081: ts->eval_times->time_point_idx = i;
4083: PetscBool is_ptime_in_sol_times = PETSC_FALSE; // If current solution has already been saved, we should not save it again
4084: 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);
4085: if (is_close && !is_ptime_in_sol_times) {
4086: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4087: ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4088: ts->eval_times->sol_idx++;
4089: ts->eval_times->time_point_idx++;
4090: }
4091: break;
4092: }
4093: }
4094: }
4096: if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));
4098: /* reset number of steps only when the step is not restarted. ARKIMEX
4099: restarts the step after an event. Resetting these counters in such case causes
4100: TSTrajectory to incorrectly save the output files
4101: */
4102: /* reset time step and iteration counters */
4103: if (!ts->steps) {
4104: ts->ksp_its = 0;
4105: ts->snes_its = 0;
4106: ts->num_snes_failures = 0;
4107: ts->reject = 0;
4108: ts->steprestart = PETSC_TRUE;
4109: ts->steprollback = PETSC_FALSE;
4110: ts->stepresize = PETSC_FALSE;
4111: ts->rhsjacobian.time = PETSC_MIN_REAL;
4112: }
4114: /* make sure initial time step does not overshoot final time or the next point in evaluation times */
4115: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
4116: PetscReal maxdt;
4117: PetscReal dt = ts->time_step;
4119: if (ts->eval_times) maxdt = ts->eval_times->time_points[ts->eval_times->time_point_idx] - ts->ptime;
4120: else maxdt = ts->max_time - ts->ptime;
4121: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
4122: }
4123: ts->reason = TS_CONVERGED_ITERATING;
4125: {
4126: PetscViewer viewer;
4127: PetscViewerFormat format;
4128: PetscBool flg;
4129: static PetscBool incall = PETSC_FALSE;
4131: if (!incall) {
4132: /* Estimate the convergence rate of the time discretization */
4133: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
4134: if (flg) {
4135: PetscConvEst conv;
4136: DM dm;
4137: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4138: PetscInt Nf;
4139: PetscBool checkTemporal = PETSC_TRUE;
4141: incall = PETSC_TRUE;
4142: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
4143: PetscCall(TSGetDM(ts, &dm));
4144: PetscCall(DMGetNumFields(dm, &Nf));
4145: PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
4146: PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
4147: PetscCall(PetscConvEstUseTS(conv, checkTemporal));
4148: PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
4149: PetscCall(PetscConvEstSetFromOptions(conv));
4150: PetscCall(PetscConvEstSetUp(conv));
4151: PetscCall(PetscConvEstGetConvRate(conv, alpha));
4152: PetscCall(PetscViewerPushFormat(viewer, format));
4153: PetscCall(PetscConvEstRateView(conv, alpha, viewer));
4154: PetscCall(PetscViewerPopFormat(viewer));
4155: PetscCall(PetscViewerDestroy(&viewer));
4156: PetscCall(PetscConvEstDestroy(&conv));
4157: PetscCall(PetscFree(alpha));
4158: incall = PETSC_FALSE;
4159: }
4160: }
4161: }
4163: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));
4165: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4166: PetscUseTypeMethod(ts, solve);
4167: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4168: ts->solvetime = ts->ptime;
4169: solution = ts->vec_sol;
4170: } else { /* Step the requested number of timesteps. */
4171: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4172: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4174: if (!ts->steps) {
4175: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4176: PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4177: }
4179: ts->start_step = ts->steps; // records starting step
4180: while (!ts->reason) {
4181: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4182: if (!ts->steprollback || (ts->stepresize && ts->resizerollback)) PetscCall(TSPreStep(ts));
4183: PetscCall(TSStep(ts));
4184: if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4185: if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4186: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4187: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4188: PetscCall(TSForwardCostIntegral(ts));
4189: if (ts->reason >= 0) ts->steps++;
4190: }
4191: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4192: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4193: PetscCall(TSForwardStep(ts));
4194: if (ts->reason >= 0) ts->steps++;
4195: }
4196: PetscCall(TSPostEvaluate(ts));
4197: 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. */
4198: if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4199: if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
4200: /* check convergence */
4201: if (!ts->reason) {
4202: if ((ts->steps - ts->start_step) >= ts->run_steps) ts->reason = TS_CONVERGED_ITS;
4203: else if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4204: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4205: }
4206: if (!ts->steprollback) {
4207: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4208: PetscCall(TSPostStep(ts));
4209: if (!ts->resizerollback) PetscCall(TSResize(ts));
4211: if (ts->eval_times && ts->eval_times->time_point_idx < ts->eval_times->num_time_points && ts->reason >= 0) {
4212: PetscCheck(ts->eval_times->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(eval_times->worktol > 0) in TSSolve()");
4213: if (PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[ts->eval_times->time_point_idx], ts->eval_times->worktol, 0)) {
4214: ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4215: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4216: ts->eval_times->sol_idx++;
4217: ts->eval_times->time_point_idx++;
4218: }
4219: }
4220: }
4221: }
4222: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4224: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4225: if (!u) u = ts->vec_sol;
4226: PetscCall(TSInterpolate(ts, ts->max_time, u));
4227: ts->solvetime = ts->max_time;
4228: solution = u;
4229: PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4230: } else {
4231: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4232: ts->solvetime = ts->ptime;
4233: solution = ts->vec_sol;
4234: }
4235: }
4237: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4238: PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4239: PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4240: if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4241: PetscFunctionReturn(PETSC_SUCCESS);
4242: }
4244: /*@
4245: TSGetTime - Gets the time of the most recently completed step.
4247: Not Collective
4249: Input Parameter:
4250: . ts - the `TS` context obtained from `TSCreate()`
4252: Output Parameter:
4253: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
4255: Level: beginner
4257: Note:
4258: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
4259: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
4261: .seealso: [](ch_ts), `TS`, `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4262: @*/
4263: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4264: {
4265: PetscFunctionBegin;
4267: PetscAssertPointer(t, 2);
4268: *t = ts->ptime;
4269: PetscFunctionReturn(PETSC_SUCCESS);
4270: }
4272: /*@
4273: TSGetPrevTime - Gets the starting time of the previously completed step.
4275: Not Collective
4277: Input Parameter:
4278: . ts - the `TS` context obtained from `TSCreate()`
4280: Output Parameter:
4281: . t - the previous time
4283: Level: beginner
4285: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4286: @*/
4287: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4288: {
4289: PetscFunctionBegin;
4291: PetscAssertPointer(t, 2);
4292: *t = ts->ptime_prev;
4293: PetscFunctionReturn(PETSC_SUCCESS);
4294: }
4296: /*@
4297: TSSetTime - Allows one to reset the time.
4299: Logically Collective
4301: Input Parameters:
4302: + ts - the `TS` context obtained from `TSCreate()`
4303: - t - the time
4305: Level: intermediate
4307: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4308: @*/
4309: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4310: {
4311: PetscFunctionBegin;
4314: ts->ptime = t;
4315: PetscFunctionReturn(PETSC_SUCCESS);
4316: }
4318: /*@
4319: TSSetOptionsPrefix - Sets the prefix used for searching for all
4320: TS options in the database.
4322: Logically Collective
4324: Input Parameters:
4325: + ts - The `TS` context
4326: - prefix - The prefix to prepend to all option names
4328: Level: advanced
4330: Note:
4331: A hyphen (-) must NOT be given at the beginning of the prefix name.
4332: The first character of all runtime options is AUTOMATICALLY the
4333: hyphen.
4335: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4336: @*/
4337: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4338: {
4339: SNES snes;
4341: PetscFunctionBegin;
4343: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4344: PetscCall(TSGetSNES(ts, &snes));
4345: PetscCall(SNESSetOptionsPrefix(snes, prefix));
4346: PetscFunctionReturn(PETSC_SUCCESS);
4347: }
4349: /*@
4350: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4351: TS options in the database.
4353: Logically Collective
4355: Input Parameters:
4356: + ts - The `TS` context
4357: - prefix - The prefix to prepend to all option names
4359: Level: advanced
4361: Note:
4362: A hyphen (-) must NOT be given at the beginning of the prefix name.
4363: The first character of all runtime options is AUTOMATICALLY the
4364: hyphen.
4366: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4367: @*/
4368: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4369: {
4370: SNES snes;
4372: PetscFunctionBegin;
4374: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4375: PetscCall(TSGetSNES(ts, &snes));
4376: PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4377: PetscFunctionReturn(PETSC_SUCCESS);
4378: }
4380: /*@
4381: TSGetOptionsPrefix - Sets the prefix used for searching for all
4382: `TS` options in the database.
4384: Not Collective
4386: Input Parameter:
4387: . ts - The `TS` context
4389: Output Parameter:
4390: . prefix - A pointer to the prefix string used
4392: Level: intermediate
4394: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4395: @*/
4396: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4397: {
4398: PetscFunctionBegin;
4400: PetscAssertPointer(prefix, 2);
4401: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4402: PetscFunctionReturn(PETSC_SUCCESS);
4403: }
4405: /*@C
4406: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4408: Not Collective, but parallel objects are returned if ts is parallel
4410: Input Parameter:
4411: . ts - The `TS` context obtained from `TSCreate()`
4413: Output Parameters:
4414: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`)
4415: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`)
4416: . func - Function to compute the Jacobian of the RHS (or `NULL`)
4417: - ctx - User-defined context for Jacobian evaluation routine (or `NULL`)
4419: Level: intermediate
4421: Note:
4422: You can pass in `NULL` for any return argument you do not need.
4424: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4426: @*/
4427: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, void **ctx)
4428: {
4429: DM dm;
4431: PetscFunctionBegin;
4432: if (Amat || Pmat) {
4433: SNES snes;
4434: PetscCall(TSGetSNES(ts, &snes));
4435: PetscCall(SNESSetUpMatrices(snes));
4436: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4437: }
4438: PetscCall(TSGetDM(ts, &dm));
4439: PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4440: PetscFunctionReturn(PETSC_SUCCESS);
4441: }
4443: /*@C
4444: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4446: Not Collective, but parallel objects are returned if ts is parallel
4448: Input Parameter:
4449: . ts - The `TS` context obtained from `TSCreate()`
4451: Output Parameters:
4452: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4453: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4454: . f - The function to compute the matrices
4455: - ctx - User-defined context for Jacobian evaluation routine
4457: Level: advanced
4459: Note:
4460: You can pass in `NULL` for any return argument you do not need.
4462: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4463: @*/
4464: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, void **ctx)
4465: {
4466: DM dm;
4468: PetscFunctionBegin;
4469: if (Amat || Pmat) {
4470: SNES snes;
4471: PetscCall(TSGetSNES(ts, &snes));
4472: PetscCall(SNESSetUpMatrices(snes));
4473: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4474: }
4475: PetscCall(TSGetDM(ts, &dm));
4476: PetscCall(DMTSGetIJacobian(dm, f, ctx));
4477: PetscFunctionReturn(PETSC_SUCCESS);
4478: }
4480: #include <petsc/private/dmimpl.h>
4481: /*@
4482: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4484: Logically Collective
4486: Input Parameters:
4487: + ts - the `TS` integrator object
4488: - dm - the dm, cannot be `NULL`
4490: Level: intermediate
4492: Notes:
4493: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4494: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4495: different problems using the same function space.
4497: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4498: @*/
4499: PetscErrorCode TSSetDM(TS ts, DM dm)
4500: {
4501: SNES snes;
4502: DMTS tsdm;
4504: PetscFunctionBegin;
4507: PetscCall(PetscObjectReference((PetscObject)dm));
4508: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4509: if (ts->dm->dmts && !dm->dmts) {
4510: PetscCall(DMCopyDMTS(ts->dm, dm));
4511: PetscCall(DMGetDMTS(ts->dm, &tsdm));
4512: /* Grant write privileges to the replacement DM */
4513: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4514: }
4515: PetscCall(DMDestroy(&ts->dm));
4516: }
4517: ts->dm = dm;
4519: PetscCall(TSGetSNES(ts, &snes));
4520: PetscCall(SNESSetDM(snes, dm));
4521: PetscFunctionReturn(PETSC_SUCCESS);
4522: }
4524: /*@
4525: TSGetDM - Gets the `DM` that may be used by some preconditioners
4527: Not Collective
4529: Input Parameter:
4530: . ts - the `TS`
4532: Output Parameter:
4533: . dm - the `DM`
4535: Level: intermediate
4537: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4538: @*/
4539: PetscErrorCode TSGetDM(TS ts, DM *dm)
4540: {
4541: PetscFunctionBegin;
4543: if (!ts->dm) {
4544: PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4545: if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4546: }
4547: *dm = ts->dm;
4548: PetscFunctionReturn(PETSC_SUCCESS);
4549: }
4551: /*@
4552: SNESTSFormFunction - Function to evaluate nonlinear residual defined by an ODE solver algorithm implemented within `TS`
4554: Logically Collective
4556: Input Parameters:
4557: + snes - nonlinear solver
4558: . U - the current state at which to evaluate the residual
4559: - ctx - user context, must be a `TS`
4561: Output Parameter:
4562: . F - the nonlinear residual
4564: Level: developer
4566: Note:
4567: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4568: It is most frequently passed to `MatFDColoringSetFunction()`.
4570: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4571: @*/
4572: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4573: {
4574: TS ts = (TS)ctx;
4576: PetscFunctionBegin;
4581: PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4582: PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4583: PetscFunctionReturn(PETSC_SUCCESS);
4584: }
4586: /*@
4587: SNESTSFormJacobian - Function to evaluate the Jacobian defined by an ODE solver algorithm implemented within `TS`
4589: Collective
4591: Input Parameters:
4592: + snes - nonlinear solver
4593: . U - the current state at which to evaluate the residual
4594: - ctx - user context, must be a `TS`
4596: Output Parameters:
4597: + A - the Jacobian
4598: - B - the matrix used to construct the preconditioner (often the same as `A`)
4600: Level: developer
4602: Note:
4603: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4605: .seealso: [](ch_ts), `SNESSetJacobian()`
4606: @*/
4607: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4608: {
4609: TS ts = (TS)ctx;
4611: PetscFunctionBegin;
4617: PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4618: PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4619: PetscFunctionReturn(PETSC_SUCCESS);
4620: }
4622: /*@C
4623: TSComputeRHSFunctionLinear - Evaluate the right-hand side via the user-provided Jacobian, for linear problems Udot = A U only
4625: Collective
4627: Input Parameters:
4628: + ts - time stepping context
4629: . t - time at which to evaluate
4630: . U - state at which to evaluate
4631: - ctx - context
4633: Output Parameter:
4634: . F - right-hand side
4636: Level: intermediate
4638: Note:
4639: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right-hand side for linear problems.
4640: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4642: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4643: @*/
4644: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4645: {
4646: Mat Arhs, Brhs;
4648: PetscFunctionBegin;
4649: PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4650: /* undo the damage caused by shifting */
4651: PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4652: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4653: PetscCall(MatMult(Arhs, U, F));
4654: PetscFunctionReturn(PETSC_SUCCESS);
4655: }
4657: /*@C
4658: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
4660: Collective
4662: Input Parameters:
4663: + ts - time stepping context
4664: . t - time at which to evaluate
4665: . U - state at which to evaluate
4666: - ctx - context
4668: Output Parameters:
4669: + A - Jacobian
4670: - B - matrix used to construct the preconditioner, often the same as `A`
4672: Level: intermediate
4674: Note:
4675: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4677: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4678: @*/
4679: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4680: {
4681: PetscFunctionBegin;
4682: PetscFunctionReturn(PETSC_SUCCESS);
4683: }
4685: /*@C
4686: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4688: Collective
4690: Input Parameters:
4691: + ts - time stepping context
4692: . t - time at which to evaluate
4693: . U - state at which to evaluate
4694: . Udot - time derivative of state vector
4695: - ctx - context
4697: Output Parameter:
4698: . F - left hand side
4700: Level: intermediate
4702: Notes:
4703: 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
4704: user is required to write their own `TSComputeIFunction()`.
4705: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4706: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4708: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4710: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4711: @*/
4712: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4713: {
4714: Mat A, B;
4716: PetscFunctionBegin;
4717: PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4718: PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4719: PetscCall(MatMult(A, Udot, F));
4720: PetscFunctionReturn(PETSC_SUCCESS);
4721: }
4723: /*@C
4724: TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE
4726: Collective
4728: Input Parameters:
4729: + ts - time stepping context
4730: . t - time at which to evaluate
4731: . U - state at which to evaluate
4732: . Udot - time derivative of state vector
4733: . shift - shift to apply
4734: - ctx - context
4736: Output Parameters:
4737: + A - pointer to operator
4738: - B - pointer to matrix from which the preconditioner is built (often `A`)
4740: Level: advanced
4742: Notes:
4743: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4745: It is only appropriate for problems of the form
4747: $$
4748: M \dot{U} = F(U,t)
4749: $$
4751: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4752: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4753: an implicit operator of the form
4755: $$
4756: shift*M + J
4757: $$
4759: 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
4760: a copy of M or reassemble it when requested.
4762: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4763: @*/
4764: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4765: {
4766: PetscFunctionBegin;
4767: PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4768: ts->ijacobian.shift = shift;
4769: PetscFunctionReturn(PETSC_SUCCESS);
4770: }
4772: /*@
4773: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4775: Not Collective
4777: Input Parameter:
4778: . ts - the `TS` context
4780: Output Parameter:
4781: . equation_type - see `TSEquationType`
4783: Level: beginner
4785: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4786: @*/
4787: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4788: {
4789: PetscFunctionBegin;
4791: PetscAssertPointer(equation_type, 2);
4792: *equation_type = ts->equation_type;
4793: PetscFunctionReturn(PETSC_SUCCESS);
4794: }
4796: /*@
4797: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4799: Not Collective
4801: Input Parameters:
4802: + ts - the `TS` context
4803: - equation_type - see `TSEquationType`
4805: Level: advanced
4807: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4808: @*/
4809: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4810: {
4811: PetscFunctionBegin;
4813: ts->equation_type = equation_type;
4814: PetscFunctionReturn(PETSC_SUCCESS);
4815: }
4817: /*@
4818: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4820: Not Collective
4822: Input Parameter:
4823: . ts - the `TS` context
4825: Output Parameter:
4826: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4827: manual pages for the individual convergence tests for complete lists
4829: Level: beginner
4831: Note:
4832: Can only be called after the call to `TSSolve()` is complete.
4834: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4835: @*/
4836: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4837: {
4838: PetscFunctionBegin;
4840: PetscAssertPointer(reason, 2);
4841: *reason = ts->reason;
4842: PetscFunctionReturn(PETSC_SUCCESS);
4843: }
4845: /*@
4846: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4848: Logically Collective; reason must contain common value
4850: Input Parameters:
4851: + ts - the `TS` context
4852: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4853: manual pages for the individual convergence tests for complete lists
4855: Level: advanced
4857: Note:
4858: Can only be called while `TSSolve()` is active.
4860: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4861: @*/
4862: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4863: {
4864: PetscFunctionBegin;
4866: ts->reason = reason;
4867: PetscFunctionReturn(PETSC_SUCCESS);
4868: }
4870: /*@
4871: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4873: Not Collective
4875: Input Parameter:
4876: . ts - the `TS` context
4878: Output Parameter:
4879: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4881: Level: beginner
4883: Note:
4884: Can only be called after the call to `TSSolve()` is complete.
4886: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4887: @*/
4888: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4889: {
4890: PetscFunctionBegin;
4892: PetscAssertPointer(ftime, 2);
4893: *ftime = ts->solvetime;
4894: PetscFunctionReturn(PETSC_SUCCESS);
4895: }
4897: /*@
4898: TSGetSNESIterations - Gets the total number of nonlinear iterations
4899: used by the time integrator.
4901: Not Collective
4903: Input Parameter:
4904: . ts - `TS` context
4906: Output Parameter:
4907: . nits - number of nonlinear iterations
4909: Level: intermediate
4911: Note:
4912: This counter is reset to zero for each successive call to `TSSolve()`.
4914: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4915: @*/
4916: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4917: {
4918: PetscFunctionBegin;
4920: PetscAssertPointer(nits, 2);
4921: *nits = ts->snes_its;
4922: PetscFunctionReturn(PETSC_SUCCESS);
4923: }
4925: /*@
4926: TSGetKSPIterations - Gets the total number of linear iterations
4927: used by the time integrator.
4929: Not Collective
4931: Input Parameter:
4932: . ts - `TS` context
4934: Output Parameter:
4935: . lits - number of linear iterations
4937: Level: intermediate
4939: Note:
4940: This counter is reset to zero for each successive call to `TSSolve()`.
4942: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`
4943: @*/
4944: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4945: {
4946: PetscFunctionBegin;
4948: PetscAssertPointer(lits, 2);
4949: *lits = ts->ksp_its;
4950: PetscFunctionReturn(PETSC_SUCCESS);
4951: }
4953: /*@
4954: TSGetStepRejections - Gets the total number of rejected steps.
4956: Not Collective
4958: Input Parameter:
4959: . ts - `TS` context
4961: Output Parameter:
4962: . rejects - number of steps rejected
4964: Level: intermediate
4966: Note:
4967: This counter is reset to zero for each successive call to `TSSolve()`.
4969: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4970: @*/
4971: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4972: {
4973: PetscFunctionBegin;
4975: PetscAssertPointer(rejects, 2);
4976: *rejects = ts->reject;
4977: PetscFunctionReturn(PETSC_SUCCESS);
4978: }
4980: /*@
4981: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
4983: Not Collective
4985: Input Parameter:
4986: . ts - `TS` context
4988: Output Parameter:
4989: . fails - number of failed nonlinear solves
4991: Level: intermediate
4993: Note:
4994: This counter is reset to zero for each successive call to `TSSolve()`.
4996: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4997: @*/
4998: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4999: {
5000: PetscFunctionBegin;
5002: PetscAssertPointer(fails, 2);
5003: *fails = ts->num_snes_failures;
5004: PetscFunctionReturn(PETSC_SUCCESS);
5005: }
5007: /*@
5008: TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails
5010: Not Collective
5012: Input Parameters:
5013: + ts - `TS` context
5014: - rejects - maximum number of rejected steps, pass `PETSC_UNLIMITED` for unlimited
5016: Options Database Key:
5017: . -ts_max_reject - Maximum number of step rejections before a step fails
5019: Level: intermediate
5021: Developer Note:
5022: The options database name is incorrect.
5024: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
5025: @*/
5026: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
5027: {
5028: PetscFunctionBegin;
5030: if (rejects == PETSC_UNLIMITED || rejects == -1) {
5031: ts->max_reject = PETSC_UNLIMITED;
5032: } else {
5033: PetscCheck(rejects >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of rejections");
5034: ts->max_reject = rejects;
5035: }
5036: PetscFunctionReturn(PETSC_SUCCESS);
5037: }
5039: /*@
5040: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves
5042: Not Collective
5044: Input Parameters:
5045: + ts - `TS` context
5046: - fails - maximum number of failed nonlinear solves, pass `PETSC_UNLIMITED` to allow any number of failures.
5048: Options Database Key:
5049: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
5051: Level: intermediate
5053: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
5054: @*/
5055: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
5056: {
5057: PetscFunctionBegin;
5059: if (fails == PETSC_UNLIMITED || fails == -1) {
5060: ts->max_snes_failures = PETSC_UNLIMITED;
5061: } else {
5062: PetscCheck(fails >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures");
5063: ts->max_snes_failures = fails;
5064: }
5065: PetscFunctionReturn(PETSC_SUCCESS);
5066: }
5068: /*@
5069: TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`
5071: Not Collective
5073: Input Parameters:
5074: + ts - `TS` context
5075: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
5077: Options Database Key:
5078: . -ts_error_if_step_fails - Error if no step succeeds
5080: Level: intermediate
5082: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
5083: @*/
5084: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
5085: {
5086: PetscFunctionBegin;
5088: ts->errorifstepfailed = err;
5089: PetscFunctionReturn(PETSC_SUCCESS);
5090: }
5092: /*@
5093: TSGetAdapt - Get the adaptive controller context for the current method
5095: Collective if controller has not yet been created
5097: Input Parameter:
5098: . ts - time stepping context
5100: Output Parameter:
5101: . adapt - adaptive controller
5103: Level: intermediate
5105: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
5106: @*/
5107: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
5108: {
5109: PetscFunctionBegin;
5111: PetscAssertPointer(adapt, 2);
5112: if (!ts->adapt) {
5113: PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
5114: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
5115: }
5116: *adapt = ts->adapt;
5117: PetscFunctionReturn(PETSC_SUCCESS);
5118: }
5120: /*@
5121: TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller
5123: Logically Collective
5125: Input Parameters:
5126: + ts - time integration context
5127: . atol - scalar absolute tolerances
5128: . vatol - vector of absolute tolerances or `NULL`, used in preference to `atol` if present
5129: . rtol - scalar relative tolerances
5130: - vrtol - vector of relative tolerances or `NULL`, used in preference to `rtol` if present
5132: Options Database Keys:
5133: + -ts_rtol <rtol> - relative tolerance for local truncation error
5134: - -ts_atol <atol> - Absolute tolerance for local truncation error
5136: Level: beginner
5138: Notes:
5139: `PETSC_CURRENT` or `PETSC_DETERMINE` may be used for `atol` or `rtol` to indicate the current value
5140: or the default value from when the object's type was set.
5142: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5143: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5144: computed only for the differential or the algebraic part then this can be done using the vector of
5145: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
5146: differential part and infinity for the algebraic part, the LTE calculation will include only the
5147: differential variables.
5149: Fortran Note:
5150: Use `PETSC_CURRENT_INTEGER` or `PETSC_DETERMINE_INTEGER`.
5152: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
5153: @*/
5154: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
5155: {
5156: PetscFunctionBegin;
5157: if (atol == (PetscReal)PETSC_DETERMINE) {
5158: ts->atol = ts->default_atol;
5159: } else if (atol != (PetscReal)PETSC_CURRENT) {
5160: PetscCheck(atol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)atol);
5161: ts->atol = atol;
5162: }
5164: if (vatol) {
5165: PetscCall(PetscObjectReference((PetscObject)vatol));
5166: PetscCall(VecDestroy(&ts->vatol));
5167: ts->vatol = vatol;
5168: }
5170: if (rtol == (PetscReal)PETSC_DETERMINE) {
5171: ts->rtol = ts->default_rtol;
5172: } else if (rtol != (PetscReal)PETSC_CURRENT) {
5173: PetscCheck(rtol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative", (double)rtol);
5174: ts->rtol = rtol;
5175: }
5177: if (vrtol) {
5178: PetscCall(PetscObjectReference((PetscObject)vrtol));
5179: PetscCall(VecDestroy(&ts->vrtol));
5180: ts->vrtol = vrtol;
5181: }
5182: PetscFunctionReturn(PETSC_SUCCESS);
5183: }
5185: /*@
5186: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
5188: Logically Collective
5190: Input Parameter:
5191: . ts - time integration context
5193: Output Parameters:
5194: + atol - scalar absolute tolerances, `NULL` to ignore
5195: . vatol - vector of absolute tolerances, `NULL` to ignore
5196: . rtol - scalar relative tolerances, `NULL` to ignore
5197: - vrtol - vector of relative tolerances, `NULL` to ignore
5199: Level: beginner
5201: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
5202: @*/
5203: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
5204: {
5205: PetscFunctionBegin;
5206: if (atol) *atol = ts->atol;
5207: if (vatol) *vatol = ts->vatol;
5208: if (rtol) *rtol = ts->rtol;
5209: if (vrtol) *vrtol = ts->vrtol;
5210: PetscFunctionReturn(PETSC_SUCCESS);
5211: }
5213: /*@
5214: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
5216: Collective
5218: Input Parameters:
5219: + ts - time stepping context
5220: . U - state vector, usually ts->vec_sol
5221: . Y - state vector to be compared to U
5222: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5224: Output Parameters:
5225: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5226: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5227: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5229: Options Database Key:
5230: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5232: Level: developer
5234: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5235: @*/
5236: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5237: {
5238: PetscInt norma_loc, norm_loc, normr_loc;
5240: PetscFunctionBegin;
5245: PetscAssertPointer(norm, 5);
5246: PetscAssertPointer(norma, 6);
5247: PetscAssertPointer(normr, 7);
5248: 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));
5249: if (wnormtype == NORM_2) {
5250: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5251: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5252: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5253: }
5254: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5255: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5256: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5257: PetscFunctionReturn(PETSC_SUCCESS);
5258: }
5260: /*@
5261: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5263: Collective
5265: Input Parameters:
5266: + ts - time stepping context
5267: . E - error vector
5268: . U - state vector, usually ts->vec_sol
5269: . Y - state vector, previous time step
5270: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5272: Output Parameters:
5273: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5274: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5275: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5277: Options Database Key:
5278: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5280: Level: developer
5282: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5283: @*/
5284: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5285: {
5286: PetscInt norma_loc, norm_loc, normr_loc;
5288: PetscFunctionBegin;
5290: 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));
5291: if (wnormtype == NORM_2) {
5292: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5293: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5294: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5295: }
5296: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5297: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5298: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5299: PetscFunctionReturn(PETSC_SUCCESS);
5300: }
5302: /*@
5303: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5305: Logically Collective
5307: Input Parameters:
5308: + ts - time stepping context
5309: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5311: Note:
5312: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5314: Level: intermediate
5316: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5317: @*/
5318: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5319: {
5320: PetscFunctionBegin;
5322: ts->cfltime_local = cfltime;
5323: ts->cfltime = -1.;
5324: PetscFunctionReturn(PETSC_SUCCESS);
5325: }
5327: /*@
5328: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5330: Collective
5332: Input Parameter:
5333: . ts - time stepping context
5335: Output Parameter:
5336: . cfltime - maximum stable time step for forward Euler
5338: Level: advanced
5340: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5341: @*/
5342: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5343: {
5344: PetscFunctionBegin;
5345: if (ts->cfltime < 0) PetscCallMPI(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5346: *cfltime = ts->cfltime;
5347: PetscFunctionReturn(PETSC_SUCCESS);
5348: }
5350: /*@
5351: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5353: Input Parameters:
5354: + ts - the `TS` context.
5355: . xl - lower bound.
5356: - xu - upper bound.
5358: Level: advanced
5360: Note:
5361: If this routine is not called then the lower and upper bounds are set to
5362: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5364: .seealso: [](ch_ts), `TS`
5365: @*/
5366: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5367: {
5368: SNES snes;
5370: PetscFunctionBegin;
5371: PetscCall(TSGetSNES(ts, &snes));
5372: PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5373: PetscFunctionReturn(PETSC_SUCCESS);
5374: }
5376: /*@
5377: TSComputeLinearStability - computes the linear stability function at a point
5379: Collective
5381: Input Parameters:
5382: + ts - the `TS` context
5383: . xr - real part of input argument
5384: - xi - imaginary part of input argument
5386: Output Parameters:
5387: + yr - real part of function value
5388: - yi - imaginary part of function value
5390: Level: developer
5392: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5393: @*/
5394: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5395: {
5396: PetscFunctionBegin;
5398: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5399: PetscFunctionReturn(PETSC_SUCCESS);
5400: }
5402: /*@
5403: TSRestartStep - Flags the solver to restart the next step
5405: Collective
5407: Input Parameter:
5408: . ts - the `TS` context obtained from `TSCreate()`
5410: Level: advanced
5412: Notes:
5413: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5414: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5415: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5416: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5417: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5418: discontinuous source terms).
5420: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5421: @*/
5422: PetscErrorCode TSRestartStep(TS ts)
5423: {
5424: PetscFunctionBegin;
5426: ts->steprestart = PETSC_TRUE;
5427: PetscFunctionReturn(PETSC_SUCCESS);
5428: }
5430: /*@
5431: TSRollBack - Rolls back one time step
5433: Collective
5435: Input Parameter:
5436: . ts - the `TS` context obtained from `TSCreate()`
5438: Level: advanced
5440: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5441: @*/
5442: PetscErrorCode TSRollBack(TS ts)
5443: {
5444: PetscFunctionBegin;
5446: PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5447: PetscTryTypeMethod(ts, rollback);
5448: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
5449: ts->time_step = ts->ptime - ts->ptime_prev;
5450: ts->ptime = ts->ptime_prev;
5451: ts->ptime_prev = ts->ptime_prev_rollback;
5452: ts->steps--;
5453: ts->steprollback = PETSC_TRUE;
5454: PetscFunctionReturn(PETSC_SUCCESS);
5455: }
5457: /*@
5458: TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step
5460: Not collective
5462: Input Parameter:
5463: . ts - the `TS` context obtained from `TSCreate()`
5465: Output Parameter:
5466: . flg - the rollback flag
5468: Level: advanced
5470: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
5471: @*/
5472: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
5473: {
5474: PetscFunctionBegin;
5476: PetscAssertPointer(flg, 2);
5477: *flg = ts->steprollback;
5478: PetscFunctionReturn(PETSC_SUCCESS);
5479: }
5481: /*@
5482: TSGetStepResize - Get the internal flag indicating if the current step is after a resize.
5484: Not collective
5486: Input Parameter:
5487: . ts - the `TS` context obtained from `TSCreate()`
5489: Output Parameter:
5490: . flg - the resize flag
5492: Level: advanced
5494: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetResize()`
5495: @*/
5496: PetscErrorCode TSGetStepResize(TS ts, PetscBool *flg)
5497: {
5498: PetscFunctionBegin;
5500: PetscAssertPointer(flg, 2);
5501: *flg = ts->stepresize;
5502: PetscFunctionReturn(PETSC_SUCCESS);
5503: }
5505: /*@
5506: TSGetStages - Get the number of stages and stage values
5508: Input Parameter:
5509: . ts - the `TS` context obtained from `TSCreate()`
5511: Output Parameters:
5512: + ns - the number of stages
5513: - Y - the current stage vectors
5515: Level: advanced
5517: Note:
5518: Both `ns` and `Y` can be `NULL`.
5520: .seealso: [](ch_ts), `TS`, `TSCreate()`
5521: @*/
5522: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5523: {
5524: PetscFunctionBegin;
5526: if (ns) PetscAssertPointer(ns, 2);
5527: if (Y) PetscAssertPointer(Y, 3);
5528: if (!ts->ops->getstages) {
5529: if (ns) *ns = 0;
5530: if (Y) *Y = NULL;
5531: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5532: PetscFunctionReturn(PETSC_SUCCESS);
5533: }
5535: /*@C
5536: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5538: Collective
5540: Input Parameters:
5541: + ts - the `TS` context
5542: . t - current timestep
5543: . U - state vector
5544: . Udot - time derivative of state vector
5545: . shift - shift to apply, see note below
5546: - ctx - an optional user context
5548: Output Parameters:
5549: + J - Jacobian matrix (not altered in this routine)
5550: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
5552: Level: intermediate
5554: Notes:
5555: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5557: dF/dU + shift*dF/dUdot
5559: Most users should not need to explicitly call this routine, as it
5560: is used internally within the nonlinear solvers.
5562: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5563: routine, then it will try to get the coloring from the matrix. This requires that the
5564: matrix have nonzero entries precomputed.
5566: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5567: @*/
5568: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5569: {
5570: SNES snes;
5571: MatFDColoring color;
5572: PetscBool hascolor, matcolor = PETSC_FALSE;
5574: PetscFunctionBegin;
5575: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5576: PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5577: if (!color) {
5578: DM dm;
5579: ISColoring iscoloring;
5581: PetscCall(TSGetDM(ts, &dm));
5582: PetscCall(DMHasColoring(dm, &hascolor));
5583: if (hascolor && !matcolor) {
5584: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5585: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5586: PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5587: PetscCall(MatFDColoringSetFromOptions(color));
5588: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5589: PetscCall(ISColoringDestroy(&iscoloring));
5590: } else {
5591: MatColoring mc;
5593: PetscCall(MatColoringCreate(B, &mc));
5594: PetscCall(MatColoringSetDistance(mc, 2));
5595: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5596: PetscCall(MatColoringSetFromOptions(mc));
5597: PetscCall(MatColoringApply(mc, &iscoloring));
5598: PetscCall(MatColoringDestroy(&mc));
5599: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5600: PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5601: PetscCall(MatFDColoringSetFromOptions(color));
5602: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5603: PetscCall(ISColoringDestroy(&iscoloring));
5604: }
5605: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5606: PetscCall(PetscObjectDereference((PetscObject)color));
5607: }
5608: PetscCall(TSGetSNES(ts, &snes));
5609: PetscCall(MatFDColoringApply(B, color, U, snes));
5610: if (J != B) {
5611: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5612: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5613: }
5614: PetscFunctionReturn(PETSC_SUCCESS);
5615: }
5617: /*@C
5618: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5620: Logically collective
5622: Input Parameters:
5623: + ts - the `TS` context
5624: - func - function called within `TSFunctionDomainError()`
5626: Calling sequence of `func`:
5627: + ts - the `TS` context
5628: . time - the current time (of the stage)
5629: . state - the state to check if it is valid
5630: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable
5632: Level: intermediate
5634: Notes:
5635: `accept` must be collectively specified.
5636: If an implicit ODE solver is being used then, in addition to providing this routine, the
5637: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5638: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5639: Use `TSGetSNES()` to obtain the `SNES` object
5641: Developer Notes:
5642: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5643: since one takes a function pointer and the other does not.
5645: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5646: @*/
5647: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5648: {
5649: PetscFunctionBegin;
5651: ts->functiondomainerror = func;
5652: PetscFunctionReturn(PETSC_SUCCESS);
5653: }
5655: /*@
5656: TSFunctionDomainError - Checks if the current state is valid
5658: Collective
5660: Input Parameters:
5661: + ts - the `TS` context
5662: . stagetime - time of the simulation
5663: - Y - state vector to check.
5665: Output Parameter:
5666: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5668: Level: developer
5670: Note:
5671: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5672: to check if the current state is valid.
5674: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5675: @*/
5676: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5677: {
5678: PetscFunctionBegin;
5682: PetscAssertPointer(accept, 4);
5683: *accept = PETSC_TRUE;
5684: if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5685: PetscFunctionReturn(PETSC_SUCCESS);
5686: }
5688: /*@
5689: TSClone - This function clones a time step `TS` object.
5691: Collective
5693: Input Parameter:
5694: . tsin - The input `TS`
5696: Output Parameter:
5697: . tsout - The output `TS` (cloned)
5699: Level: developer
5701: Notes:
5702: 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.
5703: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5705: 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
5706: .vb
5707: SNES snes_dup = NULL;
5708: TSGetSNES(ts,&snes_dup);
5709: TSSetSNES(ts,snes_dup);
5710: .ve
5712: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5713: @*/
5714: PetscErrorCode TSClone(TS tsin, TS *tsout)
5715: {
5716: TS t;
5717: SNES snes_start;
5718: DM dm;
5719: TSType type;
5721: PetscFunctionBegin;
5722: PetscAssertPointer(tsin, 1);
5723: *tsout = NULL;
5725: PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));
5727: /* General TS description */
5728: t->numbermonitors = 0;
5729: t->setupcalled = PETSC_FALSE;
5730: t->ksp_its = 0;
5731: t->snes_its = 0;
5732: t->nwork = 0;
5733: t->rhsjacobian.time = PETSC_MIN_REAL;
5734: t->rhsjacobian.scale = 1.;
5735: t->ijacobian.shift = 1.;
5737: PetscCall(TSGetSNES(tsin, &snes_start));
5738: PetscCall(TSSetSNES(t, snes_start));
5740: PetscCall(TSGetDM(tsin, &dm));
5741: PetscCall(TSSetDM(t, dm));
5743: t->adapt = tsin->adapt;
5744: PetscCall(PetscObjectReference((PetscObject)t->adapt));
5746: t->trajectory = tsin->trajectory;
5747: PetscCall(PetscObjectReference((PetscObject)t->trajectory));
5749: t->event = tsin->event;
5750: if (t->event) t->event->refct++;
5752: t->problem_type = tsin->problem_type;
5753: t->ptime = tsin->ptime;
5754: t->ptime_prev = tsin->ptime_prev;
5755: t->time_step = tsin->time_step;
5756: t->max_time = tsin->max_time;
5757: t->steps = tsin->steps;
5758: t->max_steps = tsin->max_steps;
5759: t->equation_type = tsin->equation_type;
5760: t->atol = tsin->atol;
5761: t->rtol = tsin->rtol;
5762: t->max_snes_failures = tsin->max_snes_failures;
5763: t->max_reject = tsin->max_reject;
5764: t->errorifstepfailed = tsin->errorifstepfailed;
5766: PetscCall(TSGetType(tsin, &type));
5767: PetscCall(TSSetType(t, type));
5769: t->vec_sol = NULL;
5771: t->cfltime = tsin->cfltime;
5772: t->cfltime_local = tsin->cfltime_local;
5773: t->exact_final_time = tsin->exact_final_time;
5775: t->ops[0] = tsin->ops[0];
5777: if (((PetscObject)tsin)->fortran_func_pointers) {
5778: PetscInt i;
5779: PetscCall(PetscMalloc((10) * sizeof(PetscFortranCallbackFn *), &((PetscObject)t)->fortran_func_pointers));
5780: for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5781: }
5782: *tsout = t;
5783: PetscFunctionReturn(PETSC_SUCCESS);
5784: }
5786: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5787: {
5788: TS ts = (TS)ctx;
5790: PetscFunctionBegin;
5791: PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5792: PetscFunctionReturn(PETSC_SUCCESS);
5793: }
5795: /*@
5796: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5798: Logically Collective
5800: Input Parameter:
5801: . ts - the time stepping routine
5803: Output Parameter:
5804: . flg - `PETSC_TRUE` if the multiply is likely correct
5806: Options Database Key:
5807: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5809: Level: advanced
5811: Note:
5812: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5814: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5815: @*/
5816: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5817: {
5818: Mat J, B;
5819: TSRHSJacobianFn *func;
5820: void *ctx;
5822: PetscFunctionBegin;
5823: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5824: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5825: PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5826: PetscFunctionReturn(PETSC_SUCCESS);
5827: }
5829: /*@
5830: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5832: Logically Collective
5834: Input Parameter:
5835: . ts - the time stepping routine
5837: Output Parameter:
5838: . flg - `PETSC_TRUE` if the multiply is likely correct
5840: Options Database Key:
5841: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5843: Level: advanced
5845: Notes:
5846: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5848: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5849: @*/
5850: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5851: {
5852: Mat J, B;
5853: void *ctx;
5854: TSRHSJacobianFn *func;
5856: PetscFunctionBegin;
5857: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5858: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5859: PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5860: PetscFunctionReturn(PETSC_SUCCESS);
5861: }
5863: /*@
5864: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
5866: Logically Collective
5868: Input Parameters:
5869: + ts - timestepping context
5870: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5872: Options Database Key:
5873: . -ts_use_splitrhsfunction - <true,false>
5875: Level: intermediate
5877: Note:
5878: This is only for multirate methods
5880: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5881: @*/
5882: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5883: {
5884: PetscFunctionBegin;
5886: ts->use_splitrhsfunction = use_splitrhsfunction;
5887: PetscFunctionReturn(PETSC_SUCCESS);
5888: }
5890: /*@
5891: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
5893: Not Collective
5895: Input Parameter:
5896: . ts - timestepping context
5898: Output Parameter:
5899: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5901: Level: intermediate
5903: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5904: @*/
5905: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5906: {
5907: PetscFunctionBegin;
5909: *use_splitrhsfunction = ts->use_splitrhsfunction;
5910: PetscFunctionReturn(PETSC_SUCCESS);
5911: }
5913: /*@
5914: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
5916: Logically Collective
5918: Input Parameters:
5919: + ts - the time-stepper
5920: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
5922: Level: intermediate
5924: Note:
5925: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
5927: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
5928: @*/
5929: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5930: {
5931: PetscFunctionBegin;
5933: ts->axpy_pattern = str;
5934: PetscFunctionReturn(PETSC_SUCCESS);
5935: }
5937: /*@
5938: TSSetEvaluationTimes - sets the evaluation points. The solution will be computed and stored for each time requested
5940: Collective
5942: Input Parameters:
5943: + ts - the time-stepper
5944: . n - number of the time points
5945: - time_points - array of the time points, must be increasing
5947: Options Database Key:
5948: . -ts_eval_times <t0,...tn> - Sets the evaluation times
5950: Level: intermediate
5952: Notes:
5953: The elements in `time_points` must be all increasing. They correspond to the intermediate points to be saved.
5955: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
5957: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using evaluation times may
5958: pressure the memory system when using a large number of time points.
5960: .seealso: [](ch_ts), `TS`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`, `TSSetTimeSpan()`
5961: @*/
5962: PetscErrorCode TSSetEvaluationTimes(TS ts, PetscInt n, PetscReal time_points[])
5963: {
5964: PetscBool is_sorted;
5966: PetscFunctionBegin;
5968: if (ts->eval_times) { // Reset eval_times
5969: ts->eval_times->sol_idx = 0;
5970: ts->eval_times->time_point_idx = 0;
5971: if (n != ts->eval_times->num_time_points) {
5972: PetscCall(PetscFree(ts->eval_times->time_points));
5973: PetscCall(PetscFree(ts->eval_times->sol_times));
5974: PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
5975: } else {
5976: PetscCall(PetscArrayzero(ts->eval_times->sol_times, n));
5977: for (PetscInt i = 0; i < n; i++) PetscCall(VecZeroEntries(ts->eval_times->sol_vecs[i]));
5978: }
5979: } else { // Create/initialize eval_times
5980: TSEvaluationTimes eval_times;
5981: PetscCall(PetscNew(&eval_times));
5982: PetscCall(PetscMalloc1(n, &eval_times->time_points));
5983: PetscCall(PetscMalloc1(n, &eval_times->sol_times));
5984: eval_times->reltol = 1e-6;
5985: eval_times->abstol = 10 * PETSC_MACHINE_EPSILON;
5986: eval_times->worktol = 0;
5987: ts->eval_times = eval_times;
5988: }
5989: ts->eval_times->num_time_points = n;
5990: PetscCall(PetscSortedReal(n, time_points, &is_sorted));
5991: PetscCheck(is_sorted, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "time_points array must be sorted");
5992: PetscCall(PetscArraycpy(ts->eval_times->time_points, time_points, n));
5993: // Note: ts->vec_sol not guaranteed to exist, so ts->eval_times->sol_vecs allocated at TSSolve time
5994: PetscFunctionReturn(PETSC_SUCCESS);
5995: }
5997: /*@C
5998: TSGetEvaluationTimes - gets the evaluation times set with `TSSetEvaluationTimes()`
6000: Not Collective
6002: Input Parameter:
6003: . ts - the time-stepper
6005: Output Parameters:
6006: + n - number of the time points
6007: - time_points - array of the time points
6009: Level: beginner
6011: Note:
6012: The values obtained are valid until the `TS` object is destroyed.
6014: Both `n` and `time_points` can be `NULL`.
6016: Also used to see time points set by `TSSetTimeSpan()`.
6018: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6019: @*/
6020: PetscErrorCode TSGetEvaluationTimes(TS ts, PetscInt *n, const PetscReal *time_points[])
6021: {
6022: PetscFunctionBegin;
6024: if (n) PetscAssertPointer(n, 2);
6025: if (time_points) PetscAssertPointer(time_points, 3);
6026: if (!ts->eval_times) {
6027: if (n) *n = 0;
6028: if (time_points) *time_points = NULL;
6029: } else {
6030: if (n) *n = ts->eval_times->num_time_points;
6031: if (time_points) *time_points = ts->eval_times->time_points;
6032: }
6033: PetscFunctionReturn(PETSC_SUCCESS);
6034: }
6036: /*@C
6037: TSGetEvaluationSolutions - Get the number of solutions and the solutions at the evaluation time points specified
6039: Input Parameter:
6040: . ts - the `TS` context obtained from `TSCreate()`
6042: Output Parameters:
6043: + nsol - the number of solutions
6044: . sol_times - array of solution times corresponding to the solution vectors. See note below
6045: - Sols - the solution vectors
6047: Level: intermediate
6049: Notes:
6050: Both `nsol` and `Sols` can be `NULL`.
6052: 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()`.
6053: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain evaluation times.
6055: Also used to see view solutions requested by `TSSetTimeSpan()`.
6057: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`
6058: @*/
6059: PetscErrorCode TSGetEvaluationSolutions(TS ts, PetscInt *nsol, const PetscReal *sol_times[], Vec *Sols[])
6060: {
6061: PetscFunctionBegin;
6063: if (nsol) PetscAssertPointer(nsol, 2);
6064: if (sol_times) PetscAssertPointer(sol_times, 3);
6065: if (Sols) PetscAssertPointer(Sols, 4);
6066: if (!ts->eval_times) {
6067: if (nsol) *nsol = 0;
6068: if (sol_times) *sol_times = NULL;
6069: if (Sols) *Sols = NULL;
6070: } else {
6071: if (nsol) *nsol = ts->eval_times->sol_idx;
6072: if (sol_times) *sol_times = ts->eval_times->sol_times;
6073: if (Sols) *Sols = ts->eval_times->sol_vecs;
6074: }
6075: PetscFunctionReturn(PETSC_SUCCESS);
6076: }
6078: /*@
6079: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
6081: Collective
6083: Input Parameters:
6084: + ts - the time-stepper
6085: . n - number of the time points (>=2)
6086: - 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.
6088: Options Database Key:
6089: . -ts_time_span <t0,...tf> - Sets the time span
6091: Level: intermediate
6093: Notes:
6094: 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.
6096: The elements in `span_times` must be all increasing. They correspond to the intermediate points to be saved.
6098: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
6100: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using time span may
6101: pressure the memory system when using a large number of span points.
6103: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6104: @*/
6105: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal span_times[])
6106: {
6107: PetscFunctionBegin;
6109: PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
6110: PetscCall(TSSetEvaluationTimes(ts, n, span_times));
6111: PetscCall(TSSetTime(ts, span_times[0]));
6112: PetscCall(TSSetMaxTime(ts, span_times[n - 1]));
6113: PetscFunctionReturn(PETSC_SUCCESS);
6114: }
6116: /*@
6117: TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
6119: Collective
6121: Input Parameters:
6122: + ts - the `TS` context
6123: . J - Jacobian matrix (not altered in this routine)
6124: - B - newly computed Jacobian matrix to use with preconditioner
6126: Level: intermediate
6128: Notes:
6129: This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
6130: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
6131: and multiple fields are involved.
6133: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
6134: structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
6135: usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
6136: `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.
6138: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
6139: @*/
6140: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
6141: {
6142: MatColoring mc = NULL;
6143: ISColoring iscoloring = NULL;
6144: MatFDColoring matfdcoloring = NULL;
6146: PetscFunctionBegin;
6147: /* Generate new coloring after eliminating zeros in the matrix */
6148: PetscCall(MatEliminateZeros(B, PETSC_TRUE));
6149: PetscCall(MatColoringCreate(B, &mc));
6150: PetscCall(MatColoringSetDistance(mc, 2));
6151: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
6152: PetscCall(MatColoringSetFromOptions(mc));
6153: PetscCall(MatColoringApply(mc, &iscoloring));
6154: PetscCall(MatColoringDestroy(&mc));
6155: /* Replace the old coloring with the new one */
6156: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
6157: PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
6158: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
6159: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
6160: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
6161: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
6162: PetscCall(ISColoringDestroy(&iscoloring));
6163: PetscFunctionReturn(PETSC_SUCCESS);
6164: }