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 - 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 (true|false) - 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, sizeof(typeName), &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 PetscDefined(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, PetscCtx))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 (basic|singlefile|memory|visualization) - 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`, `TSTrajectoryType`, `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: 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_Internal - Evaluates the Jacobian of the DAE
865: Collective
867: Input Parameters:
868: + ts - the `TS` context
869: . ijacobian - function to compute LHS Jacobian
870: . rhsjacobian - function to compute RHS Jacobian
871: . ctx - user context for Jacobian functions
872: . t - current timestep
873: . U - state vector
874: . Udot - time derivative of state vector
875: . shift - shift to apply, see note below
876: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSJacobian should be kept separate
878: Output Parameters:
879: + A - Jacobian matrix
880: - B - matrix from which the preconditioner is constructed; often the same as `A`
882: Level: developer
884: Notes:
885: This function exists so that a user can assemble the Jacobian pieces with functions not stores in the `DMTS`. This was necessary in `TSDISCGRAD` since two different representations of the formulation can be stored.
887: If $ F(t,U,\dot{U})=0 $ is the DAE, the required Jacobian is
888: .vb
889: dF/dU + shift*dF/dUdot
890: .ve
892: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
893: */
894: PetscErrorCode TSComputeIJacobian_Internal(TS ts, TSIJacobianFn *ijacobian, TSRHSJacobianFn *rhsjacobian, void *ctx, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
895: {
896: PetscFunctionBegin;
897: PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
899: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
900: if (ijacobian) {
901: PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
902: ts->ijacs++;
903: }
904: if (imex) {
905: if (!ijacobian) { /* system was written as Udot = G(t,U) */
906: PetscBool assembled;
907: if (rhsjacobian) {
908: Mat Arhs = NULL;
909: PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
910: if (A == Arhs) {
911: 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 */
912: ts->rhsjacobian.time = PETSC_MIN_REAL;
913: }
914: }
915: PetscCall(MatZeroEntries(A));
916: PetscCall(MatAssembled(A, &assembled));
917: if (!assembled) {
918: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
919: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
920: }
921: PetscCall(MatShift(A, shift));
922: if (A != B) {
923: PetscCall(MatZeroEntries(B));
924: PetscCall(MatAssembled(B, &assembled));
925: if (!assembled) {
926: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
927: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
928: }
929: PetscCall(MatShift(B, shift));
930: }
931: }
932: } else {
933: Mat Arhs = NULL, Brhs = NULL;
935: /* RHSJacobian needs to be converted to part of IJacobian if exists */
936: if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
937: if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
938: DM dm;
939: PetscObjectState Ustate;
940: PetscObjectId Uid;
941: TSRHSFunctionFn *rhsfunction;
943: PetscCall(TSGetDM(ts, &dm));
944: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
945: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
946: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
947: if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
948: ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */
949: PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
950: if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
951: } else {
952: PetscBool flg;
954: if (ts->rhsjacobian.reuse) { /* Undo the damage */
955: /* MatScale has a short path for this case.
956: However, this code path is taken the first time TSComputeRHSJacobian is called
957: and the matrices have not been assembled yet */
958: PetscCall(TSRecoverRHSJacobian(ts, A, B));
959: }
960: PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
961: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
962: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
963: if (!flg) {
964: PetscCall(MatScale(A, -1));
965: PetscCall(MatShift(A, shift));
966: }
967: if (A != B) {
968: PetscCall(MatScale(B, -1));
969: PetscCall(MatShift(B, shift));
970: }
971: }
972: ts->rhsjacobian.scale = -1;
973: ts->rhsjacobian.shift = shift;
974: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
975: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
976: PetscCall(MatZeroEntries(A));
977: PetscCall(MatShift(A, shift));
978: if (A != B) {
979: PetscCall(MatZeroEntries(B));
980: PetscCall(MatShift(B, shift));
981: }
982: }
983: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
984: PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
985: if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
986: }
987: }
988: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
989: PetscFunctionReturn(PETSC_SUCCESS);
990: }
992: /*@
993: TSComputeIJacobian - Evaluates the Jacobian of the DAE
995: Collective
997: Input Parameters:
998: + ts - the `TS` context
999: . t - current timestep
1000: . U - state vector
1001: . Udot - time derivative of state vector
1002: . shift - shift to apply, see note below
1003: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSJacobian should be kept separate
1005: Output Parameters:
1006: + A - Jacobian matrix
1007: - B - matrix from which the preconditioner is constructed; often the same as `A`
1009: Level: developer
1011: Notes:
1012: If $ F(t,U,\dot{U})=0 $ is the DAE, the required Jacobian is
1013: .vb
1014: dF/dU + shift*dF/dUdot
1015: .ve
1016: Most users should not need to explicitly call this routine, as it
1017: is used internally within the nonlinear solvers.
1019: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
1020: @*/
1021: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
1022: {
1023: TSIJacobianFn *ijacobian;
1024: TSRHSJacobianFn *rhsjacobian;
1025: DM dm;
1026: void *ctx;
1028: PetscFunctionBegin;
1035: PetscCall(TSGetDM(ts, &dm));
1036: PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
1037: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1038: PetscCall(TSComputeIJacobian_Internal(ts, ijacobian, rhsjacobian, ctx, t, U, Udot, shift, A, B, imex));
1039: PetscFunctionReturn(PETSC_SUCCESS);
1040: }
1042: /*@
1043: TSSetRHSFunction - Sets the routine for evaluating the function,
1044: where U_t = G(t,u).
1046: Logically Collective
1048: Input Parameters:
1049: + ts - the `TS` context obtained from `TSCreate()`
1050: . r - vector to put the computed right-hand side (or `NULL` to have it created)
1051: . f - routine for evaluating the right-hand-side function
1052: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
1054: Level: beginner
1056: Note:
1057: You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.
1059: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1060: @*/
1061: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, PetscCtx ctx)
1062: {
1063: SNES snes;
1064: Vec ralloc = NULL;
1065: DM dm;
1067: PetscFunctionBegin;
1071: PetscCall(TSGetDM(ts, &dm));
1072: PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1073: PetscCall(TSGetSNES(ts, &snes));
1074: if (!r && !ts->dm && ts->vec_sol) {
1075: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1076: r = ralloc;
1077: }
1078: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1079: PetscCall(VecDestroy(&ralloc));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: /*@
1084: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1086: Logically Collective
1088: Input Parameters:
1089: + ts - the `TS` context obtained from `TSCreate()`
1090: . f - routine for evaluating the solution
1091: - ctx - [optional] user-defined context for private data for the
1092: function evaluation routine (may be `NULL`)
1094: Options Database Keys:
1095: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1096: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
1098: Level: intermediate
1100: Notes:
1101: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1102: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1103: create closed-form solutions with non-physical forcing terms.
1105: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1107: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1108: @*/
1109: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, PetscCtx ctx)
1110: {
1111: DM dm;
1113: PetscFunctionBegin;
1115: PetscCall(TSGetDM(ts, &dm));
1116: PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: /*@
1121: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1123: Logically Collective
1125: Input Parameters:
1126: + ts - the `TS` context obtained from `TSCreate()`
1127: . func - routine for evaluating the forcing function
1128: - ctx - [optional] user-defined context for private data for the function evaluation routine
1129: (may be `NULL`)
1131: Level: intermediate
1133: Notes:
1134: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1135: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1136: definition of the problem you are solving and hence possibly introducing bugs.
1138: This replaces the ODE F(u,u_t,t) = 0 the `TS` is solving with F(u,u_t,t) - func(t) = 0
1140: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1141: parameters can be passed in the ctx variable.
1143: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1145: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
1146: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1147: @*/
1148: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, PetscCtx ctx)
1149: {
1150: DM dm;
1152: PetscFunctionBegin;
1154: PetscCall(TSGetDM(ts, &dm));
1155: PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1156: PetscFunctionReturn(PETSC_SUCCESS);
1157: }
1159: /*@
1160: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1161: where U_t = G(U,t), as well as the location to store the matrix.
1163: Logically Collective
1165: Input Parameters:
1166: + ts - the `TS` context obtained from `TSCreate()`
1167: . Amat - (approximate) location to store Jacobian matrix entries computed by `f`
1168: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1169: . f - the Jacobian evaluation routine
1170: - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1172: Level: beginner
1174: Notes:
1175: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1177: The `TS` solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f()`
1178: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1180: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
1181: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
1182: @*/
1183: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, PetscCtx ctx)
1184: {
1185: SNES snes;
1186: DM dm;
1187: TSIJacobianFn *ijacobian;
1189: PetscFunctionBegin;
1193: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1194: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1196: PetscCall(TSGetDM(ts, &dm));
1197: PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1198: PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1199: PetscCall(TSGetSNES(ts, &snes));
1200: if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1201: if (Amat) {
1202: PetscCall(PetscObjectReference((PetscObject)Amat));
1203: PetscCall(MatDestroy(&ts->Arhs));
1204: ts->Arhs = Amat;
1205: }
1206: if (Pmat) {
1207: PetscCall(PetscObjectReference((PetscObject)Pmat));
1208: PetscCall(MatDestroy(&ts->Brhs));
1209: ts->Brhs = Pmat;
1210: }
1211: PetscFunctionReturn(PETSC_SUCCESS);
1212: }
1214: /*@
1215: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1217: Logically Collective
1219: Input Parameters:
1220: + ts - the `TS` context obtained from `TSCreate()`
1221: . r - vector to hold the residual (or `NULL` to have it created internally)
1222: . f - the function evaluation routine
1223: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1225: Level: beginner
1227: Note:
1228: The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function.
1230: .seealso: [](ch_ts), `TS`, `TSIFunctionFn`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`,
1231: `TSSetIJacobian()`
1232: @*/
1233: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunctionFn *f, PetscCtx ctx)
1234: {
1235: SNES snes;
1236: Vec ralloc = NULL;
1237: DM dm;
1239: PetscFunctionBegin;
1243: PetscCall(TSGetDM(ts, &dm));
1244: PetscCall(DMTSSetIFunction(dm, f, ctx));
1246: PetscCall(TSGetSNES(ts, &snes));
1247: if (!r && !ts->dm && ts->vec_sol) {
1248: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1249: r = ralloc;
1250: }
1251: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1252: PetscCall(VecDestroy(&ralloc));
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: /*@
1257: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.
1259: Not Collective
1261: Input Parameter:
1262: . ts - the `TS` context
1264: Output Parameters:
1265: + r - vector to hold residual (or `NULL`)
1266: . func - the function to compute residual (or `NULL`)
1267: - ctx - the function context (or `NULL`)
1269: Level: advanced
1271: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1272: @*/
1273: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunctionFn **func, PetscCtxRt ctx)
1274: {
1275: SNES snes;
1276: DM dm;
1278: PetscFunctionBegin;
1280: PetscCall(TSGetSNES(ts, &snes));
1281: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1282: PetscCall(TSGetDM(ts, &dm));
1283: PetscCall(DMTSGetIFunction(dm, func, ctx));
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1287: /*@
1288: TSGetRHSFunction - Returns the vector where the right-hand side is stored and the function/context to compute it.
1290: Not Collective
1292: Input Parameter:
1293: . ts - the `TS` context
1295: Output Parameters:
1296: + r - vector to hold computed right-hand side (or `NULL`)
1297: . func - the function to compute right-hand side (or `NULL`)
1298: - ctx - the function context (or `NULL`)
1300: Level: advanced
1302: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1303: @*/
1304: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunctionFn **func, PetscCtxRt ctx)
1305: {
1306: SNES snes;
1307: DM dm;
1309: PetscFunctionBegin;
1311: PetscCall(TSGetSNES(ts, &snes));
1312: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1313: PetscCall(TSGetDM(ts, &dm));
1314: PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1315: PetscFunctionReturn(PETSC_SUCCESS);
1316: }
1318: /*@
1319: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1320: provided with `TSSetIFunction()`.
1322: Logically Collective
1324: Input Parameters:
1325: + ts - the `TS` context obtained from `TSCreate()`
1326: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
1327: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1328: . f - the Jacobian evaluation routine
1329: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1331: Level: beginner
1333: Notes:
1334: The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1336: If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
1337: space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.
1339: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1340: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1341: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1342: a and vector W depend on the integration method, step size, and past states. For example with
1343: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1344: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1346: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1348: The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f`
1349: You should not assume the values are the same in the next call to `f` as you set them in the previous call.
1351: In case `TSSetRHSJacobian()` is also used in conjunction with a fully-implicit solver,
1352: multilevel linear solvers, e.g. `PCMG`, will likely not work due to the way `TS` handles rhs matrices.
1354: .seealso: [](ch_ts), `TS`, `TSIJacobianFn`, `TSSetIFunction()`, `TSSetRHSJacobian()`,
1355: `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1356: @*/
1357: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobianFn *f, PetscCtx ctx)
1358: {
1359: SNES snes;
1360: DM dm;
1362: PetscFunctionBegin;
1366: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1367: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1369: PetscCall(TSGetDM(ts, &dm));
1370: PetscCall(DMTSSetIJacobian(dm, f, ctx));
1372: PetscCall(TSGetSNES(ts, &snes));
1373: PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1374: PetscFunctionReturn(PETSC_SUCCESS);
1375: }
1377: /*@
1378: TSRHSJacobianSetReuse - restore the RHS Jacobian before calling the user-provided `TSRHSJacobianFn` function again
1380: Logically Collective
1382: Input Parameters:
1383: + ts - `TS` context obtained from `TSCreate()`
1384: - reuse - `PETSC_TRUE` if the RHS Jacobian
1386: Level: intermediate
1388: Notes:
1389: Without this flag, `TS` will change the sign and shift the RHS Jacobian for a
1390: finite-time-step implicit solve, in which case the user function will need to recompute the
1391: entire Jacobian. The `reuse `flag must be set if the evaluation function assumes that the
1392: matrix entries have not been changed by the `TS`.
1394: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1395: @*/
1396: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1397: {
1398: PetscFunctionBegin;
1399: ts->rhsjacobian.reuse = reuse;
1400: PetscFunctionReturn(PETSC_SUCCESS);
1401: }
1403: /*@
1404: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1406: Logically Collective
1408: Input Parameters:
1409: + ts - the `TS` context obtained from `TSCreate()`
1410: . F - vector to hold the residual (or `NULL` to have it created internally)
1411: . fun - the function evaluation routine
1412: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1414: Level: beginner
1416: .seealso: [](ch_ts), `TS`, `TSI2FunctionFn`, `TSSetI2Jacobian()`, `TSSetIFunction()`,
1417: `TSCreate()`, `TSSetRHSFunction()`
1418: @*/
1419: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2FunctionFn *fun, PetscCtx ctx)
1420: {
1421: DM dm;
1423: PetscFunctionBegin;
1426: PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1427: PetscCall(TSGetDM(ts, &dm));
1428: PetscCall(DMTSSetI2Function(dm, fun, ctx));
1429: PetscFunctionReturn(PETSC_SUCCESS);
1430: }
1432: /*@
1433: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.
1435: Not Collective
1437: Input Parameter:
1438: . ts - the `TS` context
1440: Output Parameters:
1441: + r - vector to hold residual (or `NULL`)
1442: . fun - the function to compute residual (or `NULL`)
1443: - ctx - the function context (or `NULL`)
1445: Level: advanced
1447: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1448: @*/
1449: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2FunctionFn **fun, PetscCtxRt ctx)
1450: {
1451: SNES snes;
1452: DM dm;
1454: PetscFunctionBegin;
1456: PetscCall(TSGetSNES(ts, &snes));
1457: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1458: PetscCall(TSGetDM(ts, &dm));
1459: PetscCall(DMTSGetI2Function(dm, fun, ctx));
1460: PetscFunctionReturn(PETSC_SUCCESS);
1461: }
1463: /*@
1464: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1465: where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.
1467: Logically Collective
1469: Input Parameters:
1470: + ts - the `TS` context obtained from `TSCreate()`
1471: . J - matrix to hold the Jacobian values
1472: . P - matrix for constructing the preconditioner (may be same as `J`)
1473: . jac - the Jacobian evaluation routine, see `TSI2JacobianFn` for the calling sequence
1474: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1476: Level: beginner
1478: Notes:
1479: The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1481: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1482: 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.
1483: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1484: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1486: .seealso: [](ch_ts), `TS`, `TSI2JacobianFn`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1487: @*/
1488: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2JacobianFn *jac, PetscCtx ctx)
1489: {
1490: DM dm;
1492: PetscFunctionBegin;
1496: PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1497: PetscCall(TSGetDM(ts, &dm));
1498: PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1499: PetscFunctionReturn(PETSC_SUCCESS);
1500: }
1502: /*@
1503: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1505: Not Collective, but parallel objects are returned if `TS` is parallel
1507: Input Parameter:
1508: . ts - The `TS` context obtained from `TSCreate()`
1510: Output Parameters:
1511: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1512: . P - The matrix from which the preconditioner is constructed, often the same as `J`
1513: . jac - The function to compute the Jacobian matrices
1514: - ctx - User-defined context for Jacobian evaluation routine
1516: Level: advanced
1518: Note:
1519: You can pass in `NULL` for any return argument you do not need.
1521: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1522: @*/
1523: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2JacobianFn **jac, PetscCtxRt ctx)
1524: {
1525: SNES snes;
1526: DM dm;
1528: PetscFunctionBegin;
1529: PetscCall(TSGetSNES(ts, &snes));
1530: PetscCall(SNESSetUpMatrices(snes));
1531: PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1532: PetscCall(TSGetDM(ts, &dm));
1533: PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1534: PetscFunctionReturn(PETSC_SUCCESS);
1535: }
1537: /*@
1538: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1540: Collective
1542: Input Parameters:
1543: + ts - the `TS` context
1544: . t - current time
1545: . U - state vector
1546: . V - time derivative of state vector (U_t)
1547: - A - second time derivative of state vector (U_tt)
1549: Output Parameter:
1550: . F - the residual vector
1552: Level: developer
1554: Note:
1555: Most users should not need to explicitly call this routine, as it
1556: is used internally within the nonlinear solvers.
1558: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1559: @*/
1560: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1561: {
1562: DM dm;
1563: TSI2FunctionFn *I2Function;
1564: void *ctx;
1565: TSRHSFunctionFn *rhsfunction;
1567: PetscFunctionBegin;
1574: PetscCall(TSGetDM(ts, &dm));
1575: PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1576: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
1578: if (!I2Function) {
1579: PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, V, F));
1585: PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));
1587: if (rhsfunction) {
1588: Vec Frhs;
1590: PetscCall(DMGetGlobalVector(dm, &Frhs));
1591: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1592: PetscCall(VecAXPY(F, -1, Frhs));
1593: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
1594: }
1596: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, V, F));
1597: PetscFunctionReturn(PETSC_SUCCESS);
1598: }
1600: /*@
1601: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1603: Collective
1605: Input Parameters:
1606: + ts - the `TS` context
1607: . t - current timestep
1608: . U - state vector
1609: . V - time derivative of state vector
1610: . A - second time derivative of state vector
1611: . shiftV - shift to apply, see note below
1612: - shiftA - shift to apply, see note below
1614: Output Parameters:
1615: + J - Jacobian matrix
1616: - P - optional matrix used to construct the preconditioner
1618: Level: developer
1620: Notes:
1621: If $F(t,U,V,A) = 0$ is the DAE, the required Jacobian is
1623: $$
1624: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1625: $$
1627: Most users should not need to explicitly call this routine, as it
1628: is used internally within the ODE integrators.
1630: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1631: @*/
1632: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1633: {
1634: DM dm;
1635: TSI2JacobianFn *I2Jacobian;
1636: void *ctx;
1637: TSRHSJacobianFn *rhsjacobian;
1639: PetscFunctionBegin;
1647: PetscCall(TSGetDM(ts, &dm));
1648: PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1649: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1651: if (!I2Jacobian) {
1652: PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1653: PetscFunctionReturn(PETSC_SUCCESS);
1654: }
1656: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, J, P));
1657: PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1658: if (rhsjacobian) {
1659: Mat Jrhs, Prhs;
1660: PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1661: PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1662: PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1663: if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1664: }
1666: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, J, P));
1667: PetscFunctionReturn(PETSC_SUCCESS);
1668: }
1670: /*@
1671: TSSetTransientVariable - sets function to transform from state to transient variables
1673: Logically Collective
1675: Input Parameters:
1676: + ts - time stepping context on which to change the transient variable
1677: . tvar - a function that transforms to transient variables, see `TSTransientVariableFn` for the calling sequence
1678: - ctx - a context for tvar
1680: Level: advanced
1682: Notes:
1683: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1684: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1685: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1686: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1687: evaluated via the chain rule, as in
1688: .vb
1689: dF/dP + shift * dF/dCdot dC/dP.
1690: .ve
1692: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSTransientVariableFn`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1693: @*/
1694: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariableFn *tvar, PetscCtx ctx)
1695: {
1696: DM dm;
1698: PetscFunctionBegin;
1700: PetscCall(TSGetDM(ts, &dm));
1701: PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1702: PetscFunctionReturn(PETSC_SUCCESS);
1703: }
1705: /*@
1706: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1708: Logically Collective
1710: Input Parameters:
1711: + ts - TS on which to compute
1712: - U - state vector to be transformed to transient variables
1714: Output Parameter:
1715: . C - transient (conservative) variable
1717: Level: developer
1719: Developer Notes:
1720: If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed.
1721: This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are
1722: being used.
1724: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1725: @*/
1726: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1727: {
1728: DM dm;
1729: DMTS dmts;
1731: PetscFunctionBegin;
1734: PetscCall(TSGetDM(ts, &dm));
1735: PetscCall(DMGetDMTS(dm, &dmts));
1736: if (dmts->ops->transientvar) {
1738: PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1739: }
1740: PetscFunctionReturn(PETSC_SUCCESS);
1741: }
1743: /*@
1744: TSHasTransientVariable - determine whether transient variables have been set
1746: Logically Collective
1748: Input Parameter:
1749: . ts - `TS` on which to compute
1751: Output Parameter:
1752: . has - `PETSC_TRUE` if transient variables have been set
1754: Level: developer
1756: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1757: @*/
1758: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1759: {
1760: DM dm;
1761: DMTS dmts;
1763: PetscFunctionBegin;
1765: PetscCall(TSGetDM(ts, &dm));
1766: PetscCall(DMGetDMTS(dm, &dmts));
1767: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1768: PetscFunctionReturn(PETSC_SUCCESS);
1769: }
1771: /*@
1772: TS2SetSolution - Sets the initial solution and time derivative vectors
1773: for use by the `TS` routines handling second order equations.
1775: Logically Collective
1777: Input Parameters:
1778: + ts - the `TS` context obtained from `TSCreate()`
1779: . u - the solution vector
1780: - v - the time derivative vector
1782: Level: beginner
1784: .seealso: [](ch_ts), `TS`
1785: @*/
1786: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1787: {
1788: PetscFunctionBegin;
1792: PetscCall(TSSetSolution(ts, u));
1793: PetscCall(PetscObjectReference((PetscObject)v));
1794: PetscCall(VecDestroy(&ts->vec_dot));
1795: ts->vec_dot = v;
1796: PetscFunctionReturn(PETSC_SUCCESS);
1797: }
1799: /*@
1800: TS2GetSolution - Returns the solution and time derivative at the present timestep
1801: for second order equations.
1803: Not Collective
1805: Input Parameter:
1806: . ts - the `TS` context obtained from `TSCreate()`
1808: Output Parameters:
1809: + u - the vector containing the solution
1810: - v - the vector containing the time derivative
1812: Level: intermediate
1814: Notes:
1815: It is valid to call this routine inside the function
1816: that you are evaluating in order to move to the new timestep. This vector not
1817: changed until the solution at the next timestep has been calculated.
1819: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1820: @*/
1821: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1822: {
1823: PetscFunctionBegin;
1825: if (u) PetscAssertPointer(u, 2);
1826: if (v) PetscAssertPointer(v, 3);
1827: if (u) *u = ts->vec_sol;
1828: if (v) *v = ts->vec_dot;
1829: PetscFunctionReturn(PETSC_SUCCESS);
1830: }
1832: /*@
1833: TSLoad - Loads a `TS` that has been stored in binary with `TSView()`.
1835: Collective
1837: Input Parameters:
1838: + ts - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1839: some related function before a call to `TSLoad()`.
1840: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1842: Level: intermediate
1844: Note:
1845: The type is determined by the data in the file, any type set into the `TS` before this call is ignored.
1847: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1848: @*/
1849: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1850: {
1851: PetscBool isbinary;
1852: PetscInt classid;
1853: char type[256];
1854: DMTS sdm;
1855: DM dm;
1857: PetscFunctionBegin;
1860: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1861: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1863: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1864: PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1865: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1866: PetscCall(TSSetType(ts, type));
1867: PetscTryTypeMethod(ts, load, viewer);
1868: PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1869: PetscCall(DMLoad(dm, viewer));
1870: PetscCall(TSSetDM(ts, dm));
1871: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1872: PetscCall(VecLoad(ts->vec_sol, viewer));
1873: PetscCall(DMGetDMTS(ts->dm, &sdm));
1874: PetscCall(DMTSLoad(sdm, viewer));
1875: PetscFunctionReturn(PETSC_SUCCESS);
1876: }
1878: #include <petscdraw.h>
1879: #if PetscDefined(HAVE_SAWS)
1880: #include <petscviewersaws.h>
1881: #endif
1883: /*@
1884: TSViewFromOptions - View a `TS` based on values in the options database
1886: Collective
1888: Input Parameters:
1889: + ts - the `TS` context
1890: . obj - optional object that provides the prefix for the options database keys, pass `NULL` to use the options prefix of `ts`
1891: - name - command line option string to be passed by user
1893: Options Database Key:
1894: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
1896: Level: intermediate
1898: Note:
1899: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
1900: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
1902: .seealso: [](ch_ts), `TS`, `TSView()`, `PetscObjectViewFromOptions()`, `TSCreate()`, `PetscOptionsCreateViewer()`
1903: @*/
1904: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1905: {
1906: PetscFunctionBegin;
1908: PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1909: PetscFunctionReturn(PETSC_SUCCESS);
1910: }
1912: /*@
1913: TSView - Displays the `TS` data structure.
1915: Collective
1917: Input Parameters:
1918: + ts - the `TS` context obtained from `TSCreate()`
1919: - viewer - visualization context
1921: Options Database Key:
1922: . -ts_view viewer_specification - calls `TSView()` at end of `TSStep()`. See `PetscOptionsCreateViewer()` for the format of `viewer_specification`
1924: Level: beginner
1926: Notes:
1927: The available visualization contexts include
1928: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1929: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1930: output where only the first processor opens
1931: the file. All other processors send their
1932: data to the first processor to print.
1934: The user can open an alternative visualization context with
1935: `PetscViewerASCIIOpen()` - output to a specified file.
1937: In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).
1939: The "initial time step" displayed is the default time step from `TSCreate()` or that set with `TSSetTimeStep()` or `-ts_time_step`
1941: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`, `PetscOptionsCreateViewer()`, `TSViewFromOptions()`
1942: @*/
1943: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1944: {
1945: TSType type;
1946: PetscBool isascii, isstring, issundials, isbinary, isdraw;
1947: DMTS sdm;
1948: #if PetscDefined(HAVE_SAWS)
1949: PetscBool issaws;
1950: #endif
1952: PetscFunctionBegin;
1954: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1956: PetscCheckSameComm(ts, 1, viewer, 2);
1958: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1959: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1960: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1961: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1962: #if PetscDefined(HAVE_SAWS)
1963: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1964: #endif
1965: if (isascii) {
1966: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer));
1967: if (ts->ops->view) {
1968: PetscCall(PetscViewerASCIIPushTab(viewer));
1969: PetscUseTypeMethod(ts, view, viewer);
1970: PetscCall(PetscViewerASCIIPopTab(viewer));
1971: }
1972: PetscCall(PetscViewerASCIIPrintf(viewer, " initial time step=%g\n", (double)ts->initial_time_step));
1973: if (ts->max_steps < PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps));
1974: if (ts->run_steps < PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, " run steps=%" PetscInt_FMT "\n", ts->run_steps));
1975: if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time));
1976: if (ts->max_reject != PETSC_UNLIMITED) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum number of step rejections=%" PetscInt_FMT "\n", ts->max_reject));
1977: if (ts->max_snes_failures != PETSC_UNLIMITED) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum number of SNES failures allowed=%" PetscInt_FMT "\n", ts->max_snes_failures));
1978: if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs));
1979: if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs));
1980: if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs));
1981: if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs));
1982: if (ts->usessnes) {
1983: PetscBool lin;
1984: if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its));
1985: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its));
1986: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
1987: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures));
1988: }
1989: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject));
1990: if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, "));
1991: else PetscCall(PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol));
1992: if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, "using vector of absolute error tolerances\n"));
1993: else PetscCall(PetscViewerASCIIPrintf(viewer, "using absolute error tolerance of %g\n", (double)ts->atol));
1994: PetscCall(PetscViewerASCIIPushTab(viewer));
1995: PetscCall(TSAdaptView(ts->adapt, viewer));
1996: PetscCall(PetscViewerASCIIPopTab(viewer));
1997: } else if (isstring) {
1998: PetscCall(TSGetType(ts, &type));
1999: PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type));
2000: PetscTryTypeMethod(ts, view, viewer);
2001: } else if (isbinary) {
2002: PetscInt classid = TS_FILE_CLASSID;
2003: MPI_Comm comm;
2004: PetscMPIInt rank;
2005: char type[256];
2007: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
2008: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2009: if (rank == 0) {
2010: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
2011: PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
2012: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
2013: }
2014: PetscTryTypeMethod(ts, view, viewer);
2015: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
2016: PetscCall(DMView(ts->dm, viewer));
2017: PetscCall(VecView(ts->vec_sol, viewer));
2018: PetscCall(DMGetDMTS(ts->dm, &sdm));
2019: PetscCall(DMTSView(sdm, viewer));
2020: } else if (isdraw) {
2021: PetscDraw draw;
2022: char str[36];
2023: PetscReal x, y, bottom, h;
2025: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2026: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
2027: PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
2028: PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
2029: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
2030: bottom = y - h;
2031: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
2032: PetscTryTypeMethod(ts, view, viewer);
2033: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
2034: if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
2035: PetscCall(PetscDrawPopCurrentPoint(draw));
2036: #if PetscDefined(HAVE_SAWS)
2037: } else if (issaws) {
2038: PetscMPIInt rank;
2039: const char *name;
2041: PetscCall(PetscObjectGetName((PetscObject)ts, &name));
2042: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
2043: if (!((PetscObject)ts)->amsmem && rank == 0) {
2044: char dir[1024];
2046: PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
2047: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
2048: PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
2049: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
2050: PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
2051: }
2052: PetscTryTypeMethod(ts, view, viewer);
2053: #endif
2054: }
2055: if (ts->snes && ts->usessnes) {
2056: PetscCall(PetscViewerASCIIPushTab(viewer));
2057: PetscCall(SNESView(ts->snes, viewer));
2058: PetscCall(PetscViewerASCIIPopTab(viewer));
2059: }
2060: PetscCall(DMGetDMTS(ts->dm, &sdm));
2061: PetscCall(DMTSView(sdm, viewer));
2063: PetscCall(PetscViewerASCIIPushTab(viewer));
2064: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &issundials));
2065: PetscCall(PetscViewerASCIIPopTab(viewer));
2066: PetscFunctionReturn(PETSC_SUCCESS);
2067: }
2069: /*@
2070: TSSetApplicationContext - Sets an optional user-defined context for the timesteppers that may be accessed, for example inside the user provided
2071: `TS` callbacks with `TSGetApplicationContext()`
2073: Logically Collective
2075: Input Parameters:
2076: + ts - the `TS` context obtained from `TSCreate()`
2077: - ctx - application context
2079: Level: intermediate
2081: Fortran Note:
2082: This only works when `ctx` is a Fortran derived type (it cannot be a `PetscObject`), we recommend writing a Fortran interface definition for this
2083: function that tells the Fortran compiler the derived data type that is passed in as the `ctx` argument. See `TSGetApplicationContext()` for
2084: an example.
2086: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2087: @*/
2088: PetscErrorCode TSSetApplicationContext(TS ts, PetscCtx ctx)
2089: {
2090: PetscFunctionBegin;
2092: ts->ctx = ctx;
2093: PetscFunctionReturn(PETSC_SUCCESS);
2094: }
2096: /*@
2097: TSGetApplicationContext - Gets the user-defined context for the
2098: timestepper that was set with `TSSetApplicationContext()`
2100: Not Collective
2102: Input Parameter:
2103: . ts - the `TS` context obtained from `TSCreate()`
2105: Output Parameter:
2106: . ctx - a pointer to the application context
2108: Level: intermediate
2110: Fortran Notes:
2111: This only works when the context is a Fortran derived type or a `PetscObject`. Declare `ctx` with
2112: .vb
2113: type(tUsertype), pointer :: ctx
2114: .ve
2116: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2117: @*/
2118: PetscErrorCode TSGetApplicationContext(TS ts, PetscCtxRt ctx)
2119: {
2120: PetscFunctionBegin;
2122: *(void **)ctx = ts->ctx;
2123: PetscFunctionReturn(PETSC_SUCCESS);
2124: }
2126: /*@
2127: TSGetStepNumber - Gets the number of time steps completed.
2129: Not Collective
2131: Input Parameter:
2132: . ts - the `TS` context obtained from `TSCreate()`
2134: Output Parameter:
2135: . steps - number of steps completed so far
2137: Level: intermediate
2139: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2140: @*/
2141: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2142: {
2143: PetscFunctionBegin;
2145: PetscAssertPointer(steps, 2);
2146: *steps = ts->steps;
2147: PetscFunctionReturn(PETSC_SUCCESS);
2148: }
2150: /*@
2151: TSSetStepNumber - Sets the number of steps completed.
2153: Logically Collective
2155: Input Parameters:
2156: + ts - the `TS` context
2157: - steps - number of steps completed so far
2159: Level: developer
2161: Note:
2162: For most uses of the `TS` solvers the user need not explicitly call
2163: `TSSetStepNumber()`, as the step counter is appropriately updated in
2164: `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2165: reinitialize timestepping by setting the step counter to zero (and time
2166: to the initial time) to solve a similar problem with different initial
2167: conditions or parameters. Other possible use case is to continue
2168: timestepping from a previously interrupted run in such a way that `TS`
2169: monitors will be called with a initial nonzero step counter.
2171: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2172: @*/
2173: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2174: {
2175: PetscFunctionBegin;
2178: PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2179: ts->steps = steps;
2180: PetscFunctionReturn(PETSC_SUCCESS);
2181: }
2183: /*@
2184: TSSetTimeStep - Allows one to reset the timestep at any time.
2186: Logically Collective
2188: Input Parameters:
2189: + ts - the `TS` context obtained from `TSCreate()`
2190: - time_step - the size of the timestep
2192: Options Database Key:
2193: . -ts_time_step dt - provide the initial time step
2195: Level: intermediate
2197: Notes:
2198: This is only a suggestion, the actual initial time step used may differ
2200: If this is called after `TSSetUp()`, it will not change the initial time step value printed by `TSView()`
2202: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2203: @*/
2204: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2205: {
2206: PetscFunctionBegin;
2209: ts->time_step = time_step;
2210: if (ts->setupcalled == PETSC_FALSE) ts->initial_time_step = time_step;
2211: PetscFunctionReturn(PETSC_SUCCESS);
2212: }
2214: /*@
2215: TSSetExactFinalTime - Determines whether to adapt the final time step to
2216: match the exact final time, to interpolate the solution to the exact final time,
2217: or to just return at the final time `TS` computed (which may be slightly larger
2218: than the requested final time).
2220: Logically Collective
2222: Input Parameters:
2223: + ts - the time-step context
2224: - eftopt - exact final time option
2225: .vb
2226: TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded, just use it
2227: TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time if the final time is exceeded
2228: TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to ensure the computed final time exactly equals the requested final time
2229: .ve
2231: Options Database Key:
2232: . -ts_exact_final_time stepover,interpolate,matchstep - select the final step approach at runtime
2234: Level: beginner
2236: Note:
2237: If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2238: then the final time you selected.
2240: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2241: @*/
2242: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2243: {
2244: PetscFunctionBegin;
2247: ts->exact_final_time = eftopt;
2248: PetscFunctionReturn(PETSC_SUCCESS);
2249: }
2251: /*@
2252: TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`
2254: Not Collective
2256: Input Parameter:
2257: . ts - the `TS` context
2259: Output Parameter:
2260: . eftopt - exact final time option
2262: Level: beginner
2264: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2265: @*/
2266: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2267: {
2268: PetscFunctionBegin;
2270: PetscAssertPointer(eftopt, 2);
2271: *eftopt = ts->exact_final_time;
2272: PetscFunctionReturn(PETSC_SUCCESS);
2273: }
2275: /*@
2276: TSGetTimeStep - Gets the current timestep size.
2278: Not Collective
2280: Input Parameter:
2281: . ts - the `TS` context obtained from `TSCreate()`
2283: Output Parameter:
2284: . dt - the current timestep size
2286: Level: intermediate
2288: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2289: @*/
2290: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2291: {
2292: PetscFunctionBegin;
2294: PetscAssertPointer(dt, 2);
2295: *dt = ts->time_step;
2296: PetscFunctionReturn(PETSC_SUCCESS);
2297: }
2299: /*@
2300: TSGetSolution - Returns the solution at the present timestep. It
2301: is valid to call this routine inside the function that you are evaluating
2302: in order to move to the new timestep. This vector not changed until
2303: the solution at the next timestep has been calculated.
2305: Not Collective, but v returned is parallel if ts is parallel
2307: Input Parameter:
2308: . ts - the `TS` context obtained from `TSCreate()`
2310: Output Parameter:
2311: . v - the vector containing the solution
2313: Level: intermediate
2315: Note:
2316: If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
2317: final time. It returns the solution at the next timestep.
2319: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2320: @*/
2321: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2322: {
2323: PetscFunctionBegin;
2325: PetscAssertPointer(v, 2);
2326: *v = ts->vec_sol;
2327: PetscFunctionReturn(PETSC_SUCCESS);
2328: }
2330: /*@
2331: TSGetSolutionComponents - Returns any solution components at the present
2332: timestep, if available for the time integration method being used.
2333: Solution components are quantities that share the same size and
2334: structure as the solution vector.
2336: Not Collective, but v returned is parallel if ts is parallel
2338: Input Parameters:
2339: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2340: . n - If v is `NULL`, then the number of solution components is
2341: returned through n, else the n-th solution component is
2342: returned in v.
2343: - v - the vector containing the n-th solution component
2344: (may be `NULL` to use this function to find out
2345: the number of solutions components).
2347: Level: advanced
2349: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2350: @*/
2351: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2352: {
2353: PetscFunctionBegin;
2355: if (!ts->ops->getsolutioncomponents) *n = 0;
2356: else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2357: PetscFunctionReturn(PETSC_SUCCESS);
2358: }
2360: /*@
2361: TSGetAuxSolution - Returns an auxiliary solution at the present
2362: timestep, if available for the time integration method being used.
2364: Not Collective, but v returned is parallel if ts is parallel
2366: Input Parameters:
2367: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2368: - v - the vector containing the auxiliary solution
2370: Level: intermediate
2372: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2373: @*/
2374: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2375: {
2376: PetscFunctionBegin;
2378: if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2379: else PetscCall(VecZeroEntries(*v));
2380: PetscFunctionReturn(PETSC_SUCCESS);
2381: }
2383: /*@
2384: TSGetTimeError - Returns the estimated error vector, if the chosen
2385: `TSType` has an error estimation functionality and `TSSetTimeError()` was called
2387: Not Collective, but v returned is parallel if ts is parallel
2389: Input Parameters:
2390: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2391: . n - current estimate (n=0) or previous one (n=-1)
2392: - v - the vector containing the error (same size as the solution).
2394: Level: intermediate
2396: Note:
2397: MUST call after `TSSetUp()`
2399: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2400: @*/
2401: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2402: {
2403: PetscFunctionBegin;
2405: if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2406: else PetscCall(VecZeroEntries(*v));
2407: PetscFunctionReturn(PETSC_SUCCESS);
2408: }
2410: /*@
2411: TSSetTimeError - Sets the estimated error vector, if the chosen
2412: `TSType` has an error estimation functionality. This can be used
2413: to restart such a time integrator with a given error vector.
2415: Not Collective, but v returned is parallel if ts is parallel
2417: Input Parameters:
2418: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2419: - v - the vector containing the error (same size as the solution).
2421: Level: intermediate
2423: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
2424: @*/
2425: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2426: {
2427: PetscFunctionBegin;
2429: PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2430: PetscTryTypeMethod(ts, settimeerror, v);
2431: PetscFunctionReturn(PETSC_SUCCESS);
2432: }
2434: /* ----- Routines to initialize and destroy a timestepper ---- */
2435: /*@
2436: TSSetProblemType - Sets the type of problem to be solved.
2438: Not collective
2440: Input Parameters:
2441: + ts - The `TS`
2442: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2443: .vb
2444: U_t - A U = 0 (linear)
2445: U_t - A(t) U = 0 (linear)
2446: F(t,U,U_t) = 0 (nonlinear)
2447: .ve
2449: Level: beginner
2451: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2452: @*/
2453: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2454: {
2455: PetscFunctionBegin;
2457: ts->problem_type = type;
2458: if (type == TS_LINEAR) {
2459: SNES snes;
2460: PetscCall(TSGetSNES(ts, &snes));
2461: PetscCall(SNESSetType(snes, SNESKSPONLY));
2462: }
2463: PetscFunctionReturn(PETSC_SUCCESS);
2464: }
2466: /*@
2467: TSGetProblemType - Gets the type of problem to be solved.
2469: Not collective
2471: Input Parameter:
2472: . ts - The `TS`
2474: Output Parameter:
2475: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2476: .vb
2477: M U_t = A U
2478: M(t) U_t = A(t) U
2479: F(t,U,U_t)
2480: .ve
2482: Level: beginner
2484: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2485: @*/
2486: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2487: {
2488: PetscFunctionBegin;
2490: PetscAssertPointer(type, 2);
2491: *type = ts->problem_type;
2492: PetscFunctionReturn(PETSC_SUCCESS);
2493: }
2495: /*
2496: Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2497: */
2498: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2499: {
2500: PetscBool isnone;
2502: PetscFunctionBegin;
2503: PetscCall(TSGetAdapt(ts, &ts->adapt));
2504: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2506: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2507: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2508: else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2509: PetscFunctionReturn(PETSC_SUCCESS);
2510: }
2512: /*@
2513: TSSetUp - Sets up the internal data structures for the later use of a timestepper.
2515: Collective
2517: Input Parameter:
2518: . ts - the `TS` context obtained from `TSCreate()`
2520: Level: advanced
2522: Note:
2523: For basic use of the `TS` solvers the user need not explicitly call
2524: `TSSetUp()`, since these actions will automatically occur during
2525: the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this
2526: phase separately, `TSSetUp()` should be called after `TSCreate()`
2527: and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.
2529: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2530: @*/
2531: PetscErrorCode TSSetUp(TS ts)
2532: {
2533: DM dm;
2534: PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2535: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2536: TSIFunctionFn *ifun;
2537: TSIJacobianFn *ijac;
2538: TSI2JacobianFn *i2jac;
2539: TSRHSJacobianFn *rhsjac;
2541: PetscFunctionBegin;
2543: if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
2545: if (!((PetscObject)ts)->type_name) {
2546: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2547: PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2548: }
2550: if (!ts->vec_sol) {
2551: PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2552: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2553: }
2555: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2556: PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2557: ts->Jacp = ts->Jacprhs;
2558: }
2560: if (ts->quadraturets) {
2561: PetscCall(TSSetUp(ts->quadraturets));
2562: PetscCall(VecDestroy(&ts->vec_costintegrand));
2563: PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2564: }
2566: PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
2567: if (rhsjac == TSComputeRHSJacobianConstant) {
2568: Mat Amat, Pmat;
2569: SNES snes;
2570: PetscCall(TSGetSNES(ts, &snes));
2571: PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
2572: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2573: * have displaced the RHS matrix */
2574: if (Amat && Amat == ts->Arhs) {
2575: /* 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 */
2576: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
2577: PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
2578: PetscCall(MatDestroy(&Amat));
2579: }
2580: if (Pmat && Pmat == ts->Brhs) {
2581: PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
2582: PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
2583: PetscCall(MatDestroy(&Pmat));
2584: }
2585: }
2587: PetscCall(TSGetAdapt(ts, &ts->adapt));
2588: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2590: PetscTryTypeMethod(ts, setup);
2592: PetscCall(TSSetExactFinalTimeDefault(ts));
2594: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2595: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2596: */
2597: PetscCall(TSGetDM(ts, &dm));
2598: PetscCall(DMSNESGetFunction(dm, &func, NULL));
2599: if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));
2601: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2602: Otherwise, the SNES will use coloring internally to form the Jacobian.
2603: */
2604: PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2605: PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2606: PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2607: PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2608: if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));
2610: /* if time integration scheme has a starting method, call it */
2611: PetscTryTypeMethod(ts, startingmethod);
2613: ts->setupcalled = PETSC_TRUE;
2614: PetscFunctionReturn(PETSC_SUCCESS);
2615: }
2617: /*@
2618: 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
2620: Collective
2622: Input Parameter:
2623: . ts - the `TS` context obtained from `TSCreate()`
2625: Level: developer
2627: Notes:
2628: Any options set on the `TS` object, including those set with `TSSetFromOptions()` remain.
2630: See also `TSSetResize()` to change the size of the system being integrated (for example by adaptive mesh refinement) during the time integration.
2632: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSetResize()`
2633: @*/
2634: PetscErrorCode TSReset(TS ts)
2635: {
2636: TS_RHSSplitLink ilink = ts->tsrhssplit, next;
2638: PetscFunctionBegin;
2641: PetscTryTypeMethod(ts, reset);
2642: if (ts->snes) PetscCall(SNESReset(ts->snes));
2643: if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));
2645: PetscCall(MatDestroy(&ts->Arhs));
2646: PetscCall(MatDestroy(&ts->Brhs));
2647: PetscCall(VecDestroy(&ts->Frhs));
2648: PetscCall(VecDestroy(&ts->vec_sol));
2649: PetscCall(VecDestroy(&ts->vec_sol0));
2650: PetscCall(VecDestroy(&ts->vec_dot));
2651: PetscCall(VecDestroy(&ts->vatol));
2652: PetscCall(VecDestroy(&ts->vrtol));
2653: PetscCall(VecDestroyVecs(ts->nwork, &ts->work));
2655: PetscCall(MatDestroy(&ts->Jacprhs));
2656: PetscCall(MatDestroy(&ts->Jacp));
2657: if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2658: if (ts->quadraturets) {
2659: PetscCall(TSReset(ts->quadraturets));
2660: PetscCall(VecDestroy(&ts->vec_costintegrand));
2661: }
2662: while (ilink) {
2663: next = ilink->next;
2664: PetscCall(TSDestroy(&ilink->ts));
2665: PetscCall(PetscFree(ilink->splitname));
2666: PetscCall(ISDestroy(&ilink->is));
2667: PetscCall(PetscFree(ilink));
2668: ilink = next;
2669: }
2670: ts->tsrhssplit = NULL;
2671: ts->num_rhs_splits = 0;
2672: if (ts->eval_times) {
2673: PetscCall(PetscFree(ts->eval_times->time_points));
2674: PetscCall(PetscFree(ts->eval_times->sol_times));
2675: PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
2676: PetscCall(PetscFree(ts->eval_times));
2677: }
2678: ts->rhsjacobian.time = PETSC_MIN_REAL;
2679: ts->rhsjacobian.scale = 1.0;
2680: ts->ijacobian.shift = 1.0;
2681: ts->setupcalled = PETSC_FALSE;
2682: PetscFunctionReturn(PETSC_SUCCESS);
2683: }
2685: static PetscErrorCode TSResizeReset(TS);
2687: /*@
2688: TSDestroy - Destroys the timestepper context that was created
2689: with `TSCreate()`.
2691: Collective
2693: Input Parameter:
2694: . ts - the `TS` context obtained from `TSCreate()`
2696: Level: beginner
2698: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2699: @*/
2700: PetscErrorCode TSDestroy(TS *ts)
2701: {
2702: PetscFunctionBegin;
2703: if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2705: if (--((PetscObject)*ts)->refct > 0) {
2706: *ts = NULL;
2707: PetscFunctionReturn(PETSC_SUCCESS);
2708: }
2710: PetscCall(TSReset(*ts));
2711: PetscCall(TSAdjointReset(*ts));
2712: if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2713: PetscCall(TSResizeReset(*ts));
2715: /* if memory was published with SAWs then destroy it */
2716: PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2717: PetscTryTypeMethod(*ts, destroy);
2719: PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));
2721: PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2722: PetscCall(TSEventDestroy(&(*ts)->event));
2724: PetscCall(SNESDestroy(&(*ts)->snes));
2725: PetscCall(SNESDestroy(&(*ts)->snesrhssplit));
2726: PetscCall(DMDestroy(&(*ts)->dm));
2727: PetscCall(TSMonitorCancel(*ts));
2728: PetscCall(TSAdjointMonitorCancel(*ts));
2730: PetscCall(TSDestroy(&(*ts)->quadraturets));
2731: PetscCall(PetscHeaderDestroy(ts));
2732: PetscFunctionReturn(PETSC_SUCCESS);
2733: }
2735: /*@
2736: TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the `TS` timestepping context
2738: Collective
2740: Input Parameters:
2741: + ts - the `TS` context obtained from `TSCreate()`
2742: - snes - the nonlinear solver context
2744: Level: developer
2746: Note:
2747: Most users should obtain the `SNES` by calling `TSGetSNES()` rather than setting it with this function.
2749: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetKSP()`, `TSIsImplicit()`, `TSGetSNES()`
2750: @*/
2751: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2752: {
2753: PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);
2755: PetscFunctionBegin;
2758: PetscCall(PetscObjectReference((PetscObject)snes));
2759: PetscCall(SNESDestroy(&ts->snes));
2760: ts->snes = snes;
2761: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2762: PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2763: if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2764: PetscFunctionReturn(PETSC_SUCCESS);
2765: }
2767: /*@
2768: TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2769: a `TS` (timestepper) context.
2771: Not Collective, but `snes` is parallel if `ts` is parallel
2773: Input Parameter:
2774: . ts - the `TS` context obtained from `TSCreate()`
2776: Output Parameter:
2777: . snes - the nonlinear solver context
2779: Level: beginner
2781: Notes:
2782: The user can then directly manipulate the `SNES` context to set various
2783: options, etc. Likewise, the user can then extract and manipulate the
2784: `KSP`, and `PC` contexts as well.
2786: For linear problems, use `TSGetKSP()`.
2788: For integrators that do not use `SNES` (that is, explicit methods),
2789: the `snes` exists but is not used. Use `TSIsImplicit()` to determine if the
2790: method is implicit and uses `snes`.
2792: Developer Note:
2793: `TS` manages the life-cycle of the `SNES` object for all `TSType` for the life-time of the `TS` object,
2794: even explicit methods that do not use `SNES`. This is so that `SNES` options are retained between changes to the `TSType` with `TSSetType()`.
2796: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetKSP()`, `TSIsImplicit()`
2797: @*/
2798: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2799: {
2800: PetscFunctionBegin;
2802: PetscAssertPointer(snes, 2);
2803: if (!ts->snes) {
2804: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2805: PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2806: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2807: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2808: if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2809: if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2810: }
2811: *snes = ts->snes;
2812: PetscFunctionReturn(PETSC_SUCCESS);
2813: }
2815: /*@
2816: TSIsImplicit - Indicates if a `TS` represents an implicit integrator that uses `SNES`
2818: Not Collective
2820: Input Parameter:
2821: . ts - the `TS` context obtained from `TSCreate()`
2823: Output Parameter:
2824: . isimplicit - if the integrator is implicit and uses either `SNES` or `KSP`
2826: Level: beginner
2828: Note:
2829: For integrators that do not use `SNES` (that is, explicit methods), `snes` exists but is not used.
2831: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetKSP()`, `TSGetSNES()`
2832: @*/
2833: PetscErrorCode TSIsImplicit(TS ts, PetscBool *isimplicit)
2834: {
2835: PetscFunctionBegin;
2837: PetscAssertPointer(isimplicit, 2);
2838: *isimplicit = ts->usessnes;
2839: PetscFunctionReturn(PETSC_SUCCESS);
2840: }
2842: /*@
2843: TSGetKSP - Returns the `KSP` (linear solver) associated with
2844: a `TS` (timestepper) context.
2846: Not Collective, but `ksp` is parallel if `ts` is parallel
2848: Input Parameter:
2849: . ts - the `TS` context obtained from `TSCreate()`
2851: Output Parameter:
2852: . ksp - the nonlinear solver context
2854: Level: beginner
2856: Notes:
2857: The user can then directly manipulate the `KSP` context to set various
2858: options, etc. Likewise, the user can then extract and manipulate the
2859: `PC` context as well.
2861: For nonlinear problems (`TS_NONLINEAR`), use `TSGetSNES()` followed by `SNESGetKSP()`.
2863: For integrators that do not use `KSP` (that is, explicit methods),
2864: `TSGetKSP()` returns a `ksp` that is not used. Use `TSIsImplicit()` to determine if
2865: the `ksp` is actually used.
2867: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`, `TSIsImplicit()`
2868: @*/
2869: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2870: {
2871: SNES snes;
2873: PetscFunctionBegin;
2875: PetscAssertPointer(ksp, 2);
2876: PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2877: PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "For linear problems only; use TSGetSNES() then SNESGetKSP()");
2878: PetscCall(TSGetSNES(ts, &snes));
2879: PetscCall(SNESGetKSP(snes, ksp));
2880: PetscFunctionReturn(PETSC_SUCCESS);
2881: }
2883: /* ----------- Routines to set solver parameters ---------- */
2885: /*@
2886: TSSetMaxSteps - Sets the maximum number of steps to use.
2888: Logically Collective
2890: Input Parameters:
2891: + ts - the `TS` context obtained from `TSCreate()`
2892: - maxsteps - maximum number of steps to use
2894: Options Database Key:
2895: . -ts_max_steps maxsteps - Sets maxsteps
2897: Level: intermediate
2899: Note:
2900: Use `PETSC_DETERMINE` to reset the maximum number of steps to the default from when the object's type was set
2902: The default maximum number of steps is 5,000
2904: Fortran Note:
2905: Use `PETSC_DETERMINE_INTEGER`
2907: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2908: @*/
2909: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2910: {
2911: PetscFunctionBegin;
2914: if (maxsteps == PETSC_DETERMINE) {
2915: ts->max_steps = ts->default_max_steps;
2916: } else {
2917: PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2918: ts->max_steps = maxsteps;
2919: }
2920: PetscFunctionReturn(PETSC_SUCCESS);
2921: }
2923: /*@
2924: TSGetMaxSteps - Gets the maximum number of steps to use.
2926: Not Collective
2928: Input Parameter:
2929: . ts - the `TS` context obtained from `TSCreate()`
2931: Output Parameter:
2932: . maxsteps - maximum number of steps to use
2934: Level: advanced
2936: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2937: @*/
2938: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2939: {
2940: PetscFunctionBegin;
2942: PetscAssertPointer(maxsteps, 2);
2943: *maxsteps = ts->max_steps;
2944: PetscFunctionReturn(PETSC_SUCCESS);
2945: }
2947: /*@
2948: TSSetRunSteps - Sets the maximum number of steps to take in each call to `TSSolve()`.
2950: If the step count when `TSSolve()` is `start_step`, this will stop the simulation once `current_step - start_step >= run_steps`.
2951: Comparatively, `TSSetMaxSteps()` will stop if `current_step >= max_steps`.
2952: The simulation will stop when either condition is reached.
2954: Logically Collective
2956: Input Parameters:
2957: + ts - the `TS` context obtained from `TSCreate()`
2958: - runsteps - maximum number of steps to take in each call to `TSSolve()`;
2960: Options Database Key:
2961: . -ts_run_steps runsteps - Sets runsteps
2963: Level: intermediate
2965: Note:
2966: The default is `PETSC_UNLIMITED`
2968: .seealso: [](ch_ts), `TS`, `TSGetRunSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`, `TSSetMaxSteps()`
2969: @*/
2970: PetscErrorCode TSSetRunSteps(TS ts, PetscInt runsteps)
2971: {
2972: PetscFunctionBegin;
2975: if (runsteps == PETSC_DETERMINE) {
2976: ts->run_steps = PETSC_UNLIMITED;
2977: } else {
2978: 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");
2979: ts->run_steps = runsteps;
2980: }
2981: PetscFunctionReturn(PETSC_SUCCESS);
2982: }
2984: /*@
2985: TSGetRunSteps - Gets the maximum number of steps to take in each call to `TSSolve()`.
2987: Not Collective
2989: Input Parameter:
2990: . ts - the `TS` context obtained from `TSCreate()`
2992: Output Parameter:
2993: . runsteps - maximum number of steps to take in each call to `TSSolve`.
2995: Level: advanced
2997: .seealso: [](ch_ts), `TS`, `TSSetRunSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`, `TSGetMaxSteps()`
2998: @*/
2999: PetscErrorCode TSGetRunSteps(TS ts, PetscInt *runsteps)
3000: {
3001: PetscFunctionBegin;
3003: PetscAssertPointer(runsteps, 2);
3004: *runsteps = ts->run_steps;
3005: PetscFunctionReturn(PETSC_SUCCESS);
3006: }
3008: /*@
3009: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
3011: Logically Collective
3013: Input Parameters:
3014: + ts - the `TS` context obtained from `TSCreate()`
3015: - maxtime - final time to step to
3017: Options Database Key:
3018: . -ts_max_time maxtime - Sets maxtime
3020: Level: intermediate
3022: Notes:
3023: Use `PETSC_DETERMINE` to reset the maximum time to the default from when the object's type was set
3025: The default maximum time is 5.0
3027: Fortran Note:
3028: Use `PETSC_DETERMINE_REAL`
3030: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
3031: @*/
3032: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
3033: {
3034: PetscFunctionBegin;
3037: if (maxtime == PETSC_DETERMINE) {
3038: ts->max_time = ts->default_max_time;
3039: } else {
3040: ts->max_time = maxtime;
3041: }
3042: PetscFunctionReturn(PETSC_SUCCESS);
3043: }
3045: /*@
3046: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
3048: Not Collective
3050: Input Parameter:
3051: . ts - the `TS` context obtained from `TSCreate()`
3053: Output Parameter:
3054: . maxtime - final time to step to
3056: Level: advanced
3058: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
3059: @*/
3060: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
3061: {
3062: PetscFunctionBegin;
3064: PetscAssertPointer(maxtime, 2);
3065: *maxtime = ts->max_time;
3066: PetscFunctionReturn(PETSC_SUCCESS);
3067: }
3069: // PetscClangLinter pragma disable: -fdoc-*
3070: /*@
3071: TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.
3073: Level: deprecated
3075: @*/
3076: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
3077: {
3078: PetscFunctionBegin;
3080: PetscCall(TSSetTime(ts, initial_time));
3081: PetscCall(TSSetTimeStep(ts, time_step));
3082: PetscFunctionReturn(PETSC_SUCCESS);
3083: }
3085: // PetscClangLinter pragma disable: -fdoc-*
3086: /*@
3087: TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.
3089: Level: deprecated
3091: @*/
3092: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
3093: {
3094: PetscFunctionBegin;
3096: if (maxsteps) {
3097: PetscAssertPointer(maxsteps, 2);
3098: *maxsteps = ts->max_steps;
3099: }
3100: if (maxtime) {
3101: PetscAssertPointer(maxtime, 3);
3102: *maxtime = ts->max_time;
3103: }
3104: PetscFunctionReturn(PETSC_SUCCESS);
3105: }
3107: // PetscClangLinter pragma disable: -fdoc-*
3108: /*@
3109: TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.
3111: Level: deprecated
3113: @*/
3114: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
3115: {
3116: PetscFunctionBegin;
3117: if (maxsteps != PETSC_CURRENT) PetscCall(TSSetMaxSteps(ts, maxsteps));
3118: if (maxtime != (PetscReal)PETSC_CURRENT) PetscCall(TSSetMaxTime(ts, maxtime));
3119: PetscFunctionReturn(PETSC_SUCCESS);
3120: }
3122: // PetscClangLinter pragma disable: -fdoc-*
3123: /*@
3124: TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.
3126: Level: deprecated
3128: @*/
3129: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
3130: {
3131: return TSGetStepNumber(ts, steps);
3132: }
3134: // PetscClangLinter pragma disable: -fdoc-*
3135: /*@
3136: TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.
3138: Level: deprecated
3140: @*/
3141: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
3142: {
3143: return TSGetStepNumber(ts, steps);
3144: }
3146: /*@
3147: TSSetSolution - Sets the initial solution vector
3148: for use by the `TS` routines.
3150: Logically Collective
3152: Input Parameters:
3153: + ts - the `TS` context obtained from `TSCreate()`
3154: - u - the solution vector
3156: Level: beginner
3158: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
3159: @*/
3160: PetscErrorCode TSSetSolution(TS ts, Vec u)
3161: {
3162: DM dm;
3164: PetscFunctionBegin;
3167: PetscCall(PetscObjectReference((PetscObject)u));
3168: PetscCall(VecDestroy(&ts->vec_sol));
3169: ts->vec_sol = u;
3171: PetscCall(TSGetDM(ts, &dm));
3172: PetscCall(DMShellSetGlobalVector(dm, u));
3173: PetscFunctionReturn(PETSC_SUCCESS);
3174: }
3176: /*@
3177: TSSetPreStep - Sets the general-purpose function
3178: called once at the beginning of each time step.
3180: Logically Collective
3182: Input Parameters:
3183: + ts - The `TS` context obtained from `TSCreate()`
3184: - func - The function
3186: Calling sequence of `func`:
3187: . ts - the `TS` context
3189: Level: intermediate
3191: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3192: @*/
3193: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
3194: {
3195: PetscFunctionBegin;
3197: ts->prestep = func;
3198: PetscFunctionReturn(PETSC_SUCCESS);
3199: }
3201: /*@
3202: TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`
3204: Collective
3206: Input Parameter:
3207: . ts - The `TS` context obtained from `TSCreate()`
3209: Level: developer
3211: Note:
3212: `TSPreStep()` is typically used within time stepping implementations,
3213: so most users would not generally call this routine themselves.
3215: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3216: @*/
3217: PetscErrorCode TSPreStep(TS ts)
3218: {
3219: PetscFunctionBegin;
3221: if (ts->prestep) {
3222: Vec U;
3223: PetscObjectId idprev;
3224: PetscBool sameObject;
3225: PetscObjectState sprev, spost;
3227: PetscCall(TSGetSolution(ts, &U));
3228: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3229: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3230: PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3231: PetscCall(TSGetSolution(ts, &U));
3232: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3233: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3234: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3235: }
3236: PetscFunctionReturn(PETSC_SUCCESS);
3237: }
3239: /*@
3240: TSSetPreStage - Sets the general-purpose function
3241: called once at the beginning of each stage.
3243: Logically Collective
3245: Input Parameters:
3246: + ts - The `TS` context obtained from `TSCreate()`
3247: - func - The function
3249: Calling sequence of `func`:
3250: + ts - the `TS` context
3251: - stagetime - the stage time
3253: Level: intermediate
3255: Note:
3256: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3257: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3258: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3260: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3261: @*/
3262: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime))
3263: {
3264: PetscFunctionBegin;
3266: ts->prestage = func;
3267: PetscFunctionReturn(PETSC_SUCCESS);
3268: }
3270: /*@
3271: TSSetPostStage - Sets the general-purpose function
3272: called once at the end of each stage.
3274: Logically Collective
3276: Input Parameters:
3277: + ts - The `TS` context obtained from `TSCreate()`
3278: - func - The function
3280: Calling sequence of `func`:
3281: + ts - the `TS` context
3282: . stagetime - the stage time
3283: . stageindex - the stage index
3284: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3286: Level: intermediate
3288: Note:
3289: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3290: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3291: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3293: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3294: @*/
3295: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
3296: {
3297: PetscFunctionBegin;
3299: ts->poststage = func;
3300: PetscFunctionReturn(PETSC_SUCCESS);
3301: }
3303: /*@
3304: TSSetPostEvaluate - Sets the general-purpose function
3305: called at the end of each step evaluation.
3307: Logically Collective
3309: Input Parameters:
3310: + ts - The `TS` context obtained from `TSCreate()`
3311: - func - The function
3313: Calling sequence of `func`:
3314: . ts - the `TS` context
3316: Level: intermediate
3318: Note:
3319: The function set by `TSSetPostEvaluate()` is called after the solution is evaluated, or after the step rollback.
3320: Inside the `func` callback, the solution vector can be obtained with `TSGetSolution()`, and modified, if need be.
3321: The time step can be obtained with `TSGetTimeStep()`, and the time at the start of the step - via `TSGetTime()`.
3322: The potential changes to the solution vector introduced by event handling (`postevent()`) are not relevant for `TSSetPostEvaluate()`,
3323: but are relevant for `TSSetPostStep()`, according to the function call scheme in `TSSolve()`, as shown below
3324: .vb
3325: ...
3326: Step()
3327: PostEvaluate()
3328: EventHandling()
3329: step_rollback ? PostEvaluate() : PostStep()
3330: ...
3331: .ve
3332: where EventHandling() may result in one of the following three outcomes
3333: .vb
3334: (1) | successful step | solution intact
3335: (2) | successful step | solution modified by `postevent()`
3336: (3) | step_rollback | solution rolled back
3337: .ve
3339: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3340: @*/
3341: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
3342: {
3343: PetscFunctionBegin;
3345: ts->postevaluate = func;
3346: PetscFunctionReturn(PETSC_SUCCESS);
3347: }
3349: /*@
3350: TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`
3352: Collective
3354: Input Parameters:
3355: + ts - The `TS` context obtained from `TSCreate()`
3356: - stagetime - The absolute time of the current stage
3358: Level: developer
3360: Note:
3361: `TSPreStage()` is typically used within time stepping implementations,
3362: most users would not generally call this routine themselves.
3364: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3365: @*/
3366: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3367: {
3368: PetscFunctionBegin;
3370: if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3371: PetscFunctionReturn(PETSC_SUCCESS);
3372: }
3374: /*@
3375: TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`
3377: Collective
3379: Input Parameters:
3380: + ts - The `TS` context obtained from `TSCreate()`
3381: . stagetime - The absolute time of the current stage
3382: . stageindex - Stage number
3383: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3385: Level: developer
3387: Note:
3388: `TSPostStage()` is typically used within time stepping implementations,
3389: most users would not generally call this routine themselves.
3391: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3392: @*/
3393: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec Y[])
3394: {
3395: PetscFunctionBegin;
3397: if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3398: PetscFunctionReturn(PETSC_SUCCESS);
3399: }
3401: /*@
3402: TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`
3404: Collective
3406: Input Parameter:
3407: . ts - The `TS` context obtained from `TSCreate()`
3409: Level: developer
3411: Note:
3412: `TSPostEvaluate()` is typically used within time stepping implementations,
3413: most users would not generally call this routine themselves.
3415: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3416: @*/
3417: PetscErrorCode TSPostEvaluate(TS ts)
3418: {
3419: PetscFunctionBegin;
3421: if (ts->postevaluate) {
3422: Vec U;
3423: PetscObjectState sprev, spost;
3425: PetscCall(TSGetSolution(ts, &U));
3426: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3427: PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3428: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3429: if (sprev != spost) PetscCall(TSRestartStep(ts));
3430: }
3431: PetscFunctionReturn(PETSC_SUCCESS);
3432: }
3434: /*@
3435: TSSetPostStep - Sets the general-purpose function
3436: called once at the end of each successful time step.
3438: Logically Collective
3440: Input Parameters:
3441: + ts - The `TS` context obtained from `TSCreate()`
3442: - func - The function
3444: Calling sequence of `func`:
3445: . ts - the `TS` context
3447: Level: intermediate
3449: Note:
3450: The function set by `TSSetPostStep()` is called after each successful step. If the event handler locates an event at the
3451: given step, and `postevent()` modifies the solution vector, the solution vector obtained by `TSGetSolution()` inside `func` will
3452: contain the changes. To get the solution without these changes, use `TSSetPostEvaluate()` to set the appropriate callback.
3453: The scheme of the relevant function calls in `TSSolve()` is shown below
3454: .vb
3455: ...
3456: Step()
3457: PostEvaluate()
3458: EventHandling()
3459: step_rollback ? PostEvaluate() : PostStep()
3460: ...
3461: .ve
3462: where EventHandling() may result in one of the following three outcomes
3463: .vb
3464: (1) | successful step | solution intact
3465: (2) | successful step | solution modified by `postevent()`
3466: (3) | step_rollback | solution rolled back
3467: .ve
3469: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3470: @*/
3471: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
3472: {
3473: PetscFunctionBegin;
3475: ts->poststep = func;
3476: PetscFunctionReturn(PETSC_SUCCESS);
3477: }
3479: /*@
3480: TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`
3482: Collective
3484: Input Parameter:
3485: . ts - The `TS` context obtained from `TSCreate()`
3487: Note:
3488: `TSPostStep()` is typically used within time stepping implementations,
3489: so most users would not generally call this routine themselves.
3491: Level: developer
3493: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPostStep()`
3494: @*/
3495: PetscErrorCode TSPostStep(TS ts)
3496: {
3497: PetscFunctionBegin;
3499: if (ts->poststep) {
3500: Vec U;
3501: PetscObjectId idprev;
3502: PetscBool sameObject;
3503: PetscObjectState sprev, spost;
3505: PetscCall(TSGetSolution(ts, &U));
3506: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3507: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3508: PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3509: PetscCall(TSGetSolution(ts, &U));
3510: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3511: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3512: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3513: }
3514: PetscFunctionReturn(PETSC_SUCCESS);
3515: }
3517: /*@
3518: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3520: Collective
3522: Input Parameters:
3523: + ts - time stepping context
3524: - t - time to interpolate to
3526: Output Parameter:
3527: . U - state at given time
3529: Level: intermediate
3531: Developer Notes:
3532: `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3534: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3535: @*/
3536: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3537: {
3538: PetscFunctionBegin;
3541: 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);
3542: PetscUseTypeMethod(ts, interpolate, t, U);
3543: PetscFunctionReturn(PETSC_SUCCESS);
3544: }
3546: /*@
3547: TSStep - Steps one time step
3549: Collective
3551: Input Parameter:
3552: . ts - the `TS` context obtained from `TSCreate()`
3554: Level: developer
3556: Notes:
3557: The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.
3559: The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
3560: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3562: This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
3563: time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.
3565: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3566: @*/
3567: PetscErrorCode TSStep(TS ts)
3568: {
3569: static PetscBool cite = PETSC_FALSE;
3570: PetscReal ptime;
3572: PetscFunctionBegin;
3574: PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3575: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3576: " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3577: " journal = {arXiv e-preprints},\n"
3578: " eprint = {1806.01437},\n"
3579: " archivePrefix = {arXiv},\n"
3580: " year = {2018}\n}\n",
3581: &cite));
3582: PetscCall(TSSetUp(ts));
3583: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3584: if (ts->eval_times)
3585: ts->eval_times->worktol = 0; /* In each step of TSSolve() 'eval_times->worktol' will be meaningfully defined (later) only once:
3586: in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */
3588: 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>");
3589: 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()");
3590: 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");
3592: if (!ts->vec_sol0) PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sol0));
3593: PetscCall(VecCopy(ts->vec_sol, ts->vec_sol0));
3594: ts->time_step0 = ts->time_step;
3596: if (!ts->steps) ts->ptime_prev = ts->ptime;
3597: ptime = ts->ptime;
3599: ts->ptime_prev_rollback = ts->ptime_prev;
3600: ts->reason = TS_CONVERGED_ITERATING;
3602: PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3603: PetscUseTypeMethod(ts, step);
3604: PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));
3606: if (ts->reason >= 0) {
3607: ts->ptime_prev = ptime;
3608: ts->steps++;
3609: ts->steprollback = PETSC_FALSE;
3610: ts->steprestart = PETSC_FALSE;
3611: ts->stepresize = PETSC_FALSE;
3612: }
3614: if (ts->reason < 0 && ts->errorifstepfailed) {
3615: PetscCall(TSMonitorCancel(ts));
3616: if (ts->usessnes && ts->snes) PetscCall(SNESMonitorCancel(ts->snes));
3617: 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]);
3618: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3619: }
3620: PetscFunctionReturn(PETSC_SUCCESS);
3621: }
3623: /*@
3624: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3625: at the end of a time step with a given order of accuracy.
3627: Collective
3629: Input Parameters:
3630: + ts - time stepping context
3631: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3633: Input/Output Parameter:
3634: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3635: on output, the actual order of the error evaluation
3637: Output Parameter:
3638: . wlte - the weighted local truncation error norm
3640: Level: advanced
3642: Note:
3643: If the timestepper cannot evaluate the error in a particular step
3644: (eg. in the first step or restart steps after event handling),
3645: this routine returns wlte=-1.0 .
3647: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3648: @*/
3649: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3650: {
3651: PetscFunctionBegin;
3655: if (order) PetscAssertPointer(order, 3);
3657: PetscAssertPointer(wlte, 4);
3658: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3659: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3660: PetscFunctionReturn(PETSC_SUCCESS);
3661: }
3663: /*@
3664: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3666: Collective
3668: Input Parameters:
3669: + ts - time stepping context
3670: . order - desired order of accuracy
3671: - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)
3673: Output Parameter:
3674: . U - state at the end of the current step
3676: Level: advanced
3678: Notes:
3679: This function cannot be called until all stages have been evaluated.
3681: 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.
3683: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3684: @*/
3685: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3686: {
3687: PetscFunctionBegin;
3691: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3692: PetscFunctionReturn(PETSC_SUCCESS);
3693: }
3695: /*@
3696: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3698: Not collective
3700: Input Parameter:
3701: . ts - time stepping context
3703: Output Parameter:
3704: . initCondition - The function which computes an initial condition
3706: Calling sequence of `initCondition`:
3707: + ts - The timestepping context
3708: - u - The input vector in which the initial condition is stored
3710: Level: advanced
3712: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3713: @*/
3714: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3715: {
3716: PetscFunctionBegin;
3718: PetscAssertPointer(initCondition, 2);
3719: *initCondition = ts->ops->initcondition;
3720: PetscFunctionReturn(PETSC_SUCCESS);
3721: }
3723: /*@
3724: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3726: Logically collective
3728: Input Parameters:
3729: + ts - time stepping context
3730: - initCondition - The function which computes an initial condition
3732: Calling sequence of `initCondition`:
3733: + ts - The timestepping context
3734: - e - The input vector in which the initial condition is to be stored
3736: Level: advanced
3738: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3739: @*/
3740: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3741: {
3742: PetscFunctionBegin;
3745: ts->ops->initcondition = initCondition;
3746: PetscFunctionReturn(PETSC_SUCCESS);
3747: }
3749: /*@
3750: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3752: Collective
3754: Input Parameters:
3755: + ts - time stepping context
3756: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3758: Level: advanced
3760: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3761: @*/
3762: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3763: {
3764: PetscFunctionBegin;
3767: PetscTryTypeMethod(ts, initcondition, u);
3768: PetscFunctionReturn(PETSC_SUCCESS);
3769: }
3771: /*@
3772: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3774: Not collective
3776: Input Parameter:
3777: . ts - time stepping context
3779: Output Parameter:
3780: . exactError - The function which computes the solution error
3782: Calling sequence of `exactError`:
3783: + ts - The timestepping context
3784: . u - The approximate solution vector
3785: - e - The vector in which the error is stored
3787: Level: advanced
3789: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3790: @*/
3791: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3792: {
3793: PetscFunctionBegin;
3795: PetscAssertPointer(exactError, 2);
3796: *exactError = ts->ops->exacterror;
3797: PetscFunctionReturn(PETSC_SUCCESS);
3798: }
3800: /*@
3801: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3803: Logically collective
3805: Input Parameters:
3806: + ts - time stepping context
3807: - exactError - The function which computes the solution error
3809: Calling sequence of `exactError`:
3810: + ts - The timestepping context
3811: . u - The approximate solution vector
3812: - e - The vector in which the error is stored
3814: Level: advanced
3816: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3817: @*/
3818: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3819: {
3820: PetscFunctionBegin;
3823: ts->ops->exacterror = exactError;
3824: PetscFunctionReturn(PETSC_SUCCESS);
3825: }
3827: /*@
3828: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3830: Collective
3832: Input Parameters:
3833: + ts - time stepping context
3834: . u - The approximate solution
3835: - e - The `Vec` used to store the error
3837: Level: advanced
3839: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3840: @*/
3841: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3842: {
3843: PetscFunctionBegin;
3847: PetscTryTypeMethod(ts, exacterror, u, e);
3848: PetscFunctionReturn(PETSC_SUCCESS);
3849: }
3851: /*@
3852: TSSetResize - Sets the resize callbacks.
3854: Logically Collective
3856: Input Parameters:
3857: + ts - The `TS` context obtained from `TSCreate()`
3858: . rollback - Whether a resize will restart the step
3859: . setup - The setup function
3860: . transfer - The transfer function
3861: - ctx - [optional] The user-defined context
3863: Calling sequence of `setup`:
3864: + ts - the `TS` context
3865: . step - the current step
3866: . time - the current time
3867: . state - the current vector of state
3868: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3869: - ctx - user defined context
3871: Calling sequence of `transfer`:
3872: + ts - the `TS` context
3873: . nv - the number of vectors to be transferred
3874: . vecsin - array of vectors to be transferred
3875: . vecsout - array of transferred vectors
3876: - ctx - user defined context
3878: Notes:
3879: The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
3880: depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
3881: and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
3882: that the size of the discrete problem has changed.
3883: In both cases, the solver will collect the needed vectors that will be
3884: transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
3885: current solution vector, and other vectors needed by the specific solver used.
3886: For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3887: Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3888: will be automatically reset if the sizes are changed and they must be specified again by the user
3889: inside the `transfer` function.
3890: The input and output arrays passed to `transfer` are allocated by PETSc.
3891: Vectors in `vecsout` must be created by the user.
3892: Ownership of vectors in `vecsout` is transferred to PETSc.
3894: Level: advanced
3896: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3897: @*/
3898: PetscErrorCode TSSetResize(TS ts, PetscBool rollback, PetscErrorCode (*setup)(TS ts, PetscInt step, PetscReal time, Vec state, PetscBool *resize, PetscCtx ctx), PetscErrorCode (*transfer)(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], PetscCtx ctx), PetscCtx ctx)
3899: {
3900: PetscFunctionBegin;
3902: ts->resizerollback = rollback;
3903: ts->resizesetup = setup;
3904: ts->resizetransfer = transfer;
3905: ts->resizectx = ctx;
3906: PetscFunctionReturn(PETSC_SUCCESS);
3907: }
3909: /*
3910: TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.
3912: Collective
3914: Input Parameters:
3915: + ts - The `TS` context obtained from `TSCreate()`
3916: - 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.
3918: Level: developer
3920: Note:
3921: `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3922: used within time stepping implementations,
3923: so most users would not generally call this routine themselves.
3925: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3926: @*/
3927: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3928: {
3929: PetscFunctionBegin;
3931: PetscTryTypeMethod(ts, resizeregister, flg);
3932: /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3933: PetscFunctionReturn(PETSC_SUCCESS);
3934: }
3936: static PetscErrorCode TSResizeReset(TS ts)
3937: {
3938: PetscFunctionBegin;
3940: PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3941: PetscFunctionReturn(PETSC_SUCCESS);
3942: }
3944: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3945: {
3946: PetscFunctionBegin;
3949: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3950: if (ts->resizetransfer) {
3951: PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3952: PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3953: }
3954: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3955: PetscFunctionReturn(PETSC_SUCCESS);
3956: }
3958: /*@
3959: TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.
3961: Collective
3963: Input Parameters:
3964: + ts - The `TS` context obtained from `TSCreate()`
3965: . name - A string identifying the vector
3966: - vec - The vector
3968: Level: developer
3970: Note:
3971: `TSResizeRegisterVec()` is typically used within time stepping implementations,
3972: so most users would not generally call this routine themselves.
3974: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3975: @*/
3976: PetscErrorCode TSResizeRegisterVec(TS ts, const char name[], Vec vec)
3977: {
3978: PetscFunctionBegin;
3980: PetscAssertPointer(name, 2);
3982: PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3983: PetscFunctionReturn(PETSC_SUCCESS);
3984: }
3986: /*@
3987: TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.
3989: Collective
3991: Input Parameters:
3992: + ts - The `TS` context obtained from `TSCreate()`
3993: . name - A string identifying the vector
3994: - vec - The vector
3996: Level: developer
3998: Note:
3999: `TSResizeRetrieveVec()` is typically used within time stepping implementations,
4000: so most users would not generally call this routine themselves.
4002: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
4003: @*/
4004: PetscErrorCode TSResizeRetrieveVec(TS ts, const char name[], Vec *vec)
4005: {
4006: PetscFunctionBegin;
4008: PetscAssertPointer(name, 2);
4009: PetscAssertPointer(vec, 3);
4010: PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
4011: PetscFunctionReturn(PETSC_SUCCESS);
4012: }
4014: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
4015: {
4016: PetscInt cnt;
4017: PetscObjectList tmp;
4018: Vec *vecsin = NULL;
4019: const char **namesin = NULL;
4021: PetscFunctionBegin;
4022: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
4023: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
4024: if (names) PetscCall(PetscMalloc1(cnt, &namesin));
4025: if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
4026: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
4027: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
4028: if (vecs) vecsin[cnt] = (Vec)tmp->obj;
4029: if (names) namesin[cnt] = tmp->name;
4030: cnt++;
4031: }
4032: }
4033: if (nv) *nv = cnt;
4034: if (names) *names = namesin;
4035: if (vecs) *vecs = vecsin;
4036: PetscFunctionReturn(PETSC_SUCCESS);
4037: }
4039: /*@
4040: TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`
4042: Collective
4044: Input Parameter:
4045: . ts - The `TS` context obtained from `TSCreate()`
4047: Level: developer
4049: Note:
4050: `TSResize()` is typically used within time stepping implementations,
4051: so most users would not generally call this routine themselves.
4053: .seealso: [](ch_ts), `TS`, `TSSetResize()`
4054: @*/
4055: PetscErrorCode TSResize(TS ts)
4056: {
4057: PetscInt nv = 0;
4058: const char **names = NULL;
4059: Vec *vecsin = NULL;
4060: const char *solname = "ts:vec_sol";
4062: PetscFunctionBegin;
4064: if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
4065: if (ts->resizesetup) {
4066: PetscCall(VecLockReadPush(ts->vec_sol));
4067: PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &ts->stepresize, ts->resizectx));
4068: PetscCall(VecLockReadPop(ts->vec_sol));
4069: if (ts->stepresize) {
4070: if (ts->resizerollback) {
4071: PetscCall(TSRollBack(ts));
4072: ts->time_step = ts->time_step0;
4073: }
4074: PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
4075: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
4076: }
4077: }
4079: PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
4080: if (nv) {
4081: Vec *vecsout, vecsol;
4083: /* Reset internal objects */
4084: PetscCall(TSReset(ts));
4086: /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
4087: PetscCall(PetscCalloc1(nv, &vecsout));
4088: PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
4089: for (PetscInt i = 0; i < nv; i++) {
4090: const char *name;
4091: char *oname;
4093: PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
4094: PetscCall(PetscStrallocpy(name, &oname));
4095: PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
4096: if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
4097: PetscCall(PetscFree(oname));
4098: PetscCall(VecDestroy(&vecsout[i]));
4099: }
4100: PetscCall(PetscFree(vecsout));
4101: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */
4103: PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
4104: if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
4105: PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
4106: }
4108: PetscCall(PetscFree(names));
4109: PetscCall(PetscFree(vecsin));
4110: PetscCall(TSResizeReset(ts));
4111: PetscFunctionReturn(PETSC_SUCCESS);
4112: }
4114: /*@
4115: TSSolve - Steps the requested number of timesteps.
4117: Collective
4119: Input Parameters:
4120: + ts - the `TS` context obtained from `TSCreate()`
4121: - u - the solution vector (can be `NULL` if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
4122: otherwise it must contain the initial conditions and will contain the solution at the final requested time
4124: Level: beginner
4126: Note:
4127: The final time returned by this function may be different from the time of the internally
4128: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
4129: stepped over the final time.
4131: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
4132: @*/
4133: PetscErrorCode TSSolve(TS ts, Vec u)
4134: {
4135: Vec solution;
4137: PetscFunctionBegin;
4141: PetscCall(TSSetExactFinalTimeDefault(ts));
4142: 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 */
4143: if (!ts->vec_sol || u == ts->vec_sol) {
4144: PetscCall(VecDuplicate(u, &solution));
4145: PetscCall(TSSetSolution(ts, solution));
4146: PetscCall(VecDestroy(&solution)); /* grant ownership */
4147: }
4148: PetscCall(VecCopy(u, ts->vec_sol));
4149: PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
4150: } else if (u) PetscCall(TSSetSolution(ts, u));
4151: PetscCall(TSSetUp(ts));
4152: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
4154: 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>");
4155: 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()");
4156: 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");
4157: 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");
4159: if (ts->eval_times) {
4160: if (!ts->eval_times->sol_vecs) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
4161: for (PetscInt i = 0; i < ts->eval_times->num_time_points; i++) {
4162: PetscBool is_close = PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[i], ts->eval_times->reltol * ts->time_step + ts->eval_times->abstol, 0);
4163: if (ts->ptime <= ts->eval_times->time_points[i] || is_close) {
4164: ts->eval_times->time_point_idx = i;
4166: PetscBool is_ptime_in_sol_times = PETSC_FALSE; // If current solution has already been saved, we should not save it again
4167: 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);
4168: if (is_close && !is_ptime_in_sol_times) {
4169: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4170: ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4171: ts->eval_times->sol_idx++;
4172: ts->eval_times->time_point_idx++;
4173: }
4174: break;
4175: }
4176: }
4177: }
4179: if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));
4181: /* reset number of steps only when the step is not restarted. ARKIMEX
4182: restarts the step after an event. Resetting these counters in such case causes
4183: TSTrajectory to incorrectly save the output files
4184: */
4185: /* reset time step and iteration counters */
4186: if (!ts->steps) {
4187: ts->ksp_its = 0;
4188: ts->snes_its = 0;
4189: ts->num_snes_failures = 0;
4190: ts->reject = 0;
4191: ts->steprestart = PETSC_TRUE;
4192: ts->steprollback = PETSC_FALSE;
4193: ts->stepresize = PETSC_FALSE;
4194: ts->rhsjacobian.time = PETSC_MIN_REAL;
4195: }
4197: /* make sure initial time step does not overshoot final time or the next point in evaluation times */
4198: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
4199: PetscReal maxdt;
4200: PetscReal dt = ts->time_step;
4202: if (ts->eval_times) maxdt = ts->eval_times->time_points[ts->eval_times->time_point_idx] - ts->ptime;
4203: else maxdt = ts->max_time - ts->ptime;
4204: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
4205: }
4206: ts->reason = TS_CONVERGED_ITERATING;
4208: {
4209: PetscViewer viewer;
4210: PetscViewerFormat format;
4211: PetscBool flg;
4212: static PetscBool incall = PETSC_FALSE;
4214: if (!incall) {
4215: /* Estimate the convergence rate of the time discretization */
4216: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
4217: if (flg) {
4218: PetscConvEst conv;
4219: DM dm;
4220: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4221: PetscInt Nf;
4222: PetscBool checkTemporal = PETSC_TRUE;
4224: incall = PETSC_TRUE;
4225: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
4226: PetscCall(TSGetDM(ts, &dm));
4227: PetscCall(DMGetNumFields(dm, &Nf));
4228: PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
4229: PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
4230: PetscCall(PetscConvEstUseTS(conv, checkTemporal));
4231: PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
4232: PetscCall(PetscConvEstSetFromOptions(conv));
4233: PetscCall(PetscConvEstSetUp(conv));
4234: PetscCall(PetscConvEstGetConvRate(conv, alpha));
4235: PetscCall(PetscViewerPushFormat(viewer, format));
4236: PetscCall(PetscConvEstRateView(conv, alpha, viewer));
4237: PetscCall(PetscViewerPopFormat(viewer));
4238: PetscCall(PetscViewerDestroy(&viewer));
4239: PetscCall(PetscConvEstDestroy(&conv));
4240: PetscCall(PetscFree(alpha));
4241: incall = PETSC_FALSE;
4242: }
4243: }
4244: }
4246: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));
4248: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4249: PetscUseTypeMethod(ts, solve);
4250: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4251: ts->solvetime = ts->ptime;
4252: solution = ts->vec_sol;
4253: } else { /* Step the requested number of timesteps. */
4254: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4255: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4257: if (!ts->steps) {
4258: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4259: PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4260: }
4262: ts->start_step = ts->steps; // records starting step
4263: while (!ts->reason) {
4264: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4265: if (!ts->steprollback || (ts->stepresize && ts->resizerollback)) PetscCall(TSPreStep(ts));
4266: PetscCall(TSStep(ts));
4267: if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4268: if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4269: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4270: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4271: PetscCall(TSForwardCostIntegral(ts));
4272: if (ts->reason >= 0) ts->steps++;
4273: }
4274: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4275: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4276: PetscCall(TSForwardStep(ts));
4277: if (ts->reason >= 0) ts->steps++;
4278: }
4279: PetscCall(TSPostEvaluate(ts));
4280: 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. */
4281: if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4282: if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
4283: /* check convergence */
4284: if (!ts->reason) {
4285: if ((ts->steps - ts->start_step) >= ts->run_steps) ts->reason = TS_CONVERGED_ITS;
4286: else if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4287: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4288: }
4289: if (!ts->steprollback) {
4290: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4291: PetscCall(TSPostStep(ts));
4292: if (!ts->resizerollback) PetscCall(TSResize(ts));
4294: if (ts->eval_times && ts->eval_times->time_point_idx < ts->eval_times->num_time_points && ts->reason >= 0) {
4295: PetscCheck(ts->eval_times->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(eval_times->worktol > 0) in TSSolve()");
4296: if (PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[ts->eval_times->time_point_idx], ts->eval_times->worktol, 0)) {
4297: ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4298: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4299: ts->eval_times->sol_idx++;
4300: ts->eval_times->time_point_idx++;
4301: }
4302: }
4303: }
4304: }
4305: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4307: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4308: if (!u) u = ts->vec_sol;
4309: PetscCall(TSInterpolate(ts, ts->max_time, u));
4310: ts->solvetime = ts->max_time;
4311: solution = u;
4312: PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4313: } else {
4314: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4315: ts->solvetime = ts->ptime;
4316: solution = ts->vec_sol;
4317: }
4318: }
4320: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4321: PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4322: PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4323: if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4324: PetscFunctionReturn(PETSC_SUCCESS);
4325: }
4327: /*@
4328: TSGetTime - Gets the time of the most recently completed step.
4330: Not Collective
4332: Input Parameter:
4333: . ts - the `TS` context obtained from `TSCreate()`
4335: Output Parameter:
4336: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
4338: Level: beginner
4340: Note:
4341: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
4342: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
4344: .seealso: [](ch_ts), `TS`, `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4345: @*/
4346: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4347: {
4348: PetscFunctionBegin;
4350: PetscAssertPointer(t, 2);
4351: *t = ts->ptime;
4352: PetscFunctionReturn(PETSC_SUCCESS);
4353: }
4355: /*@
4356: TSGetPrevTime - Gets the starting time of the previously completed step.
4358: Not Collective
4360: Input Parameter:
4361: . ts - the `TS` context obtained from `TSCreate()`
4363: Output Parameter:
4364: . t - the previous time
4366: Level: beginner
4368: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4369: @*/
4370: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4371: {
4372: PetscFunctionBegin;
4374: PetscAssertPointer(t, 2);
4375: *t = ts->ptime_prev;
4376: PetscFunctionReturn(PETSC_SUCCESS);
4377: }
4379: /*@
4380: TSSetTime - Allows one to reset the time.
4382: Logically Collective
4384: Input Parameters:
4385: + ts - the `TS` context obtained from `TSCreate()`
4386: - t - the time
4388: Level: intermediate
4390: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4391: @*/
4392: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4393: {
4394: PetscFunctionBegin;
4397: ts->ptime = t;
4398: PetscFunctionReturn(PETSC_SUCCESS);
4399: }
4401: /*@
4402: TSSetOptionsPrefix - Sets the prefix used for searching for all
4403: TS options in the database.
4405: Logically Collective
4407: Input Parameters:
4408: + ts - The `TS` context
4409: - prefix - The prefix to prepend to all option names
4411: Level: advanced
4413: Note:
4414: A hyphen (-) must NOT be given at the beginning of the prefix name.
4415: The first character of all runtime options is AUTOMATICALLY the
4416: hyphen.
4418: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4419: @*/
4420: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4421: {
4422: SNES snes;
4424: PetscFunctionBegin;
4426: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4427: PetscCall(TSGetSNES(ts, &snes));
4428: PetscCall(SNESSetOptionsPrefix(snes, prefix));
4429: PetscFunctionReturn(PETSC_SUCCESS);
4430: }
4432: /*@
4433: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4434: TS options in the database.
4436: Logically Collective
4438: Input Parameters:
4439: + ts - The `TS` context
4440: - prefix - The prefix to prepend to all option names
4442: Level: advanced
4444: Note:
4445: A hyphen (-) must NOT be given at the beginning of the prefix name.
4446: The first character of all runtime options is AUTOMATICALLY the
4447: hyphen.
4449: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4450: @*/
4451: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4452: {
4453: SNES snes;
4455: PetscFunctionBegin;
4457: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4458: PetscCall(TSGetSNES(ts, &snes));
4459: PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4460: PetscFunctionReturn(PETSC_SUCCESS);
4461: }
4463: /*@
4464: TSGetOptionsPrefix - Sets the prefix used for searching for all
4465: `TS` options in the database.
4467: Not Collective
4469: Input Parameter:
4470: . ts - The `TS` context
4472: Output Parameter:
4473: . prefix - A pointer to the prefix string used
4475: Level: intermediate
4477: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4478: @*/
4479: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4480: {
4481: PetscFunctionBegin;
4483: PetscAssertPointer(prefix, 2);
4484: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4485: PetscFunctionReturn(PETSC_SUCCESS);
4486: }
4488: /*@
4489: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4491: Not Collective, but parallel objects are returned if ts is parallel
4493: Input Parameter:
4494: . ts - The `TS` context obtained from `TSCreate()`
4496: Output Parameters:
4497: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`)
4498: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`)
4499: . func - Function to compute the Jacobian of the RHS (or `NULL`)
4500: - ctx - User-defined context for Jacobian evaluation routine (or `NULL`)
4502: Level: intermediate
4504: Note:
4505: You can pass in `NULL` for any return argument you do not need.
4507: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4508: @*/
4509: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, PetscCtxRt ctx)
4510: {
4511: DM dm;
4513: PetscFunctionBegin;
4514: if (Amat || Pmat) {
4515: SNES snes;
4516: PetscCall(TSGetSNES(ts, &snes));
4517: PetscCall(SNESSetUpMatrices(snes));
4518: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4519: }
4520: PetscCall(TSGetDM(ts, &dm));
4521: PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4522: PetscFunctionReturn(PETSC_SUCCESS);
4523: }
4525: /*@
4526: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4528: Not Collective, but parallel objects are returned if ts is parallel
4530: Input Parameter:
4531: . ts - The `TS` context obtained from `TSCreate()`
4533: Output Parameters:
4534: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4535: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4536: . f - The function to compute the matrices
4537: - ctx - User-defined context for Jacobian evaluation routine
4539: Level: advanced
4541: Note:
4542: You can pass in `NULL` for any return argument you do not need.
4544: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4545: @*/
4546: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, PetscCtxRt ctx)
4547: {
4548: DM dm;
4550: PetscFunctionBegin;
4551: if (Amat || Pmat) {
4552: SNES snes;
4553: PetscCall(TSGetSNES(ts, &snes));
4554: PetscCall(SNESSetUpMatrices(snes));
4555: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4556: }
4557: PetscCall(TSGetDM(ts, &dm));
4558: PetscCall(DMTSGetIJacobian(dm, f, ctx));
4559: PetscFunctionReturn(PETSC_SUCCESS);
4560: }
4562: #include <petsc/private/dmimpl.h>
4563: /*@
4564: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4566: Logically Collective
4568: Input Parameters:
4569: + ts - the `TS` integrator object
4570: - dm - the dm, cannot be `NULL`
4572: Level: intermediate
4574: Notes:
4575: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4576: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4577: different problems using the same function space.
4579: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4580: @*/
4581: PetscErrorCode TSSetDM(TS ts, DM dm)
4582: {
4583: SNES snes;
4584: DMTS tsdm;
4586: PetscFunctionBegin;
4589: PetscCall(PetscObjectReference((PetscObject)dm));
4590: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4591: if (ts->dm->dmts && !dm->dmts) {
4592: PetscCall(DMCopyDMTS(ts->dm, dm));
4593: PetscCall(DMGetDMTS(ts->dm, &tsdm));
4594: /* Grant write privileges to the replacement DM */
4595: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4596: }
4597: PetscCall(DMDestroy(&ts->dm));
4598: }
4599: ts->dm = dm;
4601: PetscCall(TSGetSNES(ts, &snes));
4602: PetscCall(SNESSetDM(snes, dm));
4603: PetscFunctionReturn(PETSC_SUCCESS);
4604: }
4606: /*@
4607: TSGetDM - Gets the `DM` that may be used by some preconditioners
4609: Not Collective
4611: Input Parameter:
4612: . ts - the `TS`
4614: Output Parameter:
4615: . dm - the `DM`
4617: Level: intermediate
4619: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4620: @*/
4621: PetscErrorCode TSGetDM(TS ts, DM *dm)
4622: {
4623: PetscFunctionBegin;
4625: if (!ts->dm) {
4626: PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4627: if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4628: }
4629: *dm = ts->dm;
4630: PetscFunctionReturn(PETSC_SUCCESS);
4631: }
4633: /*@
4634: SNESTSFormFunction - Function to evaluate nonlinear residual defined by an ODE solver algorithm implemented within `TS`
4636: Logically Collective
4638: Input Parameters:
4639: + snes - nonlinear solver
4640: . U - the current state at which to evaluate the residual
4641: - ctx - application context, must be a `TS`
4643: Output Parameter:
4644: . F - the nonlinear residual
4646: Level: developer
4648: Note:
4649: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4650: It is most frequently passed to `MatFDColoringSetFunction()`.
4652: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4653: @*/
4654: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, PetscCtx ctx)
4655: {
4656: TS ts = (TS)ctx;
4658: PetscFunctionBegin;
4663: PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4664: PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4665: PetscFunctionReturn(PETSC_SUCCESS);
4666: }
4668: /*@
4669: SNESTSFormJacobian - Function to evaluate the Jacobian defined by an ODE solver algorithm implemented within `TS`
4671: Collective
4673: Input Parameters:
4674: + snes - nonlinear solver
4675: . U - the current state at which to evaluate the residual
4676: - ctx - application context, must be a `TS`
4678: Output Parameters:
4679: + A - the Jacobian
4680: - B - the matrix used to construct the preconditioner (often the same as `A`)
4682: Level: developer
4684: Note:
4685: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4687: .seealso: [](ch_ts), `SNESSetJacobian()`
4688: @*/
4689: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, PetscCtx ctx)
4690: {
4691: TS ts = (TS)ctx;
4693: PetscFunctionBegin;
4699: PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4700: PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4701: PetscFunctionReturn(PETSC_SUCCESS);
4702: }
4704: /*@
4705: TSComputeRHSFunctionLinear - Evaluate the right-hand side via the user-provided Jacobian, for linear problems Udot = A U only
4707: Collective
4709: Input Parameters:
4710: + ts - time stepping context
4711: . t - time at which to evaluate
4712: . U - state at which to evaluate
4713: - ctx - context
4715: Output Parameter:
4716: . F - right-hand side
4718: Level: intermediate
4720: Note:
4721: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right-hand side for linear problems.
4722: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4724: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4725: @*/
4726: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, PetscCtx ctx)
4727: {
4728: Mat Arhs, Brhs;
4730: PetscFunctionBegin;
4731: PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4732: /* undo the damage caused by shifting */
4733: PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4734: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4735: PetscCall(MatMult(Arhs, U, F));
4736: PetscFunctionReturn(PETSC_SUCCESS);
4737: }
4739: /*@
4740: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
4742: Collective
4744: Input Parameters:
4745: + ts - time stepping context
4746: . t - time at which to evaluate
4747: . U - state at which to evaluate
4748: - ctx - context
4750: Output Parameters:
4751: + A - Jacobian
4752: - B - matrix used to construct the preconditioner, often the same as `A`
4754: Level: intermediate
4756: Note:
4757: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4759: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4760: @*/
4761: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, PetscCtx ctx)
4762: {
4763: PetscFunctionBegin;
4764: PetscFunctionReturn(PETSC_SUCCESS);
4765: }
4767: /*@
4768: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4770: Collective
4772: Input Parameters:
4773: + ts - time stepping context
4774: . t - time at which to evaluate
4775: . U - state at which to evaluate
4776: . Udot - time derivative of state vector
4777: - ctx - context
4779: Output Parameter:
4780: . F - left hand side
4782: Level: intermediate
4784: Notes:
4785: 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
4786: user is required to write their own `TSComputeIFunction()`.
4787: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4788: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4790: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4792: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4793: @*/
4794: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
4795: {
4796: Mat A, B;
4798: PetscFunctionBegin;
4799: PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4800: PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4801: PetscCall(MatMult(A, Udot, F));
4802: PetscFunctionReturn(PETSC_SUCCESS);
4803: }
4805: /*@
4806: TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE
4808: Collective
4810: Input Parameters:
4811: + ts - time stepping context
4812: . t - time at which to evaluate
4813: . U - state at which to evaluate
4814: . Udot - time derivative of state vector
4815: . shift - shift to apply
4816: - ctx - context
4818: Output Parameters:
4819: + A - pointer to operator
4820: - B - pointer to matrix from which the preconditioner is built (often `A`)
4822: Level: advanced
4824: Notes:
4825: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4827: It is only appropriate for problems of the form
4829: $$
4830: M \dot{U} = F(U,t)
4831: $$
4833: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4834: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4835: an implicit operator of the form
4837: $$
4838: shift*M + J
4839: $$
4841: 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
4842: a copy of M or reassemble it when requested.
4844: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4845: @*/
4846: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscCtx ctx)
4847: {
4848: PetscFunctionBegin;
4849: PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4850: ts->ijacobian.shift = shift;
4851: PetscFunctionReturn(PETSC_SUCCESS);
4852: }
4854: /*@
4855: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4857: Not Collective
4859: Input Parameter:
4860: . ts - the `TS` context
4862: Output Parameter:
4863: . equation_type - see `TSEquationType`
4865: Level: beginner
4867: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4868: @*/
4869: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4870: {
4871: PetscFunctionBegin;
4873: PetscAssertPointer(equation_type, 2);
4874: *equation_type = ts->equation_type;
4875: PetscFunctionReturn(PETSC_SUCCESS);
4876: }
4878: /*@
4879: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4881: Not Collective
4883: Input Parameters:
4884: + ts - the `TS` context
4885: - equation_type - see `TSEquationType`
4887: Level: advanced
4889: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4890: @*/
4891: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4892: {
4893: PetscFunctionBegin;
4895: ts->equation_type = equation_type;
4896: PetscFunctionReturn(PETSC_SUCCESS);
4897: }
4899: /*@
4900: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4902: Not Collective
4904: Input Parameter:
4905: . ts - the `TS` context
4907: Output Parameter:
4908: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4909: manual pages for the individual convergence tests for complete lists
4911: Level: beginner
4913: Note:
4914: Can only be called after the call to `TSSolve()` is complete.
4916: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4917: @*/
4918: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4919: {
4920: PetscFunctionBegin;
4922: PetscAssertPointer(reason, 2);
4923: *reason = ts->reason;
4924: PetscFunctionReturn(PETSC_SUCCESS);
4925: }
4927: /*@
4928: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4930: Logically Collective; reason must contain common value
4932: Input Parameters:
4933: + ts - the `TS` context
4934: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4935: manual pages for the individual convergence tests for complete lists
4937: Level: advanced
4939: Note:
4940: Can only be called while `TSSolve()` is active.
4942: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4943: @*/
4944: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4945: {
4946: PetscFunctionBegin;
4948: ts->reason = reason;
4949: PetscFunctionReturn(PETSC_SUCCESS);
4950: }
4952: /*@
4953: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4955: Not Collective
4957: Input Parameter:
4958: . ts - the `TS` context
4960: Output Parameter:
4961: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4963: Level: beginner
4965: Note:
4966: Can only be called after the call to `TSSolve()` is complete.
4968: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4969: @*/
4970: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4971: {
4972: PetscFunctionBegin;
4974: PetscAssertPointer(ftime, 2);
4975: *ftime = ts->solvetime;
4976: PetscFunctionReturn(PETSC_SUCCESS);
4977: }
4979: /*@
4980: TSGetSNESIterations - Gets the total number of nonlinear iterations
4981: used by the time integrator.
4983: Not Collective
4985: Input Parameter:
4986: . ts - `TS` context
4988: Output Parameter:
4989: . nits - number of nonlinear iterations
4991: Level: intermediate
4993: Note:
4994: This counter is reset to zero for each successive call to `TSSolve()`.
4996: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4997: @*/
4998: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4999: {
5000: PetscFunctionBegin;
5002: PetscAssertPointer(nits, 2);
5003: *nits = ts->snes_its;
5004: PetscFunctionReturn(PETSC_SUCCESS);
5005: }
5007: /*@
5008: TSGetKSPIterations - Gets the total number of linear iterations
5009: used by the time integrator.
5011: Not Collective
5013: Input Parameter:
5014: . ts - `TS` context
5016: Output Parameter:
5017: . lits - number of linear iterations
5019: Level: intermediate
5021: Note:
5022: This counter is reset to zero for each successive call to `TSSolve()`.
5024: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`
5025: @*/
5026: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
5027: {
5028: PetscFunctionBegin;
5030: PetscAssertPointer(lits, 2);
5031: *lits = ts->ksp_its;
5032: PetscFunctionReturn(PETSC_SUCCESS);
5033: }
5035: /*@
5036: TSGetStepRejections - Gets the total number of rejected steps.
5038: Not Collective
5040: Input Parameter:
5041: . ts - `TS` context
5043: Output Parameter:
5044: . rejects - number of steps rejected
5046: Level: intermediate
5048: Note:
5049: This counter is reset to zero for each successive call to `TSSolve()`.
5051: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
5052: @*/
5053: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
5054: {
5055: PetscFunctionBegin;
5057: PetscAssertPointer(rejects, 2);
5058: *rejects = ts->reject;
5059: PetscFunctionReturn(PETSC_SUCCESS);
5060: }
5062: /*@
5063: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
5065: Not Collective
5067: Input Parameter:
5068: . ts - `TS` context
5070: Output Parameter:
5071: . fails - number of failed nonlinear solves
5073: Level: intermediate
5075: Note:
5076: This counter is reset to zero for each successive call to `TSSolve()`.
5078: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
5079: @*/
5080: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
5081: {
5082: PetscFunctionBegin;
5084: PetscAssertPointer(fails, 2);
5085: *fails = ts->num_snes_failures;
5086: PetscFunctionReturn(PETSC_SUCCESS);
5087: }
5089: /*@
5090: 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`
5092: Not Collective
5094: Input Parameters:
5095: + ts - `TS` context
5096: - rejects - maximum number of rejected steps, pass `PETSC_UNLIMITED` for unlimited
5098: Options Database Key:
5099: . -ts_max_step_rejections - Maximum number of step rejections before a step fails
5101: Level: intermediate
5103: Developer Note:
5104: The options database name is incorrect.
5106: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`,
5107: `TSGetConvergedReason()`, `TSSolve()`, `TS_DIVERGED_STEP_REJECTED`
5108: @*/
5109: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
5110: {
5111: PetscFunctionBegin;
5113: if (rejects == PETSC_UNLIMITED || rejects == -1) {
5114: ts->max_reject = PETSC_UNLIMITED;
5115: } else {
5116: PetscCheck(rejects >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of rejections");
5117: ts->max_reject = rejects;
5118: }
5119: PetscFunctionReturn(PETSC_SUCCESS);
5120: }
5122: /*@
5123: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves allowed before `TSSolve()` is ended with a `TSConvergedReason` of `TS_DIVERGED_NONLINEAR_SOLVE`
5125: Not Collective
5127: Input Parameters:
5128: + ts - `TS` context
5129: - fails - maximum number of failed nonlinear solves, pass `PETSC_UNLIMITED` to allow any number of failures.
5131: Options Database Key:
5132: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
5134: Level: intermediate
5136: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`,
5137: `TSGetConvergedReason()`, `TS_DIVERGED_NONLINEAR_SOLVE`, `TSConvergedReason`
5138: @*/
5139: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
5140: {
5141: PetscFunctionBegin;
5143: if (fails == PETSC_UNLIMITED || fails == -1) {
5144: ts->max_snes_failures = PETSC_UNLIMITED;
5145: } else {
5146: PetscCheck(fails >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures");
5147: ts->max_snes_failures = fails;
5148: }
5149: PetscFunctionReturn(PETSC_SUCCESS);
5150: }
5152: /*@
5153: TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`
5155: Not Collective
5157: Input Parameters:
5158: + ts - `TS` context
5159: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
5161: Options Database Key:
5162: . -ts_error_if_step_fails - Error if no step succeeds
5164: Level: intermediate
5166: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
5167: @*/
5168: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
5169: {
5170: PetscFunctionBegin;
5172: ts->errorifstepfailed = err;
5173: PetscFunctionReturn(PETSC_SUCCESS);
5174: }
5176: /*@
5177: TSGetAdapt - Get the adaptive controller context for the current method
5179: Collective if controller has not yet been created
5181: Input Parameter:
5182: . ts - time stepping context
5184: Output Parameter:
5185: . adapt - adaptive controller
5187: Level: intermediate
5189: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
5190: @*/
5191: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
5192: {
5193: PetscFunctionBegin;
5195: PetscAssertPointer(adapt, 2);
5196: if (!ts->adapt) {
5197: PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
5198: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
5199: }
5200: *adapt = ts->adapt;
5201: PetscFunctionReturn(PETSC_SUCCESS);
5202: }
5204: /*@
5205: TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller
5207: Logically Collective
5209: Input Parameters:
5210: + ts - time integration context
5211: . atol - scalar absolute tolerances
5212: . vatol - vector of absolute tolerances or `NULL`, used in preference to `atol` if present
5213: . rtol - scalar relative tolerances
5214: - vrtol - vector of relative tolerances or `NULL`, used in preference to `rtol` if present
5216: Options Database Keys:
5217: + -ts_rtol rtol - relative tolerance for local truncation error
5218: - -ts_atol atol - Absolute tolerance for local truncation error
5220: Level: beginner
5222: Notes:
5223: `PETSC_CURRENT` or `PETSC_DETERMINE` may be used for `atol` or `rtol` to indicate the current value
5224: or the default value from when the object's type was set.
5226: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5227: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5228: computed only for the differential or the algebraic part then this can be done using the vector of
5229: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
5230: differential part and infinity for the algebraic part, the LTE calculation will include only the
5231: differential variables.
5233: Fortran Note:
5234: Use `PETSC_CURRENT_INTEGER` or `PETSC_DETERMINE_INTEGER`.
5236: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
5237: @*/
5238: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
5239: {
5240: PetscFunctionBegin;
5241: if (atol == (PetscReal)PETSC_DETERMINE) {
5242: ts->atol = ts->default_atol;
5243: } else if (atol != (PetscReal)PETSC_CURRENT) {
5244: PetscCheck(atol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)atol);
5245: ts->atol = atol;
5246: }
5248: if (vatol) {
5249: PetscCall(PetscObjectReference((PetscObject)vatol));
5250: PetscCall(VecDestroy(&ts->vatol));
5251: ts->vatol = vatol;
5252: }
5254: if (rtol == (PetscReal)PETSC_DETERMINE) {
5255: ts->rtol = ts->default_rtol;
5256: } else if (rtol != (PetscReal)PETSC_CURRENT) {
5257: PetscCheck(rtol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative", (double)rtol);
5258: ts->rtol = rtol;
5259: }
5261: if (vrtol) {
5262: PetscCall(PetscObjectReference((PetscObject)vrtol));
5263: PetscCall(VecDestroy(&ts->vrtol));
5264: ts->vrtol = vrtol;
5265: }
5266: PetscFunctionReturn(PETSC_SUCCESS);
5267: }
5269: /*@
5270: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
5272: Logically Collective
5274: Input Parameter:
5275: . ts - time integration context
5277: Output Parameters:
5278: + atol - scalar absolute tolerances, `NULL` to ignore
5279: . vatol - vector of absolute tolerances, `NULL` to ignore
5280: . rtol - scalar relative tolerances, `NULL` to ignore
5281: - vrtol - vector of relative tolerances, `NULL` to ignore
5283: Level: beginner
5285: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
5286: @*/
5287: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
5288: {
5289: PetscFunctionBegin;
5290: if (atol) *atol = ts->atol;
5291: if (vatol) *vatol = ts->vatol;
5292: if (rtol) *rtol = ts->rtol;
5293: if (vrtol) *vrtol = ts->vrtol;
5294: PetscFunctionReturn(PETSC_SUCCESS);
5295: }
5297: /*@
5298: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
5300: Collective
5302: Input Parameters:
5303: + ts - time stepping context
5304: . U - state vector, usually ts->vec_sol
5305: . Y - state vector to be compared to U
5306: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5308: Output Parameters:
5309: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5310: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5311: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5313: Options Database Key:
5314: . -ts_adapt_wnormtype wnormtype - 2, INFINITY
5316: Level: developer
5318: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5319: @*/
5320: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5321: {
5322: PetscInt norma_loc, norm_loc, normr_loc;
5324: PetscFunctionBegin;
5329: PetscAssertPointer(norm, 5);
5330: PetscAssertPointer(norma, 6);
5331: PetscAssertPointer(normr, 7);
5332: 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));
5333: if (wnormtype == NORM_2) {
5334: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5335: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5336: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5337: }
5338: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5339: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5340: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5341: PetscFunctionReturn(PETSC_SUCCESS);
5342: }
5344: /*@
5345: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5347: Collective
5349: Input Parameters:
5350: + ts - time stepping context
5351: . E - error vector
5352: . U - state vector, usually ts->vec_sol
5353: . Y - state vector, previous time step
5354: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5356: Output Parameters:
5357: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5358: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5359: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5361: Options Database Key:
5362: . -ts_adapt_wnormtype wnormtype - 2, INFINITY
5364: Level: developer
5366: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5367: @*/
5368: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5369: {
5370: PetscInt norma_loc, norm_loc, normr_loc;
5372: PetscFunctionBegin;
5374: 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));
5375: if (wnormtype == NORM_2) {
5376: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5377: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5378: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5379: }
5380: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5381: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5382: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5383: PetscFunctionReturn(PETSC_SUCCESS);
5384: }
5386: /*@
5387: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5389: Logically Collective
5391: Input Parameters:
5392: + ts - time stepping context
5393: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5395: Note:
5396: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5398: Level: intermediate
5400: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5401: @*/
5402: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5403: {
5404: PetscFunctionBegin;
5406: ts->cfltime_local = cfltime;
5407: ts->cfltime = -1.;
5408: PetscFunctionReturn(PETSC_SUCCESS);
5409: }
5411: /*@
5412: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5414: Collective
5416: Input Parameter:
5417: . ts - time stepping context
5419: Output Parameter:
5420: . cfltime - maximum stable time step for forward Euler
5422: Level: advanced
5424: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5425: @*/
5426: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5427: {
5428: PetscFunctionBegin;
5429: if (ts->cfltime < 0) PetscCallMPI(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5430: *cfltime = ts->cfltime;
5431: PetscFunctionReturn(PETSC_SUCCESS);
5432: }
5434: /*@
5435: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5437: Input Parameters:
5438: + ts - the `TS` context.
5439: . xl - lower bound.
5440: - xu - upper bound.
5442: Level: advanced
5444: Note:
5445: If this routine is not called then the lower and upper bounds are set to
5446: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5448: .seealso: [](ch_ts), `TS`
5449: @*/
5450: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5451: {
5452: SNES snes;
5454: PetscFunctionBegin;
5455: PetscCall(TSGetSNES(ts, &snes));
5456: PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5457: PetscFunctionReturn(PETSC_SUCCESS);
5458: }
5460: /*@
5461: TSComputeLinearStability - computes the linear stability function at a point
5463: Collective
5465: Input Parameters:
5466: + ts - the `TS` context
5467: . xr - real part of input argument
5468: - xi - imaginary part of input argument
5470: Output Parameters:
5471: + yr - real part of function value
5472: - yi - imaginary part of function value
5474: Level: developer
5476: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5477: @*/
5478: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5479: {
5480: PetscFunctionBegin;
5482: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5483: PetscFunctionReturn(PETSC_SUCCESS);
5484: }
5486: /*@
5487: TSRestartStep - Flags the solver to restart the next step
5489: Collective
5491: Input Parameter:
5492: . ts - the `TS` context obtained from `TSCreate()`
5494: Level: advanced
5496: Notes:
5497: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5498: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5499: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5500: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5501: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5502: discontinuous source terms).
5504: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5505: @*/
5506: PetscErrorCode TSRestartStep(TS ts)
5507: {
5508: PetscFunctionBegin;
5510: ts->steprestart = PETSC_TRUE;
5511: PetscFunctionReturn(PETSC_SUCCESS);
5512: }
5514: /*@
5515: TSRollBack - Rolls back one time step
5517: Collective
5519: Input Parameter:
5520: . ts - the `TS` context obtained from `TSCreate()`
5522: Level: advanced
5524: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5525: @*/
5526: PetscErrorCode TSRollBack(TS ts)
5527: {
5528: PetscFunctionBegin;
5530: PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5531: PetscTryTypeMethod(ts, rollback);
5532: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
5533: ts->time_step = ts->ptime - ts->ptime_prev;
5534: ts->ptime = ts->ptime_prev;
5535: ts->ptime_prev = ts->ptime_prev_rollback;
5536: ts->steps--;
5537: ts->steprollback = PETSC_TRUE;
5538: PetscFunctionReturn(PETSC_SUCCESS);
5539: }
5541: /*@
5542: TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step
5544: Not collective
5546: Input Parameter:
5547: . ts - the `TS` context obtained from `TSCreate()`
5549: Output Parameter:
5550: . flg - the rollback flag
5552: Level: advanced
5554: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
5555: @*/
5556: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
5557: {
5558: PetscFunctionBegin;
5560: PetscAssertPointer(flg, 2);
5561: *flg = ts->steprollback;
5562: PetscFunctionReturn(PETSC_SUCCESS);
5563: }
5565: /*@
5566: TSGetStepResize - Get the internal flag indicating if the current step is after a resize.
5568: Not collective
5570: Input Parameter:
5571: . ts - the `TS` context obtained from `TSCreate()`
5573: Output Parameter:
5574: . flg - the resize flag
5576: Level: advanced
5578: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetResize()`
5579: @*/
5580: PetscErrorCode TSGetStepResize(TS ts, PetscBool *flg)
5581: {
5582: PetscFunctionBegin;
5584: PetscAssertPointer(flg, 2);
5585: *flg = ts->stepresize;
5586: PetscFunctionReturn(PETSC_SUCCESS);
5587: }
5589: /*@
5590: TSGetStages - Get the number of stages and stage values
5592: Input Parameter:
5593: . ts - the `TS` context obtained from `TSCreate()`
5595: Output Parameters:
5596: + ns - the number of stages
5597: - Y - the current stage vectors
5599: Level: advanced
5601: Note:
5602: Both `ns` and `Y` can be `NULL`.
5604: .seealso: [](ch_ts), `TS`, `TSCreate()`
5605: @*/
5606: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5607: {
5608: PetscFunctionBegin;
5610: if (ns) PetscAssertPointer(ns, 2);
5611: if (Y) PetscAssertPointer(Y, 3);
5612: if (!ts->ops->getstages) {
5613: if (ns) *ns = 0;
5614: if (Y) *Y = NULL;
5615: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5616: PetscFunctionReturn(PETSC_SUCCESS);
5617: }
5619: /*@
5620: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5622: Collective
5624: Input Parameters:
5625: + ts - the `TS` context
5626: . t - current timestep
5627: . U - state vector
5628: . Udot - time derivative of state vector
5629: . shift - shift to apply, see note below
5630: - ctx - an optional application context
5632: Output Parameters:
5633: + J - Jacobian matrix (not altered in this routine)
5634: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
5636: Level: intermediate
5638: Notes:
5639: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5641: dF/dU + shift*dF/dUdot
5643: Most users should not need to explicitly call this routine, as it
5644: is used internally within the nonlinear solvers.
5646: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5647: routine, then it will try to get the coloring from the matrix. This requires that the
5648: matrix have nonzero entries precomputed.
5650: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5651: @*/
5652: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, PetscCtx ctx)
5653: {
5654: SNES snes;
5655: MatFDColoring color;
5656: PetscBool hascolor, matcolor = PETSC_FALSE;
5658: PetscFunctionBegin;
5659: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5660: PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5661: if (!color) {
5662: DM dm;
5663: ISColoring iscoloring;
5665: PetscCall(TSGetDM(ts, &dm));
5666: PetscCall(DMHasColoring(dm, &hascolor));
5667: if (hascolor && !matcolor) {
5668: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5669: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5670: PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5671: PetscCall(MatFDColoringSetFromOptions(color));
5672: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5673: PetscCall(ISColoringDestroy(&iscoloring));
5674: } else {
5675: MatColoring mc;
5677: PetscCall(MatColoringCreate(B, &mc));
5678: PetscCall(MatColoringSetDistance(mc, 2));
5679: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5680: PetscCall(MatColoringSetFromOptions(mc));
5681: PetscCall(MatColoringApply(mc, &iscoloring));
5682: PetscCall(MatColoringDestroy(&mc));
5683: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5684: PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5685: PetscCall(MatFDColoringSetFromOptions(color));
5686: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5687: PetscCall(ISColoringDestroy(&iscoloring));
5688: }
5689: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5690: PetscCall(PetscObjectDereference((PetscObject)color));
5691: }
5692: PetscCall(TSGetSNES(ts, &snes));
5693: PetscCall(MatFDColoringApply(B, color, U, snes));
5694: if (J != B) {
5695: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5696: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5697: }
5698: PetscFunctionReturn(PETSC_SUCCESS);
5699: }
5701: /*@
5702: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5704: Logically collective
5706: Input Parameters:
5707: + ts - the `TS` context
5708: - func - function called within `TSFunctionDomainError()`
5710: Calling sequence of `func`:
5711: + ts - the `TS` context
5712: . time - the current time (of the stage)
5713: . state - the state to check if it is valid
5714: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable
5716: Level: intermediate
5718: Notes:
5719: `accept` must be collectively specified.
5720: If an implicit ODE solver is being used then, in addition to providing this routine, the
5721: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5722: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5723: Use `TSGetSNES()` to obtain the `SNES` object
5725: Developer Notes:
5726: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5727: since one takes a function pointer and the other does not.
5729: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5730: @*/
5731: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5732: {
5733: PetscFunctionBegin;
5735: ts->functiondomainerror = func;
5736: PetscFunctionReturn(PETSC_SUCCESS);
5737: }
5739: /*@
5740: TSFunctionDomainError - Checks if the current state is valid
5742: Collective
5744: Input Parameters:
5745: + ts - the `TS` context
5746: . stagetime - time of the simulation
5747: - Y - state vector to check.
5749: Output Parameter:
5750: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5752: Level: developer
5754: Note:
5755: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5756: to check if the current state is valid.
5758: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5759: @*/
5760: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5761: {
5762: PetscFunctionBegin;
5766: PetscAssertPointer(accept, 4);
5767: *accept = PETSC_TRUE;
5768: if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5769: PetscFunctionReturn(PETSC_SUCCESS);
5770: }
5772: /*@
5773: TSClone - This function clones a time step `TS` object.
5775: Collective
5777: Input Parameter:
5778: . tsin - The input `TS`
5780: Output Parameter:
5781: . tsout - The output `TS` (cloned)
5783: Level: developer
5785: Notes:
5786: 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.
5787: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5789: 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
5790: .vb
5791: SNES snes_dup = NULL;
5792: TSGetSNES(ts,&snes_dup);
5793: TSSetSNES(ts,snes_dup);
5794: .ve
5796: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5797: @*/
5798: PetscErrorCode TSClone(TS tsin, TS *tsout)
5799: {
5800: TS t;
5801: SNES snes_start;
5802: DM dm;
5803: TSType type;
5805: PetscFunctionBegin;
5806: PetscAssertPointer(tsin, 1);
5807: *tsout = NULL;
5809: PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));
5811: /* General TS description */
5812: t->numbermonitors = 0;
5813: t->setupcalled = PETSC_FALSE;
5814: t->ksp_its = 0;
5815: t->snes_its = 0;
5816: t->nwork = 0;
5817: t->rhsjacobian.time = PETSC_MIN_REAL;
5818: t->rhsjacobian.scale = 1.;
5819: t->ijacobian.shift = 1.;
5821: PetscCall(TSGetSNES(tsin, &snes_start));
5822: PetscCall(TSSetSNES(t, snes_start));
5824: PetscCall(TSGetDM(tsin, &dm));
5825: PetscCall(TSSetDM(t, dm));
5827: t->adapt = tsin->adapt;
5828: PetscCall(PetscObjectReference((PetscObject)t->adapt));
5830: t->trajectory = tsin->trajectory;
5831: PetscCall(PetscObjectReference((PetscObject)t->trajectory));
5833: t->event = tsin->event;
5834: if (t->event) t->event->refct++;
5836: t->problem_type = tsin->problem_type;
5837: t->ptime = tsin->ptime;
5838: t->ptime_prev = tsin->ptime_prev;
5839: t->time_step = tsin->time_step;
5840: t->max_time = tsin->max_time;
5841: t->steps = tsin->steps;
5842: t->max_steps = tsin->max_steps;
5843: t->equation_type = tsin->equation_type;
5844: t->atol = tsin->atol;
5845: t->rtol = tsin->rtol;
5846: t->max_snes_failures = tsin->max_snes_failures;
5847: t->max_reject = tsin->max_reject;
5848: t->errorifstepfailed = tsin->errorifstepfailed;
5850: PetscCall(TSGetType(tsin, &type));
5851: PetscCall(TSSetType(t, type));
5853: t->vec_sol = NULL;
5855: t->cfltime = tsin->cfltime;
5856: t->cfltime_local = tsin->cfltime_local;
5857: t->exact_final_time = tsin->exact_final_time;
5859: t->ops[0] = tsin->ops[0];
5861: if (((PetscObject)tsin)->fortran_func_pointers) {
5862: PetscCall(PetscMalloc((10) * sizeof(PetscFortranCallbackFn *), &((PetscObject)t)->fortran_func_pointers));
5863: for (PetscInt i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5864: }
5865: *tsout = t;
5866: PetscFunctionReturn(PETSC_SUCCESS);
5867: }
5869: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(PetscCtx ctx, Vec x, Vec y)
5870: {
5871: TS ts = (TS)ctx;
5873: PetscFunctionBegin;
5874: PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5875: PetscFunctionReturn(PETSC_SUCCESS);
5876: }
5878: /*@
5879: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5881: Logically Collective
5883: Input Parameter:
5884: . ts - the time stepping routine
5886: Output Parameter:
5887: . flg - `PETSC_TRUE` if the multiply is likely correct
5889: Options Database Key:
5890: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5892: Level: advanced
5894: Note:
5895: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5897: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5898: @*/
5899: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5900: {
5901: Mat J, B;
5902: TSRHSJacobianFn *func;
5903: void *ctx;
5905: PetscFunctionBegin;
5906: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5907: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5908: PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5909: PetscFunctionReturn(PETSC_SUCCESS);
5910: }
5912: /*@
5913: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5915: Logically Collective
5917: Input Parameter:
5918: . ts - the time stepping routine
5920: Output Parameter:
5921: . flg - `PETSC_TRUE` if the multiply is likely correct
5923: Options Database Key:
5924: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5926: Level: advanced
5928: Notes:
5929: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5931: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5932: @*/
5933: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5934: {
5935: Mat J, B;
5936: void *ctx;
5937: TSRHSJacobianFn *func;
5939: PetscFunctionBegin;
5940: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5941: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5942: PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5943: PetscFunctionReturn(PETSC_SUCCESS);
5944: }
5946: /*@
5947: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
5949: Logically Collective
5951: Input Parameters:
5952: + ts - timestepping context
5953: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5955: Options Database Key:
5956: . -ts_use_splitrhsfunction (true|false) - use the split RHS function for multirate solvers
5958: Level: intermediate
5960: Note:
5961: This is only for multirate methods
5963: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5964: @*/
5965: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5966: {
5967: PetscFunctionBegin;
5969: ts->use_splitrhsfunction = use_splitrhsfunction;
5970: PetscFunctionReturn(PETSC_SUCCESS);
5971: }
5973: /*@
5974: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
5976: Not Collective
5978: Input Parameter:
5979: . ts - timestepping context
5981: Output Parameter:
5982: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5984: Level: intermediate
5986: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5987: @*/
5988: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5989: {
5990: PetscFunctionBegin;
5992: *use_splitrhsfunction = ts->use_splitrhsfunction;
5993: PetscFunctionReturn(PETSC_SUCCESS);
5994: }
5996: /*@
5997: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
5999: Logically Collective
6001: Input Parameters:
6002: + ts - the time-stepper
6003: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
6005: Level: intermediate
6007: Note:
6008: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
6010: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
6011: @*/
6012: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
6013: {
6014: PetscFunctionBegin;
6016: ts->axpy_pattern = str;
6017: PetscFunctionReturn(PETSC_SUCCESS);
6018: }
6020: /*@
6021: TSSetEvaluationTimes - sets the evaluation points. The solution will be computed and stored for each time requested
6023: Collective
6025: Input Parameters:
6026: + ts - the time-stepper
6027: . n - number of the time points
6028: - time_points - array of the time points, must be increasing
6030: Options Database Key:
6031: . -ts_eval_times t0,...,tn - Sets the evaluation times
6033: Level: intermediate
6035: Notes:
6036: The elements in `time_points` must be all increasing. They correspond to the intermediate points to be saved.
6038: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
6040: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using evaluation times may
6041: pressure the memory system when using a large number of time points.
6043: .seealso: [](ch_ts), `TS`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`, `TSSetTimeSpan()`
6044: @*/
6045: PetscErrorCode TSSetEvaluationTimes(TS ts, PetscInt n, PetscReal time_points[])
6046: {
6047: PetscBool is_sorted;
6049: PetscFunctionBegin;
6051: if (ts->eval_times) { // Reset eval_times
6052: ts->eval_times->sol_idx = 0;
6053: ts->eval_times->time_point_idx = 0;
6054: if (n != ts->eval_times->num_time_points) {
6055: PetscCall(PetscFree(ts->eval_times->time_points));
6056: PetscCall(PetscFree(ts->eval_times->sol_times));
6057: PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
6058: } else {
6059: PetscCall(PetscArrayzero(ts->eval_times->sol_times, n));
6060: for (PetscInt i = 0; i < n; i++) PetscCall(VecZeroEntries(ts->eval_times->sol_vecs[i]));
6061: }
6062: } else { // Create/initialize eval_times
6063: TSEvaluationTimes eval_times;
6064: PetscCall(PetscNew(&eval_times));
6065: PetscCall(PetscMalloc1(n, &eval_times->time_points));
6066: PetscCall(PetscMalloc1(n, &eval_times->sol_times));
6067: eval_times->reltol = 1e-6;
6068: eval_times->abstol = 10 * PETSC_MACHINE_EPSILON;
6069: eval_times->worktol = 0;
6070: ts->eval_times = eval_times;
6071: }
6072: ts->eval_times->num_time_points = n;
6073: PetscCall(PetscSortedReal(n, time_points, &is_sorted));
6074: PetscCheck(is_sorted, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "time_points array must be sorted");
6075: PetscCall(PetscArraycpy(ts->eval_times->time_points, time_points, n));
6076: // Note: ts->vec_sol not guaranteed to exist, so ts->eval_times->sol_vecs allocated at TSSolve time
6077: PetscFunctionReturn(PETSC_SUCCESS);
6078: }
6080: /*@
6081: TSGetEvaluationTimes - gets the evaluation times set with `TSSetEvaluationTimes()`
6083: Not Collective
6085: Input Parameter:
6086: . ts - the time-stepper
6088: Output Parameters:
6089: + n - number of the time points
6090: - time_points - array of the time points
6092: Level: beginner
6094: Note:
6095: The values obtained are valid until the `TS` object is destroyed.
6097: Both `n` and `time_points` can be `NULL`.
6099: Also used to see time points set by `TSSetTimeSpan()`.
6101: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6102: @*/
6103: PetscErrorCode TSGetEvaluationTimes(TS ts, PetscInt *n, const PetscReal *time_points[])
6104: {
6105: PetscFunctionBegin;
6107: if (n) PetscAssertPointer(n, 2);
6108: if (time_points) PetscAssertPointer(time_points, 3);
6109: if (!ts->eval_times) {
6110: if (n) *n = 0;
6111: if (time_points) *time_points = NULL;
6112: } else {
6113: if (n) *n = ts->eval_times->num_time_points;
6114: if (time_points) *time_points = ts->eval_times->time_points;
6115: }
6116: PetscFunctionReturn(PETSC_SUCCESS);
6117: }
6119: /*@
6120: TSGetEvaluationSolutions - Get the number of solutions and the solutions at the evaluation time points specified
6122: Input Parameter:
6123: . ts - the `TS` context obtained from `TSCreate()`
6125: Output Parameters:
6126: + nsol - the number of solutions
6127: . sol_times - array of solution times corresponding to the solution vectors. See note below
6128: - Sols - the solution vectors
6130: Level: intermediate
6132: Notes:
6133: Both `nsol` and `Sols` can be `NULL`.
6135: 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()`.
6136: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain evaluation times.
6138: Also used to see view solutions requested by `TSSetTimeSpan()`.
6140: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`
6141: @*/
6142: PetscErrorCode TSGetEvaluationSolutions(TS ts, PetscInt *nsol, const PetscReal *sol_times[], Vec *Sols[])
6143: {
6144: PetscFunctionBegin;
6146: if (nsol) PetscAssertPointer(nsol, 2);
6147: if (sol_times) PetscAssertPointer(sol_times, 3);
6148: if (Sols) PetscAssertPointer(Sols, 4);
6149: if (!ts->eval_times) {
6150: if (nsol) *nsol = 0;
6151: if (sol_times) *sol_times = NULL;
6152: if (Sols) *Sols = NULL;
6153: } else {
6154: if (nsol) *nsol = ts->eval_times->sol_idx;
6155: if (sol_times) *sol_times = ts->eval_times->sol_times;
6156: if (Sols) *Sols = ts->eval_times->sol_vecs;
6157: }
6158: PetscFunctionReturn(PETSC_SUCCESS);
6159: }
6161: /*@
6162: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
6164: Collective
6166: Input Parameters:
6167: + ts - the time-stepper
6168: . n - number of the time points (>=2)
6169: - 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.
6171: Options Database Key:
6172: . -ts_time_span t0,...,tf - Sets the time span
6174: Level: intermediate
6176: Notes:
6177: 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.
6179: The elements in `span_times` must be all increasing. They correspond to the intermediate points to be saved.
6181: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
6183: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using time span may
6184: pressure the memory system when using a large number of span points.
6186: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6187: @*/
6188: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal span_times[])
6189: {
6190: PetscFunctionBegin;
6192: PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
6193: PetscCall(TSSetEvaluationTimes(ts, n, span_times));
6194: PetscCall(TSSetTime(ts, span_times[0]));
6195: PetscCall(TSSetMaxTime(ts, span_times[n - 1]));
6196: PetscFunctionReturn(PETSC_SUCCESS);
6197: }
6199: /*@
6200: TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
6202: Collective
6204: Input Parameters:
6205: + ts - the `TS` context
6206: . J - Jacobian matrix (not altered in this routine)
6207: - B - newly computed Jacobian matrix to use with preconditioner
6209: Level: intermediate
6211: Notes:
6212: This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
6213: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
6214: and multiple fields are involved.
6216: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
6217: structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
6218: usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
6219: `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.
6221: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
6222: @*/
6223: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
6224: {
6225: MatColoring mc = NULL;
6226: ISColoring iscoloring = NULL;
6227: MatFDColoring matfdcoloring = NULL;
6229: PetscFunctionBegin;
6230: /* Generate new coloring after eliminating zeros in the matrix */
6231: PetscCall(MatEliminateZeros(B, PETSC_TRUE));
6232: PetscCall(MatColoringCreate(B, &mc));
6233: PetscCall(MatColoringSetDistance(mc, 2));
6234: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
6235: PetscCall(MatColoringSetFromOptions(mc));
6236: PetscCall(MatColoringApply(mc, &iscoloring));
6237: PetscCall(MatColoringDestroy(&mc));
6238: /* Replace the old coloring with the new one */
6239: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
6240: PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
6241: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
6242: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
6243: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
6244: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
6245: PetscCall(ISColoringDestroy(&iscoloring));
6246: PetscFunctionReturn(PETSC_SUCCESS);
6247: }