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_time_step <dt> - initial time step (only a suggestion, the actual initial time step used differ)
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_step_rejections <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(PetscOptionsDeprecated("-ts_dt", "-ts_time_step", "3.25", NULL));
127: PetscCall(PetscOptionsReal("-ts_time_step", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg));
128: if (flg) PetscCall(TSSetTimeStep(ts, time_step));
129: PetscCall(PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg));
130: if (flg) PetscCall(TSSetExactFinalTime(ts, eftopt));
131: PetscCall(PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, &flg));
132: if (flg) PetscCall(TSSetMaxSNESFailures(ts, ts->max_snes_failures));
133: PetscCall(PetscOptionsDeprecated("-ts_max_reject", "-ts_max_step_rejections", "3.25", NULL));
134: PetscCall(PetscOptionsInt("-ts_max_step_rejections", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, &flg));
135: if (flg) PetscCall(TSSetMaxStepRejections(ts, ts->max_reject));
136: PetscCall(PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL));
137: PetscCall(PetscOptionsBoundedReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL, 0));
138: PetscCall(PetscOptionsBoundedReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL, 0));
140: PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL));
141: 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));
142: PetscCall(PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL));
143: #if defined(PETSC_HAVE_SAWS)
144: {
145: PetscBool set;
146: flg = PETSC_FALSE;
147: PetscCall(PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set));
148: if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)ts, flg));
149: }
150: #endif
152: /* Monitor options */
153: 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)"));
154: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL));
155: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_wall_clock_time", "Monitor wall-clock time, KSP iterations, and SNES iterations per step", "TSMonitorWallClockTime", TSMonitorWallClockTime, TSMonitorWallClockTimeSetUp));
156: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL));
157: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, TSMonitorSolutionSetup));
158: PetscCall(TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL));
159: PetscCall(PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg));
160: if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)ts, monfilename));
162: PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt));
163: if (opt) {
164: PetscInt howoften = 1;
165: DM dm;
166: PetscBool net;
168: PetscCall(PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL));
169: PetscCall(TSGetDM(ts, &dm));
170: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net));
171: if (net) {
172: TSMonitorLGCtxNetwork ctx;
173: PetscCall(TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx));
174: PetscCall(TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxNetworkDestroy));
175: PetscCall(PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL));
176: } else {
177: TSMonitorLGCtx ctx;
178: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
179: PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
180: }
181: }
183: PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt));
184: if (opt) {
185: TSMonitorLGCtx ctx;
186: PetscInt howoften = 1;
188: PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL));
189: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
190: PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
191: }
192: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL));
194: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt));
195: if (opt) {
196: TSMonitorLGCtx ctx;
197: PetscInt howoften = 1;
199: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
200: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
201: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
202: }
203: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt));
204: if (opt) {
205: TSMonitorLGCtx ctx;
206: PetscInt howoften = 1;
208: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
209: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
210: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
211: ctx->semilogy = PETSC_TRUE;
212: }
214: PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt));
215: if (opt) {
216: TSMonitorLGCtx ctx;
217: PetscInt howoften = 1;
219: PetscCall(PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL));
220: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
221: PetscCall(TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
222: }
223: PetscCall(PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt));
224: if (opt) {
225: TSMonitorLGCtx ctx;
226: PetscInt howoften = 1;
228: PetscCall(PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL));
229: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
230: PetscCall(TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
231: }
232: PetscCall(PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt));
233: if (opt) {
234: TSMonitorSPEigCtx ctx;
235: PetscInt howoften = 1;
237: PetscCall(PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL));
238: PetscCall(TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
239: PetscCall(TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscCtxDestroyFn *)TSMonitorSPEigCtxDestroy));
240: }
241: PetscCall(PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase space from the DMSwarm", "TSMonitorSPSwarm", &opt));
242: if (opt) {
243: TSMonitorSPCtx ctx;
244: PetscInt howoften = 1, retain = 0;
245: PetscBool phase = PETSC_TRUE, create = PETSC_TRUE, multispecies = PETSC_FALSE;
247: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
248: if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
249: create = PETSC_FALSE;
250: break;
251: }
252: if (create) {
253: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase space from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL));
254: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL));
255: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL));
256: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_multi_species", "Color particles by particle species", "TSMonitorSPSwarm", multispecies, &multispecies, NULL));
257: PetscCall(TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, multispecies, &ctx));
258: PetscCall(TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscCtxDestroyFn *)TSMonitorSPCtxDestroy));
259: }
260: }
261: PetscCall(PetscOptionsName("-ts_monitor_hg_swarm", "Display particle histogram from the DMSwarm", "TSMonitorHGSwarm", &opt));
262: if (opt) {
263: TSMonitorHGCtx ctx;
264: PetscInt howoften = 1, Ns = 1;
265: PetscBool velocity = PETSC_FALSE, create = PETSC_TRUE;
267: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
268: if (ts->monitor[i] == TSMonitorHGSwarmSolution) {
269: create = PETSC_FALSE;
270: break;
271: }
272: if (create) {
273: DM sw, dm;
274: PetscInt Nc, Nb;
276: PetscCall(TSGetDM(ts, &sw));
277: PetscCall(DMSwarmGetCellDM(sw, &dm));
278: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc));
279: Nb = PetscMin(20, PetscMax(10, Nc));
280: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL));
281: PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL));
282: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL));
283: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL));
284: PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx));
285: PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscCtxDestroyFn *)TSMonitorHGCtxDestroy));
286: }
287: }
288: opt = PETSC_FALSE;
289: PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt));
290: if (opt) {
291: TSMonitorDrawCtx ctx;
292: PetscInt howoften = 1;
294: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL));
295: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
296: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
297: }
298: opt = PETSC_FALSE;
299: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt));
300: if (opt) {
301: TSMonitorDrawCtx ctx;
302: PetscReal bounds[4];
303: PetscInt n = 4;
304: PetscDraw draw;
305: PetscDrawAxis axis;
307: PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL));
308: PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field");
309: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx));
310: PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw));
311: PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis));
312: PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]));
313: PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2"));
314: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
315: }
316: opt = PETSC_FALSE;
317: PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt));
318: if (opt) {
319: TSMonitorDrawCtx ctx;
320: PetscInt howoften = 1;
322: PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL));
323: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
324: PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
325: }
326: opt = PETSC_FALSE;
327: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt));
328: if (opt) {
329: TSMonitorDrawCtx ctx;
330: PetscInt howoften = 1;
332: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL));
333: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
334: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
335: }
337: opt = PETSC_FALSE;
338: 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));
339: if (flg) {
340: TSMonitorVTKCtx ctx;
342: PetscCall(TSMonitorSolutionVTKCtxCreate(monfilename, &ctx));
343: PetscCall(PetscOptionsInt("-ts_monitor_solution_vtk_interval", "Save every interval time step (-1 for last step only)", NULL, ctx->interval, &ctx->interval, NULL));
344: PetscCall(TSMonitorSet(ts, (PetscErrorCode (*)(TS, PetscInt, PetscReal, Vec, void *))TSMonitorSolutionVTK, ctx, (PetscCtxDestroyFn *)TSMonitorSolutionVTKDestroy));
345: }
347: PetscCall(PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg));
348: if (flg) {
349: TSMonitorDMDARayCtx *rayctx;
350: int ray = 0;
351: DMDirection ddir;
352: DM da;
353: PetscMPIInt rank;
355: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
356: if (dir[0] == 'x') ddir = DM_X;
357: else if (dir[0] == 'y') ddir = DM_Y;
358: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
359: sscanf(dir + 2, "%d", &ray);
361: PetscCall(PetscInfo(ts, "Displaying DMDA ray %c = %d\n", dir[0], ray));
362: PetscCall(PetscNew(&rayctx));
363: PetscCall(TSGetDM(ts, &da));
364: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
365: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank));
366: if (rank == 0) PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer));
367: rayctx->lgctx = NULL;
368: PetscCall(TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy));
369: }
370: PetscCall(PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg));
371: if (flg) {
372: TSMonitorDMDARayCtx *rayctx;
373: int ray = 0;
374: DMDirection ddir;
375: DM da;
376: PetscInt howoften = 1;
378: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
379: if (dir[0] == 'x') ddir = DM_X;
380: else if (dir[0] == 'y') ddir = DM_Y;
381: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
382: sscanf(dir + 2, "%d", &ray);
384: PetscCall(PetscInfo(ts, "Displaying LG DMDA ray %c = %d\n", dir[0], ray));
385: PetscCall(PetscNew(&rayctx));
386: PetscCall(TSGetDM(ts, &da));
387: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
388: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx));
389: PetscCall(TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy));
390: }
392: PetscCall(PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt));
393: if (opt) {
394: TSMonitorEnvelopeCtx ctx;
396: PetscCall(TSMonitorEnvelopeCtxCreate(ts, &ctx));
397: PetscCall(TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscCtxDestroyFn *)TSMonitorEnvelopeCtxDestroy));
398: }
399: flg = PETSC_FALSE;
400: PetscCall(PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt));
401: if (opt && flg) PetscCall(TSMonitorCancel(ts));
403: flg = PETSC_FALSE;
404: PetscCall(PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeIJacobianDefaultColor", flg, &flg, NULL));
405: if (flg) {
406: DM dm;
408: PetscCall(TSGetDM(ts, &dm));
409: PetscCall(DMTSUnsetIJacobianContext_Internal(dm));
410: PetscCall(TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL));
411: PetscCall(PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n"));
412: }
414: /* Handle specific TS options */
415: PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);
417: /* Handle TSAdapt options */
418: PetscCall(TSGetAdapt(ts, &ts->adapt));
419: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
420: PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject));
422: /* TS trajectory must be set after TS, since it may use some TS options above */
423: tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
424: PetscCall(PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL));
425: if (tflg) PetscCall(TSSetSaveTrajectory(ts));
427: PetscCall(TSAdjointSetFromOptions(ts, PetscOptionsObject));
429: /* process any options handlers added with PetscObjectAddOptionsHandler() */
430: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject));
431: PetscOptionsEnd();
433: if (ts->trajectory) PetscCall(TSTrajectorySetFromOptions(ts->trajectory, ts));
435: /* why do we have to do this here and not during TSSetUp? */
436: PetscCall(TSGetSNES(ts, &ts->snes));
437: if (ts->problem_type == TS_LINEAR) {
438: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
439: if (!flg) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
440: }
441: PetscCall(SNESSetFromOptions(ts->snes));
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: /*@
446: TSGetTrajectory - Gets the trajectory from a `TS` if it exists
448: Collective
450: Input Parameter:
451: . ts - the `TS` context obtained from `TSCreate()`
453: Output Parameter:
454: . tr - the `TSTrajectory` object, if it exists
456: Level: advanced
458: Note:
459: This routine should be called after all `TS` options have been set
461: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectoryCreate()`
462: @*/
463: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
464: {
465: PetscFunctionBegin;
467: *tr = ts->trajectory;
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /*@
472: TSSetSaveTrajectory - Causes the `TS` to save its solutions as it iterates forward in time in a `TSTrajectory` object
474: Collective
476: Input Parameter:
477: . ts - the `TS` context obtained from `TSCreate()`
479: Options Database Keys:
480: + -ts_save_trajectory - saves the trajectory to a file
481: - -ts_trajectory_type type - set trajectory type
483: Level: intermediate
485: Notes:
486: This routine should be called after all `TS` options have been set
488: The `TSTRAJECTORYVISUALIZATION` files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
489: MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m
491: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`
492: @*/
493: PetscErrorCode TSSetSaveTrajectory(TS ts)
494: {
495: PetscFunctionBegin;
497: if (!ts->trajectory) PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: /*@
502: TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object
504: Collective
506: Input Parameter:
507: . ts - the `TS` context obtained from `TSCreate()`
509: Level: intermediate
511: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
512: @*/
513: PetscErrorCode TSResetTrajectory(TS ts)
514: {
515: PetscFunctionBegin;
517: if (ts->trajectory) {
518: PetscCall(TSTrajectoryDestroy(&ts->trajectory));
519: PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
520: }
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from a `TS`
527: Collective
529: Input Parameter:
530: . ts - the `TS` context obtained from `TSCreate()`
532: Level: intermediate
534: .seealso: [](ch_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
535: @*/
536: PetscErrorCode TSRemoveTrajectory(TS ts)
537: {
538: PetscFunctionBegin;
540: if (ts->trajectory) PetscCall(TSTrajectoryDestroy(&ts->trajectory));
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /*@
545: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
546: set with `TSSetRHSJacobian()`.
548: Collective
550: Input Parameters:
551: + ts - the `TS` context
552: . t - current timestep
553: - U - input vector
555: Output Parameters:
556: + A - Jacobian matrix
557: - B - optional matrix used to compute the preconditioner, often the same as `A`
559: Level: developer
561: Note:
562: Most users should not need to explicitly call this routine, as it
563: is used internally within the ODE integrators.
565: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()`
566: @*/
567: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
568: {
569: PetscObjectState Ustate;
570: PetscObjectId Uid;
571: DM dm;
572: DMTS tsdm;
573: TSRHSJacobianFn *rhsjacobianfunc;
574: void *ctx;
575: TSRHSFunctionFn *rhsfunction;
577: PetscFunctionBegin;
580: PetscCheckSameComm(ts, 1, U, 3);
581: PetscCall(TSGetDM(ts, &dm));
582: PetscCall(DMGetDMTS(dm, &tsdm));
583: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
584: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx));
585: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
586: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
588: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) PetscFunctionReturn(PETSC_SUCCESS);
590: 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);
591: if (rhsjacobianfunc) {
592: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
593: PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
594: ts->rhsjacs++;
595: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
596: } else {
597: PetscCall(MatZeroEntries(A));
598: if (B && A != B) PetscCall(MatZeroEntries(B));
599: }
600: ts->rhsjacobian.time = t;
601: ts->rhsjacobian.shift = 0;
602: ts->rhsjacobian.scale = 1.;
603: PetscCall(PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid));
604: PetscCall(PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: /*@
609: TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS`
611: Collective
613: Input Parameters:
614: + ts - the `TS` context
615: . t - current time
616: - U - state vector
618: Output Parameter:
619: . y - right-hand side
621: Level: developer
623: Note:
624: Most users should not need to explicitly call this routine, as it
625: is used internally within the nonlinear solvers.
627: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
628: @*/
629: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
630: {
631: TSRHSFunctionFn *rhsfunction;
632: TSIFunctionFn *ifunction;
633: void *ctx;
634: DM dm;
636: PetscFunctionBegin;
640: PetscCall(TSGetDM(ts, &dm));
641: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx));
642: PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
644: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
646: if (rhsfunction) {
647: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, y, 0));
648: PetscCall(VecLockReadPush(U));
649: PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
650: PetscCall(VecLockReadPop(U));
651: ts->rhsfuncs++;
652: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, y, 0));
653: } else PetscCall(VecZeroEntries(y));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: /*@
658: TSComputeSolutionFunction - Evaluates the solution function.
660: Collective
662: Input Parameters:
663: + ts - the `TS` context
664: - t - current time
666: Output Parameter:
667: . U - the solution
669: Level: developer
671: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
672: @*/
673: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
674: {
675: TSSolutionFn *solutionfunction;
676: void *ctx;
677: DM dm;
679: PetscFunctionBegin;
682: PetscCall(TSGetDM(ts, &dm));
683: PetscCall(DMTSGetSolutionFunction(dm, &solutionfunction, &ctx));
684: if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
687: /*@
688: TSComputeForcingFunction - Evaluates the forcing function.
690: Collective
692: Input Parameters:
693: + ts - the `TS` context
694: - t - current time
696: Output Parameter:
697: . U - the function value
699: Level: developer
701: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
702: @*/
703: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
704: {
705: void *ctx;
706: DM dm;
707: TSForcingFn *forcing;
709: PetscFunctionBegin;
712: PetscCall(TSGetDM(ts, &dm));
713: PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));
715: if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
720: {
721: Mat A, B;
722: TSIJacobianFn *ijacobian;
724: PetscFunctionBegin;
725: if (Arhs) *Arhs = NULL;
726: if (Brhs) *Brhs = NULL;
727: PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL));
728: if (Arhs) {
729: if (!ts->Arhs) {
730: if (ijacobian) {
731: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs));
732: PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN));
733: } else {
734: ts->Arhs = A;
735: PetscCall(PetscObjectReference((PetscObject)A));
736: }
737: } else {
738: PetscBool flg;
739: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
740: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
741: if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
742: PetscCall(PetscObjectDereference((PetscObject)ts->Arhs));
743: ts->Arhs = A;
744: PetscCall(PetscObjectReference((PetscObject)A));
745: }
746: }
747: *Arhs = ts->Arhs;
748: }
749: if (Brhs) {
750: if (!ts->Brhs) {
751: if (A != B) {
752: if (ijacobian) {
753: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs));
754: } else {
755: ts->Brhs = B;
756: PetscCall(PetscObjectReference((PetscObject)B));
757: }
758: } else {
759: PetscCall(PetscObjectReference((PetscObject)ts->Arhs));
760: ts->Brhs = ts->Arhs;
761: }
762: }
763: *Brhs = ts->Brhs;
764: }
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: /*@
769: TSComputeIFunction - Evaluates the DAE residual written in the implicit form F(t,U,Udot)=0
771: Collective
773: Input Parameters:
774: + ts - the `TS` context
775: . t - current time
776: . U - state vector
777: . Udot - time derivative of state vector
778: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSFunction should be kept separate
780: Output Parameter:
781: . Y - right-hand side
783: Level: developer
785: Note:
786: Most users should not need to explicitly call this routine, as it
787: is used internally within the nonlinear solvers.
789: If the user did not write their equations in implicit form, this
790: function recasts them in implicit form.
792: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
793: @*/
794: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
795: {
796: TSIFunctionFn *ifunction;
797: TSRHSFunctionFn *rhsfunction;
798: void *ctx;
799: DM dm;
801: PetscFunctionBegin;
807: PetscCall(TSGetDM(ts, &dm));
808: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
809: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
811: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
813: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, Udot, Y));
814: if (ifunction) {
815: PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
816: ts->ifuncs++;
817: }
818: if (imex) {
819: if (!ifunction) PetscCall(VecCopy(Udot, Y));
820: } else if (rhsfunction) {
821: if (ifunction) {
822: Vec Frhs;
824: PetscCall(DMGetGlobalVector(dm, &Frhs));
825: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
826: PetscCall(VecAXPY(Y, -1, Frhs));
827: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
828: } else {
829: PetscCall(TSComputeRHSFunction(ts, t, U, Y));
830: PetscCall(VecAYPX(Y, -1, Udot));
831: }
832: }
833: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, Udot, Y));
834: PetscFunctionReturn(PETSC_SUCCESS);
835: }
837: /*
838: TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call `TSComputeRHSJacobian()` on it.
840: Note:
841: This routine is needed when one switches from `TSComputeIJacobian()` to `TSComputeRHSJacobian()` because the Jacobian matrix may be shifted or scaled in `TSComputeIJacobian()`.
843: */
844: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
845: {
846: PetscFunctionBegin;
848: PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
849: PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");
851: if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
852: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
853: if (B && B == ts->Brhs && A != B) {
854: if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
855: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
856: }
857: ts->rhsjacobian.shift = 0;
858: ts->rhsjacobian.scale = 1.;
859: PetscFunctionReturn(PETSC_SUCCESS);
860: }
862: /*@
863: TSComputeIJacobian - Evaluates the Jacobian of the DAE
865: Collective
867: Input Parameters:
868: + ts - the `TS` context
869: . t - current timestep
870: . U - state vector
871: . Udot - time derivative of state vector
872: . shift - shift to apply, see note below
873: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSJacobian should be kept separate
875: Output Parameters:
876: + A - Jacobian matrix
877: - B - matrix from which the preconditioner is constructed; often the same as `A`
879: Level: developer
881: Notes:
882: If $ F(t,U,\dot{U})=0 $ is the DAE, the required Jacobian is
883: .vb
884: dF/dU + shift*dF/dUdot
885: .ve
886: Most users should not need to explicitly call this routine, as it
887: is used internally within the nonlinear solvers.
889: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
890: @*/
891: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
892: {
893: TSIJacobianFn *ijacobian;
894: TSRHSJacobianFn *rhsjacobian;
895: DM dm;
896: void *ctx;
898: PetscFunctionBegin;
905: PetscCall(TSGetDM(ts, &dm));
906: PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
907: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
909: PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
911: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
912: if (ijacobian) {
913: PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
914: ts->ijacs++;
915: }
916: if (imex) {
917: if (!ijacobian) { /* system was written as Udot = G(t,U) */
918: PetscBool assembled;
919: if (rhsjacobian) {
920: Mat Arhs = NULL;
921: PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
922: if (A == Arhs) {
923: 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 */
924: ts->rhsjacobian.time = PETSC_MIN_REAL;
925: }
926: }
927: PetscCall(MatZeroEntries(A));
928: PetscCall(MatAssembled(A, &assembled));
929: if (!assembled) {
930: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
931: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
932: }
933: PetscCall(MatShift(A, shift));
934: if (A != B) {
935: PetscCall(MatZeroEntries(B));
936: PetscCall(MatAssembled(B, &assembled));
937: if (!assembled) {
938: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
939: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
940: }
941: PetscCall(MatShift(B, shift));
942: }
943: }
944: } else {
945: Mat Arhs = NULL, Brhs = NULL;
947: /* RHSJacobian needs to be converted to part of IJacobian if exists */
948: if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
949: if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
950: PetscObjectState Ustate;
951: PetscObjectId Uid;
952: TSRHSFunctionFn *rhsfunction;
954: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
955: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
956: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
957: if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
958: ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */
959: PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
960: if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
961: } else {
962: PetscBool flg;
964: if (ts->rhsjacobian.reuse) { /* Undo the damage */
965: /* MatScale has a short path for this case.
966: However, this code path is taken the first time TSComputeRHSJacobian is called
967: and the matrices have not been assembled yet */
968: PetscCall(TSRecoverRHSJacobian(ts, A, B));
969: }
970: PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
971: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
972: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
973: if (!flg) {
974: PetscCall(MatScale(A, -1));
975: PetscCall(MatShift(A, shift));
976: }
977: if (A != B) {
978: PetscCall(MatScale(B, -1));
979: PetscCall(MatShift(B, shift));
980: }
981: }
982: ts->rhsjacobian.scale = -1;
983: ts->rhsjacobian.shift = shift;
984: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
985: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
986: PetscCall(MatZeroEntries(A));
987: PetscCall(MatShift(A, shift));
988: if (A != B) {
989: PetscCall(MatZeroEntries(B));
990: PetscCall(MatShift(B, shift));
991: }
992: }
993: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
994: PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
995: if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
996: }
997: }
998: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
999: PetscFunctionReturn(PETSC_SUCCESS);
1000: }
1002: /*@C
1003: TSSetRHSFunction - Sets the routine for evaluating the function,
1004: where U_t = G(t,u).
1006: Logically Collective
1008: Input Parameters:
1009: + ts - the `TS` context obtained from `TSCreate()`
1010: . r - vector to put the computed right-hand side (or `NULL` to have it created)
1011: . f - routine for evaluating the right-hand-side function
1012: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
1014: Level: beginner
1016: Note:
1017: You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.
1019: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1020: @*/
1021: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, void *ctx)
1022: {
1023: SNES snes;
1024: Vec ralloc = NULL;
1025: DM dm;
1027: PetscFunctionBegin;
1031: PetscCall(TSGetDM(ts, &dm));
1032: PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1033: PetscCall(TSGetSNES(ts, &snes));
1034: if (!r && !ts->dm && ts->vec_sol) {
1035: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1036: r = ralloc;
1037: }
1038: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1039: PetscCall(VecDestroy(&ralloc));
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: /*@C
1044: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1046: Logically Collective
1048: Input Parameters:
1049: + ts - the `TS` context obtained from `TSCreate()`
1050: . f - routine for evaluating the solution
1051: - ctx - [optional] user-defined context for private data for the
1052: function evaluation routine (may be `NULL`)
1054: Options Database Keys:
1055: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1056: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
1058: Level: intermediate
1060: Notes:
1061: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1062: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1063: create closed-form solutions with non-physical forcing terms.
1065: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1067: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1068: @*/
1069: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, void *ctx)
1070: {
1071: DM dm;
1073: PetscFunctionBegin;
1075: PetscCall(TSGetDM(ts, &dm));
1076: PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1077: PetscFunctionReturn(PETSC_SUCCESS);
1078: }
1080: /*@C
1081: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1083: Logically Collective
1085: Input Parameters:
1086: + ts - the `TS` context obtained from `TSCreate()`
1087: . func - routine for evaluating the forcing function
1088: - ctx - [optional] user-defined context for private data for the function evaluation routine
1089: (may be `NULL`)
1091: Level: intermediate
1093: Notes:
1094: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1095: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1096: definition of the problem you are solving and hence possibly introducing bugs.
1098: This replaces the ODE F(u,u_t,t) = 0 the `TS` is solving with F(u,u_t,t) - func(t) = 0
1100: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1101: parameters can be passed in the ctx variable.
1103: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1105: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
1106: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1107: @*/
1108: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, void *ctx)
1109: {
1110: DM dm;
1112: PetscFunctionBegin;
1114: PetscCall(TSGetDM(ts, &dm));
1115: PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1116: PetscFunctionReturn(PETSC_SUCCESS);
1117: }
1119: /*@C
1120: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1121: where U_t = G(U,t), as well as the location to store the matrix.
1123: Logically Collective
1125: Input Parameters:
1126: + ts - the `TS` context obtained from `TSCreate()`
1127: . Amat - (approximate) location to store Jacobian matrix entries computed by `f`
1128: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1129: . f - the Jacobian evaluation routine
1130: - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1132: Level: beginner
1134: Notes:
1135: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1137: The `TS` solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f()`
1138: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1140: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
1141: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
1142: @*/
1143: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, void *ctx)
1144: {
1145: SNES snes;
1146: DM dm;
1147: TSIJacobianFn *ijacobian;
1149: PetscFunctionBegin;
1153: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1154: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1156: PetscCall(TSGetDM(ts, &dm));
1157: PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1158: PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1159: PetscCall(TSGetSNES(ts, &snes));
1160: if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1161: if (Amat) {
1162: PetscCall(PetscObjectReference((PetscObject)Amat));
1163: PetscCall(MatDestroy(&ts->Arhs));
1164: ts->Arhs = Amat;
1165: }
1166: if (Pmat) {
1167: PetscCall(PetscObjectReference((PetscObject)Pmat));
1168: PetscCall(MatDestroy(&ts->Brhs));
1169: ts->Brhs = Pmat;
1170: }
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: /*@C
1175: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1177: Logically Collective
1179: Input Parameters:
1180: + ts - the `TS` context obtained from `TSCreate()`
1181: . r - vector to hold the residual (or `NULL` to have it created internally)
1182: . f - the function evaluation routine
1183: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1185: Level: beginner
1187: Note:
1188: The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function.
1190: .seealso: [](ch_ts), `TS`, `TSIFunctionFn`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`,
1191: `TSSetIJacobian()`
1192: @*/
1193: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunctionFn *f, void *ctx)
1194: {
1195: SNES snes;
1196: Vec ralloc = NULL;
1197: DM dm;
1199: PetscFunctionBegin;
1203: PetscCall(TSGetDM(ts, &dm));
1204: PetscCall(DMTSSetIFunction(dm, f, ctx));
1206: PetscCall(TSGetSNES(ts, &snes));
1207: if (!r && !ts->dm && ts->vec_sol) {
1208: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1209: r = ralloc;
1210: }
1211: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1212: PetscCall(VecDestroy(&ralloc));
1213: PetscFunctionReturn(PETSC_SUCCESS);
1214: }
1216: /*@C
1217: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.
1219: Not Collective
1221: Input Parameter:
1222: . ts - the `TS` context
1224: Output Parameters:
1225: + r - vector to hold residual (or `NULL`)
1226: . func - the function to compute residual (or `NULL`)
1227: - ctx - the function context (or `NULL`)
1229: Level: advanced
1231: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1232: @*/
1233: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunctionFn **func, void **ctx)
1234: {
1235: SNES snes;
1236: DM dm;
1238: PetscFunctionBegin;
1240: PetscCall(TSGetSNES(ts, &snes));
1241: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1242: PetscCall(TSGetDM(ts, &dm));
1243: PetscCall(DMTSGetIFunction(dm, func, ctx));
1244: PetscFunctionReturn(PETSC_SUCCESS);
1245: }
1247: /*@C
1248: TSGetRHSFunction - Returns the vector where the right-hand side is stored and the function/context to compute it.
1250: Not Collective
1252: Input Parameter:
1253: . ts - the `TS` context
1255: Output Parameters:
1256: + r - vector to hold computed right-hand side (or `NULL`)
1257: . func - the function to compute right-hand side (or `NULL`)
1258: - ctx - the function context (or `NULL`)
1260: Level: advanced
1262: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1263: @*/
1264: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunctionFn **func, void **ctx)
1265: {
1266: SNES snes;
1267: DM dm;
1269: PetscFunctionBegin;
1271: PetscCall(TSGetSNES(ts, &snes));
1272: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1273: PetscCall(TSGetDM(ts, &dm));
1274: PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1275: PetscFunctionReturn(PETSC_SUCCESS);
1276: }
1278: /*@C
1279: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1280: provided with `TSSetIFunction()`.
1282: Logically Collective
1284: Input Parameters:
1285: + ts - the `TS` context obtained from `TSCreate()`
1286: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
1287: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1288: . f - the Jacobian evaluation routine
1289: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1291: Level: beginner
1293: Notes:
1294: The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1296: If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
1297: space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.
1299: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1300: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1301: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1302: a and vector W depend on the integration method, step size, and past states. For example with
1303: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1304: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1306: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1308: The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f`
1309: You should not assume the values are the same in the next call to `f` as you set them in the previous call.
1311: In case `TSSetRHSJacobian()` is also used in conjunction with a fully-implicit solver,
1312: multilevel linear solvers, e.g. `PCMG`, will likely not work due to the way `TS` handles rhs matrices.
1314: .seealso: [](ch_ts), `TS`, `TSIJacobianFn`, `TSSetIFunction()`, `TSSetRHSJacobian()`,
1315: `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1316: @*/
1317: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobianFn *f, void *ctx)
1318: {
1319: SNES snes;
1320: DM dm;
1322: PetscFunctionBegin;
1326: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1327: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1329: PetscCall(TSGetDM(ts, &dm));
1330: PetscCall(DMTSSetIJacobian(dm, f, ctx));
1332: PetscCall(TSGetSNES(ts, &snes));
1333: PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1334: PetscFunctionReturn(PETSC_SUCCESS);
1335: }
1337: /*@
1338: TSRHSJacobianSetReuse - restore the RHS Jacobian before calling the user-provided `TSRHSJacobianFn` function again
1340: Logically Collective
1342: Input Parameters:
1343: + ts - `TS` context obtained from `TSCreate()`
1344: - reuse - `PETSC_TRUE` if the RHS Jacobian
1346: Level: intermediate
1348: Notes:
1349: Without this flag, `TS` will change the sign and shift the RHS Jacobian for a
1350: finite-time-step implicit solve, in which case the user function will need to recompute the
1351: entire Jacobian. The `reuse `flag must be set if the evaluation function assumes that the
1352: matrix entries have not been changed by the `TS`.
1354: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1355: @*/
1356: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1357: {
1358: PetscFunctionBegin;
1359: ts->rhsjacobian.reuse = reuse;
1360: PetscFunctionReturn(PETSC_SUCCESS);
1361: }
1363: /*@C
1364: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1366: Logically Collective
1368: Input Parameters:
1369: + ts - the `TS` context obtained from `TSCreate()`
1370: . F - vector to hold the residual (or `NULL` to have it created internally)
1371: . fun - the function evaluation routine
1372: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1374: Level: beginner
1376: .seealso: [](ch_ts), `TS`, `TSI2FunctionFn`, `TSSetI2Jacobian()`, `TSSetIFunction()`,
1377: `TSCreate()`, `TSSetRHSFunction()`
1378: @*/
1379: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2FunctionFn *fun, void *ctx)
1380: {
1381: DM dm;
1383: PetscFunctionBegin;
1386: PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1387: PetscCall(TSGetDM(ts, &dm));
1388: PetscCall(DMTSSetI2Function(dm, fun, ctx));
1389: PetscFunctionReturn(PETSC_SUCCESS);
1390: }
1392: /*@C
1393: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.
1395: Not Collective
1397: Input Parameter:
1398: . ts - the `TS` context
1400: Output Parameters:
1401: + r - vector to hold residual (or `NULL`)
1402: . fun - the function to compute residual (or `NULL`)
1403: - ctx - the function context (or `NULL`)
1405: Level: advanced
1407: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1408: @*/
1409: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2FunctionFn **fun, void **ctx)
1410: {
1411: SNES snes;
1412: DM dm;
1414: PetscFunctionBegin;
1416: PetscCall(TSGetSNES(ts, &snes));
1417: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1418: PetscCall(TSGetDM(ts, &dm));
1419: PetscCall(DMTSGetI2Function(dm, fun, ctx));
1420: PetscFunctionReturn(PETSC_SUCCESS);
1421: }
1423: /*@C
1424: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1425: where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.
1427: Logically Collective
1429: Input Parameters:
1430: + ts - the `TS` context obtained from `TSCreate()`
1431: . J - matrix to hold the Jacobian values
1432: . P - matrix for constructing the preconditioner (may be same as `J`)
1433: . jac - the Jacobian evaluation routine, see `TSI2JacobianFn` for the calling sequence
1434: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1436: Level: beginner
1438: Notes:
1439: The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1441: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1442: 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.
1443: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1444: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1446: .seealso: [](ch_ts), `TS`, `TSI2JacobianFn`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1447: @*/
1448: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2JacobianFn *jac, void *ctx)
1449: {
1450: DM dm;
1452: PetscFunctionBegin;
1456: PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1457: PetscCall(TSGetDM(ts, &dm));
1458: PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1459: PetscFunctionReturn(PETSC_SUCCESS);
1460: }
1462: /*@C
1463: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1465: Not Collective, but parallel objects are returned if `TS` is parallel
1467: Input Parameter:
1468: . ts - The `TS` context obtained from `TSCreate()`
1470: Output Parameters:
1471: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1472: . P - The matrix from which the preconditioner is constructed, often the same as `J`
1473: . jac - The function to compute the Jacobian matrices
1474: - ctx - User-defined context for Jacobian evaluation routine
1476: Level: advanced
1478: Note:
1479: You can pass in `NULL` for any return argument you do not need.
1481: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1482: @*/
1483: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2JacobianFn **jac, void **ctx)
1484: {
1485: SNES snes;
1486: DM dm;
1488: PetscFunctionBegin;
1489: PetscCall(TSGetSNES(ts, &snes));
1490: PetscCall(SNESSetUpMatrices(snes));
1491: PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1492: PetscCall(TSGetDM(ts, &dm));
1493: PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1494: PetscFunctionReturn(PETSC_SUCCESS);
1495: }
1497: /*@
1498: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1500: Collective
1502: Input Parameters:
1503: + ts - the `TS` context
1504: . t - current time
1505: . U - state vector
1506: . V - time derivative of state vector (U_t)
1507: - A - second time derivative of state vector (U_tt)
1509: Output Parameter:
1510: . F - the residual vector
1512: Level: developer
1514: Note:
1515: Most users should not need to explicitly call this routine, as it
1516: is used internally within the nonlinear solvers.
1518: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1519: @*/
1520: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1521: {
1522: DM dm;
1523: TSI2FunctionFn *I2Function;
1524: void *ctx;
1525: TSRHSFunctionFn *rhsfunction;
1527: PetscFunctionBegin;
1534: PetscCall(TSGetDM(ts, &dm));
1535: PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1536: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
1538: if (!I2Function) {
1539: PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1540: PetscFunctionReturn(PETSC_SUCCESS);
1541: }
1543: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, V, F));
1545: PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));
1547: if (rhsfunction) {
1548: Vec Frhs;
1550: PetscCall(DMGetGlobalVector(dm, &Frhs));
1551: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1552: PetscCall(VecAXPY(F, -1, Frhs));
1553: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
1554: }
1556: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, V, F));
1557: PetscFunctionReturn(PETSC_SUCCESS);
1558: }
1560: /*@
1561: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1563: Collective
1565: Input Parameters:
1566: + ts - the `TS` context
1567: . t - current timestep
1568: . U - state vector
1569: . V - time derivative of state vector
1570: . A - second time derivative of state vector
1571: . shiftV - shift to apply, see note below
1572: - shiftA - shift to apply, see note below
1574: Output Parameters:
1575: + J - Jacobian matrix
1576: - P - optional matrix used to construct the preconditioner
1578: Level: developer
1580: Notes:
1581: If $F(t,U,V,A) = 0$ is the DAE, the required Jacobian is
1583: $$
1584: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1585: $$
1587: Most users should not need to explicitly call this routine, as it
1588: is used internally within the ODE integrators.
1590: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1591: @*/
1592: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1593: {
1594: DM dm;
1595: TSI2JacobianFn *I2Jacobian;
1596: void *ctx;
1597: TSRHSJacobianFn *rhsjacobian;
1599: PetscFunctionBegin;
1607: PetscCall(TSGetDM(ts, &dm));
1608: PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1609: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1611: if (!I2Jacobian) {
1612: PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, J, P));
1617: PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1618: if (rhsjacobian) {
1619: Mat Jrhs, Prhs;
1620: PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1621: PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1622: PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1623: if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1624: }
1626: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, J, P));
1627: PetscFunctionReturn(PETSC_SUCCESS);
1628: }
1630: /*@C
1631: TSSetTransientVariable - sets function to transform from state to transient variables
1633: Logically Collective
1635: Input Parameters:
1636: + ts - time stepping context on which to change the transient variable
1637: . tvar - a function that transforms to transient variables, see `TSTransientVariableFn` for the calling sequence
1638: - ctx - a context for tvar
1640: Level: advanced
1642: Notes:
1643: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1644: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1645: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1646: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1647: evaluated via the chain rule, as in
1648: .vb
1649: dF/dP + shift * dF/dCdot dC/dP.
1650: .ve
1652: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSTransientVariableFn`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1653: @*/
1654: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariableFn *tvar, void *ctx)
1655: {
1656: DM dm;
1658: PetscFunctionBegin;
1660: PetscCall(TSGetDM(ts, &dm));
1661: PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1662: PetscFunctionReturn(PETSC_SUCCESS);
1663: }
1665: /*@
1666: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1668: Logically Collective
1670: Input Parameters:
1671: + ts - TS on which to compute
1672: - U - state vector to be transformed to transient variables
1674: Output Parameter:
1675: . C - transient (conservative) variable
1677: Level: developer
1679: Developer Notes:
1680: If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed.
1681: This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are
1682: being used.
1684: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1685: @*/
1686: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1687: {
1688: DM dm;
1689: DMTS dmts;
1691: PetscFunctionBegin;
1694: PetscCall(TSGetDM(ts, &dm));
1695: PetscCall(DMGetDMTS(dm, &dmts));
1696: if (dmts->ops->transientvar) {
1698: PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1699: }
1700: PetscFunctionReturn(PETSC_SUCCESS);
1701: }
1703: /*@
1704: TSHasTransientVariable - determine whether transient variables have been set
1706: Logically Collective
1708: Input Parameter:
1709: . ts - `TS` on which to compute
1711: Output Parameter:
1712: . has - `PETSC_TRUE` if transient variables have been set
1714: Level: developer
1716: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1717: @*/
1718: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1719: {
1720: DM dm;
1721: DMTS dmts;
1723: PetscFunctionBegin;
1725: PetscCall(TSGetDM(ts, &dm));
1726: PetscCall(DMGetDMTS(dm, &dmts));
1727: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1728: PetscFunctionReturn(PETSC_SUCCESS);
1729: }
1731: /*@
1732: TS2SetSolution - Sets the initial solution and time derivative vectors
1733: for use by the `TS` routines handling second order equations.
1735: Logically Collective
1737: Input Parameters:
1738: + ts - the `TS` context obtained from `TSCreate()`
1739: . u - the solution vector
1740: - v - the time derivative vector
1742: Level: beginner
1744: .seealso: [](ch_ts), `TS`
1745: @*/
1746: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1747: {
1748: PetscFunctionBegin;
1752: PetscCall(TSSetSolution(ts, u));
1753: PetscCall(PetscObjectReference((PetscObject)v));
1754: PetscCall(VecDestroy(&ts->vec_dot));
1755: ts->vec_dot = v;
1756: PetscFunctionReturn(PETSC_SUCCESS);
1757: }
1759: /*@
1760: TS2GetSolution - Returns the solution and time derivative at the present timestep
1761: for second order equations.
1763: Not Collective
1765: Input Parameter:
1766: . ts - the `TS` context obtained from `TSCreate()`
1768: Output Parameters:
1769: + u - the vector containing the solution
1770: - v - the vector containing the time derivative
1772: Level: intermediate
1774: Notes:
1775: It is valid to call this routine inside the function
1776: that you are evaluating in order to move to the new timestep. This vector not
1777: changed until the solution at the next timestep has been calculated.
1779: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1780: @*/
1781: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1782: {
1783: PetscFunctionBegin;
1785: if (u) PetscAssertPointer(u, 2);
1786: if (v) PetscAssertPointer(v, 3);
1787: if (u) *u = ts->vec_sol;
1788: if (v) *v = ts->vec_dot;
1789: PetscFunctionReturn(PETSC_SUCCESS);
1790: }
1792: /*@
1793: TSLoad - Loads a `TS` that has been stored in binary with `TSView()`.
1795: Collective
1797: Input Parameters:
1798: + ts - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1799: some related function before a call to `TSLoad()`.
1800: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1802: Level: intermediate
1804: Note:
1805: The type is determined by the data in the file, any type set into the `TS` before this call is ignored.
1807: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1808: @*/
1809: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1810: {
1811: PetscBool isbinary;
1812: PetscInt classid;
1813: char type[256];
1814: DMTS sdm;
1815: DM dm;
1817: PetscFunctionBegin;
1820: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1821: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1823: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1824: PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1825: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1826: PetscCall(TSSetType(ts, type));
1827: PetscTryTypeMethod(ts, load, viewer);
1828: PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1829: PetscCall(DMLoad(dm, viewer));
1830: PetscCall(TSSetDM(ts, dm));
1831: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1832: PetscCall(VecLoad(ts->vec_sol, viewer));
1833: PetscCall(DMGetDMTS(ts->dm, &sdm));
1834: PetscCall(DMTSLoad(sdm, viewer));
1835: PetscFunctionReturn(PETSC_SUCCESS);
1836: }
1838: #include <petscdraw.h>
1839: #if defined(PETSC_HAVE_SAWS)
1840: #include <petscviewersaws.h>
1841: #endif
1843: /*@
1844: TSViewFromOptions - View a `TS` based on values in the options database
1846: Collective
1848: Input Parameters:
1849: + ts - the `TS` context
1850: . obj - Optional object that provides the prefix for the options database keys
1851: - name - command line option string to be passed by user
1853: Level: intermediate
1855: .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1856: @*/
1857: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1858: {
1859: PetscFunctionBegin;
1861: PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1862: PetscFunctionReturn(PETSC_SUCCESS);
1863: }
1865: /*@
1866: TSView - Prints the `TS` data structure.
1868: Collective
1870: Input Parameters:
1871: + ts - the `TS` context obtained from `TSCreate()`
1872: - viewer - visualization context
1874: Options Database Key:
1875: . -ts_view - calls `TSView()` at end of `TSStep()`
1877: Level: beginner
1879: Notes:
1880: The available visualization contexts include
1881: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1882: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1883: output where only the first processor opens
1884: the file. All other processors send their
1885: data to the first processor to print.
1887: The user can open an alternative visualization context with
1888: `PetscViewerASCIIOpen()` - output to a specified file.
1890: In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).
1892: The "initial time step" displayed is the default time step from `TSCreate()` or that set with `TSSetTimeStep()` or `-ts_time_step`
1894: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1895: @*/
1896: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1897: {
1898: TSType type;
1899: PetscBool isascii, isstring, issundials, isbinary, isdraw;
1900: DMTS sdm;
1901: #if defined(PETSC_HAVE_SAWS)
1902: PetscBool issaws;
1903: #endif
1905: PetscFunctionBegin;
1907: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1909: PetscCheckSameComm(ts, 1, viewer, 2);
1911: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1912: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1913: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1914: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1915: #if defined(PETSC_HAVE_SAWS)
1916: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1917: #endif
1918: if (isascii) {
1919: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer));
1920: if (ts->ops->view) {
1921: PetscCall(PetscViewerASCIIPushTab(viewer));
1922: PetscUseTypeMethod(ts, view, viewer);
1923: PetscCall(PetscViewerASCIIPopTab(viewer));
1924: }
1925: PetscCall(PetscViewerASCIIPrintf(viewer, " initial time step=%g\n", (double)ts->initial_time_step));
1926: if (ts->max_steps < PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps));
1927: if (ts->run_steps < PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, " run steps=%" PetscInt_FMT "\n", ts->run_steps));
1928: if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time));
1929: if (ts->max_reject != PETSC_UNLIMITED) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum number of step rejections=%" PetscInt_FMT "\n", ts->max_reject));
1930: if (ts->max_snes_failures != PETSC_UNLIMITED) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum number of SNES failures allowed=%" PetscInt_FMT "\n", ts->max_snes_failures));
1931: if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs));
1932: if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs));
1933: if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs));
1934: if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs));
1935: if (ts->usessnes) {
1936: PetscBool lin;
1937: if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its));
1938: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its));
1939: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
1940: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures));
1941: }
1942: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject));
1943: if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, "));
1944: else PetscCall(PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol));
1945: if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, "using vector of absolute error tolerances\n"));
1946: else PetscCall(PetscViewerASCIIPrintf(viewer, "using absolute error tolerance of %g\n", (double)ts->atol));
1947: PetscCall(PetscViewerASCIIPushTab(viewer));
1948: PetscCall(TSAdaptView(ts->adapt, viewer));
1949: PetscCall(PetscViewerASCIIPopTab(viewer));
1950: } else if (isstring) {
1951: PetscCall(TSGetType(ts, &type));
1952: PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type));
1953: PetscTryTypeMethod(ts, view, viewer);
1954: } else if (isbinary) {
1955: PetscInt classid = TS_FILE_CLASSID;
1956: MPI_Comm comm;
1957: PetscMPIInt rank;
1958: char type[256];
1960: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
1961: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1962: if (rank == 0) {
1963: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1964: PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
1965: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1966: }
1967: PetscTryTypeMethod(ts, view, viewer);
1968: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1969: PetscCall(DMView(ts->dm, viewer));
1970: PetscCall(VecView(ts->vec_sol, viewer));
1971: PetscCall(DMGetDMTS(ts->dm, &sdm));
1972: PetscCall(DMTSView(sdm, viewer));
1973: } else if (isdraw) {
1974: PetscDraw draw;
1975: char str[36];
1976: PetscReal x, y, bottom, h;
1978: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1979: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1980: PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
1981: PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
1982: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
1983: bottom = y - h;
1984: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1985: PetscTryTypeMethod(ts, view, viewer);
1986: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1987: if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
1988: PetscCall(PetscDrawPopCurrentPoint(draw));
1989: #if defined(PETSC_HAVE_SAWS)
1990: } else if (issaws) {
1991: PetscMPIInt rank;
1992: const char *name;
1994: PetscCall(PetscObjectGetName((PetscObject)ts, &name));
1995: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1996: if (!((PetscObject)ts)->amsmem && rank == 0) {
1997: char dir[1024];
1999: PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
2000: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
2001: PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
2002: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
2003: PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
2004: }
2005: PetscTryTypeMethod(ts, view, viewer);
2006: #endif
2007: }
2008: if (ts->snes && ts->usessnes) {
2009: PetscCall(PetscViewerASCIIPushTab(viewer));
2010: PetscCall(SNESView(ts->snes, viewer));
2011: PetscCall(PetscViewerASCIIPopTab(viewer));
2012: }
2013: PetscCall(DMGetDMTS(ts->dm, &sdm));
2014: PetscCall(DMTSView(sdm, viewer));
2016: PetscCall(PetscViewerASCIIPushTab(viewer));
2017: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &issundials));
2018: PetscCall(PetscViewerASCIIPopTab(viewer));
2019: PetscFunctionReturn(PETSC_SUCCESS);
2020: }
2022: /*@
2023: TSSetApplicationContext - Sets an optional user-defined context for the timesteppers that may be accessed, for example inside the user provided
2024: `TS` callbacks with `TSGetApplicationContext()`
2026: Logically Collective
2028: Input Parameters:
2029: + ts - the `TS` context obtained from `TSCreate()`
2030: - ctx - user context
2032: Level: intermediate
2034: Fortran Note:
2035: This only works when `ctx` is a Fortran derived type (it cannot be a `PetscObject`), we recommend writing a Fortran interface definition for this
2036: function that tells the Fortran compiler the derived data type that is passed in as the `ctx` argument. See `TSGetApplicationContext()` for
2037: an example.
2039: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2040: @*/
2041: PetscErrorCode TSSetApplicationContext(TS ts, PeCtx ctx)
2042: {
2043: PetscFunctionBegin;
2045: ts->ctx = ctx;
2046: PetscFunctionReturn(PETSC_SUCCESS);
2047: }
2049: /*@
2050: TSGetApplicationContext - Gets the user-defined context for the
2051: timestepper that was set with `TSSetApplicationContext()`
2053: Not Collective
2055: Input Parameter:
2056: . ts - the `TS` context obtained from `TSCreate()`
2058: Output Parameter:
2059: . ctx - a pointer to the user context
2061: Level: intermediate
2063: Fortran Notes:
2064: 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
2065: function that tells the Fortran compiler the derived data type that is returned as the `ctx` argument. For example,
2066: .vb
2067: Interface TSGetApplicationContext
2068: Subroutine TSGetApplicationContext(ts,ctx,ierr)
2069: #include <petsc/finclude/petscts.h>
2070: use petscts
2071: TS ts
2072: type(tUsertype), pointer :: ctx
2073: PetscErrorCode ierr
2074: End Subroutine
2075: End Interface TSGetApplicationContext
2076: .ve
2078: The prototype for `ctx` must be
2079: .vb
2080: type(tUsertype), pointer :: ctx
2081: .ve
2083: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2084: @*/
2085: PetscErrorCode TSGetApplicationContext(TS ts, void *ctx)
2086: {
2087: PetscFunctionBegin;
2089: *(void **)ctx = ts->ctx;
2090: PetscFunctionReturn(PETSC_SUCCESS);
2091: }
2093: /*@
2094: TSGetStepNumber - Gets the number of time steps completed.
2096: Not Collective
2098: Input Parameter:
2099: . ts - the `TS` context obtained from `TSCreate()`
2101: Output Parameter:
2102: . steps - number of steps completed so far
2104: Level: intermediate
2106: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2107: @*/
2108: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2109: {
2110: PetscFunctionBegin;
2112: PetscAssertPointer(steps, 2);
2113: *steps = ts->steps;
2114: PetscFunctionReturn(PETSC_SUCCESS);
2115: }
2117: /*@
2118: TSSetStepNumber - Sets the number of steps completed.
2120: Logically Collective
2122: Input Parameters:
2123: + ts - the `TS` context
2124: - steps - number of steps completed so far
2126: Level: developer
2128: Note:
2129: For most uses of the `TS` solvers the user need not explicitly call
2130: `TSSetStepNumber()`, as the step counter is appropriately updated in
2131: `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2132: reinitialize timestepping by setting the step counter to zero (and time
2133: to the initial time) to solve a similar problem with different initial
2134: conditions or parameters. Other possible use case is to continue
2135: timestepping from a previously interrupted run in such a way that `TS`
2136: monitors will be called with a initial nonzero step counter.
2138: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2139: @*/
2140: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2141: {
2142: PetscFunctionBegin;
2145: PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2146: ts->steps = steps;
2147: PetscFunctionReturn(PETSC_SUCCESS);
2148: }
2150: /*@
2151: TSSetTimeStep - Allows one to reset the timestep at any time.
2153: Logically Collective
2155: Input Parameters:
2156: + ts - the `TS` context obtained from `TSCreate()`
2157: - time_step - the size of the timestep
2159: Options Database Key:
2160: . -ts_time_step <dt> - provide the initial time step
2162: Level: intermediate
2164: Notes:
2165: This is only a suggestion, the actual initial time step used may differ
2167: If this is called after `TSSetUp()`, it will not change the initial time step value printed by `TSView()`
2169: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2170: @*/
2171: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2172: {
2173: PetscFunctionBegin;
2176: ts->time_step = time_step;
2177: if (ts->setupcalled == PETSC_FALSE) ts->initial_time_step = time_step;
2178: PetscFunctionReturn(PETSC_SUCCESS);
2179: }
2181: /*@
2182: TSSetExactFinalTime - Determines whether to adapt the final time step to
2183: match the exact final time, to interpolate the solution to the exact final time,
2184: or to just return at the final time `TS` computed (which may be slightly larger
2185: than the requested final time).
2187: Logically Collective
2189: Input Parameters:
2190: + ts - the time-step context
2191: - eftopt - exact final time option
2192: .vb
2193: TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded, just use it
2194: TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time if the final time is exceeded
2195: TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to ensure the computed final time exactly equals the requested final time
2196: .ve
2198: Options Database Key:
2199: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step approach at runtime
2201: Level: beginner
2203: Note:
2204: If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2205: then the final time you selected.
2207: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2208: @*/
2209: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2210: {
2211: PetscFunctionBegin;
2214: ts->exact_final_time = eftopt;
2215: PetscFunctionReturn(PETSC_SUCCESS);
2216: }
2218: /*@
2219: TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`
2221: Not Collective
2223: Input Parameter:
2224: . ts - the `TS` context
2226: Output Parameter:
2227: . eftopt - exact final time option
2229: Level: beginner
2231: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2232: @*/
2233: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2234: {
2235: PetscFunctionBegin;
2237: PetscAssertPointer(eftopt, 2);
2238: *eftopt = ts->exact_final_time;
2239: PetscFunctionReturn(PETSC_SUCCESS);
2240: }
2242: /*@
2243: TSGetTimeStep - Gets the current timestep size.
2245: Not Collective
2247: Input Parameter:
2248: . ts - the `TS` context obtained from `TSCreate()`
2250: Output Parameter:
2251: . dt - the current timestep size
2253: Level: intermediate
2255: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2256: @*/
2257: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2258: {
2259: PetscFunctionBegin;
2261: PetscAssertPointer(dt, 2);
2262: *dt = ts->time_step;
2263: PetscFunctionReturn(PETSC_SUCCESS);
2264: }
2266: /*@
2267: TSGetSolution - Returns the solution at the present timestep. It
2268: is valid to call this routine inside the function that you are evaluating
2269: in order to move to the new timestep. This vector not changed until
2270: the solution at the next timestep has been calculated.
2272: Not Collective, but v returned is parallel if ts is parallel
2274: Input Parameter:
2275: . ts - the `TS` context obtained from `TSCreate()`
2277: Output Parameter:
2278: . v - the vector containing the solution
2280: Level: intermediate
2282: Note:
2283: If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
2284: final time. It returns the solution at the next timestep.
2286: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2287: @*/
2288: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2289: {
2290: PetscFunctionBegin;
2292: PetscAssertPointer(v, 2);
2293: *v = ts->vec_sol;
2294: PetscFunctionReturn(PETSC_SUCCESS);
2295: }
2297: /*@
2298: TSGetSolutionComponents - Returns any solution components at the present
2299: timestep, if available for the time integration method being used.
2300: Solution components are quantities that share the same size and
2301: structure as the solution vector.
2303: Not Collective, but v returned is parallel if ts is parallel
2305: Input Parameters:
2306: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2307: . n - If v is `NULL`, then the number of solution components is
2308: returned through n, else the n-th solution component is
2309: returned in v.
2310: - v - the vector containing the n-th solution component
2311: (may be `NULL` to use this function to find out
2312: the number of solutions components).
2314: Level: advanced
2316: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2317: @*/
2318: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2319: {
2320: PetscFunctionBegin;
2322: if (!ts->ops->getsolutioncomponents) *n = 0;
2323: else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2324: PetscFunctionReturn(PETSC_SUCCESS);
2325: }
2327: /*@
2328: TSGetAuxSolution - Returns an auxiliary solution at the present
2329: timestep, if available for the time integration method being used.
2331: Not Collective, but v returned is parallel if ts is parallel
2333: Input Parameters:
2334: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2335: - v - the vector containing the auxiliary solution
2337: Level: intermediate
2339: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2340: @*/
2341: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2342: {
2343: PetscFunctionBegin;
2345: if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2346: else PetscCall(VecZeroEntries(*v));
2347: PetscFunctionReturn(PETSC_SUCCESS);
2348: }
2350: /*@
2351: TSGetTimeError - Returns the estimated error vector, if the chosen
2352: `TSType` has an error estimation functionality and `TSSetTimeError()` was called
2354: Not Collective, but v returned is parallel if ts is parallel
2356: Input Parameters:
2357: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2358: . n - current estimate (n=0) or previous one (n=-1)
2359: - v - the vector containing the error (same size as the solution).
2361: Level: intermediate
2363: Note:
2364: MUST call after `TSSetUp()`
2366: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2367: @*/
2368: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2369: {
2370: PetscFunctionBegin;
2372: if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2373: else PetscCall(VecZeroEntries(*v));
2374: PetscFunctionReturn(PETSC_SUCCESS);
2375: }
2377: /*@
2378: TSSetTimeError - Sets the estimated error vector, if the chosen
2379: `TSType` has an error estimation functionality. This can be used
2380: to restart such a time integrator with a given error vector.
2382: Not Collective, but v returned is parallel if ts is parallel
2384: Input Parameters:
2385: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2386: - v - the vector containing the error (same size as the solution).
2388: Level: intermediate
2390: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
2391: @*/
2392: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2393: {
2394: PetscFunctionBegin;
2396: PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2397: PetscTryTypeMethod(ts, settimeerror, v);
2398: PetscFunctionReturn(PETSC_SUCCESS);
2399: }
2401: /* ----- Routines to initialize and destroy a timestepper ---- */
2402: /*@
2403: TSSetProblemType - Sets the type of problem to be solved.
2405: Not collective
2407: Input Parameters:
2408: + ts - The `TS`
2409: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2410: .vb
2411: U_t - A U = 0 (linear)
2412: U_t - A(t) U = 0 (linear)
2413: F(t,U,U_t) = 0 (nonlinear)
2414: .ve
2416: Level: beginner
2418: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2419: @*/
2420: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2421: {
2422: PetscFunctionBegin;
2424: ts->problem_type = type;
2425: if (type == TS_LINEAR) {
2426: SNES snes;
2427: PetscCall(TSGetSNES(ts, &snes));
2428: PetscCall(SNESSetType(snes, SNESKSPONLY));
2429: }
2430: PetscFunctionReturn(PETSC_SUCCESS);
2431: }
2433: /*@
2434: TSGetProblemType - Gets the type of problem to be solved.
2436: Not collective
2438: Input Parameter:
2439: . ts - The `TS`
2441: Output Parameter:
2442: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2443: .vb
2444: M U_t = A U
2445: M(t) U_t = A(t) U
2446: F(t,U,U_t)
2447: .ve
2449: Level: beginner
2451: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2452: @*/
2453: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2454: {
2455: PetscFunctionBegin;
2457: PetscAssertPointer(type, 2);
2458: *type = ts->problem_type;
2459: PetscFunctionReturn(PETSC_SUCCESS);
2460: }
2462: /*
2463: Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2464: */
2465: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2466: {
2467: PetscBool isnone;
2469: PetscFunctionBegin;
2470: PetscCall(TSGetAdapt(ts, &ts->adapt));
2471: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2473: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2474: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2475: else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2476: PetscFunctionReturn(PETSC_SUCCESS);
2477: }
2479: /*@
2480: TSSetUp - Sets up the internal data structures for the later use of a timestepper.
2482: Collective
2484: Input Parameter:
2485: . ts - the `TS` context obtained from `TSCreate()`
2487: Level: advanced
2489: Note:
2490: For basic use of the `TS` solvers the user need not explicitly call
2491: `TSSetUp()`, since these actions will automatically occur during
2492: the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this
2493: phase separately, `TSSetUp()` should be called after `TSCreate()`
2494: and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.
2496: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2497: @*/
2498: PetscErrorCode TSSetUp(TS ts)
2499: {
2500: DM dm;
2501: PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2502: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2503: TSIFunctionFn *ifun;
2504: TSIJacobianFn *ijac;
2505: TSI2JacobianFn *i2jac;
2506: TSRHSJacobianFn *rhsjac;
2508: PetscFunctionBegin;
2510: if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
2512: if (!((PetscObject)ts)->type_name) {
2513: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2514: PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2515: }
2517: if (!ts->vec_sol) {
2518: PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2519: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2520: }
2522: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2523: PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2524: ts->Jacp = ts->Jacprhs;
2525: }
2527: if (ts->quadraturets) {
2528: PetscCall(TSSetUp(ts->quadraturets));
2529: PetscCall(VecDestroy(&ts->vec_costintegrand));
2530: PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2531: }
2533: PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
2534: if (rhsjac == TSComputeRHSJacobianConstant) {
2535: Mat Amat, Pmat;
2536: SNES snes;
2537: PetscCall(TSGetSNES(ts, &snes));
2538: PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
2539: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2540: * have displaced the RHS matrix */
2541: if (Amat && Amat == ts->Arhs) {
2542: /* 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 */
2543: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
2544: PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
2545: PetscCall(MatDestroy(&Amat));
2546: }
2547: if (Pmat && Pmat == ts->Brhs) {
2548: PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
2549: PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
2550: PetscCall(MatDestroy(&Pmat));
2551: }
2552: }
2554: PetscCall(TSGetAdapt(ts, &ts->adapt));
2555: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2557: PetscTryTypeMethod(ts, setup);
2559: PetscCall(TSSetExactFinalTimeDefault(ts));
2561: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2562: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2563: */
2564: PetscCall(TSGetDM(ts, &dm));
2565: PetscCall(DMSNESGetFunction(dm, &func, NULL));
2566: if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));
2568: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2569: Otherwise, the SNES will use coloring internally to form the Jacobian.
2570: */
2571: PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2572: PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2573: PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2574: PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2575: if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));
2577: /* if time integration scheme has a starting method, call it */
2578: PetscTryTypeMethod(ts, startingmethod);
2580: ts->setupcalled = PETSC_TRUE;
2581: PetscFunctionReturn(PETSC_SUCCESS);
2582: }
2584: /*@
2585: 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
2587: Collective
2589: Input Parameter:
2590: . ts - the `TS` context obtained from `TSCreate()`
2592: Level: developer
2594: Notes:
2595: Any options set on the `TS` object, including those set with `TSSetFromOptions()` remain.
2597: See also `TSSetResize()` to change the size of the system being integrated (for example by adaptive mesh refinement) during the time integration.
2599: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSetResize()`
2600: @*/
2601: PetscErrorCode TSReset(TS ts)
2602: {
2603: TS_RHSSplitLink ilink = ts->tsrhssplit, next;
2605: PetscFunctionBegin;
2608: PetscTryTypeMethod(ts, reset);
2609: if (ts->snes) PetscCall(SNESReset(ts->snes));
2610: if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));
2612: PetscCall(MatDestroy(&ts->Arhs));
2613: PetscCall(MatDestroy(&ts->Brhs));
2614: PetscCall(VecDestroy(&ts->Frhs));
2615: PetscCall(VecDestroy(&ts->vec_sol));
2616: PetscCall(VecDestroy(&ts->vec_sol0));
2617: PetscCall(VecDestroy(&ts->vec_dot));
2618: PetscCall(VecDestroy(&ts->vatol));
2619: PetscCall(VecDestroy(&ts->vrtol));
2620: PetscCall(VecDestroyVecs(ts->nwork, &ts->work));
2622: PetscCall(MatDestroy(&ts->Jacprhs));
2623: PetscCall(MatDestroy(&ts->Jacp));
2624: if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2625: if (ts->quadraturets) {
2626: PetscCall(TSReset(ts->quadraturets));
2627: PetscCall(VecDestroy(&ts->vec_costintegrand));
2628: }
2629: while (ilink) {
2630: next = ilink->next;
2631: PetscCall(TSDestroy(&ilink->ts));
2632: PetscCall(PetscFree(ilink->splitname));
2633: PetscCall(ISDestroy(&ilink->is));
2634: PetscCall(PetscFree(ilink));
2635: ilink = next;
2636: }
2637: ts->tsrhssplit = NULL;
2638: ts->num_rhs_splits = 0;
2639: if (ts->eval_times) {
2640: PetscCall(PetscFree(ts->eval_times->time_points));
2641: PetscCall(PetscFree(ts->eval_times->sol_times));
2642: PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
2643: PetscCall(PetscFree(ts->eval_times));
2644: }
2645: ts->rhsjacobian.time = PETSC_MIN_REAL;
2646: ts->rhsjacobian.scale = 1.0;
2647: ts->ijacobian.shift = 1.0;
2648: ts->setupcalled = PETSC_FALSE;
2649: PetscFunctionReturn(PETSC_SUCCESS);
2650: }
2652: static PetscErrorCode TSResizeReset(TS);
2654: /*@
2655: TSDestroy - Destroys the timestepper context that was created
2656: with `TSCreate()`.
2658: Collective
2660: Input Parameter:
2661: . ts - the `TS` context obtained from `TSCreate()`
2663: Level: beginner
2665: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2666: @*/
2667: PetscErrorCode TSDestroy(TS *ts)
2668: {
2669: PetscFunctionBegin;
2670: if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2672: if (--((PetscObject)*ts)->refct > 0) {
2673: *ts = NULL;
2674: PetscFunctionReturn(PETSC_SUCCESS);
2675: }
2677: PetscCall(TSReset(*ts));
2678: PetscCall(TSAdjointReset(*ts));
2679: if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2680: PetscCall(TSResizeReset(*ts));
2682: /* if memory was published with SAWs then destroy it */
2683: PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2684: PetscTryTypeMethod(*ts, destroy);
2686: PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));
2688: PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2689: PetscCall(TSEventDestroy(&(*ts)->event));
2691: PetscCall(SNESDestroy(&(*ts)->snes));
2692: PetscCall(SNESDestroy(&(*ts)->snesrhssplit));
2693: PetscCall(DMDestroy(&(*ts)->dm));
2694: PetscCall(TSMonitorCancel(*ts));
2695: PetscCall(TSAdjointMonitorCancel(*ts));
2697: PetscCall(TSDestroy(&(*ts)->quadraturets));
2698: PetscCall(PetscHeaderDestroy(ts));
2699: PetscFunctionReturn(PETSC_SUCCESS);
2700: }
2702: /*@
2703: TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2704: a `TS` (timestepper) context. Valid only for nonlinear problems.
2706: Not Collective, but snes is parallel if ts is parallel
2708: Input Parameter:
2709: . ts - the `TS` context obtained from `TSCreate()`
2711: Output Parameter:
2712: . snes - the nonlinear solver context
2714: Level: beginner
2716: Notes:
2717: The user can then directly manipulate the `SNES` context to set various
2718: options, etc. Likewise, the user can then extract and manipulate the
2719: `KSP`, and `PC` contexts as well.
2721: `TSGetSNES()` does not work for integrators that do not use `SNES`; in
2722: this case `TSGetSNES()` returns `NULL` in `snes`.
2724: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2725: @*/
2726: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2727: {
2728: PetscFunctionBegin;
2730: PetscAssertPointer(snes, 2);
2731: if (!ts->snes) {
2732: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2733: PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2734: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2735: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2736: if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2737: if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2738: }
2739: *snes = ts->snes;
2740: PetscFunctionReturn(PETSC_SUCCESS);
2741: }
2743: /*@
2744: TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the `TS` timestepping context
2746: Collective
2748: Input Parameters:
2749: + ts - the `TS` context obtained from `TSCreate()`
2750: - snes - the nonlinear solver context
2752: Level: developer
2754: Note:
2755: Most users should have the `TS` created by calling `TSGetSNES()`
2757: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2758: @*/
2759: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2760: {
2761: PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);
2763: PetscFunctionBegin;
2766: PetscCall(PetscObjectReference((PetscObject)snes));
2767: PetscCall(SNESDestroy(&ts->snes));
2769: ts->snes = snes;
2771: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2772: PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2773: if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2774: PetscFunctionReturn(PETSC_SUCCESS);
2775: }
2777: /*@
2778: TSGetKSP - Returns the `KSP` (linear solver) associated with
2779: a `TS` (timestepper) context.
2781: Not Collective, but `ksp` is parallel if `ts` is parallel
2783: Input Parameter:
2784: . ts - the `TS` context obtained from `TSCreate()`
2786: Output Parameter:
2787: . ksp - the nonlinear solver context
2789: Level: beginner
2791: Notes:
2792: The user can then directly manipulate the `KSP` context to set various
2793: options, etc. Likewise, the user can then extract and manipulate the
2794: `PC` context as well.
2796: `TSGetKSP()` does not work for integrators that do not use `KSP`;
2797: in this case `TSGetKSP()` returns `NULL` in `ksp`.
2799: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2800: @*/
2801: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2802: {
2803: SNES snes;
2805: PetscFunctionBegin;
2807: PetscAssertPointer(ksp, 2);
2808: PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2809: PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
2810: PetscCall(TSGetSNES(ts, &snes));
2811: PetscCall(SNESGetKSP(snes, ksp));
2812: PetscFunctionReturn(PETSC_SUCCESS);
2813: }
2815: /* ----------- Routines to set solver parameters ---------- */
2817: /*@
2818: TSSetMaxSteps - Sets the maximum number of steps to use.
2820: Logically Collective
2822: Input Parameters:
2823: + ts - the `TS` context obtained from `TSCreate()`
2824: - maxsteps - maximum number of steps to use
2826: Options Database Key:
2827: . -ts_max_steps <maxsteps> - Sets maxsteps
2829: Level: intermediate
2831: Note:
2832: Use `PETSC_DETERMINE` to reset the maximum number of steps to the default from when the object's type was set
2834: The default maximum number of steps is 5,000
2836: Fortran Note:
2837: Use `PETSC_DETERMINE_INTEGER`
2839: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2840: @*/
2841: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2842: {
2843: PetscFunctionBegin;
2846: if (maxsteps == PETSC_DETERMINE) {
2847: ts->max_steps = ts->default_max_steps;
2848: } else {
2849: PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2850: ts->max_steps = maxsteps;
2851: }
2852: PetscFunctionReturn(PETSC_SUCCESS);
2853: }
2855: /*@
2856: TSGetMaxSteps - Gets the maximum number of steps to use.
2858: Not Collective
2860: Input Parameter:
2861: . ts - the `TS` context obtained from `TSCreate()`
2863: Output Parameter:
2864: . maxsteps - maximum number of steps to use
2866: Level: advanced
2868: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2869: @*/
2870: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2871: {
2872: PetscFunctionBegin;
2874: PetscAssertPointer(maxsteps, 2);
2875: *maxsteps = ts->max_steps;
2876: PetscFunctionReturn(PETSC_SUCCESS);
2877: }
2879: /*@
2880: TSSetRunSteps - Sets the maximum number of steps to take in each call to `TSSolve()`.
2882: If the step count when `TSSolve()` is `start_step`, this will stop the simulation once `current_step - start_step >= run_steps`.
2883: Comparatively, `TSSetMaxSteps()` will stop if `current_step >= max_steps`.
2884: The simulation will stop when either condition is reached.
2886: Logically Collective
2888: Input Parameters:
2889: + ts - the `TS` context obtained from `TSCreate()`
2890: - runsteps - maximum number of steps to take in each call to `TSSolve()`;
2892: Options Database Key:
2893: . -ts_run_steps <runsteps> - Sets runsteps
2895: Level: intermediate
2897: Note:
2898: The default is `PETSC_UNLIMITED`
2900: .seealso: [](ch_ts), `TS`, `TSGetRunSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`, `TSSetMaxSteps()`
2901: @*/
2902: PetscErrorCode TSSetRunSteps(TS ts, PetscInt runsteps)
2903: {
2904: PetscFunctionBegin;
2907: if (runsteps == PETSC_DETERMINE) {
2908: ts->run_steps = PETSC_UNLIMITED;
2909: } else {
2910: 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");
2911: ts->run_steps = runsteps;
2912: }
2913: PetscFunctionReturn(PETSC_SUCCESS);
2914: }
2916: /*@
2917: TSGetRunSteps - Gets the maximum number of steps to take in each call to `TSSolve()`.
2919: Not Collective
2921: Input Parameter:
2922: . ts - the `TS` context obtained from `TSCreate()`
2924: Output Parameter:
2925: . runsteps - maximum number of steps to take in each call to `TSSolve`.
2927: Level: advanced
2929: .seealso: [](ch_ts), `TS`, `TSSetRunSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`, `TSGetMaxSteps()`
2930: @*/
2931: PetscErrorCode TSGetRunSteps(TS ts, PetscInt *runsteps)
2932: {
2933: PetscFunctionBegin;
2935: PetscAssertPointer(runsteps, 2);
2936: *runsteps = ts->run_steps;
2937: PetscFunctionReturn(PETSC_SUCCESS);
2938: }
2940: /*@
2941: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
2943: Logically Collective
2945: Input Parameters:
2946: + ts - the `TS` context obtained from `TSCreate()`
2947: - maxtime - final time to step to
2949: Options Database Key:
2950: . -ts_max_time <maxtime> - Sets maxtime
2952: Level: intermediate
2954: Notes:
2955: Use `PETSC_DETERMINE` to reset the maximum time to the default from when the object's type was set
2957: The default maximum time is 5.0
2959: Fortran Note:
2960: Use `PETSC_DETERMINE_REAL`
2962: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2963: @*/
2964: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2965: {
2966: PetscFunctionBegin;
2969: if (maxtime == PETSC_DETERMINE) {
2970: ts->max_time = ts->default_max_time;
2971: } else {
2972: ts->max_time = maxtime;
2973: }
2974: PetscFunctionReturn(PETSC_SUCCESS);
2975: }
2977: /*@
2978: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
2980: Not Collective
2982: Input Parameter:
2983: . ts - the `TS` context obtained from `TSCreate()`
2985: Output Parameter:
2986: . maxtime - final time to step to
2988: Level: advanced
2990: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2991: @*/
2992: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2993: {
2994: PetscFunctionBegin;
2996: PetscAssertPointer(maxtime, 2);
2997: *maxtime = ts->max_time;
2998: PetscFunctionReturn(PETSC_SUCCESS);
2999: }
3001: // PetscClangLinter pragma disable: -fdoc-*
3002: /*@
3003: TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.
3005: Level: deprecated
3007: @*/
3008: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
3009: {
3010: PetscFunctionBegin;
3012: PetscCall(TSSetTime(ts, initial_time));
3013: PetscCall(TSSetTimeStep(ts, time_step));
3014: PetscFunctionReturn(PETSC_SUCCESS);
3015: }
3017: // PetscClangLinter pragma disable: -fdoc-*
3018: /*@
3019: TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.
3021: Level: deprecated
3023: @*/
3024: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
3025: {
3026: PetscFunctionBegin;
3028: if (maxsteps) {
3029: PetscAssertPointer(maxsteps, 2);
3030: *maxsteps = ts->max_steps;
3031: }
3032: if (maxtime) {
3033: PetscAssertPointer(maxtime, 3);
3034: *maxtime = ts->max_time;
3035: }
3036: PetscFunctionReturn(PETSC_SUCCESS);
3037: }
3039: // PetscClangLinter pragma disable: -fdoc-*
3040: /*@
3041: TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.
3043: Level: deprecated
3045: @*/
3046: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
3047: {
3048: PetscFunctionBegin;
3049: if (maxsteps != PETSC_CURRENT) PetscCall(TSSetMaxSteps(ts, maxsteps));
3050: if (maxtime != (PetscReal)PETSC_CURRENT) PetscCall(TSSetMaxTime(ts, maxtime));
3051: PetscFunctionReturn(PETSC_SUCCESS);
3052: }
3054: // PetscClangLinter pragma disable: -fdoc-*
3055: /*@
3056: TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.
3058: Level: deprecated
3060: @*/
3061: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
3062: {
3063: return TSGetStepNumber(ts, steps);
3064: }
3066: // PetscClangLinter pragma disable: -fdoc-*
3067: /*@
3068: TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.
3070: Level: deprecated
3072: @*/
3073: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
3074: {
3075: return TSGetStepNumber(ts, steps);
3076: }
3078: /*@
3079: TSSetSolution - Sets the initial solution vector
3080: for use by the `TS` routines.
3082: Logically Collective
3084: Input Parameters:
3085: + ts - the `TS` context obtained from `TSCreate()`
3086: - u - the solution vector
3088: Level: beginner
3090: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
3091: @*/
3092: PetscErrorCode TSSetSolution(TS ts, Vec u)
3093: {
3094: DM dm;
3096: PetscFunctionBegin;
3099: PetscCall(PetscObjectReference((PetscObject)u));
3100: PetscCall(VecDestroy(&ts->vec_sol));
3101: ts->vec_sol = u;
3103: PetscCall(TSGetDM(ts, &dm));
3104: PetscCall(DMShellSetGlobalVector(dm, u));
3105: PetscFunctionReturn(PETSC_SUCCESS);
3106: }
3108: /*@C
3109: TSSetPreStep - Sets the general-purpose function
3110: called once at the beginning of each time step.
3112: Logically Collective
3114: Input Parameters:
3115: + ts - The `TS` context obtained from `TSCreate()`
3116: - func - The function
3118: Calling sequence of `func`:
3119: . ts - the `TS` context
3121: Level: intermediate
3123: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3124: @*/
3125: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
3126: {
3127: PetscFunctionBegin;
3129: ts->prestep = func;
3130: PetscFunctionReturn(PETSC_SUCCESS);
3131: }
3133: /*@
3134: TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`
3136: Collective
3138: Input Parameter:
3139: . ts - The `TS` context obtained from `TSCreate()`
3141: Level: developer
3143: Note:
3144: `TSPreStep()` is typically used within time stepping implementations,
3145: so most users would not generally call this routine themselves.
3147: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3148: @*/
3149: PetscErrorCode TSPreStep(TS ts)
3150: {
3151: PetscFunctionBegin;
3153: if (ts->prestep) {
3154: Vec U;
3155: PetscObjectId idprev;
3156: PetscBool sameObject;
3157: PetscObjectState sprev, spost;
3159: PetscCall(TSGetSolution(ts, &U));
3160: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3161: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3162: PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3163: PetscCall(TSGetSolution(ts, &U));
3164: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3165: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3166: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3167: }
3168: PetscFunctionReturn(PETSC_SUCCESS);
3169: }
3171: /*@C
3172: TSSetPreStage - Sets the general-purpose function
3173: called once at the beginning of each stage.
3175: Logically Collective
3177: Input Parameters:
3178: + ts - The `TS` context obtained from `TSCreate()`
3179: - func - The function
3181: Calling sequence of `func`:
3182: + ts - the `TS` context
3183: - stagetime - the stage time
3185: Level: intermediate
3187: Note:
3188: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3189: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3190: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3192: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3193: @*/
3194: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime))
3195: {
3196: PetscFunctionBegin;
3198: ts->prestage = func;
3199: PetscFunctionReturn(PETSC_SUCCESS);
3200: }
3202: /*@C
3203: TSSetPostStage - Sets the general-purpose function
3204: called once at the end of each stage.
3206: Logically Collective
3208: Input Parameters:
3209: + ts - The `TS` context obtained from `TSCreate()`
3210: - func - The function
3212: Calling sequence of `func`:
3213: + ts - the `TS` context
3214: . stagetime - the stage time
3215: . stageindex - the stage index
3216: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3218: Level: intermediate
3220: Note:
3221: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3222: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3223: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3225: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3226: @*/
3227: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
3228: {
3229: PetscFunctionBegin;
3231: ts->poststage = func;
3232: PetscFunctionReturn(PETSC_SUCCESS);
3233: }
3235: /*@C
3236: TSSetPostEvaluate - Sets the general-purpose function
3237: called at the end of each step evaluation.
3239: Logically Collective
3241: Input Parameters:
3242: + ts - The `TS` context obtained from `TSCreate()`
3243: - func - The function
3245: Calling sequence of `func`:
3246: . ts - the `TS` context
3248: Level: intermediate
3250: Note:
3251: The function set by `TSSetPostEvaluate()` is called after the solution is evaluated, or after the step rollback.
3252: Inside the `func` callback, the solution vector can be obtained with `TSGetSolution()`, and modified, if need be.
3253: The time step can be obtained with `TSGetTimeStep()`, and the time at the start of the step - via `TSGetTime()`.
3254: The potential changes to the solution vector introduced by event handling (`postevent()`) are not relevant for `TSSetPostEvaluate()`,
3255: but are relevant for `TSSetPostStep()`, according to the function call scheme in `TSSolve()`, as shown below
3256: .vb
3257: ...
3258: Step()
3259: PostEvaluate()
3260: EventHandling()
3261: step_rollback ? PostEvaluate() : PostStep()
3262: ...
3263: .ve
3264: where EventHandling() may result in one of the following three outcomes
3265: .vb
3266: (1) | successful step | solution intact
3267: (2) | successful step | solution modified by `postevent()`
3268: (3) | step_rollback | solution rolled back
3269: .ve
3271: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3272: @*/
3273: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
3274: {
3275: PetscFunctionBegin;
3277: ts->postevaluate = func;
3278: PetscFunctionReturn(PETSC_SUCCESS);
3279: }
3281: /*@
3282: TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`
3284: Collective
3286: Input Parameters:
3287: + ts - The `TS` context obtained from `TSCreate()`
3288: - stagetime - The absolute time of the current stage
3290: Level: developer
3292: Note:
3293: `TSPreStage()` is typically used within time stepping implementations,
3294: most users would not generally call this routine themselves.
3296: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3297: @*/
3298: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3299: {
3300: PetscFunctionBegin;
3302: if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3303: PetscFunctionReturn(PETSC_SUCCESS);
3304: }
3306: /*@
3307: TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`
3309: Collective
3311: Input Parameters:
3312: + ts - The `TS` context obtained from `TSCreate()`
3313: . stagetime - The absolute time of the current stage
3314: . stageindex - Stage number
3315: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3317: Level: developer
3319: Note:
3320: `TSPostStage()` is typically used within time stepping implementations,
3321: most users would not generally call this routine themselves.
3323: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3324: @*/
3325: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec Y[])
3326: {
3327: PetscFunctionBegin;
3329: if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3330: PetscFunctionReturn(PETSC_SUCCESS);
3331: }
3333: /*@
3334: TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`
3336: Collective
3338: Input Parameter:
3339: . ts - The `TS` context obtained from `TSCreate()`
3341: Level: developer
3343: Note:
3344: `TSPostEvaluate()` is typically used within time stepping implementations,
3345: most users would not generally call this routine themselves.
3347: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3348: @*/
3349: PetscErrorCode TSPostEvaluate(TS ts)
3350: {
3351: PetscFunctionBegin;
3353: if (ts->postevaluate) {
3354: Vec U;
3355: PetscObjectState sprev, spost;
3357: PetscCall(TSGetSolution(ts, &U));
3358: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3359: PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3360: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3361: if (sprev != spost) PetscCall(TSRestartStep(ts));
3362: }
3363: PetscFunctionReturn(PETSC_SUCCESS);
3364: }
3366: /*@C
3367: TSSetPostStep - Sets the general-purpose function
3368: called once at the end of each successful time step.
3370: Logically Collective
3372: Input Parameters:
3373: + ts - The `TS` context obtained from `TSCreate()`
3374: - func - The function
3376: Calling sequence of `func`:
3377: . ts - the `TS` context
3379: Level: intermediate
3381: Note:
3382: The function set by `TSSetPostStep()` is called after each successful step. If the event handler locates an event at the
3383: given step, and `postevent()` modifies the solution vector, the solution vector obtained by `TSGetSolution()` inside `func` will
3384: contain the changes. To get the solution without these changes, use `TSSetPostEvaluate()` to set the appropriate callback.
3385: The scheme of the relevant function calls in `TSSolve()` is shown below
3386: .vb
3387: ...
3388: Step()
3389: PostEvaluate()
3390: EventHandling()
3391: step_rollback ? PostEvaluate() : PostStep()
3392: ...
3393: .ve
3394: where EventHandling() may result in one of the following three outcomes
3395: .vb
3396: (1) | successful step | solution intact
3397: (2) | successful step | solution modified by `postevent()`
3398: (3) | step_rollback | solution rolled back
3399: .ve
3401: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3402: @*/
3403: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
3404: {
3405: PetscFunctionBegin;
3407: ts->poststep = func;
3408: PetscFunctionReturn(PETSC_SUCCESS);
3409: }
3411: /*@
3412: TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`
3414: Collective
3416: Input Parameter:
3417: . ts - The `TS` context obtained from `TSCreate()`
3419: Note:
3420: `TSPostStep()` is typically used within time stepping implementations,
3421: so most users would not generally call this routine themselves.
3423: Level: developer
3425: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPostStep()`
3426: @*/
3427: PetscErrorCode TSPostStep(TS ts)
3428: {
3429: PetscFunctionBegin;
3431: if (ts->poststep) {
3432: Vec U;
3433: PetscObjectId idprev;
3434: PetscBool sameObject;
3435: PetscObjectState sprev, spost;
3437: PetscCall(TSGetSolution(ts, &U));
3438: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3439: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3440: PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3441: PetscCall(TSGetSolution(ts, &U));
3442: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3443: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3444: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3445: }
3446: PetscFunctionReturn(PETSC_SUCCESS);
3447: }
3449: /*@
3450: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3452: Collective
3454: Input Parameters:
3455: + ts - time stepping context
3456: - t - time to interpolate to
3458: Output Parameter:
3459: . U - state at given time
3461: Level: intermediate
3463: Developer Notes:
3464: `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3466: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3467: @*/
3468: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3469: {
3470: PetscFunctionBegin;
3473: 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);
3474: PetscUseTypeMethod(ts, interpolate, t, U);
3475: PetscFunctionReturn(PETSC_SUCCESS);
3476: }
3478: /*@
3479: TSStep - Steps one time step
3481: Collective
3483: Input Parameter:
3484: . ts - the `TS` context obtained from `TSCreate()`
3486: Level: developer
3488: Notes:
3489: The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.
3491: The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
3492: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3494: This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
3495: time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.
3497: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3498: @*/
3499: PetscErrorCode TSStep(TS ts)
3500: {
3501: static PetscBool cite = PETSC_FALSE;
3502: PetscReal ptime;
3504: PetscFunctionBegin;
3506: PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3507: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3508: " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3509: " journal = {arXiv e-preprints},\n"
3510: " eprint = {1806.01437},\n"
3511: " archivePrefix = {arXiv},\n"
3512: " year = {2018}\n}\n",
3513: &cite));
3514: PetscCall(TSSetUp(ts));
3515: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3516: if (ts->eval_times)
3517: ts->eval_times->worktol = 0; /* In each step of TSSolve() 'eval_times->worktol' will be meaningfully defined (later) only once:
3518: in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */
3520: 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>");
3521: 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()");
3522: 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");
3524: if (!ts->vec_sol0) PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sol0));
3525: PetscCall(VecCopy(ts->vec_sol, ts->vec_sol0));
3526: ts->time_step0 = ts->time_step;
3528: if (!ts->steps) ts->ptime_prev = ts->ptime;
3529: ptime = ts->ptime;
3531: ts->ptime_prev_rollback = ts->ptime_prev;
3532: ts->reason = TS_CONVERGED_ITERATING;
3534: PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3535: PetscUseTypeMethod(ts, step);
3536: PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));
3538: if (ts->reason >= 0) {
3539: ts->ptime_prev = ptime;
3540: ts->steps++;
3541: ts->steprollback = PETSC_FALSE;
3542: ts->steprestart = PETSC_FALSE;
3543: ts->stepresize = PETSC_FALSE;
3544: }
3546: if (ts->reason < 0 && ts->errorifstepfailed) {
3547: PetscCall(TSMonitorCancel(ts));
3548: if (ts->usessnes && ts->snes) PetscCall(SNESMonitorCancel(ts->snes));
3549: 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]);
3550: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3551: }
3552: PetscFunctionReturn(PETSC_SUCCESS);
3553: }
3555: /*@
3556: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3557: at the end of a time step with a given order of accuracy.
3559: Collective
3561: Input Parameters:
3562: + ts - time stepping context
3563: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3565: Input/Output Parameter:
3566: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3567: on output, the actual order of the error evaluation
3569: Output Parameter:
3570: . wlte - the weighted local truncation error norm
3572: Level: advanced
3574: Note:
3575: If the timestepper cannot evaluate the error in a particular step
3576: (eg. in the first step or restart steps after event handling),
3577: this routine returns wlte=-1.0 .
3579: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3580: @*/
3581: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3582: {
3583: PetscFunctionBegin;
3587: if (order) PetscAssertPointer(order, 3);
3589: PetscAssertPointer(wlte, 4);
3590: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3591: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3592: PetscFunctionReturn(PETSC_SUCCESS);
3593: }
3595: /*@
3596: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3598: Collective
3600: Input Parameters:
3601: + ts - time stepping context
3602: . order - desired order of accuracy
3603: - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)
3605: Output Parameter:
3606: . U - state at the end of the current step
3608: Level: advanced
3610: Notes:
3611: This function cannot be called until all stages have been evaluated.
3613: 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.
3615: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3616: @*/
3617: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3618: {
3619: PetscFunctionBegin;
3623: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3624: PetscFunctionReturn(PETSC_SUCCESS);
3625: }
3627: /*@C
3628: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3630: Not collective
3632: Input Parameter:
3633: . ts - time stepping context
3635: Output Parameter:
3636: . initCondition - The function which computes an initial condition
3638: Calling sequence of `initCondition`:
3639: + ts - The timestepping context
3640: - u - The input vector in which the initial condition is stored
3642: Level: advanced
3644: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3645: @*/
3646: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3647: {
3648: PetscFunctionBegin;
3650: PetscAssertPointer(initCondition, 2);
3651: *initCondition = ts->ops->initcondition;
3652: PetscFunctionReturn(PETSC_SUCCESS);
3653: }
3655: /*@C
3656: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3658: Logically collective
3660: Input Parameters:
3661: + ts - time stepping context
3662: - initCondition - The function which computes an initial condition
3664: Calling sequence of `initCondition`:
3665: + ts - The timestepping context
3666: - e - The input vector in which the initial condition is to be stored
3668: Level: advanced
3670: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3671: @*/
3672: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3673: {
3674: PetscFunctionBegin;
3677: ts->ops->initcondition = initCondition;
3678: PetscFunctionReturn(PETSC_SUCCESS);
3679: }
3681: /*@
3682: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3684: Collective
3686: Input Parameters:
3687: + ts - time stepping context
3688: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3690: Level: advanced
3692: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3693: @*/
3694: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3695: {
3696: PetscFunctionBegin;
3699: PetscTryTypeMethod(ts, initcondition, u);
3700: PetscFunctionReturn(PETSC_SUCCESS);
3701: }
3703: /*@C
3704: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3706: Not collective
3708: Input Parameter:
3709: . ts - time stepping context
3711: Output Parameter:
3712: . exactError - The function which computes the solution error
3714: Calling sequence of `exactError`:
3715: + ts - The timestepping context
3716: . u - The approximate solution vector
3717: - e - The vector in which the error is stored
3719: Level: advanced
3721: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3722: @*/
3723: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3724: {
3725: PetscFunctionBegin;
3727: PetscAssertPointer(exactError, 2);
3728: *exactError = ts->ops->exacterror;
3729: PetscFunctionReturn(PETSC_SUCCESS);
3730: }
3732: /*@C
3733: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3735: Logically collective
3737: Input Parameters:
3738: + ts - time stepping context
3739: - exactError - The function which computes the solution error
3741: Calling sequence of `exactError`:
3742: + ts - The timestepping context
3743: . u - The approximate solution vector
3744: - e - The vector in which the error is stored
3746: Level: advanced
3748: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3749: @*/
3750: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3751: {
3752: PetscFunctionBegin;
3755: ts->ops->exacterror = exactError;
3756: PetscFunctionReturn(PETSC_SUCCESS);
3757: }
3759: /*@
3760: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3762: Collective
3764: Input Parameters:
3765: + ts - time stepping context
3766: . u - The approximate solution
3767: - e - The `Vec` used to store the error
3769: Level: advanced
3771: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3772: @*/
3773: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3774: {
3775: PetscFunctionBegin;
3779: PetscTryTypeMethod(ts, exacterror, u, e);
3780: PetscFunctionReturn(PETSC_SUCCESS);
3781: }
3783: /*@C
3784: TSSetResize - Sets the resize callbacks.
3786: Logically Collective
3788: Input Parameters:
3789: + ts - The `TS` context obtained from `TSCreate()`
3790: . rollback - Whether a resize will restart the step
3791: . setup - The setup function
3792: . transfer - The transfer function
3793: - ctx - [optional] The user-defined context
3795: Calling sequence of `setup`:
3796: + ts - the `TS` context
3797: . step - the current step
3798: . time - the current time
3799: . state - the current vector of state
3800: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3801: - ctx - user defined context
3803: Calling sequence of `transfer`:
3804: + ts - the `TS` context
3805: . nv - the number of vectors to be transferred
3806: . vecsin - array of vectors to be transferred
3807: . vecsout - array of transferred vectors
3808: - ctx - user defined context
3810: Notes:
3811: The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
3812: depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
3813: and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
3814: that the size of the discrete problem has changed.
3815: In both cases, the solver will collect the needed vectors that will be
3816: transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
3817: current solution vector, and other vectors needed by the specific solver used.
3818: For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3819: Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3820: will be automatically reset if the sizes are changed and they must be specified again by the user
3821: inside the `transfer` function.
3822: The input and output arrays passed to `transfer` are allocated by PETSc.
3823: Vectors in `vecsout` must be created by the user.
3824: Ownership of vectors in `vecsout` is transferred to PETSc.
3826: Level: advanced
3828: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3829: @*/
3830: 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)
3831: {
3832: PetscFunctionBegin;
3834: ts->resizerollback = rollback;
3835: ts->resizesetup = setup;
3836: ts->resizetransfer = transfer;
3837: ts->resizectx = ctx;
3838: PetscFunctionReturn(PETSC_SUCCESS);
3839: }
3841: /*
3842: TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.
3844: Collective
3846: Input Parameters:
3847: + ts - The `TS` context obtained from `TSCreate()`
3848: - 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.
3850: Level: developer
3852: Note:
3853: `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3854: used within time stepping implementations,
3855: so most users would not generally call this routine themselves.
3857: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3858: @*/
3859: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3860: {
3861: PetscFunctionBegin;
3863: PetscTryTypeMethod(ts, resizeregister, flg);
3864: /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3865: PetscFunctionReturn(PETSC_SUCCESS);
3866: }
3868: static PetscErrorCode TSResizeReset(TS ts)
3869: {
3870: PetscFunctionBegin;
3872: PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3873: PetscFunctionReturn(PETSC_SUCCESS);
3874: }
3876: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3877: {
3878: PetscFunctionBegin;
3881: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3882: if (ts->resizetransfer) {
3883: PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3884: PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3885: }
3886: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3887: PetscFunctionReturn(PETSC_SUCCESS);
3888: }
3890: /*@C
3891: TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.
3893: Collective
3895: Input Parameters:
3896: + ts - The `TS` context obtained from `TSCreate()`
3897: . name - A string identifying the vector
3898: - vec - The vector
3900: Level: developer
3902: Note:
3903: `TSResizeRegisterVec()` is typically used within time stepping implementations,
3904: so most users would not generally call this routine themselves.
3906: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3907: @*/
3908: PetscErrorCode TSResizeRegisterVec(TS ts, const char name[], Vec vec)
3909: {
3910: PetscFunctionBegin;
3912: PetscAssertPointer(name, 2);
3914: PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3915: PetscFunctionReturn(PETSC_SUCCESS);
3916: }
3918: /*@C
3919: TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.
3921: Collective
3923: Input Parameters:
3924: + ts - The `TS` context obtained from `TSCreate()`
3925: . name - A string identifying the vector
3926: - vec - The vector
3928: Level: developer
3930: Note:
3931: `TSResizeRetrieveVec()` is typically used within time stepping implementations,
3932: so most users would not generally call this routine themselves.
3934: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
3935: @*/
3936: PetscErrorCode TSResizeRetrieveVec(TS ts, const char name[], Vec *vec)
3937: {
3938: PetscFunctionBegin;
3940: PetscAssertPointer(name, 2);
3941: PetscAssertPointer(vec, 3);
3942: PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
3943: PetscFunctionReturn(PETSC_SUCCESS);
3944: }
3946: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
3947: {
3948: PetscInt cnt;
3949: PetscObjectList tmp;
3950: Vec *vecsin = NULL;
3951: const char **namesin = NULL;
3953: PetscFunctionBegin;
3954: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
3955: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
3956: if (names) PetscCall(PetscMalloc1(cnt, &namesin));
3957: if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
3958: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
3959: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
3960: if (vecs) vecsin[cnt] = (Vec)tmp->obj;
3961: if (names) namesin[cnt] = tmp->name;
3962: cnt++;
3963: }
3964: }
3965: if (nv) *nv = cnt;
3966: if (names) *names = namesin;
3967: if (vecs) *vecs = vecsin;
3968: PetscFunctionReturn(PETSC_SUCCESS);
3969: }
3971: /*@
3972: TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`
3974: Collective
3976: Input Parameter:
3977: . ts - The `TS` context obtained from `TSCreate()`
3979: Level: developer
3981: Note:
3982: `TSResize()` is typically used within time stepping implementations,
3983: so most users would not generally call this routine themselves.
3985: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3986: @*/
3987: PetscErrorCode TSResize(TS ts)
3988: {
3989: PetscInt nv = 0;
3990: const char **names = NULL;
3991: Vec *vecsin = NULL;
3992: const char *solname = "ts:vec_sol";
3994: PetscFunctionBegin;
3996: if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
3997: if (ts->resizesetup) {
3998: PetscCall(VecLockReadPush(ts->vec_sol));
3999: PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &ts->stepresize, ts->resizectx));
4000: PetscCall(VecLockReadPop(ts->vec_sol));
4001: if (ts->stepresize) {
4002: if (ts->resizerollback) {
4003: PetscCall(TSRollBack(ts));
4004: ts->time_step = ts->time_step0;
4005: }
4006: PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
4007: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
4008: }
4009: }
4011: PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
4012: if (nv) {
4013: Vec *vecsout, vecsol;
4015: /* Reset internal objects */
4016: PetscCall(TSReset(ts));
4018: /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
4019: PetscCall(PetscCalloc1(nv, &vecsout));
4020: PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
4021: for (PetscInt i = 0; i < nv; i++) {
4022: const char *name;
4023: char *oname;
4025: PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
4026: PetscCall(PetscStrallocpy(name, &oname));
4027: PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
4028: if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
4029: PetscCall(PetscFree(oname));
4030: PetscCall(VecDestroy(&vecsout[i]));
4031: }
4032: PetscCall(PetscFree(vecsout));
4033: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */
4035: PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
4036: if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
4037: PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
4038: }
4040: PetscCall(PetscFree(names));
4041: PetscCall(PetscFree(vecsin));
4042: PetscCall(TSResizeReset(ts));
4043: PetscFunctionReturn(PETSC_SUCCESS);
4044: }
4046: /*@
4047: TSSolve - Steps the requested number of timesteps.
4049: Collective
4051: Input Parameters:
4052: + ts - the `TS` context obtained from `TSCreate()`
4053: - u - the solution vector (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
4054: otherwise it must contain the initial conditions and will contain the solution at the final requested time
4056: Level: beginner
4058: Notes:
4059: The final time returned by this function may be different from the time of the internally
4060: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
4061: stepped over the final time.
4063: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
4064: @*/
4065: PetscErrorCode TSSolve(TS ts, Vec u)
4066: {
4067: Vec solution;
4069: PetscFunctionBegin;
4073: PetscCall(TSSetExactFinalTimeDefault(ts));
4074: 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 */
4075: if (!ts->vec_sol || u == ts->vec_sol) {
4076: PetscCall(VecDuplicate(u, &solution));
4077: PetscCall(TSSetSolution(ts, solution));
4078: PetscCall(VecDestroy(&solution)); /* grant ownership */
4079: }
4080: PetscCall(VecCopy(u, ts->vec_sol));
4081: PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
4082: } else if (u) PetscCall(TSSetSolution(ts, u));
4083: PetscCall(TSSetUp(ts));
4084: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
4086: 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>");
4087: 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()");
4088: 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");
4089: 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");
4091: if (ts->eval_times) {
4092: if (!ts->eval_times->sol_vecs) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
4093: for (PetscInt i = 0; i < ts->eval_times->num_time_points; i++) {
4094: PetscBool is_close = PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[i], ts->eval_times->reltol * ts->time_step + ts->eval_times->abstol, 0);
4095: if (ts->ptime <= ts->eval_times->time_points[i] || is_close) {
4096: ts->eval_times->time_point_idx = i;
4098: PetscBool is_ptime_in_sol_times = PETSC_FALSE; // If current solution has already been saved, we should not save it again
4099: 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);
4100: if (is_close && !is_ptime_in_sol_times) {
4101: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4102: ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4103: ts->eval_times->sol_idx++;
4104: ts->eval_times->time_point_idx++;
4105: }
4106: break;
4107: }
4108: }
4109: }
4111: if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));
4113: /* reset number of steps only when the step is not restarted. ARKIMEX
4114: restarts the step after an event. Resetting these counters in such case causes
4115: TSTrajectory to incorrectly save the output files
4116: */
4117: /* reset time step and iteration counters */
4118: if (!ts->steps) {
4119: ts->ksp_its = 0;
4120: ts->snes_its = 0;
4121: ts->num_snes_failures = 0;
4122: ts->reject = 0;
4123: ts->steprestart = PETSC_TRUE;
4124: ts->steprollback = PETSC_FALSE;
4125: ts->stepresize = PETSC_FALSE;
4126: ts->rhsjacobian.time = PETSC_MIN_REAL;
4127: }
4129: /* make sure initial time step does not overshoot final time or the next point in evaluation times */
4130: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
4131: PetscReal maxdt;
4132: PetscReal dt = ts->time_step;
4134: if (ts->eval_times) maxdt = ts->eval_times->time_points[ts->eval_times->time_point_idx] - ts->ptime;
4135: else maxdt = ts->max_time - ts->ptime;
4136: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
4137: }
4138: ts->reason = TS_CONVERGED_ITERATING;
4140: {
4141: PetscViewer viewer;
4142: PetscViewerFormat format;
4143: PetscBool flg;
4144: static PetscBool incall = PETSC_FALSE;
4146: if (!incall) {
4147: /* Estimate the convergence rate of the time discretization */
4148: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
4149: if (flg) {
4150: PetscConvEst conv;
4151: DM dm;
4152: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4153: PetscInt Nf;
4154: PetscBool checkTemporal = PETSC_TRUE;
4156: incall = PETSC_TRUE;
4157: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
4158: PetscCall(TSGetDM(ts, &dm));
4159: PetscCall(DMGetNumFields(dm, &Nf));
4160: PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
4161: PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
4162: PetscCall(PetscConvEstUseTS(conv, checkTemporal));
4163: PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
4164: PetscCall(PetscConvEstSetFromOptions(conv));
4165: PetscCall(PetscConvEstSetUp(conv));
4166: PetscCall(PetscConvEstGetConvRate(conv, alpha));
4167: PetscCall(PetscViewerPushFormat(viewer, format));
4168: PetscCall(PetscConvEstRateView(conv, alpha, viewer));
4169: PetscCall(PetscViewerPopFormat(viewer));
4170: PetscCall(PetscViewerDestroy(&viewer));
4171: PetscCall(PetscConvEstDestroy(&conv));
4172: PetscCall(PetscFree(alpha));
4173: incall = PETSC_FALSE;
4174: }
4175: }
4176: }
4178: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));
4180: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4181: PetscUseTypeMethod(ts, solve);
4182: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4183: ts->solvetime = ts->ptime;
4184: solution = ts->vec_sol;
4185: } else { /* Step the requested number of timesteps. */
4186: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4187: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4189: if (!ts->steps) {
4190: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4191: PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4192: }
4194: ts->start_step = ts->steps; // records starting step
4195: while (!ts->reason) {
4196: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4197: if (!ts->steprollback || (ts->stepresize && ts->resizerollback)) PetscCall(TSPreStep(ts));
4198: PetscCall(TSStep(ts));
4199: if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4200: if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4201: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4202: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4203: PetscCall(TSForwardCostIntegral(ts));
4204: if (ts->reason >= 0) ts->steps++;
4205: }
4206: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4207: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4208: PetscCall(TSForwardStep(ts));
4209: if (ts->reason >= 0) ts->steps++;
4210: }
4211: PetscCall(TSPostEvaluate(ts));
4212: 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. */
4213: if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4214: if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
4215: /* check convergence */
4216: if (!ts->reason) {
4217: if ((ts->steps - ts->start_step) >= ts->run_steps) ts->reason = TS_CONVERGED_ITS;
4218: else if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4219: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4220: }
4221: if (!ts->steprollback) {
4222: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4223: PetscCall(TSPostStep(ts));
4224: if (!ts->resizerollback) PetscCall(TSResize(ts));
4226: if (ts->eval_times && ts->eval_times->time_point_idx < ts->eval_times->num_time_points && ts->reason >= 0) {
4227: PetscCheck(ts->eval_times->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(eval_times->worktol > 0) in TSSolve()");
4228: if (PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[ts->eval_times->time_point_idx], ts->eval_times->worktol, 0)) {
4229: ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4230: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4231: ts->eval_times->sol_idx++;
4232: ts->eval_times->time_point_idx++;
4233: }
4234: }
4235: }
4236: }
4237: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4239: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4240: if (!u) u = ts->vec_sol;
4241: PetscCall(TSInterpolate(ts, ts->max_time, u));
4242: ts->solvetime = ts->max_time;
4243: solution = u;
4244: PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4245: } else {
4246: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4247: ts->solvetime = ts->ptime;
4248: solution = ts->vec_sol;
4249: }
4250: }
4252: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4253: PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4254: PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4255: if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4256: PetscFunctionReturn(PETSC_SUCCESS);
4257: }
4259: /*@
4260: TSGetTime - Gets the time of the most recently completed step.
4262: Not Collective
4264: Input Parameter:
4265: . ts - the `TS` context obtained from `TSCreate()`
4267: Output Parameter:
4268: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
4270: Level: beginner
4272: Note:
4273: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
4274: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
4276: .seealso: [](ch_ts), `TS`, `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4277: @*/
4278: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4279: {
4280: PetscFunctionBegin;
4282: PetscAssertPointer(t, 2);
4283: *t = ts->ptime;
4284: PetscFunctionReturn(PETSC_SUCCESS);
4285: }
4287: /*@
4288: TSGetPrevTime - Gets the starting time of the previously completed step.
4290: Not Collective
4292: Input Parameter:
4293: . ts - the `TS` context obtained from `TSCreate()`
4295: Output Parameter:
4296: . t - the previous time
4298: Level: beginner
4300: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4301: @*/
4302: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4303: {
4304: PetscFunctionBegin;
4306: PetscAssertPointer(t, 2);
4307: *t = ts->ptime_prev;
4308: PetscFunctionReturn(PETSC_SUCCESS);
4309: }
4311: /*@
4312: TSSetTime - Allows one to reset the time.
4314: Logically Collective
4316: Input Parameters:
4317: + ts - the `TS` context obtained from `TSCreate()`
4318: - t - the time
4320: Level: intermediate
4322: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4323: @*/
4324: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4325: {
4326: PetscFunctionBegin;
4329: ts->ptime = t;
4330: PetscFunctionReturn(PETSC_SUCCESS);
4331: }
4333: /*@
4334: TSSetOptionsPrefix - Sets the prefix used for searching for all
4335: TS options in the database.
4337: Logically Collective
4339: Input Parameters:
4340: + ts - The `TS` context
4341: - prefix - The prefix to prepend to all option names
4343: Level: advanced
4345: Note:
4346: A hyphen (-) must NOT be given at the beginning of the prefix name.
4347: The first character of all runtime options is AUTOMATICALLY the
4348: hyphen.
4350: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4351: @*/
4352: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4353: {
4354: SNES snes;
4356: PetscFunctionBegin;
4358: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4359: PetscCall(TSGetSNES(ts, &snes));
4360: PetscCall(SNESSetOptionsPrefix(snes, prefix));
4361: PetscFunctionReturn(PETSC_SUCCESS);
4362: }
4364: /*@
4365: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4366: TS options in the database.
4368: Logically Collective
4370: Input Parameters:
4371: + ts - The `TS` context
4372: - prefix - The prefix to prepend to all option names
4374: Level: advanced
4376: Note:
4377: A hyphen (-) must NOT be given at the beginning of the prefix name.
4378: The first character of all runtime options is AUTOMATICALLY the
4379: hyphen.
4381: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4382: @*/
4383: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4384: {
4385: SNES snes;
4387: PetscFunctionBegin;
4389: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4390: PetscCall(TSGetSNES(ts, &snes));
4391: PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4392: PetscFunctionReturn(PETSC_SUCCESS);
4393: }
4395: /*@
4396: TSGetOptionsPrefix - Sets the prefix used for searching for all
4397: `TS` options in the database.
4399: Not Collective
4401: Input Parameter:
4402: . ts - The `TS` context
4404: Output Parameter:
4405: . prefix - A pointer to the prefix string used
4407: Level: intermediate
4409: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4410: @*/
4411: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4412: {
4413: PetscFunctionBegin;
4415: PetscAssertPointer(prefix, 2);
4416: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4417: PetscFunctionReturn(PETSC_SUCCESS);
4418: }
4420: /*@C
4421: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4423: Not Collective, but parallel objects are returned if ts is parallel
4425: Input Parameter:
4426: . ts - The `TS` context obtained from `TSCreate()`
4428: Output Parameters:
4429: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`)
4430: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`)
4431: . func - Function to compute the Jacobian of the RHS (or `NULL`)
4432: - ctx - User-defined context for Jacobian evaluation routine (or `NULL`)
4434: Level: intermediate
4436: Note:
4437: You can pass in `NULL` for any return argument you do not need.
4439: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4441: @*/
4442: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, void **ctx)
4443: {
4444: DM dm;
4446: PetscFunctionBegin;
4447: if (Amat || Pmat) {
4448: SNES snes;
4449: PetscCall(TSGetSNES(ts, &snes));
4450: PetscCall(SNESSetUpMatrices(snes));
4451: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4452: }
4453: PetscCall(TSGetDM(ts, &dm));
4454: PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4455: PetscFunctionReturn(PETSC_SUCCESS);
4456: }
4458: /*@C
4459: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4461: Not Collective, but parallel objects are returned if ts is parallel
4463: Input Parameter:
4464: . ts - The `TS` context obtained from `TSCreate()`
4466: Output Parameters:
4467: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4468: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4469: . f - The function to compute the matrices
4470: - ctx - User-defined context for Jacobian evaluation routine
4472: Level: advanced
4474: Note:
4475: You can pass in `NULL` for any return argument you do not need.
4477: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4478: @*/
4479: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, void **ctx)
4480: {
4481: DM dm;
4483: PetscFunctionBegin;
4484: if (Amat || Pmat) {
4485: SNES snes;
4486: PetscCall(TSGetSNES(ts, &snes));
4487: PetscCall(SNESSetUpMatrices(snes));
4488: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4489: }
4490: PetscCall(TSGetDM(ts, &dm));
4491: PetscCall(DMTSGetIJacobian(dm, f, ctx));
4492: PetscFunctionReturn(PETSC_SUCCESS);
4493: }
4495: #include <petsc/private/dmimpl.h>
4496: /*@
4497: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4499: Logically Collective
4501: Input Parameters:
4502: + ts - the `TS` integrator object
4503: - dm - the dm, cannot be `NULL`
4505: Level: intermediate
4507: Notes:
4508: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4509: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4510: different problems using the same function space.
4512: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4513: @*/
4514: PetscErrorCode TSSetDM(TS ts, DM dm)
4515: {
4516: SNES snes;
4517: DMTS tsdm;
4519: PetscFunctionBegin;
4522: PetscCall(PetscObjectReference((PetscObject)dm));
4523: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4524: if (ts->dm->dmts && !dm->dmts) {
4525: PetscCall(DMCopyDMTS(ts->dm, dm));
4526: PetscCall(DMGetDMTS(ts->dm, &tsdm));
4527: /* Grant write privileges to the replacement DM */
4528: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4529: }
4530: PetscCall(DMDestroy(&ts->dm));
4531: }
4532: ts->dm = dm;
4534: PetscCall(TSGetSNES(ts, &snes));
4535: PetscCall(SNESSetDM(snes, dm));
4536: PetscFunctionReturn(PETSC_SUCCESS);
4537: }
4539: /*@
4540: TSGetDM - Gets the `DM` that may be used by some preconditioners
4542: Not Collective
4544: Input Parameter:
4545: . ts - the `TS`
4547: Output Parameter:
4548: . dm - the `DM`
4550: Level: intermediate
4552: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4553: @*/
4554: PetscErrorCode TSGetDM(TS ts, DM *dm)
4555: {
4556: PetscFunctionBegin;
4558: if (!ts->dm) {
4559: PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4560: if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4561: }
4562: *dm = ts->dm;
4563: PetscFunctionReturn(PETSC_SUCCESS);
4564: }
4566: /*@
4567: SNESTSFormFunction - Function to evaluate nonlinear residual defined by an ODE solver algorithm implemented within `TS`
4569: Logically Collective
4571: Input Parameters:
4572: + snes - nonlinear solver
4573: . U - the current state at which to evaluate the residual
4574: - ctx - user context, must be a `TS`
4576: Output Parameter:
4577: . F - the nonlinear residual
4579: Level: developer
4581: Note:
4582: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4583: It is most frequently passed to `MatFDColoringSetFunction()`.
4585: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4586: @*/
4587: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4588: {
4589: TS ts = (TS)ctx;
4591: PetscFunctionBegin;
4596: PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4597: PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4598: PetscFunctionReturn(PETSC_SUCCESS);
4599: }
4601: /*@
4602: SNESTSFormJacobian - Function to evaluate the Jacobian defined by an ODE solver algorithm implemented within `TS`
4604: Collective
4606: Input Parameters:
4607: + snes - nonlinear solver
4608: . U - the current state at which to evaluate the residual
4609: - ctx - user context, must be a `TS`
4611: Output Parameters:
4612: + A - the Jacobian
4613: - B - the matrix used to construct the preconditioner (often the same as `A`)
4615: Level: developer
4617: Note:
4618: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4620: .seealso: [](ch_ts), `SNESSetJacobian()`
4621: @*/
4622: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4623: {
4624: TS ts = (TS)ctx;
4626: PetscFunctionBegin;
4632: PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4633: PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4634: PetscFunctionReturn(PETSC_SUCCESS);
4635: }
4637: /*@C
4638: TSComputeRHSFunctionLinear - Evaluate the right-hand side via the user-provided Jacobian, for linear problems Udot = A U only
4640: Collective
4642: Input Parameters:
4643: + ts - time stepping context
4644: . t - time at which to evaluate
4645: . U - state at which to evaluate
4646: - ctx - context
4648: Output Parameter:
4649: . F - right-hand side
4651: Level: intermediate
4653: Note:
4654: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right-hand side for linear problems.
4655: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4657: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4658: @*/
4659: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4660: {
4661: Mat Arhs, Brhs;
4663: PetscFunctionBegin;
4664: PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4665: /* undo the damage caused by shifting */
4666: PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4667: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4668: PetscCall(MatMult(Arhs, U, F));
4669: PetscFunctionReturn(PETSC_SUCCESS);
4670: }
4672: /*@C
4673: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
4675: Collective
4677: Input Parameters:
4678: + ts - time stepping context
4679: . t - time at which to evaluate
4680: . U - state at which to evaluate
4681: - ctx - context
4683: Output Parameters:
4684: + A - Jacobian
4685: - B - matrix used to construct the preconditioner, often the same as `A`
4687: Level: intermediate
4689: Note:
4690: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4692: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4693: @*/
4694: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4695: {
4696: PetscFunctionBegin;
4697: PetscFunctionReturn(PETSC_SUCCESS);
4698: }
4700: /*@C
4701: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4703: Collective
4705: Input Parameters:
4706: + ts - time stepping context
4707: . t - time at which to evaluate
4708: . U - state at which to evaluate
4709: . Udot - time derivative of state vector
4710: - ctx - context
4712: Output Parameter:
4713: . F - left hand side
4715: Level: intermediate
4717: Notes:
4718: 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
4719: user is required to write their own `TSComputeIFunction()`.
4720: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4721: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4723: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4725: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4726: @*/
4727: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4728: {
4729: Mat A, B;
4731: PetscFunctionBegin;
4732: PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4733: PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4734: PetscCall(MatMult(A, Udot, F));
4735: PetscFunctionReturn(PETSC_SUCCESS);
4736: }
4738: /*@C
4739: TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE
4741: Collective
4743: Input Parameters:
4744: + ts - time stepping context
4745: . t - time at which to evaluate
4746: . U - state at which to evaluate
4747: . Udot - time derivative of state vector
4748: . shift - shift to apply
4749: - ctx - context
4751: Output Parameters:
4752: + A - pointer to operator
4753: - B - pointer to matrix from which the preconditioner is built (often `A`)
4755: Level: advanced
4757: Notes:
4758: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4760: It is only appropriate for problems of the form
4762: $$
4763: M \dot{U} = F(U,t)
4764: $$
4766: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4767: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4768: an implicit operator of the form
4770: $$
4771: shift*M + J
4772: $$
4774: 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
4775: a copy of M or reassemble it when requested.
4777: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4778: @*/
4779: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4780: {
4781: PetscFunctionBegin;
4782: PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4783: ts->ijacobian.shift = shift;
4784: PetscFunctionReturn(PETSC_SUCCESS);
4785: }
4787: /*@
4788: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4790: Not Collective
4792: Input Parameter:
4793: . ts - the `TS` context
4795: Output Parameter:
4796: . equation_type - see `TSEquationType`
4798: Level: beginner
4800: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4801: @*/
4802: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4803: {
4804: PetscFunctionBegin;
4806: PetscAssertPointer(equation_type, 2);
4807: *equation_type = ts->equation_type;
4808: PetscFunctionReturn(PETSC_SUCCESS);
4809: }
4811: /*@
4812: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4814: Not Collective
4816: Input Parameters:
4817: + ts - the `TS` context
4818: - equation_type - see `TSEquationType`
4820: Level: advanced
4822: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4823: @*/
4824: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4825: {
4826: PetscFunctionBegin;
4828: ts->equation_type = equation_type;
4829: PetscFunctionReturn(PETSC_SUCCESS);
4830: }
4832: /*@
4833: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4835: Not Collective
4837: Input Parameter:
4838: . ts - the `TS` context
4840: Output Parameter:
4841: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4842: manual pages for the individual convergence tests for complete lists
4844: Level: beginner
4846: Note:
4847: Can only be called after the call to `TSSolve()` is complete.
4849: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4850: @*/
4851: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4852: {
4853: PetscFunctionBegin;
4855: PetscAssertPointer(reason, 2);
4856: *reason = ts->reason;
4857: PetscFunctionReturn(PETSC_SUCCESS);
4858: }
4860: /*@
4861: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4863: Logically Collective; reason must contain common value
4865: Input Parameters:
4866: + ts - the `TS` context
4867: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4868: manual pages for the individual convergence tests for complete lists
4870: Level: advanced
4872: Note:
4873: Can only be called while `TSSolve()` is active.
4875: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4876: @*/
4877: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4878: {
4879: PetscFunctionBegin;
4881: ts->reason = reason;
4882: PetscFunctionReturn(PETSC_SUCCESS);
4883: }
4885: /*@
4886: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4888: Not Collective
4890: Input Parameter:
4891: . ts - the `TS` context
4893: Output Parameter:
4894: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4896: Level: beginner
4898: Note:
4899: Can only be called after the call to `TSSolve()` is complete.
4901: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4902: @*/
4903: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4904: {
4905: PetscFunctionBegin;
4907: PetscAssertPointer(ftime, 2);
4908: *ftime = ts->solvetime;
4909: PetscFunctionReturn(PETSC_SUCCESS);
4910: }
4912: /*@
4913: TSGetSNESIterations - Gets the total number of nonlinear iterations
4914: used by the time integrator.
4916: Not Collective
4918: Input Parameter:
4919: . ts - `TS` context
4921: Output Parameter:
4922: . nits - number of nonlinear iterations
4924: Level: intermediate
4926: Note:
4927: This counter is reset to zero for each successive call to `TSSolve()`.
4929: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4930: @*/
4931: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4932: {
4933: PetscFunctionBegin;
4935: PetscAssertPointer(nits, 2);
4936: *nits = ts->snes_its;
4937: PetscFunctionReturn(PETSC_SUCCESS);
4938: }
4940: /*@
4941: TSGetKSPIterations - Gets the total number of linear iterations
4942: used by the time integrator.
4944: Not Collective
4946: Input Parameter:
4947: . ts - `TS` context
4949: Output Parameter:
4950: . lits - number of linear iterations
4952: Level: intermediate
4954: Note:
4955: This counter is reset to zero for each successive call to `TSSolve()`.
4957: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`
4958: @*/
4959: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4960: {
4961: PetscFunctionBegin;
4963: PetscAssertPointer(lits, 2);
4964: *lits = ts->ksp_its;
4965: PetscFunctionReturn(PETSC_SUCCESS);
4966: }
4968: /*@
4969: TSGetStepRejections - Gets the total number of rejected steps.
4971: Not Collective
4973: Input Parameter:
4974: . ts - `TS` context
4976: Output Parameter:
4977: . rejects - number of steps rejected
4979: Level: intermediate
4981: Note:
4982: This counter is reset to zero for each successive call to `TSSolve()`.
4984: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4985: @*/
4986: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4987: {
4988: PetscFunctionBegin;
4990: PetscAssertPointer(rejects, 2);
4991: *rejects = ts->reject;
4992: PetscFunctionReturn(PETSC_SUCCESS);
4993: }
4995: /*@
4996: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
4998: Not Collective
5000: Input Parameter:
5001: . ts - `TS` context
5003: Output Parameter:
5004: . fails - number of failed nonlinear solves
5006: Level: intermediate
5008: Note:
5009: This counter is reset to zero for each successive call to `TSSolve()`.
5011: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
5012: @*/
5013: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
5014: {
5015: PetscFunctionBegin;
5017: PetscAssertPointer(fails, 2);
5018: *fails = ts->num_snes_failures;
5019: PetscFunctionReturn(PETSC_SUCCESS);
5020: }
5022: /*@
5023: TSSetMaxStepRejections - Sets the maximum number of step rejections allowed in a single time-step attempt before a time step fails in `TSSolve()` with `TS_DIVERGED_STEP_REJECTED`
5025: Not Collective
5027: Input Parameters:
5028: + ts - `TS` context
5029: - rejects - maximum number of rejected steps, pass `PETSC_UNLIMITED` for unlimited
5031: Options Database Key:
5032: . -ts_max_step_rejections - Maximum number of step rejections before a step fails
5034: Level: intermediate
5036: Developer Note:
5037: The options database name is incorrect.
5039: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`,
5040: `TSGetConvergedReason()`, `TSSolve()`, `TS_DIVERGED_STEP_REJECTED`
5041: @*/
5042: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
5043: {
5044: PetscFunctionBegin;
5046: if (rejects == PETSC_UNLIMITED || rejects == -1) {
5047: ts->max_reject = PETSC_UNLIMITED;
5048: } else {
5049: PetscCheck(rejects >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of rejections");
5050: ts->max_reject = rejects;
5051: }
5052: PetscFunctionReturn(PETSC_SUCCESS);
5053: }
5055: /*@
5056: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves allowed before `TSSolve()` is ended with a `TSConvergedReason` of `TS_DIVERGED_NONLINEAR_SOLVE`
5058: Not Collective
5060: Input Parameters:
5061: + ts - `TS` context
5062: - fails - maximum number of failed nonlinear solves, pass `PETSC_UNLIMITED` to allow any number of failures.
5064: Options Database Key:
5065: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
5067: Level: intermediate
5069: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`,
5070: `TSGetConvergedReason()`, `TS_DIVERGED_NONLINEAR_SOLVE`, `TSConvergedReason`
5071: @*/
5072: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
5073: {
5074: PetscFunctionBegin;
5076: if (fails == PETSC_UNLIMITED || fails == -1) {
5077: ts->max_snes_failures = PETSC_UNLIMITED;
5078: } else {
5079: PetscCheck(fails >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures");
5080: ts->max_snes_failures = fails;
5081: }
5082: PetscFunctionReturn(PETSC_SUCCESS);
5083: }
5085: /*@
5086: TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`
5088: Not Collective
5090: Input Parameters:
5091: + ts - `TS` context
5092: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
5094: Options Database Key:
5095: . -ts_error_if_step_fails - Error if no step succeeds
5097: Level: intermediate
5099: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
5100: @*/
5101: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
5102: {
5103: PetscFunctionBegin;
5105: ts->errorifstepfailed = err;
5106: PetscFunctionReturn(PETSC_SUCCESS);
5107: }
5109: /*@
5110: TSGetAdapt - Get the adaptive controller context for the current method
5112: Collective if controller has not yet been created
5114: Input Parameter:
5115: . ts - time stepping context
5117: Output Parameter:
5118: . adapt - adaptive controller
5120: Level: intermediate
5122: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
5123: @*/
5124: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
5125: {
5126: PetscFunctionBegin;
5128: PetscAssertPointer(adapt, 2);
5129: if (!ts->adapt) {
5130: PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
5131: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
5132: }
5133: *adapt = ts->adapt;
5134: PetscFunctionReturn(PETSC_SUCCESS);
5135: }
5137: /*@
5138: TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller
5140: Logically Collective
5142: Input Parameters:
5143: + ts - time integration context
5144: . atol - scalar absolute tolerances
5145: . vatol - vector of absolute tolerances or `NULL`, used in preference to `atol` if present
5146: . rtol - scalar relative tolerances
5147: - vrtol - vector of relative tolerances or `NULL`, used in preference to `rtol` if present
5149: Options Database Keys:
5150: + -ts_rtol <rtol> - relative tolerance for local truncation error
5151: - -ts_atol <atol> - Absolute tolerance for local truncation error
5153: Level: beginner
5155: Notes:
5156: `PETSC_CURRENT` or `PETSC_DETERMINE` may be used for `atol` or `rtol` to indicate the current value
5157: or the default value from when the object's type was set.
5159: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5160: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5161: computed only for the differential or the algebraic part then this can be done using the vector of
5162: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
5163: differential part and infinity for the algebraic part, the LTE calculation will include only the
5164: differential variables.
5166: Fortran Note:
5167: Use `PETSC_CURRENT_INTEGER` or `PETSC_DETERMINE_INTEGER`.
5169: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
5170: @*/
5171: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
5172: {
5173: PetscFunctionBegin;
5174: if (atol == (PetscReal)PETSC_DETERMINE) {
5175: ts->atol = ts->default_atol;
5176: } else if (atol != (PetscReal)PETSC_CURRENT) {
5177: PetscCheck(atol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)atol);
5178: ts->atol = atol;
5179: }
5181: if (vatol) {
5182: PetscCall(PetscObjectReference((PetscObject)vatol));
5183: PetscCall(VecDestroy(&ts->vatol));
5184: ts->vatol = vatol;
5185: }
5187: if (rtol == (PetscReal)PETSC_DETERMINE) {
5188: ts->rtol = ts->default_rtol;
5189: } else if (rtol != (PetscReal)PETSC_CURRENT) {
5190: PetscCheck(rtol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative", (double)rtol);
5191: ts->rtol = rtol;
5192: }
5194: if (vrtol) {
5195: PetscCall(PetscObjectReference((PetscObject)vrtol));
5196: PetscCall(VecDestroy(&ts->vrtol));
5197: ts->vrtol = vrtol;
5198: }
5199: PetscFunctionReturn(PETSC_SUCCESS);
5200: }
5202: /*@
5203: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
5205: Logically Collective
5207: Input Parameter:
5208: . ts - time integration context
5210: Output Parameters:
5211: + atol - scalar absolute tolerances, `NULL` to ignore
5212: . vatol - vector of absolute tolerances, `NULL` to ignore
5213: . rtol - scalar relative tolerances, `NULL` to ignore
5214: - vrtol - vector of relative tolerances, `NULL` to ignore
5216: Level: beginner
5218: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
5219: @*/
5220: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
5221: {
5222: PetscFunctionBegin;
5223: if (atol) *atol = ts->atol;
5224: if (vatol) *vatol = ts->vatol;
5225: if (rtol) *rtol = ts->rtol;
5226: if (vrtol) *vrtol = ts->vrtol;
5227: PetscFunctionReturn(PETSC_SUCCESS);
5228: }
5230: /*@
5231: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
5233: Collective
5235: Input Parameters:
5236: + ts - time stepping context
5237: . U - state vector, usually ts->vec_sol
5238: . Y - state vector to be compared to U
5239: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5241: Output Parameters:
5242: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5243: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5244: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5246: Options Database Key:
5247: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5249: Level: developer
5251: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5252: @*/
5253: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5254: {
5255: PetscInt norma_loc, norm_loc, normr_loc;
5257: PetscFunctionBegin;
5262: PetscAssertPointer(norm, 5);
5263: PetscAssertPointer(norma, 6);
5264: PetscAssertPointer(normr, 7);
5265: 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));
5266: if (wnormtype == NORM_2) {
5267: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5268: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5269: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5270: }
5271: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5272: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5273: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5274: PetscFunctionReturn(PETSC_SUCCESS);
5275: }
5277: /*@
5278: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5280: Collective
5282: Input Parameters:
5283: + ts - time stepping context
5284: . E - error vector
5285: . U - state vector, usually ts->vec_sol
5286: . Y - state vector, previous time step
5287: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5289: Output Parameters:
5290: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5291: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5292: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5294: Options Database Key:
5295: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5297: Level: developer
5299: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5300: @*/
5301: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5302: {
5303: PetscInt norma_loc, norm_loc, normr_loc;
5305: PetscFunctionBegin;
5307: 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));
5308: if (wnormtype == NORM_2) {
5309: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5310: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5311: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5312: }
5313: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5314: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5315: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5316: PetscFunctionReturn(PETSC_SUCCESS);
5317: }
5319: /*@
5320: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5322: Logically Collective
5324: Input Parameters:
5325: + ts - time stepping context
5326: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5328: Note:
5329: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5331: Level: intermediate
5333: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5334: @*/
5335: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5336: {
5337: PetscFunctionBegin;
5339: ts->cfltime_local = cfltime;
5340: ts->cfltime = -1.;
5341: PetscFunctionReturn(PETSC_SUCCESS);
5342: }
5344: /*@
5345: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5347: Collective
5349: Input Parameter:
5350: . ts - time stepping context
5352: Output Parameter:
5353: . cfltime - maximum stable time step for forward Euler
5355: Level: advanced
5357: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5358: @*/
5359: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5360: {
5361: PetscFunctionBegin;
5362: if (ts->cfltime < 0) PetscCallMPI(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5363: *cfltime = ts->cfltime;
5364: PetscFunctionReturn(PETSC_SUCCESS);
5365: }
5367: /*@
5368: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5370: Input Parameters:
5371: + ts - the `TS` context.
5372: . xl - lower bound.
5373: - xu - upper bound.
5375: Level: advanced
5377: Note:
5378: If this routine is not called then the lower and upper bounds are set to
5379: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5381: .seealso: [](ch_ts), `TS`
5382: @*/
5383: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5384: {
5385: SNES snes;
5387: PetscFunctionBegin;
5388: PetscCall(TSGetSNES(ts, &snes));
5389: PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5390: PetscFunctionReturn(PETSC_SUCCESS);
5391: }
5393: /*@
5394: TSComputeLinearStability - computes the linear stability function at a point
5396: Collective
5398: Input Parameters:
5399: + ts - the `TS` context
5400: . xr - real part of input argument
5401: - xi - imaginary part of input argument
5403: Output Parameters:
5404: + yr - real part of function value
5405: - yi - imaginary part of function value
5407: Level: developer
5409: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5410: @*/
5411: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5412: {
5413: PetscFunctionBegin;
5415: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5416: PetscFunctionReturn(PETSC_SUCCESS);
5417: }
5419: /*@
5420: TSRestartStep - Flags the solver to restart the next step
5422: Collective
5424: Input Parameter:
5425: . ts - the `TS` context obtained from `TSCreate()`
5427: Level: advanced
5429: Notes:
5430: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5431: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5432: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5433: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5434: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5435: discontinuous source terms).
5437: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5438: @*/
5439: PetscErrorCode TSRestartStep(TS ts)
5440: {
5441: PetscFunctionBegin;
5443: ts->steprestart = PETSC_TRUE;
5444: PetscFunctionReturn(PETSC_SUCCESS);
5445: }
5447: /*@
5448: TSRollBack - Rolls back one time step
5450: Collective
5452: Input Parameter:
5453: . ts - the `TS` context obtained from `TSCreate()`
5455: Level: advanced
5457: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5458: @*/
5459: PetscErrorCode TSRollBack(TS ts)
5460: {
5461: PetscFunctionBegin;
5463: PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5464: PetscTryTypeMethod(ts, rollback);
5465: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
5466: ts->time_step = ts->ptime - ts->ptime_prev;
5467: ts->ptime = ts->ptime_prev;
5468: ts->ptime_prev = ts->ptime_prev_rollback;
5469: ts->steps--;
5470: ts->steprollback = PETSC_TRUE;
5471: PetscFunctionReturn(PETSC_SUCCESS);
5472: }
5474: /*@
5475: TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step
5477: Not collective
5479: Input Parameter:
5480: . ts - the `TS` context obtained from `TSCreate()`
5482: Output Parameter:
5483: . flg - the rollback flag
5485: Level: advanced
5487: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
5488: @*/
5489: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
5490: {
5491: PetscFunctionBegin;
5493: PetscAssertPointer(flg, 2);
5494: *flg = ts->steprollback;
5495: PetscFunctionReturn(PETSC_SUCCESS);
5496: }
5498: /*@
5499: TSGetStepResize - Get the internal flag indicating if the current step is after a resize.
5501: Not collective
5503: Input Parameter:
5504: . ts - the `TS` context obtained from `TSCreate()`
5506: Output Parameter:
5507: . flg - the resize flag
5509: Level: advanced
5511: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetResize()`
5512: @*/
5513: PetscErrorCode TSGetStepResize(TS ts, PetscBool *flg)
5514: {
5515: PetscFunctionBegin;
5517: PetscAssertPointer(flg, 2);
5518: *flg = ts->stepresize;
5519: PetscFunctionReturn(PETSC_SUCCESS);
5520: }
5522: /*@
5523: TSGetStages - Get the number of stages and stage values
5525: Input Parameter:
5526: . ts - the `TS` context obtained from `TSCreate()`
5528: Output Parameters:
5529: + ns - the number of stages
5530: - Y - the current stage vectors
5532: Level: advanced
5534: Note:
5535: Both `ns` and `Y` can be `NULL`.
5537: .seealso: [](ch_ts), `TS`, `TSCreate()`
5538: @*/
5539: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5540: {
5541: PetscFunctionBegin;
5543: if (ns) PetscAssertPointer(ns, 2);
5544: if (Y) PetscAssertPointer(Y, 3);
5545: if (!ts->ops->getstages) {
5546: if (ns) *ns = 0;
5547: if (Y) *Y = NULL;
5548: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5549: PetscFunctionReturn(PETSC_SUCCESS);
5550: }
5552: /*@C
5553: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5555: Collective
5557: Input Parameters:
5558: + ts - the `TS` context
5559: . t - current timestep
5560: . U - state vector
5561: . Udot - time derivative of state vector
5562: . shift - shift to apply, see note below
5563: - ctx - an optional user context
5565: Output Parameters:
5566: + J - Jacobian matrix (not altered in this routine)
5567: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
5569: Level: intermediate
5571: Notes:
5572: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5574: dF/dU + shift*dF/dUdot
5576: Most users should not need to explicitly call this routine, as it
5577: is used internally within the nonlinear solvers.
5579: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5580: routine, then it will try to get the coloring from the matrix. This requires that the
5581: matrix have nonzero entries precomputed.
5583: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5584: @*/
5585: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5586: {
5587: SNES snes;
5588: MatFDColoring color;
5589: PetscBool hascolor, matcolor = PETSC_FALSE;
5591: PetscFunctionBegin;
5592: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5593: PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5594: if (!color) {
5595: DM dm;
5596: ISColoring iscoloring;
5598: PetscCall(TSGetDM(ts, &dm));
5599: PetscCall(DMHasColoring(dm, &hascolor));
5600: if (hascolor && !matcolor) {
5601: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5602: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5603: PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5604: PetscCall(MatFDColoringSetFromOptions(color));
5605: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5606: PetscCall(ISColoringDestroy(&iscoloring));
5607: } else {
5608: MatColoring mc;
5610: PetscCall(MatColoringCreate(B, &mc));
5611: PetscCall(MatColoringSetDistance(mc, 2));
5612: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5613: PetscCall(MatColoringSetFromOptions(mc));
5614: PetscCall(MatColoringApply(mc, &iscoloring));
5615: PetscCall(MatColoringDestroy(&mc));
5616: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5617: PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5618: PetscCall(MatFDColoringSetFromOptions(color));
5619: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5620: PetscCall(ISColoringDestroy(&iscoloring));
5621: }
5622: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5623: PetscCall(PetscObjectDereference((PetscObject)color));
5624: }
5625: PetscCall(TSGetSNES(ts, &snes));
5626: PetscCall(MatFDColoringApply(B, color, U, snes));
5627: if (J != B) {
5628: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5629: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5630: }
5631: PetscFunctionReturn(PETSC_SUCCESS);
5632: }
5634: /*@C
5635: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5637: Logically collective
5639: Input Parameters:
5640: + ts - the `TS` context
5641: - func - function called within `TSFunctionDomainError()`
5643: Calling sequence of `func`:
5644: + ts - the `TS` context
5645: . time - the current time (of the stage)
5646: . state - the state to check if it is valid
5647: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable
5649: Level: intermediate
5651: Notes:
5652: `accept` must be collectively specified.
5653: If an implicit ODE solver is being used then, in addition to providing this routine, the
5654: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5655: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5656: Use `TSGetSNES()` to obtain the `SNES` object
5658: Developer Notes:
5659: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5660: since one takes a function pointer and the other does not.
5662: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5663: @*/
5664: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5665: {
5666: PetscFunctionBegin;
5668: ts->functiondomainerror = func;
5669: PetscFunctionReturn(PETSC_SUCCESS);
5670: }
5672: /*@
5673: TSFunctionDomainError - Checks if the current state is valid
5675: Collective
5677: Input Parameters:
5678: + ts - the `TS` context
5679: . stagetime - time of the simulation
5680: - Y - state vector to check.
5682: Output Parameter:
5683: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5685: Level: developer
5687: Note:
5688: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5689: to check if the current state is valid.
5691: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5692: @*/
5693: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5694: {
5695: PetscFunctionBegin;
5699: PetscAssertPointer(accept, 4);
5700: *accept = PETSC_TRUE;
5701: if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5702: PetscFunctionReturn(PETSC_SUCCESS);
5703: }
5705: /*@
5706: TSClone - This function clones a time step `TS` object.
5708: Collective
5710: Input Parameter:
5711: . tsin - The input `TS`
5713: Output Parameter:
5714: . tsout - The output `TS` (cloned)
5716: Level: developer
5718: Notes:
5719: 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.
5720: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5722: 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
5723: .vb
5724: SNES snes_dup = NULL;
5725: TSGetSNES(ts,&snes_dup);
5726: TSSetSNES(ts,snes_dup);
5727: .ve
5729: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5730: @*/
5731: PetscErrorCode TSClone(TS tsin, TS *tsout)
5732: {
5733: TS t;
5734: SNES snes_start;
5735: DM dm;
5736: TSType type;
5738: PetscFunctionBegin;
5739: PetscAssertPointer(tsin, 1);
5740: *tsout = NULL;
5742: PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));
5744: /* General TS description */
5745: t->numbermonitors = 0;
5746: t->setupcalled = PETSC_FALSE;
5747: t->ksp_its = 0;
5748: t->snes_its = 0;
5749: t->nwork = 0;
5750: t->rhsjacobian.time = PETSC_MIN_REAL;
5751: t->rhsjacobian.scale = 1.;
5752: t->ijacobian.shift = 1.;
5754: PetscCall(TSGetSNES(tsin, &snes_start));
5755: PetscCall(TSSetSNES(t, snes_start));
5757: PetscCall(TSGetDM(tsin, &dm));
5758: PetscCall(TSSetDM(t, dm));
5760: t->adapt = tsin->adapt;
5761: PetscCall(PetscObjectReference((PetscObject)t->adapt));
5763: t->trajectory = tsin->trajectory;
5764: PetscCall(PetscObjectReference((PetscObject)t->trajectory));
5766: t->event = tsin->event;
5767: if (t->event) t->event->refct++;
5769: t->problem_type = tsin->problem_type;
5770: t->ptime = tsin->ptime;
5771: t->ptime_prev = tsin->ptime_prev;
5772: t->time_step = tsin->time_step;
5773: t->max_time = tsin->max_time;
5774: t->steps = tsin->steps;
5775: t->max_steps = tsin->max_steps;
5776: t->equation_type = tsin->equation_type;
5777: t->atol = tsin->atol;
5778: t->rtol = tsin->rtol;
5779: t->max_snes_failures = tsin->max_snes_failures;
5780: t->max_reject = tsin->max_reject;
5781: t->errorifstepfailed = tsin->errorifstepfailed;
5783: PetscCall(TSGetType(tsin, &type));
5784: PetscCall(TSSetType(t, type));
5786: t->vec_sol = NULL;
5788: t->cfltime = tsin->cfltime;
5789: t->cfltime_local = tsin->cfltime_local;
5790: t->exact_final_time = tsin->exact_final_time;
5792: t->ops[0] = tsin->ops[0];
5794: if (((PetscObject)tsin)->fortran_func_pointers) {
5795: PetscInt i;
5796: PetscCall(PetscMalloc((10) * sizeof(PetscFortranCallbackFn *), &((PetscObject)t)->fortran_func_pointers));
5797: for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5798: }
5799: *tsout = t;
5800: PetscFunctionReturn(PETSC_SUCCESS);
5801: }
5803: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5804: {
5805: TS ts = (TS)ctx;
5807: PetscFunctionBegin;
5808: PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5809: PetscFunctionReturn(PETSC_SUCCESS);
5810: }
5812: /*@
5813: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5815: Logically Collective
5817: Input Parameter:
5818: . ts - the time stepping routine
5820: Output Parameter:
5821: . flg - `PETSC_TRUE` if the multiply is likely correct
5823: Options Database Key:
5824: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5826: Level: advanced
5828: Note:
5829: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5831: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5832: @*/
5833: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5834: {
5835: Mat J, B;
5836: TSRHSJacobianFn *func;
5837: void *ctx;
5839: PetscFunctionBegin;
5840: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5841: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5842: PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5843: PetscFunctionReturn(PETSC_SUCCESS);
5844: }
5846: /*@
5847: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5849: Logically Collective
5851: Input Parameter:
5852: . ts - the time stepping routine
5854: Output Parameter:
5855: . flg - `PETSC_TRUE` if the multiply is likely correct
5857: Options Database Key:
5858: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5860: Level: advanced
5862: Notes:
5863: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5865: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5866: @*/
5867: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5868: {
5869: Mat J, B;
5870: void *ctx;
5871: TSRHSJacobianFn *func;
5873: PetscFunctionBegin;
5874: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5875: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5876: PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5877: PetscFunctionReturn(PETSC_SUCCESS);
5878: }
5880: /*@
5881: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
5883: Logically Collective
5885: Input Parameters:
5886: + ts - timestepping context
5887: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5889: Options Database Key:
5890: . -ts_use_splitrhsfunction - <true,false>
5892: Level: intermediate
5894: Note:
5895: This is only for multirate methods
5897: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5898: @*/
5899: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5900: {
5901: PetscFunctionBegin;
5903: ts->use_splitrhsfunction = use_splitrhsfunction;
5904: PetscFunctionReturn(PETSC_SUCCESS);
5905: }
5907: /*@
5908: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
5910: Not Collective
5912: Input Parameter:
5913: . ts - timestepping context
5915: Output Parameter:
5916: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5918: Level: intermediate
5920: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5921: @*/
5922: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5923: {
5924: PetscFunctionBegin;
5926: *use_splitrhsfunction = ts->use_splitrhsfunction;
5927: PetscFunctionReturn(PETSC_SUCCESS);
5928: }
5930: /*@
5931: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
5933: Logically Collective
5935: Input Parameters:
5936: + ts - the time-stepper
5937: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
5939: Level: intermediate
5941: Note:
5942: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
5944: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
5945: @*/
5946: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5947: {
5948: PetscFunctionBegin;
5950: ts->axpy_pattern = str;
5951: PetscFunctionReturn(PETSC_SUCCESS);
5952: }
5954: /*@
5955: TSSetEvaluationTimes - sets the evaluation points. The solution will be computed and stored for each time requested
5957: Collective
5959: Input Parameters:
5960: + ts - the time-stepper
5961: . n - number of the time points
5962: - time_points - array of the time points, must be increasing
5964: Options Database Key:
5965: . -ts_eval_times <t0,...tn> - Sets the evaluation times
5967: Level: intermediate
5969: Notes:
5970: The elements in `time_points` must be all increasing. They correspond to the intermediate points to be saved.
5972: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
5974: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using evaluation times may
5975: pressure the memory system when using a large number of time points.
5977: .seealso: [](ch_ts), `TS`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`, `TSSetTimeSpan()`
5978: @*/
5979: PetscErrorCode TSSetEvaluationTimes(TS ts, PetscInt n, PetscReal time_points[])
5980: {
5981: PetscBool is_sorted;
5983: PetscFunctionBegin;
5985: if (ts->eval_times) { // Reset eval_times
5986: ts->eval_times->sol_idx = 0;
5987: ts->eval_times->time_point_idx = 0;
5988: if (n != ts->eval_times->num_time_points) {
5989: PetscCall(PetscFree(ts->eval_times->time_points));
5990: PetscCall(PetscFree(ts->eval_times->sol_times));
5991: PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
5992: } else {
5993: PetscCall(PetscArrayzero(ts->eval_times->sol_times, n));
5994: for (PetscInt i = 0; i < n; i++) PetscCall(VecZeroEntries(ts->eval_times->sol_vecs[i]));
5995: }
5996: } else { // Create/initialize eval_times
5997: TSEvaluationTimes eval_times;
5998: PetscCall(PetscNew(&eval_times));
5999: PetscCall(PetscMalloc1(n, &eval_times->time_points));
6000: PetscCall(PetscMalloc1(n, &eval_times->sol_times));
6001: eval_times->reltol = 1e-6;
6002: eval_times->abstol = 10 * PETSC_MACHINE_EPSILON;
6003: eval_times->worktol = 0;
6004: ts->eval_times = eval_times;
6005: }
6006: ts->eval_times->num_time_points = n;
6007: PetscCall(PetscSortedReal(n, time_points, &is_sorted));
6008: PetscCheck(is_sorted, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "time_points array must be sorted");
6009: PetscCall(PetscArraycpy(ts->eval_times->time_points, time_points, n));
6010: // Note: ts->vec_sol not guaranteed to exist, so ts->eval_times->sol_vecs allocated at TSSolve time
6011: PetscFunctionReturn(PETSC_SUCCESS);
6012: }
6014: /*@C
6015: TSGetEvaluationTimes - gets the evaluation times set with `TSSetEvaluationTimes()`
6017: Not Collective
6019: Input Parameter:
6020: . ts - the time-stepper
6022: Output Parameters:
6023: + n - number of the time points
6024: - time_points - array of the time points
6026: Level: beginner
6028: Note:
6029: The values obtained are valid until the `TS` object is destroyed.
6031: Both `n` and `time_points` can be `NULL`.
6033: Also used to see time points set by `TSSetTimeSpan()`.
6035: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6036: @*/
6037: PetscErrorCode TSGetEvaluationTimes(TS ts, PetscInt *n, const PetscReal *time_points[])
6038: {
6039: PetscFunctionBegin;
6041: if (n) PetscAssertPointer(n, 2);
6042: if (time_points) PetscAssertPointer(time_points, 3);
6043: if (!ts->eval_times) {
6044: if (n) *n = 0;
6045: if (time_points) *time_points = NULL;
6046: } else {
6047: if (n) *n = ts->eval_times->num_time_points;
6048: if (time_points) *time_points = ts->eval_times->time_points;
6049: }
6050: PetscFunctionReturn(PETSC_SUCCESS);
6051: }
6053: /*@C
6054: TSGetEvaluationSolutions - Get the number of solutions and the solutions at the evaluation time points specified
6056: Input Parameter:
6057: . ts - the `TS` context obtained from `TSCreate()`
6059: Output Parameters:
6060: + nsol - the number of solutions
6061: . sol_times - array of solution times corresponding to the solution vectors. See note below
6062: - Sols - the solution vectors
6064: Level: intermediate
6066: Notes:
6067: Both `nsol` and `Sols` can be `NULL`.
6069: 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()`.
6070: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain evaluation times.
6072: Also used to see view solutions requested by `TSSetTimeSpan()`.
6074: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`
6075: @*/
6076: PetscErrorCode TSGetEvaluationSolutions(TS ts, PetscInt *nsol, const PetscReal *sol_times[], Vec *Sols[])
6077: {
6078: PetscFunctionBegin;
6080: if (nsol) PetscAssertPointer(nsol, 2);
6081: if (sol_times) PetscAssertPointer(sol_times, 3);
6082: if (Sols) PetscAssertPointer(Sols, 4);
6083: if (!ts->eval_times) {
6084: if (nsol) *nsol = 0;
6085: if (sol_times) *sol_times = NULL;
6086: if (Sols) *Sols = NULL;
6087: } else {
6088: if (nsol) *nsol = ts->eval_times->sol_idx;
6089: if (sol_times) *sol_times = ts->eval_times->sol_times;
6090: if (Sols) *Sols = ts->eval_times->sol_vecs;
6091: }
6092: PetscFunctionReturn(PETSC_SUCCESS);
6093: }
6095: /*@
6096: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
6098: Collective
6100: Input Parameters:
6101: + ts - the time-stepper
6102: . n - number of the time points (>=2)
6103: - 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.
6105: Options Database Key:
6106: . -ts_time_span <t0,...tf> - Sets the time span
6108: Level: intermediate
6110: Notes:
6111: 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.
6113: The elements in `span_times` must be all increasing. They correspond to the intermediate points to be saved.
6115: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
6117: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using time span may
6118: pressure the memory system when using a large number of span points.
6120: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6121: @*/
6122: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal span_times[])
6123: {
6124: PetscFunctionBegin;
6126: PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
6127: PetscCall(TSSetEvaluationTimes(ts, n, span_times));
6128: PetscCall(TSSetTime(ts, span_times[0]));
6129: PetscCall(TSSetMaxTime(ts, span_times[n - 1]));
6130: PetscFunctionReturn(PETSC_SUCCESS);
6131: }
6133: /*@
6134: TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
6136: Collective
6138: Input Parameters:
6139: + ts - the `TS` context
6140: . J - Jacobian matrix (not altered in this routine)
6141: - B - newly computed Jacobian matrix to use with preconditioner
6143: Level: intermediate
6145: Notes:
6146: This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
6147: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
6148: and multiple fields are involved.
6150: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
6151: structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
6152: usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
6153: `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.
6155: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
6156: @*/
6157: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
6158: {
6159: MatColoring mc = NULL;
6160: ISColoring iscoloring = NULL;
6161: MatFDColoring matfdcoloring = NULL;
6163: PetscFunctionBegin;
6164: /* Generate new coloring after eliminating zeros in the matrix */
6165: PetscCall(MatEliminateZeros(B, PETSC_TRUE));
6166: PetscCall(MatColoringCreate(B, &mc));
6167: PetscCall(MatColoringSetDistance(mc, 2));
6168: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
6169: PetscCall(MatColoringSetFromOptions(mc));
6170: PetscCall(MatColoringApply(mc, &iscoloring));
6171: PetscCall(MatColoringDestroy(&mc));
6172: /* Replace the old coloring with the new one */
6173: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
6174: PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
6175: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
6176: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
6177: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
6178: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
6179: PetscCall(ISColoringDestroy(&iscoloring));
6180: PetscFunctionReturn(PETSC_SUCCESS);
6181: }