Actual source code: ts.c
1: #include <petsc/private/tsimpl.h>
2: #include <petscdmda.h>
3: #include <petscdmshell.h>
4: #include <petscdmplex.h>
5: #include <petscdmswarm.h>
6: #include <petscviewer.h>
7: #include <petscdraw.h>
8: #include <petscconvest.h>
10: /* Logging support */
11: PetscClassId TS_CLASSID, DMTS_CLASSID;
12: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;
14: const char *const TSExactFinalTimeOptions[] = {"UNSPECIFIED", "STEPOVER", "INTERPOLATE", "MATCHSTEP", "TSExactFinalTimeOption", "TS_EXACTFINALTIME_", NULL};
16: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt, TSAdaptType default_type)
17: {
18: PetscFunctionBegin;
20: PetscAssertPointer(default_type, 2);
21: if (!((PetscObject)adapt)->type_name) PetscCall(TSAdaptSetType(adapt, default_type));
22: PetscFunctionReturn(PETSC_SUCCESS);
23: }
25: /*@
26: TSSetFromOptions - Sets various `TS` parameters from the options database
28: Collective
30: Input Parameter:
31: . ts - the `TS` context obtained from `TSCreate()`
33: Options Database Keys:
34: + -ts_type <type> - EULER, BEULER, SUNDIALS, PSEUDO, CN, RK, THETA, ALPHA, GLLE, SSP, GLEE, BSYMP, IRK, see `TSType`
35: . -ts_save_trajectory - checkpoint the solution at each time-step
36: . -ts_max_time <time> - maximum time to compute to
37: . -ts_time_span <t0,...tf> - sets the time span, solutions are computed and stored for each indicated time, init_time and max_time are set
38: . -ts_eval_times <t0,...tn> - time points where solutions are computed and stored for each indicated time
39: . -ts_max_steps <steps> - maximum number of time-steps to take
40: . -ts_init_time <time> - initial time to start computation
41: . -ts_final_time <time> - final time to compute to (deprecated: use `-ts_max_time`)
42: . -ts_dt <dt> - initial time step
43: . -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
44: . -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
45: . -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
46: . -ts_error_if_step_fails <true,false> - Error if no step succeeds
47: . -ts_rtol <rtol> - relative tolerance for local truncation error
48: . -ts_atol <atol> - Absolute tolerance for local truncation error
49: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function
50: . -ts_rhs_jacobian_test_mult_transpose - test the Jacobian at each iteration against finite difference with RHS function
51: . -ts_adjoint_solve <yes,no> - After solving the ODE/DAE solve the adjoint problem (requires `-ts_save_trajectory`)
52: . -ts_fd_color - Use finite differences with coloring to compute IJacobian
53: . -ts_monitor - print information at each timestep
54: . -ts_monitor_cancel - Cancel all monitors
55: . -ts_monitor_lg_solution - Monitor solution graphically
56: . -ts_monitor_lg_error - Monitor error graphically
57: . -ts_monitor_error - Monitors norm of error
58: . -ts_monitor_lg_timestep - Monitor timestep size graphically
59: . -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
60: . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
61: . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
62: . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
63: . -ts_monitor_draw_solution - Monitor solution graphically
64: . -ts_monitor_draw_solution_phase <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
65: . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
66: . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
67: . -ts_monitor_solution_interval <interval> - output once every interval (default=1) time steps. Use -1 to only output at the end of the simulation
68: . -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)
69: . -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
70: - -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
72: Level: beginner
74: Notes:
75: See `SNESSetFromOptions()` and `KSPSetFromOptions()` for how to control the nonlinear and linear solves used by the time-stepper.
77: Certain `SNES` options get reset for each new nonlinear solver, for example `-snes_lag_jacobian its` and `-snes_lag_preconditioner its`, in order
78: to retain them over the multiple nonlinear solves that `TS` uses you must also provide `-snes_lag_jacobian_persists true` and
79: `-snes_lag_preconditioner_persists true`
81: Developer Notes:
82: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
84: .seealso: [](ch_ts), `TS`, `TSGetType()`
85: @*/
86: PetscErrorCode TSSetFromOptions(TS ts)
87: {
88: PetscBool opt, flg, tflg;
89: char monfilename[PETSC_MAX_PATH_LEN];
90: PetscReal time_step, eval_times[100];
91: PetscInt num_eval_times = PETSC_STATIC_ARRAY_LENGTH(eval_times);
92: TSExactFinalTimeOption eftopt;
93: char dir[16];
94: TSIFunctionFn *ifun;
95: const char *defaultType;
96: char typeName[256];
98: PetscFunctionBegin;
101: PetscCall(TSRegisterAll());
102: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
104: PetscObjectOptionsBegin((PetscObject)ts);
105: if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
106: else defaultType = ifun ? TSBEULER : TSEULER;
107: PetscCall(PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, 256, &opt));
108: if (opt) PetscCall(TSSetType(ts, typeName));
109: else PetscCall(TSSetType(ts, defaultType));
111: /* Handle generic TS options */
112: PetscCall(PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL));
113: PetscCall(PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL));
114: PetscCall(PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", eval_times, &num_eval_times, &flg));
115: if (flg) PetscCall(TSSetTimeSpan(ts, num_eval_times, eval_times));
116: num_eval_times = PETSC_STATIC_ARRAY_LENGTH(eval_times);
117: PetscCall(PetscOptionsRealArray("-ts_eval_times", "Evaluation time points", "TSSetEvaluationTimes", eval_times, &num_eval_times, &opt));
118: PetscCheck(flg != opt || (!flg && !opt), PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "May not provide -ts_time_span and -ts_eval_times simultaneously");
119: if (opt) PetscCall(TSSetEvaluationTimes(ts, num_eval_times, eval_times));
120: PetscCall(PetscOptionsInt("-ts_max_steps", "Maximum number of time steps", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL));
121: PetscCall(PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL));
122: PetscCall(PetscOptionsReal("-ts_dt", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg));
123: if (flg) PetscCall(TSSetTimeStep(ts, time_step));
124: PetscCall(PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg));
125: if (flg) PetscCall(TSSetExactFinalTime(ts, eftopt));
126: PetscCall(PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, &flg));
127: if (flg) PetscCall(TSSetMaxSNESFailures(ts, ts->max_snes_failures));
128: PetscCall(PetscOptionsInt("-ts_max_reject", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, &flg));
129: if (flg) PetscCall(TSSetMaxStepRejections(ts, ts->max_reject));
130: PetscCall(PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL));
131: PetscCall(PetscOptionsBoundedReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL, 0));
132: PetscCall(PetscOptionsBoundedReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL, 0));
134: PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL));
135: 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));
136: PetscCall(PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL));
137: #if defined(PETSC_HAVE_SAWS)
138: {
139: PetscBool set;
140: flg = PETSC_FALSE;
141: PetscCall(PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set));
142: if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)ts, flg));
143: }
144: #endif
146: /* Monitor options */
147: PetscCall(PetscOptionsInt("-ts_monitor_frequency", "Number of time steps between monitor output", "TSMonitorSetFrequency", ts->monitorFrequency, &ts->monitorFrequency, NULL));
148: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL));
149: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL));
150: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, NULL));
151: PetscCall(TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL));
153: PetscCall(PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg));
154: if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)ts, monfilename));
156: PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt));
157: if (opt) {
158: PetscInt howoften = 1;
159: DM dm;
160: PetscBool net;
162: PetscCall(PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL));
163: PetscCall(TSGetDM(ts, &dm));
164: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net));
165: if (net) {
166: TSMonitorLGCtxNetwork ctx;
167: PetscCall(TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx));
168: PetscCall(TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxNetworkDestroy));
169: PetscCall(PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL));
170: } else {
171: TSMonitorLGCtx ctx;
172: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
173: PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
174: }
175: }
177: PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt));
178: if (opt) {
179: TSMonitorLGCtx ctx;
180: PetscInt howoften = 1;
182: PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL));
183: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
184: PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
185: }
186: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL));
188: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt));
189: if (opt) {
190: TSMonitorLGCtx ctx;
191: PetscInt howoften = 1;
193: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
194: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
195: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
196: }
197: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt));
198: if (opt) {
199: TSMonitorLGCtx ctx;
200: PetscInt howoften = 1;
202: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
203: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
204: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
205: ctx->semilogy = PETSC_TRUE;
206: }
208: PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt));
209: if (opt) {
210: TSMonitorLGCtx ctx;
211: PetscInt howoften = 1;
213: PetscCall(PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL));
214: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
215: PetscCall(TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
216: }
217: PetscCall(PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt));
218: if (opt) {
219: TSMonitorLGCtx ctx;
220: PetscInt howoften = 1;
222: PetscCall(PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL));
223: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
224: PetscCall(TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
225: }
226: PetscCall(PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt));
227: if (opt) {
228: TSMonitorSPEigCtx ctx;
229: PetscInt howoften = 1;
231: PetscCall(PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL));
232: PetscCall(TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
233: PetscCall(TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscCtxDestroyFn *)TSMonitorSPEigCtxDestroy));
234: }
235: PetscCall(PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase space from the DMSwarm", "TSMonitorSPSwarm", &opt));
236: if (opt) {
237: TSMonitorSPCtx ctx;
238: PetscInt howoften = 1, retain = 0;
239: PetscBool phase = PETSC_TRUE, create = PETSC_TRUE, multispecies = PETSC_FALSE;
241: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
242: if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
243: create = PETSC_FALSE;
244: break;
245: }
246: if (create) {
247: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase space from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL));
248: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL));
249: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL));
250: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_multi_species", "Color particles by particle species", "TSMonitorSPSwarm", multispecies, &multispecies, NULL));
251: PetscCall(TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, multispecies, &ctx));
252: PetscCall(TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscCtxDestroyFn *)TSMonitorSPCtxDestroy));
253: }
254: }
255: PetscCall(PetscOptionsName("-ts_monitor_hg_swarm", "Display particle histogram from the DMSwarm", "TSMonitorHGSwarm", &opt));
256: if (opt) {
257: TSMonitorHGCtx ctx;
258: PetscInt howoften = 1, Ns = 1;
259: PetscBool velocity = PETSC_FALSE, create = PETSC_TRUE;
261: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
262: if (ts->monitor[i] == TSMonitorHGSwarmSolution) {
263: create = PETSC_FALSE;
264: break;
265: }
266: if (create) {
267: DM sw, dm;
268: PetscInt Nc, Nb;
270: PetscCall(TSGetDM(ts, &sw));
271: PetscCall(DMSwarmGetCellDM(sw, &dm));
272: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc));
273: Nb = PetscMin(20, PetscMax(10, Nc));
274: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL));
275: PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL));
276: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL));
277: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL));
278: PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx));
279: PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscCtxDestroyFn *)TSMonitorHGCtxDestroy));
280: }
281: }
282: opt = PETSC_FALSE;
283: PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt));
284: if (opt) {
285: TSMonitorDrawCtx ctx;
286: PetscInt howoften = 1;
288: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL));
289: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
290: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
291: }
292: opt = PETSC_FALSE;
293: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt));
294: if (opt) {
295: TSMonitorDrawCtx ctx;
296: PetscReal bounds[4];
297: PetscInt n = 4;
298: PetscDraw draw;
299: PetscDrawAxis axis;
301: PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL));
302: PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field");
303: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx));
304: PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw));
305: PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis));
306: PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]));
307: PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2"));
308: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
309: }
310: opt = PETSC_FALSE;
311: PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt));
312: if (opt) {
313: TSMonitorDrawCtx ctx;
314: PetscInt howoften = 1;
316: PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL));
317: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
318: PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
319: }
320: opt = PETSC_FALSE;
321: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt));
322: if (opt) {
323: TSMonitorDrawCtx ctx;
324: PetscInt howoften = 1;
326: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL));
327: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
328: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
329: }
331: opt = PETSC_FALSE;
332: 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));
333: if (flg) {
334: TSMonitorVTKCtx ctx;
336: PetscCall(TSMonitorSolutionVTKCtxCreate(monfilename, &ctx));
337: PetscCall(PetscOptionsInt("-ts_monitor_solution_vtk_interval", "Save every interval time step (-1 for last step only)", NULL, ctx->interval, &ctx->interval, NULL));
338: PetscCall(TSMonitorSet(ts, (PetscErrorCode (*)(TS, PetscInt, PetscReal, Vec, void *))TSMonitorSolutionVTK, ctx, (PetscCtxDestroyFn *)TSMonitorSolutionVTKDestroy));
339: }
341: PetscCall(PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg));
342: if (flg) {
343: TSMonitorDMDARayCtx *rayctx;
344: int ray = 0;
345: DMDirection ddir;
346: DM da;
347: PetscMPIInt rank;
349: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
350: if (dir[0] == 'x') ddir = DM_X;
351: else if (dir[0] == 'y') ddir = DM_Y;
352: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
353: sscanf(dir + 2, "%d", &ray);
355: PetscCall(PetscInfo(ts, "Displaying DMDA ray %c = %d\n", dir[0], ray));
356: PetscCall(PetscNew(&rayctx));
357: PetscCall(TSGetDM(ts, &da));
358: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
359: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank));
360: if (rank == 0) PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer));
361: rayctx->lgctx = NULL;
362: PetscCall(TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy));
363: }
364: PetscCall(PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg));
365: if (flg) {
366: TSMonitorDMDARayCtx *rayctx;
367: int ray = 0;
368: DMDirection ddir;
369: DM da;
370: PetscInt howoften = 1;
372: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
373: if (dir[0] == 'x') ddir = DM_X;
374: else if (dir[0] == 'y') ddir = DM_Y;
375: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
376: sscanf(dir + 2, "%d", &ray);
378: PetscCall(PetscInfo(ts, "Displaying LG DMDA ray %c = %d\n", dir[0], ray));
379: PetscCall(PetscNew(&rayctx));
380: PetscCall(TSGetDM(ts, &da));
381: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
382: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx));
383: PetscCall(TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy));
384: }
386: PetscCall(PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt));
387: if (opt) {
388: TSMonitorEnvelopeCtx ctx;
390: PetscCall(TSMonitorEnvelopeCtxCreate(ts, &ctx));
391: PetscCall(TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscCtxDestroyFn *)TSMonitorEnvelopeCtxDestroy));
392: }
393: flg = PETSC_FALSE;
394: PetscCall(PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt));
395: if (opt && flg) PetscCall(TSMonitorCancel(ts));
397: flg = PETSC_FALSE;
398: PetscCall(PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeIJacobianDefaultColor", flg, &flg, NULL));
399: if (flg) {
400: DM dm;
402: PetscCall(TSGetDM(ts, &dm));
403: PetscCall(DMTSUnsetIJacobianContext_Internal(dm));
404: PetscCall(TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL));
405: PetscCall(PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n"));
406: }
408: /* Handle specific TS options */
409: PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);
411: /* Handle TSAdapt options */
412: PetscCall(TSGetAdapt(ts, &ts->adapt));
413: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
414: PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject));
416: /* TS trajectory must be set after TS, since it may use some TS options above */
417: tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
418: PetscCall(PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL));
419: if (tflg) PetscCall(TSSetSaveTrajectory(ts));
421: PetscCall(TSAdjointSetFromOptions(ts, PetscOptionsObject));
423: /* process any options handlers added with PetscObjectAddOptionsHandler() */
424: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject));
425: PetscOptionsEnd();
427: if (ts->trajectory) PetscCall(TSTrajectorySetFromOptions(ts->trajectory, ts));
429: /* why do we have to do this here and not during TSSetUp? */
430: PetscCall(TSGetSNES(ts, &ts->snes));
431: if (ts->problem_type == TS_LINEAR) {
432: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
433: if (!flg) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
434: }
435: PetscCall(SNESSetFromOptions(ts->snes));
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: /*@
440: TSGetTrajectory - Gets the trajectory from a `TS` if it exists
442: Collective
444: Input Parameter:
445: . ts - the `TS` context obtained from `TSCreate()`
447: Output Parameter:
448: . tr - the `TSTrajectory` object, if it exists
450: Level: advanced
452: Note:
453: This routine should be called after all `TS` options have been set
455: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectoryCreate()`
456: @*/
457: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
458: {
459: PetscFunctionBegin;
461: *tr = ts->trajectory;
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: /*@
466: TSSetSaveTrajectory - Causes the `TS` to save its solutions as it iterates forward in time in a `TSTrajectory` object
468: Collective
470: Input Parameter:
471: . ts - the `TS` context obtained from `TSCreate()`
473: Options Database Keys:
474: + -ts_save_trajectory - saves the trajectory to a file
475: - -ts_trajectory_type type - set trajectory type
477: Level: intermediate
479: Notes:
480: This routine should be called after all `TS` options have been set
482: The `TSTRAJECTORYVISUALIZATION` files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
483: MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m
485: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`
486: @*/
487: PetscErrorCode TSSetSaveTrajectory(TS ts)
488: {
489: PetscFunctionBegin;
491: if (!ts->trajectory) PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@
496: TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object
498: Collective
500: Input Parameter:
501: . ts - the `TS` context obtained from `TSCreate()`
503: Level: intermediate
505: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
506: @*/
507: PetscErrorCode TSResetTrajectory(TS ts)
508: {
509: PetscFunctionBegin;
511: if (ts->trajectory) {
512: PetscCall(TSTrajectoryDestroy(&ts->trajectory));
513: PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
514: }
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: /*@
519: TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from a `TS`
521: Collective
523: Input Parameter:
524: . ts - the `TS` context obtained from `TSCreate()`
526: Level: intermediate
528: .seealso: [](ch_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
529: @*/
530: PetscErrorCode TSRemoveTrajectory(TS ts)
531: {
532: PetscFunctionBegin;
534: if (ts->trajectory) PetscCall(TSTrajectoryDestroy(&ts->trajectory));
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: /*@
539: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
540: set with `TSSetRHSJacobian()`.
542: Collective
544: Input Parameters:
545: + ts - the `TS` context
546: . t - current timestep
547: - U - input vector
549: Output Parameters:
550: + A - Jacobian matrix
551: - B - optional matrix used to compute the preconditioner, often the same as `A`
553: Level: developer
555: Note:
556: Most users should not need to explicitly call this routine, as it
557: is used internally within the ODE integrators.
559: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()`
560: @*/
561: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
562: {
563: PetscObjectState Ustate;
564: PetscObjectId Uid;
565: DM dm;
566: DMTS tsdm;
567: TSRHSJacobianFn *rhsjacobianfunc;
568: void *ctx;
569: TSRHSFunctionFn *rhsfunction;
571: PetscFunctionBegin;
574: PetscCheckSameComm(ts, 1, U, 3);
575: PetscCall(TSGetDM(ts, &dm));
576: PetscCall(DMGetDMTS(dm, &tsdm));
577: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
578: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx));
579: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
580: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
582: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) PetscFunctionReturn(PETSC_SUCCESS);
584: PetscCheck(ts->rhsjacobian.shift == 0.0 || !ts->rhsjacobian.reuse, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Should not call TSComputeRHSJacobian() on a shifted matrix (shift=%lf) when RHSJacobian is reusable.", (double)ts->rhsjacobian.shift);
585: if (rhsjacobianfunc) {
586: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
587: PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
588: ts->rhsjacs++;
589: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
590: } else {
591: PetscCall(MatZeroEntries(A));
592: if (B && A != B) PetscCall(MatZeroEntries(B));
593: }
594: ts->rhsjacobian.time = t;
595: ts->rhsjacobian.shift = 0;
596: ts->rhsjacobian.scale = 1.;
597: PetscCall(PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid));
598: PetscCall(PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate));
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: /*@
603: TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS`
605: Collective
607: Input Parameters:
608: + ts - the `TS` context
609: . t - current time
610: - U - state vector
612: Output Parameter:
613: . y - right-hand side
615: Level: developer
617: Note:
618: Most users should not need to explicitly call this routine, as it
619: is used internally within the nonlinear solvers.
621: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
622: @*/
623: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
624: {
625: TSRHSFunctionFn *rhsfunction;
626: TSIFunctionFn *ifunction;
627: void *ctx;
628: DM dm;
630: PetscFunctionBegin;
634: PetscCall(TSGetDM(ts, &dm));
635: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx));
636: PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
638: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
640: if (rhsfunction) {
641: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, y, 0));
642: PetscCall(VecLockReadPush(U));
643: PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
644: PetscCall(VecLockReadPop(U));
645: ts->rhsfuncs++;
646: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, y, 0));
647: } else PetscCall(VecZeroEntries(y));
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
651: /*@
652: TSComputeSolutionFunction - Evaluates the solution function.
654: Collective
656: Input Parameters:
657: + ts - the `TS` context
658: - t - current time
660: Output Parameter:
661: . U - the solution
663: Level: developer
665: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
666: @*/
667: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
668: {
669: TSSolutionFn *solutionfunction;
670: void *ctx;
671: DM dm;
673: PetscFunctionBegin;
676: PetscCall(TSGetDM(ts, &dm));
677: PetscCall(DMTSGetSolutionFunction(dm, &solutionfunction, &ctx));
678: if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
681: /*@
682: TSComputeForcingFunction - Evaluates the forcing function.
684: Collective
686: Input Parameters:
687: + ts - the `TS` context
688: - t - current time
690: Output Parameter:
691: . U - the function value
693: Level: developer
695: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
696: @*/
697: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
698: {
699: void *ctx;
700: DM dm;
701: TSForcingFn *forcing;
703: PetscFunctionBegin;
706: PetscCall(TSGetDM(ts, &dm));
707: PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));
709: if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
714: {
715: Mat A, B;
716: TSIJacobianFn *ijacobian;
718: PetscFunctionBegin;
719: if (Arhs) *Arhs = NULL;
720: if (Brhs) *Brhs = NULL;
721: PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL));
722: if (Arhs) {
723: if (!ts->Arhs) {
724: if (ijacobian) {
725: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs));
726: PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN));
727: } else {
728: ts->Arhs = A;
729: PetscCall(PetscObjectReference((PetscObject)A));
730: }
731: } else {
732: PetscBool flg;
733: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
734: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
735: if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
736: PetscCall(PetscObjectDereference((PetscObject)ts->Arhs));
737: ts->Arhs = A;
738: PetscCall(PetscObjectReference((PetscObject)A));
739: }
740: }
741: *Arhs = ts->Arhs;
742: }
743: if (Brhs) {
744: if (!ts->Brhs) {
745: if (A != B) {
746: if (ijacobian) {
747: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs));
748: } else {
749: ts->Brhs = B;
750: PetscCall(PetscObjectReference((PetscObject)B));
751: }
752: } else {
753: PetscCall(PetscObjectReference((PetscObject)ts->Arhs));
754: ts->Brhs = ts->Arhs;
755: }
756: }
757: *Brhs = ts->Brhs;
758: }
759: PetscFunctionReturn(PETSC_SUCCESS);
760: }
762: /*@
763: TSComputeIFunction - Evaluates the DAE residual written in the implicit form F(t,U,Udot)=0
765: Collective
767: Input Parameters:
768: + ts - the `TS` context
769: . t - current time
770: . U - state vector
771: . Udot - time derivative of state vector
772: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSFunction should be kept separate
774: Output Parameter:
775: . Y - right-hand side
777: Level: developer
779: Note:
780: Most users should not need to explicitly call this routine, as it
781: is used internally within the nonlinear solvers.
783: If the user did not write their equations in implicit form, this
784: function recasts them in implicit form.
786: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
787: @*/
788: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
789: {
790: TSIFunctionFn *ifunction;
791: TSRHSFunctionFn *rhsfunction;
792: void *ctx;
793: DM dm;
795: PetscFunctionBegin;
801: PetscCall(TSGetDM(ts, &dm));
802: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
803: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
805: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
807: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, Udot, Y));
808: if (ifunction) {
809: PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
810: ts->ifuncs++;
811: }
812: if (imex) {
813: if (!ifunction) PetscCall(VecCopy(Udot, Y));
814: } else if (rhsfunction) {
815: if (ifunction) {
816: Vec Frhs;
818: PetscCall(DMGetGlobalVector(dm, &Frhs));
819: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
820: PetscCall(VecAXPY(Y, -1, Frhs));
821: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
822: } else {
823: PetscCall(TSComputeRHSFunction(ts, t, U, Y));
824: PetscCall(VecAYPX(Y, -1, Udot));
825: }
826: }
827: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, Udot, Y));
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }
831: /*
832: TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call `TSComputeRHSJacobian()` on it.
834: Note:
835: This routine is needed when one switches from `TSComputeIJacobian()` to `TSComputeRHSJacobian()` because the Jacobian matrix may be shifted or scaled in `TSComputeIJacobian()`.
837: */
838: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
839: {
840: PetscFunctionBegin;
842: PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
843: PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");
845: if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
846: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
847: if (B && B == ts->Brhs && A != B) {
848: if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
849: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
850: }
851: ts->rhsjacobian.shift = 0;
852: ts->rhsjacobian.scale = 1.;
853: PetscFunctionReturn(PETSC_SUCCESS);
854: }
856: /*@
857: TSComputeIJacobian - Evaluates the Jacobian of the DAE
859: Collective
861: Input Parameters:
862: + ts - the `TS` context
863: . t - current timestep
864: . U - state vector
865: . Udot - time derivative of state vector
866: . shift - shift to apply, see note below
867: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSJacobian should be kept separate
869: Output Parameters:
870: + A - Jacobian matrix
871: - B - matrix from which the preconditioner is constructed; often the same as `A`
873: Level: developer
875: Notes:
876: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
877: .vb
878: dF/dU + shift*dF/dUdot
879: .ve
880: Most users should not need to explicitly call this routine, as it
881: is used internally within the nonlinear solvers.
883: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
884: @*/
885: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
886: {
887: TSIJacobianFn *ijacobian;
888: TSRHSJacobianFn *rhsjacobian;
889: DM dm;
890: void *ctx;
892: PetscFunctionBegin;
899: PetscCall(TSGetDM(ts, &dm));
900: PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
901: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
903: PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
905: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
906: if (ijacobian) {
907: PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
908: ts->ijacs++;
909: }
910: if (imex) {
911: if (!ijacobian) { /* system was written as Udot = G(t,U) */
912: PetscBool assembled;
913: if (rhsjacobian) {
914: Mat Arhs = NULL;
915: PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
916: if (A == Arhs) {
917: 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 */
918: ts->rhsjacobian.time = PETSC_MIN_REAL;
919: }
920: }
921: PetscCall(MatZeroEntries(A));
922: PetscCall(MatAssembled(A, &assembled));
923: if (!assembled) {
924: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
925: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
926: }
927: PetscCall(MatShift(A, shift));
928: if (A != B) {
929: PetscCall(MatZeroEntries(B));
930: PetscCall(MatAssembled(B, &assembled));
931: if (!assembled) {
932: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
933: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
934: }
935: PetscCall(MatShift(B, shift));
936: }
937: }
938: } else {
939: Mat Arhs = NULL, Brhs = NULL;
941: /* RHSJacobian needs to be converted to part of IJacobian if exists */
942: if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
943: if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
944: PetscObjectState Ustate;
945: PetscObjectId Uid;
946: TSRHSFunctionFn *rhsfunction;
948: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
949: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
950: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
951: if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
952: ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */
953: PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
954: if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
955: } else {
956: PetscBool flg;
958: if (ts->rhsjacobian.reuse) { /* Undo the damage */
959: /* MatScale has a short path for this case.
960: However, this code path is taken the first time TSComputeRHSJacobian is called
961: and the matrices have not been assembled yet */
962: PetscCall(TSRecoverRHSJacobian(ts, A, B));
963: }
964: PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
965: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
966: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
967: if (!flg) {
968: PetscCall(MatScale(A, -1));
969: PetscCall(MatShift(A, shift));
970: }
971: if (A != B) {
972: PetscCall(MatScale(B, -1));
973: PetscCall(MatShift(B, shift));
974: }
975: }
976: ts->rhsjacobian.scale = -1;
977: ts->rhsjacobian.shift = shift;
978: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
979: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
980: PetscCall(MatZeroEntries(A));
981: PetscCall(MatShift(A, shift));
982: if (A != B) {
983: PetscCall(MatZeroEntries(B));
984: PetscCall(MatShift(B, shift));
985: }
986: }
987: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
988: PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
989: if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
990: }
991: }
992: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
993: PetscFunctionReturn(PETSC_SUCCESS);
994: }
996: /*@C
997: TSSetRHSFunction - Sets the routine for evaluating the function,
998: where U_t = G(t,u).
1000: Logically Collective
1002: Input Parameters:
1003: + ts - the `TS` context obtained from `TSCreate()`
1004: . r - vector to put the computed right-hand side (or `NULL` to have it created)
1005: . f - routine for evaluating the right-hand-side function
1006: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
1008: Level: beginner
1010: Note:
1011: You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.
1013: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1014: @*/
1015: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, void *ctx)
1016: {
1017: SNES snes;
1018: Vec ralloc = NULL;
1019: DM dm;
1021: PetscFunctionBegin;
1025: PetscCall(TSGetDM(ts, &dm));
1026: PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1027: PetscCall(TSGetSNES(ts, &snes));
1028: if (!r && !ts->dm && ts->vec_sol) {
1029: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1030: r = ralloc;
1031: }
1032: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1033: PetscCall(VecDestroy(&ralloc));
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: /*@C
1038: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1040: Logically Collective
1042: Input Parameters:
1043: + ts - the `TS` context obtained from `TSCreate()`
1044: . f - routine for evaluating the solution
1045: - ctx - [optional] user-defined context for private data for the
1046: function evaluation routine (may be `NULL`)
1048: Options Database Keys:
1049: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1050: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
1052: Level: intermediate
1054: Notes:
1055: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1056: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1057: create closed-form solutions with non-physical forcing terms.
1059: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1061: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1062: @*/
1063: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, void *ctx)
1064: {
1065: DM dm;
1067: PetscFunctionBegin;
1069: PetscCall(TSGetDM(ts, &dm));
1070: PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: /*@C
1075: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1077: Logically Collective
1079: Input Parameters:
1080: + ts - the `TS` context obtained from `TSCreate()`
1081: . func - routine for evaluating the forcing function
1082: - ctx - [optional] user-defined context for private data for the function evaluation routine
1083: (may be `NULL`)
1085: Level: intermediate
1087: Notes:
1088: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1089: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1090: definition of the problem you are solving and hence possibly introducing bugs.
1092: This replaces the ODE F(u,u_t,t) = 0 the `TS` is solving with F(u,u_t,t) - func(t) = 0
1094: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1095: parameters can be passed in the ctx variable.
1097: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1099: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
1100: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1101: @*/
1102: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, void *ctx)
1103: {
1104: DM dm;
1106: PetscFunctionBegin;
1108: PetscCall(TSGetDM(ts, &dm));
1109: PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1110: PetscFunctionReturn(PETSC_SUCCESS);
1111: }
1113: /*@C
1114: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1115: where U_t = G(U,t), as well as the location to store the matrix.
1117: Logically Collective
1119: Input Parameters:
1120: + ts - the `TS` context obtained from `TSCreate()`
1121: . Amat - (approximate) location to store Jacobian matrix entries computed by `f`
1122: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1123: . f - the Jacobian evaluation routine
1124: - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1126: Level: beginner
1128: Notes:
1129: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1131: The `TS` solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f()`
1132: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1134: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
1135: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
1136: @*/
1137: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, void *ctx)
1138: {
1139: SNES snes;
1140: DM dm;
1141: TSIJacobianFn *ijacobian;
1143: PetscFunctionBegin;
1147: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1148: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1150: PetscCall(TSGetDM(ts, &dm));
1151: PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1152: PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1153: PetscCall(TSGetSNES(ts, &snes));
1154: if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1155: if (Amat) {
1156: PetscCall(PetscObjectReference((PetscObject)Amat));
1157: PetscCall(MatDestroy(&ts->Arhs));
1158: ts->Arhs = Amat;
1159: }
1160: if (Pmat) {
1161: PetscCall(PetscObjectReference((PetscObject)Pmat));
1162: PetscCall(MatDestroy(&ts->Brhs));
1163: ts->Brhs = Pmat;
1164: }
1165: PetscFunctionReturn(PETSC_SUCCESS);
1166: }
1168: /*@C
1169: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1171: Logically Collective
1173: Input Parameters:
1174: + ts - the `TS` context obtained from `TSCreate()`
1175: . r - vector to hold the residual (or `NULL` to have it created internally)
1176: . f - the function evaluation routine
1177: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1179: Level: beginner
1181: Note:
1182: The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function.
1184: .seealso: [](ch_ts), `TS`, `TSIFunctionFn`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`,
1185: `TSSetIJacobian()`
1186: @*/
1187: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunctionFn *f, void *ctx)
1188: {
1189: SNES snes;
1190: Vec ralloc = NULL;
1191: DM dm;
1193: PetscFunctionBegin;
1197: PetscCall(TSGetDM(ts, &dm));
1198: PetscCall(DMTSSetIFunction(dm, f, ctx));
1200: PetscCall(TSGetSNES(ts, &snes));
1201: if (!r && !ts->dm && ts->vec_sol) {
1202: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1203: r = ralloc;
1204: }
1205: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1206: PetscCall(VecDestroy(&ralloc));
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: /*@C
1211: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.
1213: Not Collective
1215: Input Parameter:
1216: . ts - the `TS` context
1218: Output Parameters:
1219: + r - vector to hold residual (or `NULL`)
1220: . func - the function to compute residual (or `NULL`)
1221: - ctx - the function context (or `NULL`)
1223: Level: advanced
1225: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1226: @*/
1227: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunctionFn **func, void **ctx)
1228: {
1229: SNES snes;
1230: DM dm;
1232: PetscFunctionBegin;
1234: PetscCall(TSGetSNES(ts, &snes));
1235: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1236: PetscCall(TSGetDM(ts, &dm));
1237: PetscCall(DMTSGetIFunction(dm, func, ctx));
1238: PetscFunctionReturn(PETSC_SUCCESS);
1239: }
1241: /*@C
1242: TSGetRHSFunction - Returns the vector where the right-hand side is stored and the function/context to compute it.
1244: Not Collective
1246: Input Parameter:
1247: . ts - the `TS` context
1249: Output Parameters:
1250: + r - vector to hold computed right-hand side (or `NULL`)
1251: . func - the function to compute right-hand side (or `NULL`)
1252: - ctx - the function context (or `NULL`)
1254: Level: advanced
1256: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1257: @*/
1258: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunctionFn **func, void **ctx)
1259: {
1260: SNES snes;
1261: DM dm;
1263: PetscFunctionBegin;
1265: PetscCall(TSGetSNES(ts, &snes));
1266: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1267: PetscCall(TSGetDM(ts, &dm));
1268: PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1269: PetscFunctionReturn(PETSC_SUCCESS);
1270: }
1272: /*@C
1273: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1274: provided with `TSSetIFunction()`.
1276: Logically Collective
1278: Input Parameters:
1279: + ts - the `TS` context obtained from `TSCreate()`
1280: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
1281: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1282: . f - the Jacobian evaluation routine
1283: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1285: Level: beginner
1287: Notes:
1288: The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1290: If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
1291: space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.
1293: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1294: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1295: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1296: a and vector W depend on the integration method, step size, and past states. For example with
1297: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1298: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1300: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1302: The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f`
1303: You should not assume the values are the same in the next call to `f` as you set them in the previous call.
1305: In case `TSSetRHSJacobian()` is also used in conjunction with a fully-implicit solver,
1306: multilevel linear solvers, e.g. `PCMG`, will likely not work due to the way `TS` handles rhs matrices.
1308: .seealso: [](ch_ts), `TS`, `TSIJacobianFn`, `TSSetIFunction()`, `TSSetRHSJacobian()`,
1309: `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1310: @*/
1311: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobianFn *f, void *ctx)
1312: {
1313: SNES snes;
1314: DM dm;
1316: PetscFunctionBegin;
1320: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1321: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1323: PetscCall(TSGetDM(ts, &dm));
1324: PetscCall(DMTSSetIJacobian(dm, f, ctx));
1326: PetscCall(TSGetSNES(ts, &snes));
1327: PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1328: PetscFunctionReturn(PETSC_SUCCESS);
1329: }
1331: /*@
1332: TSRHSJacobianSetReuse - restore the RHS Jacobian before calling the user-provided `TSRHSJacobianFn` function again
1334: Logically Collective
1336: Input Parameters:
1337: + ts - `TS` context obtained from `TSCreate()`
1338: - reuse - `PETSC_TRUE` if the RHS Jacobian
1340: Level: intermediate
1342: Notes:
1343: Without this flag, `TS` will change the sign and shift the RHS Jacobian for a
1344: finite-time-step implicit solve, in which case the user function will need to recompute the
1345: entire Jacobian. The `reuse `flag must be set if the evaluation function assumes that the
1346: matrix entries have not been changed by the `TS`.
1348: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1349: @*/
1350: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1351: {
1352: PetscFunctionBegin;
1353: ts->rhsjacobian.reuse = reuse;
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: /*@C
1358: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1360: Logically Collective
1362: Input Parameters:
1363: + ts - the `TS` context obtained from `TSCreate()`
1364: . F - vector to hold the residual (or `NULL` to have it created internally)
1365: . fun - the function evaluation routine
1366: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1368: Level: beginner
1370: .seealso: [](ch_ts), `TS`, `TSI2FunctionFn`, `TSSetI2Jacobian()`, `TSSetIFunction()`,
1371: `TSCreate()`, `TSSetRHSFunction()`
1372: @*/
1373: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2FunctionFn *fun, void *ctx)
1374: {
1375: DM dm;
1377: PetscFunctionBegin;
1380: PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1381: PetscCall(TSGetDM(ts, &dm));
1382: PetscCall(DMTSSetI2Function(dm, fun, ctx));
1383: PetscFunctionReturn(PETSC_SUCCESS);
1384: }
1386: /*@C
1387: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.
1389: Not Collective
1391: Input Parameter:
1392: . ts - the `TS` context
1394: Output Parameters:
1395: + r - vector to hold residual (or `NULL`)
1396: . fun - the function to compute residual (or `NULL`)
1397: - ctx - the function context (or `NULL`)
1399: Level: advanced
1401: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1402: @*/
1403: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2FunctionFn **fun, void **ctx)
1404: {
1405: SNES snes;
1406: DM dm;
1408: PetscFunctionBegin;
1410: PetscCall(TSGetSNES(ts, &snes));
1411: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1412: PetscCall(TSGetDM(ts, &dm));
1413: PetscCall(DMTSGetI2Function(dm, fun, ctx));
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: /*@C
1418: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1419: where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.
1421: Logically Collective
1423: Input Parameters:
1424: + ts - the `TS` context obtained from `TSCreate()`
1425: . J - matrix to hold the Jacobian values
1426: . P - matrix for constructing the preconditioner (may be same as `J`)
1427: . jac - the Jacobian evaluation routine, see `TSI2JacobianFn` for the calling sequence
1428: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1430: Level: beginner
1432: Notes:
1433: The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1435: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1436: 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.
1437: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1438: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1440: .seealso: [](ch_ts), `TS`, `TSI2JacobianFn`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1441: @*/
1442: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2JacobianFn *jac, void *ctx)
1443: {
1444: DM dm;
1446: PetscFunctionBegin;
1450: PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1451: PetscCall(TSGetDM(ts, &dm));
1452: PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1453: PetscFunctionReturn(PETSC_SUCCESS);
1454: }
1456: /*@C
1457: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1459: Not Collective, but parallel objects are returned if `TS` is parallel
1461: Input Parameter:
1462: . ts - The `TS` context obtained from `TSCreate()`
1464: Output Parameters:
1465: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1466: . P - The matrix from which the preconditioner is constructed, often the same as `J`
1467: . jac - The function to compute the Jacobian matrices
1468: - ctx - User-defined context for Jacobian evaluation routine
1470: Level: advanced
1472: Note:
1473: You can pass in `NULL` for any return argument you do not need.
1475: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1476: @*/
1477: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2JacobianFn **jac, void **ctx)
1478: {
1479: SNES snes;
1480: DM dm;
1482: PetscFunctionBegin;
1483: PetscCall(TSGetSNES(ts, &snes));
1484: PetscCall(SNESSetUpMatrices(snes));
1485: PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1486: PetscCall(TSGetDM(ts, &dm));
1487: PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1488: PetscFunctionReturn(PETSC_SUCCESS);
1489: }
1491: /*@
1492: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1494: Collective
1496: Input Parameters:
1497: + ts - the `TS` context
1498: . t - current time
1499: . U - state vector
1500: . V - time derivative of state vector (U_t)
1501: - A - second time derivative of state vector (U_tt)
1503: Output Parameter:
1504: . F - the residual vector
1506: Level: developer
1508: Note:
1509: Most users should not need to explicitly call this routine, as it
1510: is used internally within the nonlinear solvers.
1512: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1513: @*/
1514: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1515: {
1516: DM dm;
1517: TSI2FunctionFn *I2Function;
1518: void *ctx;
1519: TSRHSFunctionFn *rhsfunction;
1521: PetscFunctionBegin;
1528: PetscCall(TSGetDM(ts, &dm));
1529: PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1530: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
1532: if (!I2Function) {
1533: PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1534: PetscFunctionReturn(PETSC_SUCCESS);
1535: }
1537: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, V, F));
1539: PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));
1541: if (rhsfunction) {
1542: Vec Frhs;
1544: PetscCall(DMGetGlobalVector(dm, &Frhs));
1545: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1546: PetscCall(VecAXPY(F, -1, Frhs));
1547: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
1548: }
1550: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, V, F));
1551: PetscFunctionReturn(PETSC_SUCCESS);
1552: }
1554: /*@
1555: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1557: Collective
1559: Input Parameters:
1560: + ts - the `TS` context
1561: . t - current timestep
1562: . U - state vector
1563: . V - time derivative of state vector
1564: . A - second time derivative of state vector
1565: . shiftV - shift to apply, see note below
1566: - shiftA - shift to apply, see note below
1568: Output Parameters:
1569: + J - Jacobian matrix
1570: - P - optional matrix used to construct the preconditioner
1572: Level: developer
1574: Notes:
1575: If $F(t,U,V,A) = 0$ is the DAE, the required Jacobian is
1577: $$
1578: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1579: $$
1581: Most users should not need to explicitly call this routine, as it
1582: is used internally within the ODE integrators.
1584: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1585: @*/
1586: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1587: {
1588: DM dm;
1589: TSI2JacobianFn *I2Jacobian;
1590: void *ctx;
1591: TSRHSJacobianFn *rhsjacobian;
1593: PetscFunctionBegin;
1601: PetscCall(TSGetDM(ts, &dm));
1602: PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1603: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1605: if (!I2Jacobian) {
1606: PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1607: PetscFunctionReturn(PETSC_SUCCESS);
1608: }
1610: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, J, P));
1611: PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1612: if (rhsjacobian) {
1613: Mat Jrhs, Prhs;
1614: PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1615: PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1616: PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1617: if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1618: }
1620: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, J, P));
1621: PetscFunctionReturn(PETSC_SUCCESS);
1622: }
1624: /*@C
1625: TSSetTransientVariable - sets function to transform from state to transient variables
1627: Logically Collective
1629: Input Parameters:
1630: + ts - time stepping context on which to change the transient variable
1631: . tvar - a function that transforms to transient variables, see `TSTransientVariableFn` for the calling sequence
1632: - ctx - a context for tvar
1634: Level: advanced
1636: Notes:
1637: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1638: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1639: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1640: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1641: evaluated via the chain rule, as in
1642: .vb
1643: dF/dP + shift * dF/dCdot dC/dP.
1644: .ve
1646: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSTransientVariableFn`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1647: @*/
1648: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariableFn *tvar, void *ctx)
1649: {
1650: DM dm;
1652: PetscFunctionBegin;
1654: PetscCall(TSGetDM(ts, &dm));
1655: PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1656: PetscFunctionReturn(PETSC_SUCCESS);
1657: }
1659: /*@
1660: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1662: Logically Collective
1664: Input Parameters:
1665: + ts - TS on which to compute
1666: - U - state vector to be transformed to transient variables
1668: Output Parameter:
1669: . C - transient (conservative) variable
1671: Level: developer
1673: Developer Notes:
1674: If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed.
1675: This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are
1676: being used.
1678: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1679: @*/
1680: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1681: {
1682: DM dm;
1683: DMTS dmts;
1685: PetscFunctionBegin;
1688: PetscCall(TSGetDM(ts, &dm));
1689: PetscCall(DMGetDMTS(dm, &dmts));
1690: if (dmts->ops->transientvar) {
1692: PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1693: }
1694: PetscFunctionReturn(PETSC_SUCCESS);
1695: }
1697: /*@
1698: TSHasTransientVariable - determine whether transient variables have been set
1700: Logically Collective
1702: Input Parameter:
1703: . ts - `TS` on which to compute
1705: Output Parameter:
1706: . has - `PETSC_TRUE` if transient variables have been set
1708: Level: developer
1710: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1711: @*/
1712: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1713: {
1714: DM dm;
1715: DMTS dmts;
1717: PetscFunctionBegin;
1719: PetscCall(TSGetDM(ts, &dm));
1720: PetscCall(DMGetDMTS(dm, &dmts));
1721: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1722: PetscFunctionReturn(PETSC_SUCCESS);
1723: }
1725: /*@
1726: TS2SetSolution - Sets the initial solution and time derivative vectors
1727: for use by the `TS` routines handling second order equations.
1729: Logically Collective
1731: Input Parameters:
1732: + ts - the `TS` context obtained from `TSCreate()`
1733: . u - the solution vector
1734: - v - the time derivative vector
1736: Level: beginner
1738: .seealso: [](ch_ts), `TS`
1739: @*/
1740: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1741: {
1742: PetscFunctionBegin;
1746: PetscCall(TSSetSolution(ts, u));
1747: PetscCall(PetscObjectReference((PetscObject)v));
1748: PetscCall(VecDestroy(&ts->vec_dot));
1749: ts->vec_dot = v;
1750: PetscFunctionReturn(PETSC_SUCCESS);
1751: }
1753: /*@
1754: TS2GetSolution - Returns the solution and time derivative at the present timestep
1755: for second order equations.
1757: Not Collective
1759: Input Parameter:
1760: . ts - the `TS` context obtained from `TSCreate()`
1762: Output Parameters:
1763: + u - the vector containing the solution
1764: - v - the vector containing the time derivative
1766: Level: intermediate
1768: Notes:
1769: It is valid to call this routine inside the function
1770: that you are evaluating in order to move to the new timestep. This vector not
1771: changed until the solution at the next timestep has been calculated.
1773: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1774: @*/
1775: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1776: {
1777: PetscFunctionBegin;
1779: if (u) PetscAssertPointer(u, 2);
1780: if (v) PetscAssertPointer(v, 3);
1781: if (u) *u = ts->vec_sol;
1782: if (v) *v = ts->vec_dot;
1783: PetscFunctionReturn(PETSC_SUCCESS);
1784: }
1786: /*@
1787: TSLoad - Loads a `TS` that has been stored in binary with `TSView()`.
1789: Collective
1791: Input Parameters:
1792: + ts - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1793: some related function before a call to `TSLoad()`.
1794: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1796: Level: intermediate
1798: Note:
1799: The type is determined by the data in the file, any type set into the `TS` before this call is ignored.
1801: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1802: @*/
1803: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1804: {
1805: PetscBool isbinary;
1806: PetscInt classid;
1807: char type[256];
1808: DMTS sdm;
1809: DM dm;
1811: PetscFunctionBegin;
1814: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1815: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1817: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1818: PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1819: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1820: PetscCall(TSSetType(ts, type));
1821: PetscTryTypeMethod(ts, load, viewer);
1822: PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1823: PetscCall(DMLoad(dm, viewer));
1824: PetscCall(TSSetDM(ts, dm));
1825: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1826: PetscCall(VecLoad(ts->vec_sol, viewer));
1827: PetscCall(DMGetDMTS(ts->dm, &sdm));
1828: PetscCall(DMTSLoad(sdm, viewer));
1829: PetscFunctionReturn(PETSC_SUCCESS);
1830: }
1832: #include <petscdraw.h>
1833: #if defined(PETSC_HAVE_SAWS)
1834: #include <petscviewersaws.h>
1835: #endif
1837: /*@
1838: TSViewFromOptions - View a `TS` based on values in the options database
1840: Collective
1842: Input Parameters:
1843: + ts - the `TS` context
1844: . obj - Optional object that provides the prefix for the options database keys
1845: - name - command line option string to be passed by user
1847: Level: intermediate
1849: .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1850: @*/
1851: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1852: {
1853: PetscFunctionBegin;
1855: PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1856: PetscFunctionReturn(PETSC_SUCCESS);
1857: }
1859: /*@
1860: TSView - Prints the `TS` data structure.
1862: Collective
1864: Input Parameters:
1865: + ts - the `TS` context obtained from `TSCreate()`
1866: - viewer - visualization context
1868: Options Database Key:
1869: . -ts_view - calls `TSView()` at end of `TSStep()`
1871: Level: beginner
1873: Notes:
1874: The available visualization contexts include
1875: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1876: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1877: output where only the first processor opens
1878: the file. All other processors send their
1879: data to the first processor to print.
1881: The user can open an alternative visualization context with
1882: `PetscViewerASCIIOpen()` - output to a specified file.
1884: In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).
1886: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1887: @*/
1888: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1889: {
1890: TSType type;
1891: PetscBool iascii, isstring, isundials, isbinary, isdraw;
1892: DMTS sdm;
1893: #if defined(PETSC_HAVE_SAWS)
1894: PetscBool issaws;
1895: #endif
1897: PetscFunctionBegin;
1899: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1901: PetscCheckSameComm(ts, 1, viewer, 2);
1903: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1904: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1905: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1906: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1907: #if defined(PETSC_HAVE_SAWS)
1908: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1909: #endif
1910: if (iascii) {
1911: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer));
1912: if (ts->ops->view) {
1913: PetscCall(PetscViewerASCIIPushTab(viewer));
1914: PetscUseTypeMethod(ts, view, viewer);
1915: PetscCall(PetscViewerASCIIPopTab(viewer));
1916: }
1917: if (ts->max_steps < PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps));
1918: if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time));
1919: if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs));
1920: if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs));
1921: if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs));
1922: if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs));
1923: if (ts->usessnes) {
1924: PetscBool lin;
1925: if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its));
1926: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its));
1927: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
1928: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures));
1929: }
1930: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject));
1931: if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, "));
1932: else PetscCall(PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol));
1933: if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of absolute error tolerances\n"));
1934: else PetscCall(PetscViewerASCIIPrintf(viewer, " using absolute error tolerance of %g\n", (double)ts->atol));
1935: PetscCall(PetscViewerASCIIPushTab(viewer));
1936: PetscCall(TSAdaptView(ts->adapt, viewer));
1937: PetscCall(PetscViewerASCIIPopTab(viewer));
1938: } else if (isstring) {
1939: PetscCall(TSGetType(ts, &type));
1940: PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type));
1941: PetscTryTypeMethod(ts, view, viewer);
1942: } else if (isbinary) {
1943: PetscInt classid = TS_FILE_CLASSID;
1944: MPI_Comm comm;
1945: PetscMPIInt rank;
1946: char type[256];
1948: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
1949: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1950: if (rank == 0) {
1951: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1952: PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
1953: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1954: }
1955: PetscTryTypeMethod(ts, view, viewer);
1956: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1957: PetscCall(DMView(ts->dm, viewer));
1958: PetscCall(VecView(ts->vec_sol, viewer));
1959: PetscCall(DMGetDMTS(ts->dm, &sdm));
1960: PetscCall(DMTSView(sdm, viewer));
1961: } else if (isdraw) {
1962: PetscDraw draw;
1963: char str[36];
1964: PetscReal x, y, bottom, h;
1966: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1967: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1968: PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
1969: PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
1970: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
1971: bottom = y - h;
1972: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1973: PetscTryTypeMethod(ts, view, viewer);
1974: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1975: if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
1976: PetscCall(PetscDrawPopCurrentPoint(draw));
1977: #if defined(PETSC_HAVE_SAWS)
1978: } else if (issaws) {
1979: PetscMPIInt rank;
1980: const char *name;
1982: PetscCall(PetscObjectGetName((PetscObject)ts, &name));
1983: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1984: if (!((PetscObject)ts)->amsmem && rank == 0) {
1985: char dir[1024];
1987: PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
1988: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
1989: PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
1990: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
1991: PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
1992: }
1993: PetscTryTypeMethod(ts, view, viewer);
1994: #endif
1995: }
1996: if (ts->snes && ts->usessnes) {
1997: PetscCall(PetscViewerASCIIPushTab(viewer));
1998: PetscCall(SNESView(ts->snes, viewer));
1999: PetscCall(PetscViewerASCIIPopTab(viewer));
2000: }
2001: PetscCall(DMGetDMTS(ts->dm, &sdm));
2002: PetscCall(DMTSView(sdm, viewer));
2004: PetscCall(PetscViewerASCIIPushTab(viewer));
2005: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials));
2006: PetscCall(PetscViewerASCIIPopTab(viewer));
2007: PetscFunctionReturn(PETSC_SUCCESS);
2008: }
2010: /*@
2011: TSSetApplicationContext - Sets an optional user-defined context for the timesteppers that may be accessed, for example inside the user provided
2012: `TS` callbacks with `TSGetApplicationContext()`
2014: Logically Collective
2016: Input Parameters:
2017: + ts - the `TS` context obtained from `TSCreate()`
2018: - ctx - user context
2020: Level: intermediate
2022: Fortran Note:
2023: This only works when `ctx` is a Fortran derived type (it cannot be a `PetscObject`), we recommend writing a Fortran interface definition for this
2024: function that tells the Fortran compiler the derived data type that is passed in as the `ctx` argument. See `TSGetApplicationContext()` for
2025: an example.
2027: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2028: @*/
2029: PetscErrorCode TSSetApplicationContext(TS ts, PeCtx ctx)
2030: {
2031: PetscFunctionBegin;
2033: ts->ctx = ctx;
2034: PetscFunctionReturn(PETSC_SUCCESS);
2035: }
2037: /*@
2038: TSGetApplicationContext - Gets the user-defined context for the
2039: timestepper that was set with `TSSetApplicationContext()`
2041: Not Collective
2043: Input Parameter:
2044: . ts - the `TS` context obtained from `TSCreate()`
2046: Output Parameter:
2047: . ctx - a pointer to the user context
2049: Level: intermediate
2051: Fortran Notes:
2052: This only works when the context is a Fortran derived type (it cannot be a `PetscObject`) and you **must** write a Fortran interface definition for this
2053: function that tells the Fortran compiler the derived data type that is returned as the `ctx` argument. For example,
2054: .vb
2055: Interface TSGetApplicationContext
2056: Subroutine TSGetApplicationContext(ts,ctx,ierr)
2057: #include <petsc/finclude/petscts.h>
2058: use petscts
2059: TS ts
2060: type(tUsertype), pointer :: ctx
2061: PetscErrorCode ierr
2062: End Subroutine
2063: End Interface TSGetApplicationContext
2064: .ve
2066: The prototype for `ctx` must be
2067: .vb
2068: type(tUsertype), pointer :: ctx
2069: .ve
2071: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2072: @*/
2073: PetscErrorCode TSGetApplicationContext(TS ts, void *ctx)
2074: {
2075: PetscFunctionBegin;
2077: *(void **)ctx = ts->ctx;
2078: PetscFunctionReturn(PETSC_SUCCESS);
2079: }
2081: /*@
2082: TSGetStepNumber - Gets the number of time steps completed.
2084: Not Collective
2086: Input Parameter:
2087: . ts - the `TS` context obtained from `TSCreate()`
2089: Output Parameter:
2090: . steps - number of steps completed so far
2092: Level: intermediate
2094: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2095: @*/
2096: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2097: {
2098: PetscFunctionBegin;
2100: PetscAssertPointer(steps, 2);
2101: *steps = ts->steps;
2102: PetscFunctionReturn(PETSC_SUCCESS);
2103: }
2105: /*@
2106: TSSetStepNumber - Sets the number of steps completed.
2108: Logically Collective
2110: Input Parameters:
2111: + ts - the `TS` context
2112: - steps - number of steps completed so far
2114: Level: developer
2116: Note:
2117: For most uses of the `TS` solvers the user need not explicitly call
2118: `TSSetStepNumber()`, as the step counter is appropriately updated in
2119: `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2120: reinitialize timestepping by setting the step counter to zero (and time
2121: to the initial time) to solve a similar problem with different initial
2122: conditions or parameters. Other possible use case is to continue
2123: timestepping from a previously interrupted run in such a way that `TS`
2124: monitors will be called with a initial nonzero step counter.
2126: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2127: @*/
2128: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2129: {
2130: PetscFunctionBegin;
2133: PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2134: ts->steps = steps;
2135: PetscFunctionReturn(PETSC_SUCCESS);
2136: }
2138: /*@
2139: TSSetTimeStep - Allows one to reset the timestep at any time,
2140: useful for simple pseudo-timestepping codes.
2142: Logically Collective
2144: Input Parameters:
2145: + ts - the `TS` context obtained from `TSCreate()`
2146: - time_step - the size of the timestep
2148: Level: intermediate
2150: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2151: @*/
2152: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2153: {
2154: PetscFunctionBegin;
2157: ts->time_step = time_step;
2158: PetscFunctionReturn(PETSC_SUCCESS);
2159: }
2161: /*@
2162: TSSetExactFinalTime - Determines whether to adapt the final time step to
2163: match the exact final time, to interpolate the solution to the exact final time,
2164: or to just return at the final time `TS` computed (which may be slightly larger
2165: than the requested final time).
2167: Logically Collective
2169: Input Parameters:
2170: + ts - the time-step context
2171: - eftopt - exact final time option
2172: .vb
2173: TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded, just use it
2174: TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time if the final time is exceeded
2175: TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to ensure the computed final time exactly equals the requested final time
2176: .ve
2178: Options Database Key:
2179: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step approach at runtime
2181: Level: beginner
2183: Note:
2184: If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2185: then the final time you selected.
2187: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2188: @*/
2189: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2190: {
2191: PetscFunctionBegin;
2194: ts->exact_final_time = eftopt;
2195: PetscFunctionReturn(PETSC_SUCCESS);
2196: }
2198: /*@
2199: TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`
2201: Not Collective
2203: Input Parameter:
2204: . ts - the `TS` context
2206: Output Parameter:
2207: . eftopt - exact final time option
2209: Level: beginner
2211: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2212: @*/
2213: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2214: {
2215: PetscFunctionBegin;
2217: PetscAssertPointer(eftopt, 2);
2218: *eftopt = ts->exact_final_time;
2219: PetscFunctionReturn(PETSC_SUCCESS);
2220: }
2222: /*@
2223: TSGetTimeStep - Gets the current timestep size.
2225: Not Collective
2227: Input Parameter:
2228: . ts - the `TS` context obtained from `TSCreate()`
2230: Output Parameter:
2231: . dt - the current timestep size
2233: Level: intermediate
2235: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2236: @*/
2237: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2238: {
2239: PetscFunctionBegin;
2241: PetscAssertPointer(dt, 2);
2242: *dt = ts->time_step;
2243: PetscFunctionReturn(PETSC_SUCCESS);
2244: }
2246: /*@
2247: TSGetSolution - Returns the solution at the present timestep. It
2248: is valid to call this routine inside the function that you are evaluating
2249: in order to move to the new timestep. This vector not changed until
2250: the solution at the next timestep has been calculated.
2252: Not Collective, but v returned is parallel if ts is parallel
2254: Input Parameter:
2255: . ts - the `TS` context obtained from `TSCreate()`
2257: Output Parameter:
2258: . v - the vector containing the solution
2260: Level: intermediate
2262: Note:
2263: If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
2264: final time. It returns the solution at the next timestep.
2266: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2267: @*/
2268: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2269: {
2270: PetscFunctionBegin;
2272: PetscAssertPointer(v, 2);
2273: *v = ts->vec_sol;
2274: PetscFunctionReturn(PETSC_SUCCESS);
2275: }
2277: /*@
2278: TSGetSolutionComponents - Returns any solution components at the present
2279: timestep, if available for the time integration method being used.
2280: Solution components are quantities that share the same size and
2281: structure as the solution vector.
2283: Not Collective, but v returned is parallel if ts is parallel
2285: Input Parameters:
2286: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2287: . n - If v is `NULL`, then the number of solution components is
2288: returned through n, else the n-th solution component is
2289: returned in v.
2290: - v - the vector containing the n-th solution component
2291: (may be `NULL` to use this function to find out
2292: the number of solutions components).
2294: Level: advanced
2296: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2297: @*/
2298: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2299: {
2300: PetscFunctionBegin;
2302: if (!ts->ops->getsolutioncomponents) *n = 0;
2303: else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2304: PetscFunctionReturn(PETSC_SUCCESS);
2305: }
2307: /*@
2308: TSGetAuxSolution - Returns an auxiliary solution at the present
2309: timestep, if available for the time integration method being used.
2311: Not Collective, but v returned is parallel if ts is parallel
2313: Input Parameters:
2314: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2315: - v - the vector containing the auxiliary solution
2317: Level: intermediate
2319: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2320: @*/
2321: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2322: {
2323: PetscFunctionBegin;
2325: if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2326: else PetscCall(VecZeroEntries(*v));
2327: PetscFunctionReturn(PETSC_SUCCESS);
2328: }
2330: /*@
2331: TSGetTimeError - Returns the estimated error vector, if the chosen
2332: `TSType` has an error estimation functionality and `TSSetTimeError()` was called
2334: Not Collective, but v returned is parallel if ts is parallel
2336: Input Parameters:
2337: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2338: . n - current estimate (n=0) or previous one (n=-1)
2339: - v - the vector containing the error (same size as the solution).
2341: Level: intermediate
2343: Note:
2344: MUST call after `TSSetUp()`
2346: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2347: @*/
2348: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2349: {
2350: PetscFunctionBegin;
2352: if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2353: else PetscCall(VecZeroEntries(*v));
2354: PetscFunctionReturn(PETSC_SUCCESS);
2355: }
2357: /*@
2358: TSSetTimeError - Sets the estimated error vector, if the chosen
2359: `TSType` has an error estimation functionality. This can be used
2360: to restart such a time integrator with a given error vector.
2362: Not Collective, but v returned is parallel if ts is parallel
2364: Input Parameters:
2365: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2366: - v - the vector containing the error (same size as the solution).
2368: Level: intermediate
2370: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
2371: @*/
2372: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2373: {
2374: PetscFunctionBegin;
2376: PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2377: PetscTryTypeMethod(ts, settimeerror, v);
2378: PetscFunctionReturn(PETSC_SUCCESS);
2379: }
2381: /* ----- Routines to initialize and destroy a timestepper ---- */
2382: /*@
2383: TSSetProblemType - Sets the type of problem to be solved.
2385: Not collective
2387: Input Parameters:
2388: + ts - The `TS`
2389: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2390: .vb
2391: U_t - A U = 0 (linear)
2392: U_t - A(t) U = 0 (linear)
2393: F(t,U,U_t) = 0 (nonlinear)
2394: .ve
2396: Level: beginner
2398: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2399: @*/
2400: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2401: {
2402: PetscFunctionBegin;
2404: ts->problem_type = type;
2405: if (type == TS_LINEAR) {
2406: SNES snes;
2407: PetscCall(TSGetSNES(ts, &snes));
2408: PetscCall(SNESSetType(snes, SNESKSPONLY));
2409: }
2410: PetscFunctionReturn(PETSC_SUCCESS);
2411: }
2413: /*@
2414: TSGetProblemType - Gets the type of problem to be solved.
2416: Not collective
2418: Input Parameter:
2419: . ts - The `TS`
2421: Output Parameter:
2422: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2423: .vb
2424: M U_t = A U
2425: M(t) U_t = A(t) U
2426: F(t,U,U_t)
2427: .ve
2429: Level: beginner
2431: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2432: @*/
2433: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2434: {
2435: PetscFunctionBegin;
2437: PetscAssertPointer(type, 2);
2438: *type = ts->problem_type;
2439: PetscFunctionReturn(PETSC_SUCCESS);
2440: }
2442: /*
2443: Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2444: */
2445: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2446: {
2447: PetscBool isnone;
2449: PetscFunctionBegin;
2450: PetscCall(TSGetAdapt(ts, &ts->adapt));
2451: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2453: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2454: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2455: else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2456: PetscFunctionReturn(PETSC_SUCCESS);
2457: }
2459: /*@
2460: TSSetUp - Sets up the internal data structures for the later use of a timestepper.
2462: Collective
2464: Input Parameter:
2465: . ts - the `TS` context obtained from `TSCreate()`
2467: Level: advanced
2469: Note:
2470: For basic use of the `TS` solvers the user need not explicitly call
2471: `TSSetUp()`, since these actions will automatically occur during
2472: the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this
2473: phase separately, `TSSetUp()` should be called after `TSCreate()`
2474: and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.
2476: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2477: @*/
2478: PetscErrorCode TSSetUp(TS ts)
2479: {
2480: DM dm;
2481: PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2482: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2483: TSIFunctionFn *ifun;
2484: TSIJacobianFn *ijac;
2485: TSI2JacobianFn *i2jac;
2486: TSRHSJacobianFn *rhsjac;
2488: PetscFunctionBegin;
2490: if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
2492: if (!((PetscObject)ts)->type_name) {
2493: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2494: PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2495: }
2497: if (!ts->vec_sol) {
2498: PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2499: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2500: }
2502: if (ts->eval_times) {
2503: if (!ts->eval_times->sol_vecs) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
2504: if (!ts->eval_times->sol_times) PetscCall(PetscMalloc1(ts->eval_times->num_time_points, &ts->eval_times->sol_times));
2505: }
2506: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2507: PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2508: ts->Jacp = ts->Jacprhs;
2509: }
2511: if (ts->quadraturets) {
2512: PetscCall(TSSetUp(ts->quadraturets));
2513: PetscCall(VecDestroy(&ts->vec_costintegrand));
2514: PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2515: }
2517: PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
2518: if (rhsjac == TSComputeRHSJacobianConstant) {
2519: Mat Amat, Pmat;
2520: SNES snes;
2521: PetscCall(TSGetSNES(ts, &snes));
2522: PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
2523: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2524: * have displaced the RHS matrix */
2525: if (Amat && Amat == ts->Arhs) {
2526: /* 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 */
2527: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
2528: PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
2529: PetscCall(MatDestroy(&Amat));
2530: }
2531: if (Pmat && Pmat == ts->Brhs) {
2532: PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
2533: PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
2534: PetscCall(MatDestroy(&Pmat));
2535: }
2536: }
2538: PetscCall(TSGetAdapt(ts, &ts->adapt));
2539: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2541: PetscTryTypeMethod(ts, setup);
2543: PetscCall(TSSetExactFinalTimeDefault(ts));
2545: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2546: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2547: */
2548: PetscCall(TSGetDM(ts, &dm));
2549: PetscCall(DMSNESGetFunction(dm, &func, NULL));
2550: if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));
2552: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2553: Otherwise, the SNES will use coloring internally to form the Jacobian.
2554: */
2555: PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2556: PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2557: PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2558: PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2559: if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));
2561: /* if time integration scheme has a starting method, call it */
2562: PetscTryTypeMethod(ts, startingmethod);
2564: ts->setupcalled = PETSC_TRUE;
2565: PetscFunctionReturn(PETSC_SUCCESS);
2566: }
2568: /*@
2569: 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
2571: Collective
2573: Input Parameter:
2574: . ts - the `TS` context obtained from `TSCreate()`
2576: Level: developer
2578: Notes:
2579: Any options set on the `TS` object, including those set with `TSSetFromOptions()` remain.
2581: See also `TSSetResize()` to change the size of the system being integrated (for example by adaptive mesh refinement) during the time integration.
2583: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSetResize()`
2584: @*/
2585: PetscErrorCode TSReset(TS ts)
2586: {
2587: TS_RHSSplitLink ilink = ts->tsrhssplit, next;
2589: PetscFunctionBegin;
2592: PetscTryTypeMethod(ts, reset);
2593: if (ts->snes) PetscCall(SNESReset(ts->snes));
2594: if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));
2596: PetscCall(MatDestroy(&ts->Arhs));
2597: PetscCall(MatDestroy(&ts->Brhs));
2598: PetscCall(VecDestroy(&ts->Frhs));
2599: PetscCall(VecDestroy(&ts->vec_sol));
2600: PetscCall(VecDestroy(&ts->vec_sol0));
2601: PetscCall(VecDestroy(&ts->vec_dot));
2602: PetscCall(VecDestroy(&ts->vatol));
2603: PetscCall(VecDestroy(&ts->vrtol));
2604: PetscCall(VecDestroyVecs(ts->nwork, &ts->work));
2606: PetscCall(MatDestroy(&ts->Jacprhs));
2607: PetscCall(MatDestroy(&ts->Jacp));
2608: if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2609: if (ts->quadraturets) {
2610: PetscCall(TSReset(ts->quadraturets));
2611: PetscCall(VecDestroy(&ts->vec_costintegrand));
2612: }
2613: while (ilink) {
2614: next = ilink->next;
2615: PetscCall(TSDestroy(&ilink->ts));
2616: PetscCall(PetscFree(ilink->splitname));
2617: PetscCall(ISDestroy(&ilink->is));
2618: PetscCall(PetscFree(ilink));
2619: ilink = next;
2620: }
2621: ts->tsrhssplit = NULL;
2622: ts->num_rhs_splits = 0;
2623: if (ts->eval_times) {
2624: PetscCall(PetscFree(ts->eval_times->time_points));
2625: PetscCall(PetscFree(ts->eval_times->sol_times));
2626: PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
2627: PetscCall(PetscFree(ts->eval_times));
2628: }
2629: ts->rhsjacobian.time = PETSC_MIN_REAL;
2630: ts->rhsjacobian.scale = 1.0;
2631: ts->ijacobian.shift = 1.0;
2632: ts->setupcalled = PETSC_FALSE;
2633: PetscFunctionReturn(PETSC_SUCCESS);
2634: }
2636: static PetscErrorCode TSResizeReset(TS);
2638: /*@
2639: TSDestroy - Destroys the timestepper context that was created
2640: with `TSCreate()`.
2642: Collective
2644: Input Parameter:
2645: . ts - the `TS` context obtained from `TSCreate()`
2647: Level: beginner
2649: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2650: @*/
2651: PetscErrorCode TSDestroy(TS *ts)
2652: {
2653: PetscFunctionBegin;
2654: if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2656: if (--((PetscObject)*ts)->refct > 0) {
2657: *ts = NULL;
2658: PetscFunctionReturn(PETSC_SUCCESS);
2659: }
2661: PetscCall(TSReset(*ts));
2662: PetscCall(TSAdjointReset(*ts));
2663: if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2664: PetscCall(TSResizeReset(*ts));
2666: /* if memory was published with SAWs then destroy it */
2667: PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2668: PetscTryTypeMethod(*ts, destroy);
2670: PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));
2672: PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2673: PetscCall(TSEventDestroy(&(*ts)->event));
2675: PetscCall(SNESDestroy(&(*ts)->snes));
2676: PetscCall(SNESDestroy(&(*ts)->snesrhssplit));
2677: PetscCall(DMDestroy(&(*ts)->dm));
2678: PetscCall(TSMonitorCancel(*ts));
2679: PetscCall(TSAdjointMonitorCancel(*ts));
2681: PetscCall(TSDestroy(&(*ts)->quadraturets));
2682: PetscCall(PetscHeaderDestroy(ts));
2683: PetscFunctionReturn(PETSC_SUCCESS);
2684: }
2686: /*@
2687: TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2688: a `TS` (timestepper) context. Valid only for nonlinear problems.
2690: Not Collective, but snes is parallel if ts is parallel
2692: Input Parameter:
2693: . ts - the `TS` context obtained from `TSCreate()`
2695: Output Parameter:
2696: . snes - the nonlinear solver context
2698: Level: beginner
2700: Notes:
2701: The user can then directly manipulate the `SNES` context to set various
2702: options, etc. Likewise, the user can then extract and manipulate the
2703: `KSP`, and `PC` contexts as well.
2705: `TSGetSNES()` does not work for integrators that do not use `SNES`; in
2706: this case `TSGetSNES()` returns `NULL` in `snes`.
2708: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2709: @*/
2710: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2711: {
2712: PetscFunctionBegin;
2714: PetscAssertPointer(snes, 2);
2715: if (!ts->snes) {
2716: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2717: PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2718: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2719: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2720: if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2721: if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2722: }
2723: *snes = ts->snes;
2724: PetscFunctionReturn(PETSC_SUCCESS);
2725: }
2727: /*@
2728: TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the `TS` timestepping context
2730: Collective
2732: Input Parameters:
2733: + ts - the `TS` context obtained from `TSCreate()`
2734: - snes - the nonlinear solver context
2736: Level: developer
2738: Note:
2739: Most users should have the `TS` created by calling `TSGetSNES()`
2741: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2742: @*/
2743: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2744: {
2745: PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);
2747: PetscFunctionBegin;
2750: PetscCall(PetscObjectReference((PetscObject)snes));
2751: PetscCall(SNESDestroy(&ts->snes));
2753: ts->snes = snes;
2755: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2756: PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2757: if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2758: PetscFunctionReturn(PETSC_SUCCESS);
2759: }
2761: /*@
2762: TSGetKSP - Returns the `KSP` (linear solver) associated with
2763: a `TS` (timestepper) context.
2765: Not Collective, but `ksp` is parallel if `ts` is parallel
2767: Input Parameter:
2768: . ts - the `TS` context obtained from `TSCreate()`
2770: Output Parameter:
2771: . ksp - the nonlinear solver context
2773: Level: beginner
2775: Notes:
2776: The user can then directly manipulate the `KSP` context to set various
2777: options, etc. Likewise, the user can then extract and manipulate the
2778: `PC` context as well.
2780: `TSGetKSP()` does not work for integrators that do not use `KSP`;
2781: in this case `TSGetKSP()` returns `NULL` in `ksp`.
2783: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2784: @*/
2785: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2786: {
2787: SNES snes;
2789: PetscFunctionBegin;
2791: PetscAssertPointer(ksp, 2);
2792: PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2793: PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
2794: PetscCall(TSGetSNES(ts, &snes));
2795: PetscCall(SNESGetKSP(snes, ksp));
2796: PetscFunctionReturn(PETSC_SUCCESS);
2797: }
2799: /* ----------- Routines to set solver parameters ---------- */
2801: /*@
2802: TSSetMaxSteps - Sets the maximum number of steps to use.
2804: Logically Collective
2806: Input Parameters:
2807: + ts - the `TS` context obtained from `TSCreate()`
2808: - maxsteps - maximum number of steps to use
2810: Options Database Key:
2811: . -ts_max_steps <maxsteps> - Sets maxsteps
2813: Level: intermediate
2815: Note:
2816: Use `PETSC_DETERMINE` to reset the maximum number of steps to the default from when the object's type was set
2818: The default maximum number of steps is 5,000
2820: Fortran Note:
2821: Use `PETSC_DETERMINE_INTEGER`
2823: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2824: @*/
2825: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2826: {
2827: PetscFunctionBegin;
2830: if (maxsteps == PETSC_DETERMINE) {
2831: ts->max_steps = ts->default_max_steps;
2832: } else {
2833: PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2834: ts->max_steps = maxsteps;
2835: }
2836: PetscFunctionReturn(PETSC_SUCCESS);
2837: }
2839: /*@
2840: TSGetMaxSteps - Gets the maximum number of steps to use.
2842: Not Collective
2844: Input Parameter:
2845: . ts - the `TS` context obtained from `TSCreate()`
2847: Output Parameter:
2848: . maxsteps - maximum number of steps to use
2850: Level: advanced
2852: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2853: @*/
2854: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2855: {
2856: PetscFunctionBegin;
2858: PetscAssertPointer(maxsteps, 2);
2859: *maxsteps = ts->max_steps;
2860: PetscFunctionReturn(PETSC_SUCCESS);
2861: }
2863: /*@
2864: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
2866: Logically Collective
2868: Input Parameters:
2869: + ts - the `TS` context obtained from `TSCreate()`
2870: - maxtime - final time to step to
2872: Options Database Key:
2873: . -ts_max_time <maxtime> - Sets maxtime
2875: Level: intermediate
2877: Notes:
2878: Use `PETSC_DETERMINE` to reset the maximum time to the default from when the object's type was set
2880: The default maximum time is 5.0
2882: Fortran Note:
2883: Use `PETSC_DETERMINE_REAL`
2885: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2886: @*/
2887: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2888: {
2889: PetscFunctionBegin;
2892: if (maxtime == PETSC_DETERMINE) {
2893: ts->max_time = ts->default_max_time;
2894: } else {
2895: ts->max_time = maxtime;
2896: }
2897: PetscFunctionReturn(PETSC_SUCCESS);
2898: }
2900: /*@
2901: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
2903: Not Collective
2905: Input Parameter:
2906: . ts - the `TS` context obtained from `TSCreate()`
2908: Output Parameter:
2909: . maxtime - final time to step to
2911: Level: advanced
2913: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2914: @*/
2915: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2916: {
2917: PetscFunctionBegin;
2919: PetscAssertPointer(maxtime, 2);
2920: *maxtime = ts->max_time;
2921: PetscFunctionReturn(PETSC_SUCCESS);
2922: }
2924: // PetscClangLinter pragma disable: -fdoc-*
2925: /*@
2926: TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.
2928: Level: deprecated
2930: @*/
2931: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2932: {
2933: PetscFunctionBegin;
2935: PetscCall(TSSetTime(ts, initial_time));
2936: PetscCall(TSSetTimeStep(ts, time_step));
2937: PetscFunctionReturn(PETSC_SUCCESS);
2938: }
2940: // PetscClangLinter pragma disable: -fdoc-*
2941: /*@
2942: TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.
2944: Level: deprecated
2946: @*/
2947: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2948: {
2949: PetscFunctionBegin;
2951: if (maxsteps) {
2952: PetscAssertPointer(maxsteps, 2);
2953: *maxsteps = ts->max_steps;
2954: }
2955: if (maxtime) {
2956: PetscAssertPointer(maxtime, 3);
2957: *maxtime = ts->max_time;
2958: }
2959: PetscFunctionReturn(PETSC_SUCCESS);
2960: }
2962: // PetscClangLinter pragma disable: -fdoc-*
2963: /*@
2964: TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.
2966: Level: deprecated
2968: @*/
2969: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2970: {
2971: PetscFunctionBegin;
2972: if (maxsteps != PETSC_CURRENT) PetscCall(TSSetMaxSteps(ts, maxsteps));
2973: if (maxtime != (PetscReal)PETSC_CURRENT) PetscCall(TSSetMaxTime(ts, maxtime));
2974: PetscFunctionReturn(PETSC_SUCCESS);
2975: }
2977: // PetscClangLinter pragma disable: -fdoc-*
2978: /*@
2979: TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.
2981: Level: deprecated
2983: @*/
2984: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
2985: {
2986: return TSGetStepNumber(ts, steps);
2987: }
2989: // PetscClangLinter pragma disable: -fdoc-*
2990: /*@
2991: TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.
2993: Level: deprecated
2995: @*/
2996: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
2997: {
2998: return TSGetStepNumber(ts, steps);
2999: }
3001: /*@
3002: TSSetSolution - Sets the initial solution vector
3003: for use by the `TS` routines.
3005: Logically Collective
3007: Input Parameters:
3008: + ts - the `TS` context obtained from `TSCreate()`
3009: - u - the solution vector
3011: Level: beginner
3013: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
3014: @*/
3015: PetscErrorCode TSSetSolution(TS ts, Vec u)
3016: {
3017: DM dm;
3019: PetscFunctionBegin;
3022: PetscCall(PetscObjectReference((PetscObject)u));
3023: PetscCall(VecDestroy(&ts->vec_sol));
3024: ts->vec_sol = u;
3026: PetscCall(TSGetDM(ts, &dm));
3027: PetscCall(DMShellSetGlobalVector(dm, u));
3028: PetscFunctionReturn(PETSC_SUCCESS);
3029: }
3031: /*@C
3032: TSSetPreStep - Sets the general-purpose function
3033: called once at the beginning of each time step.
3035: Logically Collective
3037: Input Parameters:
3038: + ts - The `TS` context obtained from `TSCreate()`
3039: - func - The function
3041: Calling sequence of `func`:
3042: . ts - the `TS` context
3044: Level: intermediate
3046: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3047: @*/
3048: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
3049: {
3050: PetscFunctionBegin;
3052: ts->prestep = func;
3053: PetscFunctionReturn(PETSC_SUCCESS);
3054: }
3056: /*@
3057: TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`
3059: Collective
3061: Input Parameter:
3062: . ts - The `TS` context obtained from `TSCreate()`
3064: Level: developer
3066: Note:
3067: `TSPreStep()` is typically used within time stepping implementations,
3068: so most users would not generally call this routine themselves.
3070: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3071: @*/
3072: PetscErrorCode TSPreStep(TS ts)
3073: {
3074: PetscFunctionBegin;
3076: if (ts->prestep) {
3077: Vec U;
3078: PetscObjectId idprev;
3079: PetscBool sameObject;
3080: PetscObjectState sprev, spost;
3082: PetscCall(TSGetSolution(ts, &U));
3083: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3084: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3085: PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3086: PetscCall(TSGetSolution(ts, &U));
3087: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3088: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3089: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3090: }
3091: PetscFunctionReturn(PETSC_SUCCESS);
3092: }
3094: /*@C
3095: TSSetPreStage - Sets the general-purpose function
3096: called once at the beginning of each stage.
3098: Logically Collective
3100: Input Parameters:
3101: + ts - The `TS` context obtained from `TSCreate()`
3102: - func - The function
3104: Calling sequence of `func`:
3105: + ts - the `TS` context
3106: - stagetime - the stage time
3108: Level: intermediate
3110: Note:
3111: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3112: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3113: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3115: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3116: @*/
3117: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime))
3118: {
3119: PetscFunctionBegin;
3121: ts->prestage = func;
3122: PetscFunctionReturn(PETSC_SUCCESS);
3123: }
3125: /*@C
3126: TSSetPostStage - Sets the general-purpose function
3127: called once at the end of each stage.
3129: Logically Collective
3131: Input Parameters:
3132: + ts - The `TS` context obtained from `TSCreate()`
3133: - func - The function
3135: Calling sequence of `func`:
3136: + ts - the `TS` context
3137: . stagetime - the stage time
3138: . stageindex - the stage index
3139: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3141: Level: intermediate
3143: Note:
3144: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3145: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3146: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3148: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3149: @*/
3150: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
3151: {
3152: PetscFunctionBegin;
3154: ts->poststage = func;
3155: PetscFunctionReturn(PETSC_SUCCESS);
3156: }
3158: /*@C
3159: TSSetPostEvaluate - Sets the general-purpose function
3160: called at the end of each step evaluation.
3162: Logically Collective
3164: Input Parameters:
3165: + ts - The `TS` context obtained from `TSCreate()`
3166: - func - The function
3168: Calling sequence of `func`:
3169: . ts - the `TS` context
3171: Level: intermediate
3173: Note:
3174: The function set by `TSSetPostEvaluate()` is called after the solution is evaluated, or after the step rollback.
3175: Inside the `func` callback, the solution vector can be obtained with `TSGetSolution()`, and modified, if need be.
3176: The time step can be obtained with `TSGetTimeStep()`, and the time at the start of the step - via `TSGetTime()`.
3177: The potential changes to the solution vector introduced by event handling (`postevent()`) are not relevant for `TSSetPostEvaluate()`,
3178: but are relevant for `TSSetPostStep()`, according to the function call scheme in `TSSolve()`, as shown below
3179: .vb
3180: ...
3181: Step()
3182: PostEvaluate()
3183: EventHandling()
3184: step_rollback ? PostEvaluate() : PostStep()
3185: ...
3186: .ve
3187: where EventHandling() may result in one of the following three outcomes
3188: .vb
3189: (1) | successful step | solution intact
3190: (2) | successful step | solution modified by `postevent()`
3191: (3) | step_rollback | solution rolled back
3192: .ve
3194: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3195: @*/
3196: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
3197: {
3198: PetscFunctionBegin;
3200: ts->postevaluate = func;
3201: PetscFunctionReturn(PETSC_SUCCESS);
3202: }
3204: /*@
3205: TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`
3207: Collective
3209: Input Parameters:
3210: + ts - The `TS` context obtained from `TSCreate()`
3211: - stagetime - The absolute time of the current stage
3213: Level: developer
3215: Note:
3216: `TSPreStage()` is typically used within time stepping implementations,
3217: most users would not generally call this routine themselves.
3219: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3220: @*/
3221: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3222: {
3223: PetscFunctionBegin;
3225: if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3226: PetscFunctionReturn(PETSC_SUCCESS);
3227: }
3229: /*@
3230: TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`
3232: Collective
3234: Input Parameters:
3235: + ts - The `TS` context obtained from `TSCreate()`
3236: . stagetime - The absolute time of the current stage
3237: . stageindex - Stage number
3238: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3240: Level: developer
3242: Note:
3243: `TSPostStage()` is typically used within time stepping implementations,
3244: most users would not generally call this routine themselves.
3246: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3247: @*/
3248: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec Y[])
3249: {
3250: PetscFunctionBegin;
3252: if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3253: PetscFunctionReturn(PETSC_SUCCESS);
3254: }
3256: /*@
3257: TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`
3259: Collective
3261: Input Parameter:
3262: . ts - The `TS` context obtained from `TSCreate()`
3264: Level: developer
3266: Note:
3267: `TSPostEvaluate()` is typically used within time stepping implementations,
3268: most users would not generally call this routine themselves.
3270: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3271: @*/
3272: PetscErrorCode TSPostEvaluate(TS ts)
3273: {
3274: PetscFunctionBegin;
3276: if (ts->postevaluate) {
3277: Vec U;
3278: PetscObjectState sprev, spost;
3280: PetscCall(TSGetSolution(ts, &U));
3281: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3282: PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3283: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3284: if (sprev != spost) PetscCall(TSRestartStep(ts));
3285: }
3286: PetscFunctionReturn(PETSC_SUCCESS);
3287: }
3289: /*@C
3290: TSSetPostStep - Sets the general-purpose function
3291: called once at the end of each successful time step.
3293: Logically Collective
3295: Input Parameters:
3296: + ts - The `TS` context obtained from `TSCreate()`
3297: - func - The function
3299: Calling sequence of `func`:
3300: . ts - the `TS` context
3302: Level: intermediate
3304: Note:
3305: The function set by `TSSetPostStep()` is called after each successful step. If the event handler locates an event at the
3306: given step, and `postevent()` modifies the solution vector, the solution vector obtained by `TSGetSolution()` inside `func` will
3307: contain the changes. To get the solution without these changes, use `TSSetPostEvaluate()` to set the appropriate callback.
3308: The scheme of the relevant function calls in `TSSolve()` is shown below
3309: .vb
3310: ...
3311: Step()
3312: PostEvaluate()
3313: EventHandling()
3314: step_rollback ? PostEvaluate() : PostStep()
3315: ...
3316: .ve
3317: where EventHandling() may result in one of the following three outcomes
3318: .vb
3319: (1) | successful step | solution intact
3320: (2) | successful step | solution modified by `postevent()`
3321: (3) | step_rollback | solution rolled back
3322: .ve
3324: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3325: @*/
3326: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
3327: {
3328: PetscFunctionBegin;
3330: ts->poststep = func;
3331: PetscFunctionReturn(PETSC_SUCCESS);
3332: }
3334: /*@
3335: TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`
3337: Collective
3339: Input Parameter:
3340: . ts - The `TS` context obtained from `TSCreate()`
3342: Note:
3343: `TSPostStep()` is typically used within time stepping implementations,
3344: so most users would not generally call this routine themselves.
3346: Level: developer
3348: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPostStep()`
3349: @*/
3350: PetscErrorCode TSPostStep(TS ts)
3351: {
3352: PetscFunctionBegin;
3354: if (ts->poststep) {
3355: Vec U;
3356: PetscObjectId idprev;
3357: PetscBool sameObject;
3358: PetscObjectState sprev, spost;
3360: PetscCall(TSGetSolution(ts, &U));
3361: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3362: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3363: PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3364: PetscCall(TSGetSolution(ts, &U));
3365: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3366: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3367: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3368: }
3369: PetscFunctionReturn(PETSC_SUCCESS);
3370: }
3372: /*@
3373: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3375: Collective
3377: Input Parameters:
3378: + ts - time stepping context
3379: - t - time to interpolate to
3381: Output Parameter:
3382: . U - state at given time
3384: Level: intermediate
3386: Developer Notes:
3387: `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3389: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3390: @*/
3391: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3392: {
3393: PetscFunctionBegin;
3396: 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);
3397: PetscUseTypeMethod(ts, interpolate, t, U);
3398: PetscFunctionReturn(PETSC_SUCCESS);
3399: }
3401: /*@
3402: TSStep - Steps one time step
3404: Collective
3406: Input Parameter:
3407: . ts - the `TS` context obtained from `TSCreate()`
3409: Level: developer
3411: Notes:
3412: The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.
3414: The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
3415: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3417: This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
3418: time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.
3420: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3421: @*/
3422: PetscErrorCode TSStep(TS ts)
3423: {
3424: static PetscBool cite = PETSC_FALSE;
3425: PetscReal ptime;
3427: PetscFunctionBegin;
3429: PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3430: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3431: " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3432: " journal = {arXiv e-preprints},\n"
3433: " eprint = {1806.01437},\n"
3434: " archivePrefix = {arXiv},\n"
3435: " year = {2018}\n}\n",
3436: &cite));
3437: PetscCall(TSSetUp(ts));
3438: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3439: if (ts->eval_times)
3440: ts->eval_times->worktol = 0; /* In each step of TSSolve() 'eval_times->worktol' will be meaningfully defined (later) only once:
3441: in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */
3443: PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_INT_MAX, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
3444: 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()");
3445: 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");
3447: if (!ts->vec_sol0) PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sol0));
3448: PetscCall(VecCopy(ts->vec_sol, ts->vec_sol0));
3449: ts->time_step0 = ts->time_step;
3451: if (!ts->steps) ts->ptime_prev = ts->ptime;
3452: ptime = ts->ptime;
3454: ts->ptime_prev_rollback = ts->ptime_prev;
3455: ts->reason = TS_CONVERGED_ITERATING;
3457: PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3458: PetscUseTypeMethod(ts, step);
3459: PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));
3461: if (ts->reason >= 0) {
3462: ts->ptime_prev = ptime;
3463: ts->steps++;
3464: ts->steprollback = PETSC_FALSE;
3465: ts->steprestart = PETSC_FALSE;
3466: ts->stepresize = PETSC_FALSE;
3467: }
3469: if (ts->reason < 0 && ts->errorifstepfailed) {
3470: PetscCall(TSMonitorCancel(ts));
3471: 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]);
3472: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3473: }
3474: PetscFunctionReturn(PETSC_SUCCESS);
3475: }
3477: /*@
3478: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3479: at the end of a time step with a given order of accuracy.
3481: Collective
3483: Input Parameters:
3484: + ts - time stepping context
3485: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3487: Input/Output Parameter:
3488: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3489: on output, the actual order of the error evaluation
3491: Output Parameter:
3492: . wlte - the weighted local truncation error norm
3494: Level: advanced
3496: Note:
3497: If the timestepper cannot evaluate the error in a particular step
3498: (eg. in the first step or restart steps after event handling),
3499: this routine returns wlte=-1.0 .
3501: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3502: @*/
3503: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3504: {
3505: PetscFunctionBegin;
3509: if (order) PetscAssertPointer(order, 3);
3511: PetscAssertPointer(wlte, 4);
3512: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3513: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3514: PetscFunctionReturn(PETSC_SUCCESS);
3515: }
3517: /*@
3518: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3520: Collective
3522: Input Parameters:
3523: + ts - time stepping context
3524: . order - desired order of accuracy
3525: - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)
3527: Output Parameter:
3528: . U - state at the end of the current step
3530: Level: advanced
3532: Notes:
3533: This function cannot be called until all stages have been evaluated.
3535: 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.
3537: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3538: @*/
3539: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3540: {
3541: PetscFunctionBegin;
3545: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3546: PetscFunctionReturn(PETSC_SUCCESS);
3547: }
3549: /*@C
3550: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3552: Not collective
3554: Input Parameter:
3555: . ts - time stepping context
3557: Output Parameter:
3558: . initCondition - The function which computes an initial condition
3560: Calling sequence of `initCondition`:
3561: + ts - The timestepping context
3562: - u - The input vector in which the initial condition is stored
3564: Level: advanced
3566: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3567: @*/
3568: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3569: {
3570: PetscFunctionBegin;
3572: PetscAssertPointer(initCondition, 2);
3573: *initCondition = ts->ops->initcondition;
3574: PetscFunctionReturn(PETSC_SUCCESS);
3575: }
3577: /*@C
3578: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3580: Logically collective
3582: Input Parameters:
3583: + ts - time stepping context
3584: - initCondition - The function which computes an initial condition
3586: Calling sequence of `initCondition`:
3587: + ts - The timestepping context
3588: - e - The input vector in which the initial condition is to be stored
3590: Level: advanced
3592: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3593: @*/
3594: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3595: {
3596: PetscFunctionBegin;
3599: ts->ops->initcondition = initCondition;
3600: PetscFunctionReturn(PETSC_SUCCESS);
3601: }
3603: /*@
3604: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3606: Collective
3608: Input Parameters:
3609: + ts - time stepping context
3610: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3612: Level: advanced
3614: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3615: @*/
3616: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3617: {
3618: PetscFunctionBegin;
3621: PetscTryTypeMethod(ts, initcondition, u);
3622: PetscFunctionReturn(PETSC_SUCCESS);
3623: }
3625: /*@C
3626: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3628: Not collective
3630: Input Parameter:
3631: . ts - time stepping context
3633: Output Parameter:
3634: . exactError - The function which computes the solution error
3636: Calling sequence of `exactError`:
3637: + ts - The timestepping context
3638: . u - The approximate solution vector
3639: - e - The vector in which the error is stored
3641: Level: advanced
3643: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3644: @*/
3645: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3646: {
3647: PetscFunctionBegin;
3649: PetscAssertPointer(exactError, 2);
3650: *exactError = ts->ops->exacterror;
3651: PetscFunctionReturn(PETSC_SUCCESS);
3652: }
3654: /*@C
3655: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3657: Logically collective
3659: Input Parameters:
3660: + ts - time stepping context
3661: - exactError - The function which computes the solution error
3663: Calling sequence of `exactError`:
3664: + ts - The timestepping context
3665: . u - The approximate solution vector
3666: - e - The vector in which the error is stored
3668: Level: advanced
3670: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3671: @*/
3672: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3673: {
3674: PetscFunctionBegin;
3677: ts->ops->exacterror = exactError;
3678: PetscFunctionReturn(PETSC_SUCCESS);
3679: }
3681: /*@
3682: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3684: Collective
3686: Input Parameters:
3687: + ts - time stepping context
3688: . u - The approximate solution
3689: - e - The `Vec` used to store the error
3691: Level: advanced
3693: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3694: @*/
3695: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3696: {
3697: PetscFunctionBegin;
3701: PetscTryTypeMethod(ts, exacterror, u, e);
3702: PetscFunctionReturn(PETSC_SUCCESS);
3703: }
3705: /*@C
3706: TSSetResize - Sets the resize callbacks.
3708: Logically Collective
3710: Input Parameters:
3711: + ts - The `TS` context obtained from `TSCreate()`
3712: . rollback - Whether a resize will restart the step
3713: . setup - The setup function
3714: . transfer - The transfer function
3715: - ctx - [optional] The user-defined context
3717: Calling sequence of `setup`:
3718: + ts - the `TS` context
3719: . step - the current step
3720: . time - the current time
3721: . state - the current vector of state
3722: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3723: - ctx - user defined context
3725: Calling sequence of `transfer`:
3726: + ts - the `TS` context
3727: . nv - the number of vectors to be transferred
3728: . vecsin - array of vectors to be transferred
3729: . vecsout - array of transferred vectors
3730: - ctx - user defined context
3732: Notes:
3733: The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
3734: depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
3735: and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
3736: that the size of the discrete problem has changed.
3737: In both cases, the solver will collect the needed vectors that will be
3738: transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
3739: current solution vector, and other vectors needed by the specific solver used.
3740: For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3741: Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3742: will be automatically reset if the sizes are changed and they must be specified again by the user
3743: inside the `transfer` function.
3744: The input and output arrays passed to `transfer` are allocated by PETSc.
3745: Vectors in `vecsout` must be created by the user.
3746: Ownership of vectors in `vecsout` is transferred to PETSc.
3748: Level: advanced
3750: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3751: @*/
3752: PetscErrorCode TSSetResize(TS ts, PetscBool rollback, PetscErrorCode (*setup)(TS ts, PetscInt step, PetscReal time, Vec state, PetscBool *resize, void *ctx), PetscErrorCode (*transfer)(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx), void *ctx)
3753: {
3754: PetscFunctionBegin;
3756: ts->resizerollback = rollback;
3757: ts->resizesetup = setup;
3758: ts->resizetransfer = transfer;
3759: ts->resizectx = ctx;
3760: PetscFunctionReturn(PETSC_SUCCESS);
3761: }
3763: /*
3764: TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.
3766: Collective
3768: Input Parameters:
3769: + ts - The `TS` context obtained from `TSCreate()`
3770: - 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.
3772: Level: developer
3774: Note:
3775: `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3776: used within time stepping implementations,
3777: so most users would not generally call this routine themselves.
3779: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3780: @*/
3781: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3782: {
3783: PetscFunctionBegin;
3785: PetscTryTypeMethod(ts, resizeregister, flg);
3786: /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3787: PetscFunctionReturn(PETSC_SUCCESS);
3788: }
3790: static PetscErrorCode TSResizeReset(TS ts)
3791: {
3792: PetscFunctionBegin;
3794: PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3795: PetscFunctionReturn(PETSC_SUCCESS);
3796: }
3798: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3799: {
3800: PetscFunctionBegin;
3803: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3804: if (ts->resizetransfer) {
3805: PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3806: PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3807: }
3808: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3809: PetscFunctionReturn(PETSC_SUCCESS);
3810: }
3812: /*@C
3813: TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.
3815: Collective
3817: Input Parameters:
3818: + ts - The `TS` context obtained from `TSCreate()`
3819: . name - A string identifying the vector
3820: - vec - The vector
3822: Level: developer
3824: Note:
3825: `TSResizeRegisterVec()` is typically used within time stepping implementations,
3826: so most users would not generally call this routine themselves.
3828: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3829: @*/
3830: PetscErrorCode TSResizeRegisterVec(TS ts, const char name[], Vec vec)
3831: {
3832: PetscFunctionBegin;
3834: PetscAssertPointer(name, 2);
3836: PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3837: PetscFunctionReturn(PETSC_SUCCESS);
3838: }
3840: /*@C
3841: TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.
3843: Collective
3845: Input Parameters:
3846: + ts - The `TS` context obtained from `TSCreate()`
3847: . name - A string identifying the vector
3848: - vec - The vector
3850: Level: developer
3852: Note:
3853: `TSResizeRetrieveVec()` is typically used within time stepping implementations,
3854: so most users would not generally call this routine themselves.
3856: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
3857: @*/
3858: PetscErrorCode TSResizeRetrieveVec(TS ts, const char name[], Vec *vec)
3859: {
3860: PetscFunctionBegin;
3862: PetscAssertPointer(name, 2);
3863: PetscAssertPointer(vec, 3);
3864: PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
3865: PetscFunctionReturn(PETSC_SUCCESS);
3866: }
3868: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
3869: {
3870: PetscInt cnt;
3871: PetscObjectList tmp;
3872: Vec *vecsin = NULL;
3873: const char **namesin = NULL;
3875: PetscFunctionBegin;
3876: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
3877: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
3878: if (names) PetscCall(PetscMalloc1(cnt, &namesin));
3879: if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
3880: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
3881: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
3882: if (vecs) vecsin[cnt] = (Vec)tmp->obj;
3883: if (names) namesin[cnt] = tmp->name;
3884: cnt++;
3885: }
3886: }
3887: if (nv) *nv = cnt;
3888: if (names) *names = namesin;
3889: if (vecs) *vecs = vecsin;
3890: PetscFunctionReturn(PETSC_SUCCESS);
3891: }
3893: /*@
3894: TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`
3896: Collective
3898: Input Parameter:
3899: . ts - The `TS` context obtained from `TSCreate()`
3901: Level: developer
3903: Note:
3904: `TSResize()` is typically used within time stepping implementations,
3905: so most users would not generally call this routine themselves.
3907: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3908: @*/
3909: PetscErrorCode TSResize(TS ts)
3910: {
3911: PetscInt nv = 0;
3912: const char **names = NULL;
3913: Vec *vecsin = NULL;
3914: const char *solname = "ts:vec_sol";
3916: PetscFunctionBegin;
3918: if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
3919: if (ts->resizesetup) {
3920: PetscCall(VecLockReadPush(ts->vec_sol));
3921: PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &ts->stepresize, ts->resizectx));
3922: PetscCall(VecLockReadPop(ts->vec_sol));
3923: if (ts->stepresize) {
3924: if (ts->resizerollback) {
3925: PetscCall(TSRollBack(ts));
3926: ts->time_step = ts->time_step0;
3927: }
3928: PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
3929: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
3930: }
3931: }
3933: PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
3934: if (nv) {
3935: Vec *vecsout, vecsol;
3937: /* Reset internal objects */
3938: PetscCall(TSReset(ts));
3940: /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
3941: PetscCall(PetscCalloc1(nv, &vecsout));
3942: PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
3943: for (PetscInt i = 0; i < nv; i++) {
3944: const char *name;
3945: char *oname;
3947: PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
3948: PetscCall(PetscStrallocpy(name, &oname));
3949: PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
3950: if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
3951: PetscCall(PetscFree(oname));
3952: PetscCall(VecDestroy(&vecsout[i]));
3953: }
3954: PetscCall(PetscFree(vecsout));
3955: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */
3957: PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
3958: if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
3959: PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
3960: }
3962: PetscCall(PetscFree(names));
3963: PetscCall(PetscFree(vecsin));
3964: PetscCall(TSResizeReset(ts));
3965: PetscFunctionReturn(PETSC_SUCCESS);
3966: }
3968: /*@
3969: TSSolve - Steps the requested number of timesteps.
3971: Collective
3973: Input Parameters:
3974: + ts - the `TS` context obtained from `TSCreate()`
3975: - u - the solution vector (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
3976: otherwise it must contain the initial conditions and will contain the solution at the final requested time
3978: Level: beginner
3980: Notes:
3981: The final time returned by this function may be different from the time of the internally
3982: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
3983: stepped over the final time.
3985: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3986: @*/
3987: PetscErrorCode TSSolve(TS ts, Vec u)
3988: {
3989: Vec solution;
3991: PetscFunctionBegin;
3995: PetscCall(TSSetExactFinalTimeDefault(ts));
3996: 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 */
3997: if (!ts->vec_sol || u == ts->vec_sol) {
3998: PetscCall(VecDuplicate(u, &solution));
3999: PetscCall(TSSetSolution(ts, solution));
4000: PetscCall(VecDestroy(&solution)); /* grant ownership */
4001: }
4002: PetscCall(VecCopy(u, ts->vec_sol));
4003: PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
4004: } else if (u) PetscCall(TSSetSolution(ts, u));
4005: PetscCall(TSSetUp(ts));
4006: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
4008: PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_INT_MAX, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
4009: 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()");
4010: 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");
4011: 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");
4013: if (ts->eval_times) {
4014: for (PetscInt i = 0; i < ts->eval_times->num_time_points; i++) {
4015: PetscBool is_close = PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[i], ts->eval_times->reltol * ts->time_step + ts->eval_times->abstol, 0);
4016: if (ts->ptime <= ts->eval_times->time_points[i] || is_close) {
4017: ts->eval_times->time_point_idx = i;
4018: if (is_close) { /* starting point in evaluation times */
4019: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_ctr]));
4020: ts->eval_times->sol_times[ts->eval_times->sol_ctr] = ts->ptime;
4021: ts->eval_times->sol_ctr++;
4022: ts->eval_times->time_point_idx++;
4023: }
4024: break;
4025: }
4026: }
4027: }
4029: if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));
4031: /* reset number of steps only when the step is not restarted. ARKIMEX
4032: restarts the step after an event. Resetting these counters in such case causes
4033: TSTrajectory to incorrectly save the output files
4034: */
4035: /* reset time step and iteration counters */
4036: if (!ts->steps) {
4037: ts->ksp_its = 0;
4038: ts->snes_its = 0;
4039: ts->num_snes_failures = 0;
4040: ts->reject = 0;
4041: ts->steprestart = PETSC_TRUE;
4042: ts->steprollback = PETSC_FALSE;
4043: ts->stepresize = PETSC_FALSE;
4044: ts->rhsjacobian.time = PETSC_MIN_REAL;
4045: }
4047: /* make sure initial time step does not overshoot final time or the next point in evaluation times */
4048: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
4049: PetscReal maxdt;
4050: PetscReal dt = ts->time_step;
4052: if (ts->eval_times) maxdt = ts->eval_times->time_points[ts->eval_times->time_point_idx] - ts->ptime;
4053: else maxdt = ts->max_time - ts->ptime;
4054: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
4055: }
4056: ts->reason = TS_CONVERGED_ITERATING;
4058: {
4059: PetscViewer viewer;
4060: PetscViewerFormat format;
4061: PetscBool flg;
4062: static PetscBool incall = PETSC_FALSE;
4064: if (!incall) {
4065: /* Estimate the convergence rate of the time discretization */
4066: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
4067: if (flg) {
4068: PetscConvEst conv;
4069: DM dm;
4070: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4071: PetscInt Nf;
4072: PetscBool checkTemporal = PETSC_TRUE;
4074: incall = PETSC_TRUE;
4075: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
4076: PetscCall(TSGetDM(ts, &dm));
4077: PetscCall(DMGetNumFields(dm, &Nf));
4078: PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
4079: PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
4080: PetscCall(PetscConvEstUseTS(conv, checkTemporal));
4081: PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
4082: PetscCall(PetscConvEstSetFromOptions(conv));
4083: PetscCall(PetscConvEstSetUp(conv));
4084: PetscCall(PetscConvEstGetConvRate(conv, alpha));
4085: PetscCall(PetscViewerPushFormat(viewer, format));
4086: PetscCall(PetscConvEstRateView(conv, alpha, viewer));
4087: PetscCall(PetscViewerPopFormat(viewer));
4088: PetscCall(PetscViewerDestroy(&viewer));
4089: PetscCall(PetscConvEstDestroy(&conv));
4090: PetscCall(PetscFree(alpha));
4091: incall = PETSC_FALSE;
4092: }
4093: }
4094: }
4096: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));
4098: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4099: PetscUseTypeMethod(ts, solve);
4100: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4101: ts->solvetime = ts->ptime;
4102: solution = ts->vec_sol;
4103: } else { /* Step the requested number of timesteps. */
4104: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4105: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4107: if (!ts->steps) {
4108: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4109: PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4110: }
4112: while (!ts->reason) {
4113: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4114: if (!ts->steprollback || (ts->stepresize && ts->resizerollback)) PetscCall(TSPreStep(ts));
4115: PetscCall(TSStep(ts));
4116: if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4117: if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4118: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4119: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4120: PetscCall(TSForwardCostIntegral(ts));
4121: if (ts->reason >= 0) ts->steps++;
4122: }
4123: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4124: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4125: PetscCall(TSForwardStep(ts));
4126: if (ts->reason >= 0) ts->steps++;
4127: }
4128: PetscCall(TSPostEvaluate(ts));
4129: 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. */
4130: if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4131: if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
4132: /* check convergence */
4133: if (!ts->reason) {
4134: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4135: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4136: }
4137: if (!ts->steprollback) {
4138: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4139: PetscCall(TSPostStep(ts));
4140: if (!ts->resizerollback) PetscCall(TSResize(ts));
4142: if (ts->eval_times && ts->eval_times->time_point_idx < ts->eval_times->num_time_points && ts->reason >= 0) {
4143: PetscCheck(ts->eval_times->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(eval_times->worktol > 0) in TSSolve()");
4144: if (PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[ts->eval_times->time_point_idx], ts->eval_times->worktol, 0)) {
4145: ts->eval_times->sol_times[ts->eval_times->sol_ctr] = ts->ptime;
4146: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_ctr]));
4147: ts->eval_times->sol_ctr++;
4148: ts->eval_times->time_point_idx++;
4149: }
4150: }
4151: }
4152: }
4153: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4155: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4156: if (!u) u = ts->vec_sol;
4157: PetscCall(TSInterpolate(ts, ts->max_time, u));
4158: ts->solvetime = ts->max_time;
4159: solution = u;
4160: PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4161: } else {
4162: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4163: ts->solvetime = ts->ptime;
4164: solution = ts->vec_sol;
4165: }
4166: }
4168: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4169: PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4170: PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4171: if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4172: PetscFunctionReturn(PETSC_SUCCESS);
4173: }
4175: /*@
4176: TSGetTime - Gets the time of the most recently completed step.
4178: Not Collective
4180: Input Parameter:
4181: . ts - the `TS` context obtained from `TSCreate()`
4183: Output Parameter:
4184: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
4186: Level: beginner
4188: Note:
4189: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
4190: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
4192: .seealso: [](ch_ts), `TS`, `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4193: @*/
4194: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4195: {
4196: PetscFunctionBegin;
4198: PetscAssertPointer(t, 2);
4199: *t = ts->ptime;
4200: PetscFunctionReturn(PETSC_SUCCESS);
4201: }
4203: /*@
4204: TSGetPrevTime - Gets the starting time of the previously completed step.
4206: Not Collective
4208: Input Parameter:
4209: . ts - the `TS` context obtained from `TSCreate()`
4211: Output Parameter:
4212: . t - the previous time
4214: Level: beginner
4216: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4217: @*/
4218: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4219: {
4220: PetscFunctionBegin;
4222: PetscAssertPointer(t, 2);
4223: *t = ts->ptime_prev;
4224: PetscFunctionReturn(PETSC_SUCCESS);
4225: }
4227: /*@
4228: TSSetTime - Allows one to reset the time.
4230: Logically Collective
4232: Input Parameters:
4233: + ts - the `TS` context obtained from `TSCreate()`
4234: - t - the time
4236: Level: intermediate
4238: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4239: @*/
4240: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4241: {
4242: PetscFunctionBegin;
4245: ts->ptime = t;
4246: PetscFunctionReturn(PETSC_SUCCESS);
4247: }
4249: /*@
4250: TSSetOptionsPrefix - Sets the prefix used for searching for all
4251: TS options in the database.
4253: Logically Collective
4255: Input Parameters:
4256: + ts - The `TS` context
4257: - prefix - The prefix to prepend to all option names
4259: Level: advanced
4261: Note:
4262: A hyphen (-) must NOT be given at the beginning of the prefix name.
4263: The first character of all runtime options is AUTOMATICALLY the
4264: hyphen.
4266: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4267: @*/
4268: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4269: {
4270: SNES snes;
4272: PetscFunctionBegin;
4274: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4275: PetscCall(TSGetSNES(ts, &snes));
4276: PetscCall(SNESSetOptionsPrefix(snes, prefix));
4277: PetscFunctionReturn(PETSC_SUCCESS);
4278: }
4280: /*@
4281: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4282: TS options in the database.
4284: Logically Collective
4286: Input Parameters:
4287: + ts - The `TS` context
4288: - prefix - The prefix to prepend to all option names
4290: Level: advanced
4292: Note:
4293: A hyphen (-) must NOT be given at the beginning of the prefix name.
4294: The first character of all runtime options is AUTOMATICALLY the
4295: hyphen.
4297: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4298: @*/
4299: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4300: {
4301: SNES snes;
4303: PetscFunctionBegin;
4305: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4306: PetscCall(TSGetSNES(ts, &snes));
4307: PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4308: PetscFunctionReturn(PETSC_SUCCESS);
4309: }
4311: /*@
4312: TSGetOptionsPrefix - Sets the prefix used for searching for all
4313: `TS` options in the database.
4315: Not Collective
4317: Input Parameter:
4318: . ts - The `TS` context
4320: Output Parameter:
4321: . prefix - A pointer to the prefix string used
4323: Level: intermediate
4325: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4326: @*/
4327: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4328: {
4329: PetscFunctionBegin;
4331: PetscAssertPointer(prefix, 2);
4332: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4333: PetscFunctionReturn(PETSC_SUCCESS);
4334: }
4336: /*@C
4337: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4339: Not Collective, but parallel objects are returned if ts is parallel
4341: Input Parameter:
4342: . ts - The `TS` context obtained from `TSCreate()`
4344: Output Parameters:
4345: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`)
4346: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`)
4347: . func - Function to compute the Jacobian of the RHS (or `NULL`)
4348: - ctx - User-defined context for Jacobian evaluation routine (or `NULL`)
4350: Level: intermediate
4352: Note:
4353: You can pass in `NULL` for any return argument you do not need.
4355: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4357: @*/
4358: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, void **ctx)
4359: {
4360: DM dm;
4362: PetscFunctionBegin;
4363: if (Amat || Pmat) {
4364: SNES snes;
4365: PetscCall(TSGetSNES(ts, &snes));
4366: PetscCall(SNESSetUpMatrices(snes));
4367: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4368: }
4369: PetscCall(TSGetDM(ts, &dm));
4370: PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4371: PetscFunctionReturn(PETSC_SUCCESS);
4372: }
4374: /*@C
4375: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4377: Not Collective, but parallel objects are returned if ts is parallel
4379: Input Parameter:
4380: . ts - The `TS` context obtained from `TSCreate()`
4382: Output Parameters:
4383: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4384: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4385: . f - The function to compute the matrices
4386: - ctx - User-defined context for Jacobian evaluation routine
4388: Level: advanced
4390: Note:
4391: You can pass in `NULL` for any return argument you do not need.
4393: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4394: @*/
4395: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, void **ctx)
4396: {
4397: DM dm;
4399: PetscFunctionBegin;
4400: if (Amat || Pmat) {
4401: SNES snes;
4402: PetscCall(TSGetSNES(ts, &snes));
4403: PetscCall(SNESSetUpMatrices(snes));
4404: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4405: }
4406: PetscCall(TSGetDM(ts, &dm));
4407: PetscCall(DMTSGetIJacobian(dm, f, ctx));
4408: PetscFunctionReturn(PETSC_SUCCESS);
4409: }
4411: #include <petsc/private/dmimpl.h>
4412: /*@
4413: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4415: Logically Collective
4417: Input Parameters:
4418: + ts - the `TS` integrator object
4419: - dm - the dm, cannot be `NULL`
4421: Level: intermediate
4423: Notes:
4424: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4425: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4426: different problems using the same function space.
4428: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4429: @*/
4430: PetscErrorCode TSSetDM(TS ts, DM dm)
4431: {
4432: SNES snes;
4433: DMTS tsdm;
4435: PetscFunctionBegin;
4438: PetscCall(PetscObjectReference((PetscObject)dm));
4439: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4440: if (ts->dm->dmts && !dm->dmts) {
4441: PetscCall(DMCopyDMTS(ts->dm, dm));
4442: PetscCall(DMGetDMTS(ts->dm, &tsdm));
4443: /* Grant write privileges to the replacement DM */
4444: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4445: }
4446: PetscCall(DMDestroy(&ts->dm));
4447: }
4448: ts->dm = dm;
4450: PetscCall(TSGetSNES(ts, &snes));
4451: PetscCall(SNESSetDM(snes, dm));
4452: PetscFunctionReturn(PETSC_SUCCESS);
4453: }
4455: /*@
4456: TSGetDM - Gets the `DM` that may be used by some preconditioners
4458: Not Collective
4460: Input Parameter:
4461: . ts - the `TS`
4463: Output Parameter:
4464: . dm - the `DM`
4466: Level: intermediate
4468: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4469: @*/
4470: PetscErrorCode TSGetDM(TS ts, DM *dm)
4471: {
4472: PetscFunctionBegin;
4474: if (!ts->dm) {
4475: PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4476: if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4477: }
4478: *dm = ts->dm;
4479: PetscFunctionReturn(PETSC_SUCCESS);
4480: }
4482: /*@
4483: SNESTSFormFunction - Function to evaluate nonlinear residual
4485: Logically Collective
4487: Input Parameters:
4488: + snes - nonlinear solver
4489: . U - the current state at which to evaluate the residual
4490: - ctx - user context, must be a TS
4492: Output Parameter:
4493: . F - the nonlinear residual
4495: Level: advanced
4497: Note:
4498: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4499: It is most frequently passed to `MatFDColoringSetFunction()`.
4501: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4502: @*/
4503: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4504: {
4505: TS ts = (TS)ctx;
4507: PetscFunctionBegin;
4512: PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4513: PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4514: PetscFunctionReturn(PETSC_SUCCESS);
4515: }
4517: /*@
4518: SNESTSFormJacobian - Function to evaluate the Jacobian
4520: Collective
4522: Input Parameters:
4523: + snes - nonlinear solver
4524: . U - the current state at which to evaluate the residual
4525: - ctx - user context, must be a `TS`
4527: Output Parameters:
4528: + A - the Jacobian
4529: - B - the matrix used to construct the preconditioner (often the same as `A`)
4531: Level: developer
4533: Note:
4534: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4536: .seealso: [](ch_ts), `SNESSetJacobian()`
4537: @*/
4538: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4539: {
4540: TS ts = (TS)ctx;
4542: PetscFunctionBegin;
4548: PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4549: PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4550: PetscFunctionReturn(PETSC_SUCCESS);
4551: }
4553: /*@C
4554: TSComputeRHSFunctionLinear - Evaluate the right-hand side via the user-provided Jacobian, for linear problems Udot = A U only
4556: Collective
4558: Input Parameters:
4559: + ts - time stepping context
4560: . t - time at which to evaluate
4561: . U - state at which to evaluate
4562: - ctx - context
4564: Output Parameter:
4565: . F - right-hand side
4567: Level: intermediate
4569: Note:
4570: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right-hand side for linear problems.
4571: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4573: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4574: @*/
4575: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4576: {
4577: Mat Arhs, Brhs;
4579: PetscFunctionBegin;
4580: PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4581: /* undo the damage caused by shifting */
4582: PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4583: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4584: PetscCall(MatMult(Arhs, U, F));
4585: PetscFunctionReturn(PETSC_SUCCESS);
4586: }
4588: /*@C
4589: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
4591: Collective
4593: Input Parameters:
4594: + ts - time stepping context
4595: . t - time at which to evaluate
4596: . U - state at which to evaluate
4597: - ctx - context
4599: Output Parameters:
4600: + A - Jacobian
4601: - B - matrix used to construct the preconditioner, often the same as `A`
4603: Level: intermediate
4605: Note:
4606: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4608: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4609: @*/
4610: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4611: {
4612: PetscFunctionBegin;
4613: PetscFunctionReturn(PETSC_SUCCESS);
4614: }
4616: /*@C
4617: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4619: Collective
4621: Input Parameters:
4622: + ts - time stepping context
4623: . t - time at which to evaluate
4624: . U - state at which to evaluate
4625: . Udot - time derivative of state vector
4626: - ctx - context
4628: Output Parameter:
4629: . F - left hand side
4631: Level: intermediate
4633: Notes:
4634: 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
4635: user is required to write their own `TSComputeIFunction()`.
4636: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4637: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4639: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4641: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4642: @*/
4643: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4644: {
4645: Mat A, B;
4647: PetscFunctionBegin;
4648: PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4649: PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4650: PetscCall(MatMult(A, Udot, F));
4651: PetscFunctionReturn(PETSC_SUCCESS);
4652: }
4654: /*@C
4655: TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE
4657: Collective
4659: Input Parameters:
4660: + ts - time stepping context
4661: . t - time at which to evaluate
4662: . U - state at which to evaluate
4663: . Udot - time derivative of state vector
4664: . shift - shift to apply
4665: - ctx - context
4667: Output Parameters:
4668: + A - pointer to operator
4669: - B - pointer to matrix from which the preconditioner is built (often `A`)
4671: Level: advanced
4673: Notes:
4674: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4676: It is only appropriate for problems of the form
4678: $$
4679: M \dot{U} = F(U,t)
4680: $$
4682: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4683: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4684: an implicit operator of the form
4686: $$
4687: shift*M + J
4688: $$
4690: 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
4691: a copy of M or reassemble it when requested.
4693: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4694: @*/
4695: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4696: {
4697: PetscFunctionBegin;
4698: PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4699: ts->ijacobian.shift = shift;
4700: PetscFunctionReturn(PETSC_SUCCESS);
4701: }
4703: /*@
4704: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4706: Not Collective
4708: Input Parameter:
4709: . ts - the `TS` context
4711: Output Parameter:
4712: . equation_type - see `TSEquationType`
4714: Level: beginner
4716: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4717: @*/
4718: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4719: {
4720: PetscFunctionBegin;
4722: PetscAssertPointer(equation_type, 2);
4723: *equation_type = ts->equation_type;
4724: PetscFunctionReturn(PETSC_SUCCESS);
4725: }
4727: /*@
4728: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4730: Not Collective
4732: Input Parameters:
4733: + ts - the `TS` context
4734: - equation_type - see `TSEquationType`
4736: Level: advanced
4738: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4739: @*/
4740: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4741: {
4742: PetscFunctionBegin;
4744: ts->equation_type = equation_type;
4745: PetscFunctionReturn(PETSC_SUCCESS);
4746: }
4748: /*@
4749: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4751: Not Collective
4753: Input Parameter:
4754: . ts - the `TS` context
4756: Output Parameter:
4757: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4758: manual pages for the individual convergence tests for complete lists
4760: Level: beginner
4762: Note:
4763: Can only be called after the call to `TSSolve()` is complete.
4765: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4766: @*/
4767: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4768: {
4769: PetscFunctionBegin;
4771: PetscAssertPointer(reason, 2);
4772: *reason = ts->reason;
4773: PetscFunctionReturn(PETSC_SUCCESS);
4774: }
4776: /*@
4777: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4779: Logically Collective; reason must contain common value
4781: Input Parameters:
4782: + ts - the `TS` context
4783: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4784: manual pages for the individual convergence tests for complete lists
4786: Level: advanced
4788: Note:
4789: Can only be called while `TSSolve()` is active.
4791: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4792: @*/
4793: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4794: {
4795: PetscFunctionBegin;
4797: ts->reason = reason;
4798: PetscFunctionReturn(PETSC_SUCCESS);
4799: }
4801: /*@
4802: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4804: Not Collective
4806: Input Parameter:
4807: . ts - the `TS` context
4809: Output Parameter:
4810: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4812: Level: beginner
4814: Note:
4815: Can only be called after the call to `TSSolve()` is complete.
4817: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4818: @*/
4819: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4820: {
4821: PetscFunctionBegin;
4823: PetscAssertPointer(ftime, 2);
4824: *ftime = ts->solvetime;
4825: PetscFunctionReturn(PETSC_SUCCESS);
4826: }
4828: /*@
4829: TSGetSNESIterations - Gets the total number of nonlinear iterations
4830: used by the time integrator.
4832: Not Collective
4834: Input Parameter:
4835: . ts - `TS` context
4837: Output Parameter:
4838: . nits - number of nonlinear iterations
4840: Level: intermediate
4842: Note:
4843: This counter is reset to zero for each successive call to `TSSolve()`.
4845: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4846: @*/
4847: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4848: {
4849: PetscFunctionBegin;
4851: PetscAssertPointer(nits, 2);
4852: *nits = ts->snes_its;
4853: PetscFunctionReturn(PETSC_SUCCESS);
4854: }
4856: /*@
4857: TSGetKSPIterations - Gets the total number of linear iterations
4858: used by the time integrator.
4860: Not Collective
4862: Input Parameter:
4863: . ts - `TS` context
4865: Output Parameter:
4866: . lits - number of linear iterations
4868: Level: intermediate
4870: Note:
4871: This counter is reset to zero for each successive call to `TSSolve()`.
4873: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`
4874: @*/
4875: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4876: {
4877: PetscFunctionBegin;
4879: PetscAssertPointer(lits, 2);
4880: *lits = ts->ksp_its;
4881: PetscFunctionReturn(PETSC_SUCCESS);
4882: }
4884: /*@
4885: TSGetStepRejections - Gets the total number of rejected steps.
4887: Not Collective
4889: Input Parameter:
4890: . ts - `TS` context
4892: Output Parameter:
4893: . rejects - number of steps rejected
4895: Level: intermediate
4897: Note:
4898: This counter is reset to zero for each successive call to `TSSolve()`.
4900: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4901: @*/
4902: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4903: {
4904: PetscFunctionBegin;
4906: PetscAssertPointer(rejects, 2);
4907: *rejects = ts->reject;
4908: PetscFunctionReturn(PETSC_SUCCESS);
4909: }
4911: /*@
4912: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
4914: Not Collective
4916: Input Parameter:
4917: . ts - `TS` context
4919: Output Parameter:
4920: . fails - number of failed nonlinear solves
4922: Level: intermediate
4924: Note:
4925: This counter is reset to zero for each successive call to `TSSolve()`.
4927: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4928: @*/
4929: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4930: {
4931: PetscFunctionBegin;
4933: PetscAssertPointer(fails, 2);
4934: *fails = ts->num_snes_failures;
4935: PetscFunctionReturn(PETSC_SUCCESS);
4936: }
4938: /*@
4939: TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails
4941: Not Collective
4943: Input Parameters:
4944: + ts - `TS` context
4945: - rejects - maximum number of rejected steps, pass `PETSC_UNLIMITED` for unlimited
4947: Options Database Key:
4948: . -ts_max_reject - Maximum number of step rejections before a step fails
4950: Level: intermediate
4952: Developer Note:
4953: The options database name is incorrect.
4955: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4956: @*/
4957: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4958: {
4959: PetscFunctionBegin;
4961: if (rejects == PETSC_UNLIMITED || rejects == -1) {
4962: ts->max_reject = PETSC_UNLIMITED;
4963: } else {
4964: PetscCheck(rejects >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of rejections");
4965: ts->max_reject = rejects;
4966: }
4967: PetscFunctionReturn(PETSC_SUCCESS);
4968: }
4970: /*@
4971: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves
4973: Not Collective
4975: Input Parameters:
4976: + ts - `TS` context
4977: - fails - maximum number of failed nonlinear solves, pass `PETSC_UNLIMITED` to allow any number of failures.
4979: Options Database Key:
4980: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
4982: Level: intermediate
4984: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4985: @*/
4986: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4987: {
4988: PetscFunctionBegin;
4990: if (fails == PETSC_UNLIMITED || fails == -1) {
4991: ts->max_snes_failures = PETSC_UNLIMITED;
4992: } else {
4993: PetscCheck(fails >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures");
4994: ts->max_snes_failures = fails;
4995: }
4996: PetscFunctionReturn(PETSC_SUCCESS);
4997: }
4999: /*@
5000: TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`
5002: Not Collective
5004: Input Parameters:
5005: + ts - `TS` context
5006: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
5008: Options Database Key:
5009: . -ts_error_if_step_fails - Error if no step succeeds
5011: Level: intermediate
5013: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
5014: @*/
5015: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
5016: {
5017: PetscFunctionBegin;
5019: ts->errorifstepfailed = err;
5020: PetscFunctionReturn(PETSC_SUCCESS);
5021: }
5023: /*@
5024: TSGetAdapt - Get the adaptive controller context for the current method
5026: Collective if controller has not yet been created
5028: Input Parameter:
5029: . ts - time stepping context
5031: Output Parameter:
5032: . adapt - adaptive controller
5034: Level: intermediate
5036: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
5037: @*/
5038: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
5039: {
5040: PetscFunctionBegin;
5042: PetscAssertPointer(adapt, 2);
5043: if (!ts->adapt) {
5044: PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
5045: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
5046: }
5047: *adapt = ts->adapt;
5048: PetscFunctionReturn(PETSC_SUCCESS);
5049: }
5051: /*@
5052: TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller
5054: Logically Collective
5056: Input Parameters:
5057: + ts - time integration context
5058: . atol - scalar absolute tolerances
5059: . vatol - vector of absolute tolerances or `NULL`, used in preference to `atol` if present
5060: . rtol - scalar relative tolerances
5061: - vrtol - vector of relative tolerances or `NULL`, used in preference to `rtol` if present
5063: Options Database Keys:
5064: + -ts_rtol <rtol> - relative tolerance for local truncation error
5065: - -ts_atol <atol> - Absolute tolerance for local truncation error
5067: Level: beginner
5069: Notes:
5070: `PETSC_CURRENT` or `PETSC_DETERMINE` may be used for `atol` or `rtol` to indicate the current value
5071: or the default value from when the object's type was set.
5073: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5074: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5075: computed only for the differential or the algebraic part then this can be done using the vector of
5076: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
5077: differential part and infinity for the algebraic part, the LTE calculation will include only the
5078: differential variables.
5080: Fortran Note:
5081: Use `PETSC_CURRENT_INTEGER` or `PETSC_DETERMINE_INTEGER`.
5083: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
5084: @*/
5085: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
5086: {
5087: PetscFunctionBegin;
5088: if (atol == (PetscReal)PETSC_DETERMINE) {
5089: ts->atol = ts->default_atol;
5090: } else if (atol != (PetscReal)PETSC_CURRENT) {
5091: PetscCheck(atol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)atol);
5092: ts->atol = atol;
5093: }
5095: if (vatol) {
5096: PetscCall(PetscObjectReference((PetscObject)vatol));
5097: PetscCall(VecDestroy(&ts->vatol));
5098: ts->vatol = vatol;
5099: }
5101: if (rtol == (PetscReal)PETSC_DETERMINE) {
5102: ts->rtol = ts->default_rtol;
5103: } else if (rtol != (PetscReal)PETSC_CURRENT) {
5104: PetscCheck(rtol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative", (double)rtol);
5105: ts->rtol = rtol;
5106: }
5108: if (vrtol) {
5109: PetscCall(PetscObjectReference((PetscObject)vrtol));
5110: PetscCall(VecDestroy(&ts->vrtol));
5111: ts->vrtol = vrtol;
5112: }
5113: PetscFunctionReturn(PETSC_SUCCESS);
5114: }
5116: /*@
5117: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
5119: Logically Collective
5121: Input Parameter:
5122: . ts - time integration context
5124: Output Parameters:
5125: + atol - scalar absolute tolerances, `NULL` to ignore
5126: . vatol - vector of absolute tolerances, `NULL` to ignore
5127: . rtol - scalar relative tolerances, `NULL` to ignore
5128: - vrtol - vector of relative tolerances, `NULL` to ignore
5130: Level: beginner
5132: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
5133: @*/
5134: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
5135: {
5136: PetscFunctionBegin;
5137: if (atol) *atol = ts->atol;
5138: if (vatol) *vatol = ts->vatol;
5139: if (rtol) *rtol = ts->rtol;
5140: if (vrtol) *vrtol = ts->vrtol;
5141: PetscFunctionReturn(PETSC_SUCCESS);
5142: }
5144: /*@
5145: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
5147: Collective
5149: Input Parameters:
5150: + ts - time stepping context
5151: . U - state vector, usually ts->vec_sol
5152: . Y - state vector to be compared to U
5153: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5155: Output Parameters:
5156: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5157: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5158: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5160: Options Database Key:
5161: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5163: Level: developer
5165: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5166: @*/
5167: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5168: {
5169: PetscInt norma_loc, norm_loc, normr_loc;
5171: PetscFunctionBegin;
5173: 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));
5174: if (wnormtype == NORM_2) {
5175: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5176: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5177: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5178: }
5179: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5180: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5181: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5182: PetscFunctionReturn(PETSC_SUCCESS);
5183: }
5185: /*@
5186: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5188: Collective
5190: Input Parameters:
5191: + ts - time stepping context
5192: . E - error vector
5193: . U - state vector, usually ts->vec_sol
5194: . Y - state vector, previous time step
5195: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5197: Output Parameters:
5198: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5199: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5200: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5202: Options Database Key:
5203: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5205: Level: developer
5207: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5208: @*/
5209: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5210: {
5211: PetscInt norma_loc, norm_loc, normr_loc;
5213: PetscFunctionBegin;
5215: 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));
5216: if (wnormtype == NORM_2) {
5217: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5218: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5219: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5220: }
5221: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5222: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5223: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5224: PetscFunctionReturn(PETSC_SUCCESS);
5225: }
5227: /*@
5228: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5230: Logically Collective
5232: Input Parameters:
5233: + ts - time stepping context
5234: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5236: Note:
5237: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5239: Level: intermediate
5241: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5242: @*/
5243: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5244: {
5245: PetscFunctionBegin;
5247: ts->cfltime_local = cfltime;
5248: ts->cfltime = -1.;
5249: PetscFunctionReturn(PETSC_SUCCESS);
5250: }
5252: /*@
5253: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5255: Collective
5257: Input Parameter:
5258: . ts - time stepping context
5260: Output Parameter:
5261: . cfltime - maximum stable time step for forward Euler
5263: Level: advanced
5265: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5266: @*/
5267: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5268: {
5269: PetscFunctionBegin;
5270: if (ts->cfltime < 0) PetscCallMPI(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5271: *cfltime = ts->cfltime;
5272: PetscFunctionReturn(PETSC_SUCCESS);
5273: }
5275: /*@
5276: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5278: Input Parameters:
5279: + ts - the `TS` context.
5280: . xl - lower bound.
5281: - xu - upper bound.
5283: Level: advanced
5285: Note:
5286: If this routine is not called then the lower and upper bounds are set to
5287: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5289: .seealso: [](ch_ts), `TS`
5290: @*/
5291: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5292: {
5293: SNES snes;
5295: PetscFunctionBegin;
5296: PetscCall(TSGetSNES(ts, &snes));
5297: PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5298: PetscFunctionReturn(PETSC_SUCCESS);
5299: }
5301: /*@
5302: TSComputeLinearStability - computes the linear stability function at a point
5304: Collective
5306: Input Parameters:
5307: + ts - the `TS` context
5308: . xr - real part of input argument
5309: - xi - imaginary part of input argument
5311: Output Parameters:
5312: + yr - real part of function value
5313: - yi - imaginary part of function value
5315: Level: developer
5317: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5318: @*/
5319: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5320: {
5321: PetscFunctionBegin;
5323: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5324: PetscFunctionReturn(PETSC_SUCCESS);
5325: }
5327: /*@
5328: TSRestartStep - Flags the solver to restart the next step
5330: Collective
5332: Input Parameter:
5333: . ts - the `TS` context obtained from `TSCreate()`
5335: Level: advanced
5337: Notes:
5338: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5339: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5340: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5341: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5342: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5343: discontinuous source terms).
5345: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5346: @*/
5347: PetscErrorCode TSRestartStep(TS ts)
5348: {
5349: PetscFunctionBegin;
5351: ts->steprestart = PETSC_TRUE;
5352: PetscFunctionReturn(PETSC_SUCCESS);
5353: }
5355: /*@
5356: TSRollBack - Rolls back one time step
5358: Collective
5360: Input Parameter:
5361: . ts - the `TS` context obtained from `TSCreate()`
5363: Level: advanced
5365: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5366: @*/
5367: PetscErrorCode TSRollBack(TS ts)
5368: {
5369: PetscFunctionBegin;
5371: PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5372: PetscTryTypeMethod(ts, rollback);
5373: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
5374: ts->time_step = ts->ptime - ts->ptime_prev;
5375: ts->ptime = ts->ptime_prev;
5376: ts->ptime_prev = ts->ptime_prev_rollback;
5377: ts->steps--;
5378: ts->steprollback = PETSC_TRUE;
5379: PetscFunctionReturn(PETSC_SUCCESS);
5380: }
5382: /*@
5383: TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step
5385: Not collective
5387: Input Parameter:
5388: . ts - the `TS` context obtained from `TSCreate()`
5390: Output Parameter:
5391: . flg - the rollback flag
5393: Level: advanced
5395: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
5396: @*/
5397: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
5398: {
5399: PetscFunctionBegin;
5401: PetscAssertPointer(flg, 2);
5402: *flg = ts->steprollback;
5403: PetscFunctionReturn(PETSC_SUCCESS);
5404: }
5406: /*@
5407: TSGetStepResize - Get the internal flag indicating if the current step is after a resize.
5409: Not collective
5411: Input Parameter:
5412: . ts - the `TS` context obtained from `TSCreate()`
5414: Output Parameter:
5415: . flg - the resize flag
5417: Level: advanced
5419: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetResize()`
5420: @*/
5421: PetscErrorCode TSGetStepResize(TS ts, PetscBool *flg)
5422: {
5423: PetscFunctionBegin;
5425: PetscAssertPointer(flg, 2);
5426: *flg = ts->stepresize;
5427: PetscFunctionReturn(PETSC_SUCCESS);
5428: }
5430: /*@
5431: TSGetStages - Get the number of stages and stage values
5433: Input Parameter:
5434: . ts - the `TS` context obtained from `TSCreate()`
5436: Output Parameters:
5437: + ns - the number of stages
5438: - Y - the current stage vectors
5440: Level: advanced
5442: Note:
5443: Both `ns` and `Y` can be `NULL`.
5445: .seealso: [](ch_ts), `TS`, `TSCreate()`
5446: @*/
5447: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5448: {
5449: PetscFunctionBegin;
5451: if (ns) PetscAssertPointer(ns, 2);
5452: if (Y) PetscAssertPointer(Y, 3);
5453: if (!ts->ops->getstages) {
5454: if (ns) *ns = 0;
5455: if (Y) *Y = NULL;
5456: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5457: PetscFunctionReturn(PETSC_SUCCESS);
5458: }
5460: /*@C
5461: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5463: Collective
5465: Input Parameters:
5466: + ts - the `TS` context
5467: . t - current timestep
5468: . U - state vector
5469: . Udot - time derivative of state vector
5470: . shift - shift to apply, see note below
5471: - ctx - an optional user context
5473: Output Parameters:
5474: + J - Jacobian matrix (not altered in this routine)
5475: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
5477: Level: intermediate
5479: Notes:
5480: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5482: dF/dU + shift*dF/dUdot
5484: Most users should not need to explicitly call this routine, as it
5485: is used internally within the nonlinear solvers.
5487: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5488: routine, then it will try to get the coloring from the matrix. This requires that the
5489: matrix have nonzero entries precomputed.
5491: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5492: @*/
5493: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5494: {
5495: SNES snes;
5496: MatFDColoring color;
5497: PetscBool hascolor, matcolor = PETSC_FALSE;
5499: PetscFunctionBegin;
5500: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5501: PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5502: if (!color) {
5503: DM dm;
5504: ISColoring iscoloring;
5506: PetscCall(TSGetDM(ts, &dm));
5507: PetscCall(DMHasColoring(dm, &hascolor));
5508: if (hascolor && !matcolor) {
5509: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5510: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5511: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode (*)(void))SNESTSFormFunction, (void *)ts));
5512: PetscCall(MatFDColoringSetFromOptions(color));
5513: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5514: PetscCall(ISColoringDestroy(&iscoloring));
5515: } else {
5516: MatColoring mc;
5518: PetscCall(MatColoringCreate(B, &mc));
5519: PetscCall(MatColoringSetDistance(mc, 2));
5520: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5521: PetscCall(MatColoringSetFromOptions(mc));
5522: PetscCall(MatColoringApply(mc, &iscoloring));
5523: PetscCall(MatColoringDestroy(&mc));
5524: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5525: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode (*)(void))SNESTSFormFunction, (void *)ts));
5526: PetscCall(MatFDColoringSetFromOptions(color));
5527: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5528: PetscCall(ISColoringDestroy(&iscoloring));
5529: }
5530: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5531: PetscCall(PetscObjectDereference((PetscObject)color));
5532: }
5533: PetscCall(TSGetSNES(ts, &snes));
5534: PetscCall(MatFDColoringApply(B, color, U, snes));
5535: if (J != B) {
5536: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5537: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5538: }
5539: PetscFunctionReturn(PETSC_SUCCESS);
5540: }
5542: /*@C
5543: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5545: Input Parameters:
5546: + ts - the `TS` context
5547: - func - function called within `TSFunctionDomainError()`
5549: Calling sequence of `func`:
5550: + ts - the `TS` context
5551: . time - the current time (of the stage)
5552: . state - the state to check if it is valid
5553: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable
5555: Level: intermediate
5557: Notes:
5558: If an implicit ODE solver is being used then, in addition to providing this routine, the
5559: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5560: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5561: Use `TSGetSNES()` to obtain the `SNES` object
5563: Developer Notes:
5564: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5565: since one takes a function pointer and the other does not.
5567: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5568: @*/
5569: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5570: {
5571: PetscFunctionBegin;
5573: ts->functiondomainerror = func;
5574: PetscFunctionReturn(PETSC_SUCCESS);
5575: }
5577: /*@
5578: TSFunctionDomainError - Checks if the current state is valid
5580: Input Parameters:
5581: + ts - the `TS` context
5582: . stagetime - time of the simulation
5583: - Y - state vector to check.
5585: Output Parameter:
5586: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5588: Level: developer
5590: Note:
5591: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5592: to check if the current state is valid.
5594: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5595: @*/
5596: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5597: {
5598: PetscFunctionBegin;
5600: *accept = PETSC_TRUE;
5601: if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5602: PetscFunctionReturn(PETSC_SUCCESS);
5603: }
5605: /*@
5606: TSClone - This function clones a time step `TS` object.
5608: Collective
5610: Input Parameter:
5611: . tsin - The input `TS`
5613: Output Parameter:
5614: . tsout - The output `TS` (cloned)
5616: Level: developer
5618: Notes:
5619: 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.
5620: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5622: 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
5623: .vb
5624: SNES snes_dup = NULL;
5625: TSGetSNES(ts,&snes_dup);
5626: TSSetSNES(ts,snes_dup);
5627: .ve
5629: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5630: @*/
5631: PetscErrorCode TSClone(TS tsin, TS *tsout)
5632: {
5633: TS t;
5634: SNES snes_start;
5635: DM dm;
5636: TSType type;
5638: PetscFunctionBegin;
5639: PetscAssertPointer(tsin, 1);
5640: *tsout = NULL;
5642: PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));
5644: /* General TS description */
5645: t->numbermonitors = 0;
5646: t->monitorFrequency = 1;
5647: t->setupcalled = 0;
5648: t->ksp_its = 0;
5649: t->snes_its = 0;
5650: t->nwork = 0;
5651: t->rhsjacobian.time = PETSC_MIN_REAL;
5652: t->rhsjacobian.scale = 1.;
5653: t->ijacobian.shift = 1.;
5655: PetscCall(TSGetSNES(tsin, &snes_start));
5656: PetscCall(TSSetSNES(t, snes_start));
5658: PetscCall(TSGetDM(tsin, &dm));
5659: PetscCall(TSSetDM(t, dm));
5661: t->adapt = tsin->adapt;
5662: PetscCall(PetscObjectReference((PetscObject)t->adapt));
5664: t->trajectory = tsin->trajectory;
5665: PetscCall(PetscObjectReference((PetscObject)t->trajectory));
5667: t->event = tsin->event;
5668: if (t->event) t->event->refct++;
5670: t->problem_type = tsin->problem_type;
5671: t->ptime = tsin->ptime;
5672: t->ptime_prev = tsin->ptime_prev;
5673: t->time_step = tsin->time_step;
5674: t->max_time = tsin->max_time;
5675: t->steps = tsin->steps;
5676: t->max_steps = tsin->max_steps;
5677: t->equation_type = tsin->equation_type;
5678: t->atol = tsin->atol;
5679: t->rtol = tsin->rtol;
5680: t->max_snes_failures = tsin->max_snes_failures;
5681: t->max_reject = tsin->max_reject;
5682: t->errorifstepfailed = tsin->errorifstepfailed;
5684: PetscCall(TSGetType(tsin, &type));
5685: PetscCall(TSSetType(t, type));
5687: t->vec_sol = NULL;
5689: t->cfltime = tsin->cfltime;
5690: t->cfltime_local = tsin->cfltime_local;
5691: t->exact_final_time = tsin->exact_final_time;
5693: t->ops[0] = tsin->ops[0];
5695: if (((PetscObject)tsin)->fortran_func_pointers) {
5696: PetscInt i;
5697: PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers));
5698: for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5699: }
5700: *tsout = t;
5701: PetscFunctionReturn(PETSC_SUCCESS);
5702: }
5704: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5705: {
5706: TS ts = (TS)ctx;
5708: PetscFunctionBegin;
5709: PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5710: PetscFunctionReturn(PETSC_SUCCESS);
5711: }
5713: /*@
5714: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5716: Logically Collective
5718: Input Parameter:
5719: . ts - the time stepping routine
5721: Output Parameter:
5722: . flg - `PETSC_TRUE` if the multiply is likely correct
5724: Options Database Key:
5725: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5727: Level: advanced
5729: Note:
5730: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5732: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5733: @*/
5734: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5735: {
5736: Mat J, B;
5737: TSRHSJacobianFn *func;
5738: void *ctx;
5740: PetscFunctionBegin;
5741: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5742: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5743: PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5744: PetscFunctionReturn(PETSC_SUCCESS);
5745: }
5747: /*@
5748: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5750: Logically Collective
5752: Input Parameter:
5753: . ts - the time stepping routine
5755: Output Parameter:
5756: . flg - `PETSC_TRUE` if the multiply is likely correct
5758: Options Database Key:
5759: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5761: Level: advanced
5763: Notes:
5764: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5766: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5767: @*/
5768: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5769: {
5770: Mat J, B;
5771: void *ctx;
5772: TSRHSJacobianFn *func;
5774: PetscFunctionBegin;
5775: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5776: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5777: PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5778: PetscFunctionReturn(PETSC_SUCCESS);
5779: }
5781: /*@
5782: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
5784: Logically Collective
5786: Input Parameters:
5787: + ts - timestepping context
5788: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5790: Options Database Key:
5791: . -ts_use_splitrhsfunction - <true,false>
5793: Level: intermediate
5795: Note:
5796: This is only for multirate methods
5798: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5799: @*/
5800: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5801: {
5802: PetscFunctionBegin;
5804: ts->use_splitrhsfunction = use_splitrhsfunction;
5805: PetscFunctionReturn(PETSC_SUCCESS);
5806: }
5808: /*@
5809: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
5811: Not Collective
5813: Input Parameter:
5814: . ts - timestepping context
5816: Output Parameter:
5817: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5819: Level: intermediate
5821: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5822: @*/
5823: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5824: {
5825: PetscFunctionBegin;
5827: *use_splitrhsfunction = ts->use_splitrhsfunction;
5828: PetscFunctionReturn(PETSC_SUCCESS);
5829: }
5831: /*@
5832: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
5834: Logically Collective
5836: Input Parameters:
5837: + ts - the time-stepper
5838: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
5840: Level: intermediate
5842: Note:
5843: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
5845: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
5846: @*/
5847: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5848: {
5849: PetscFunctionBegin;
5851: ts->axpy_pattern = str;
5852: PetscFunctionReturn(PETSC_SUCCESS);
5853: }
5855: /*@
5856: TSSetEvaluationTimes - sets the evaluation points. The solution will be computed and stored for each time requested
5858: Collective
5860: Input Parameters:
5861: + ts - the time-stepper
5862: . n - number of the time points
5863: - time_points - array of the time points
5865: Options Database Key:
5866: . -ts_eval_times <t0,...tn> - Sets the evaluation times
5868: Level: intermediate
5870: Notes:
5871: The elements in `time_points` must be all increasing. They correspond to the intermediate points for time integration.
5873: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
5875: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using evaluation times may
5876: pressure the memory system when using a large number of time points.
5878: .seealso: [](ch_ts), `TS`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`, `TSSetTimeSpan()`
5879: @*/
5880: PetscErrorCode TSSetEvaluationTimes(TS ts, PetscInt n, PetscReal time_points[])
5881: {
5882: PetscBool is_sorted;
5884: PetscFunctionBegin;
5886: if (ts->eval_times && n != ts->eval_times->num_time_points) {
5887: PetscCall(PetscFree(ts->eval_times->time_points));
5888: PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
5889: PetscCall(PetscMalloc1(n, &ts->eval_times->time_points));
5890: }
5891: if (!ts->eval_times) {
5892: TSEvaluationTimes eval_times;
5893: PetscCall(PetscNew(&eval_times));
5894: PetscCall(PetscMalloc1(n, &eval_times->time_points));
5895: eval_times->reltol = 1e-6;
5896: eval_times->abstol = 10 * PETSC_MACHINE_EPSILON;
5897: eval_times->worktol = 0;
5898: ts->eval_times = eval_times;
5899: }
5900: ts->eval_times->num_time_points = n;
5901: PetscCall(PetscSortedReal(n, time_points, &is_sorted));
5902: PetscCheck(is_sorted, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "time_points array must be sorted");
5903: PetscCall(PetscArraycpy(ts->eval_times->time_points, time_points, n));
5904: PetscFunctionReturn(PETSC_SUCCESS);
5905: }
5907: /*@C
5908: TSGetEvaluationTimes - gets the evaluation times set with `TSSetEvaluationTimes()`
5910: Not Collective
5912: Input Parameter:
5913: . ts - the time-stepper
5915: Output Parameters:
5916: + n - number of the time points
5917: - time_points - array of the time points
5919: Level: beginner
5921: Note:
5922: The values obtained are valid until the `TS` object is destroyed.
5924: Both `n` and `time_points` can be `NULL`.
5926: Also used to see time points set by `TSSetTimeSpan()`.
5928: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationSolutions()`
5929: @*/
5930: PetscErrorCode TSGetEvaluationTimes(TS ts, PetscInt *n, const PetscReal *time_points[])
5931: {
5932: PetscFunctionBegin;
5934: if (n) PetscAssertPointer(n, 2);
5935: if (time_points) PetscAssertPointer(time_points, 3);
5936: if (!ts->eval_times) {
5937: if (n) *n = 0;
5938: if (time_points) *time_points = NULL;
5939: } else {
5940: if (n) *n = ts->eval_times->num_time_points;
5941: if (time_points) *time_points = ts->eval_times->time_points;
5942: }
5943: PetscFunctionReturn(PETSC_SUCCESS);
5944: }
5946: /*@C
5947: TSGetEvaluationSolutions - Get the number of solutions and the solutions at the evaluation time points specified
5949: Input Parameter:
5950: . ts - the `TS` context obtained from `TSCreate()`
5952: Output Parameters:
5953: + nsol - the number of solutions
5954: . sol_times - array of solution times corresponding to the solution vectors. See note below
5955: - Sols - the solution vectors
5957: Level: intermediate
5959: Notes:
5960: Both `nsol` and `Sols` can be `NULL`.
5962: 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()`.
5963: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain evaluation times.
5965: Also used to see view solutions requested by `TSSetTimeSpan()`.
5967: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`
5968: @*/
5969: PetscErrorCode TSGetEvaluationSolutions(TS ts, PetscInt *nsol, const PetscReal *sol_times[], Vec *Sols[])
5970: {
5971: PetscFunctionBegin;
5973: if (nsol) PetscAssertPointer(nsol, 2);
5974: if (sol_times) PetscAssertPointer(sol_times, 3);
5975: if (Sols) PetscAssertPointer(Sols, 4);
5976: if (!ts->eval_times) {
5977: if (nsol) *nsol = 0;
5978: if (sol_times) *sol_times = NULL;
5979: if (Sols) *Sols = NULL;
5980: } else {
5981: if (nsol) *nsol = ts->eval_times->sol_ctr;
5982: if (sol_times) *sol_times = ts->eval_times->sol_times;
5983: if (Sols) *Sols = ts->eval_times->sol_vecs;
5984: }
5985: PetscFunctionReturn(PETSC_SUCCESS);
5986: }
5988: /*@
5989: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
5991: Collective
5993: Input Parameters:
5994: + ts - the time-stepper
5995: . n - number of the time points (>=2)
5996: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5998: Options Database Key:
5999: . -ts_time_span <t0,...tf> - Sets the time span
6001: Level: intermediate
6003: Notes:
6004: 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.
6006: The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
6008: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
6010: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using time span may
6011: pressure the memory system when using a large number of span points.
6013: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6014: @*/
6015: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal span_times[])
6016: {
6017: PetscFunctionBegin;
6019: PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
6020: PetscCall(TSSetEvaluationTimes(ts, n, span_times));
6021: PetscCall(TSSetTime(ts, span_times[0]));
6022: PetscCall(TSSetMaxTime(ts, span_times[n - 1]));
6023: PetscFunctionReturn(PETSC_SUCCESS);
6024: }
6026: /*@
6027: TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
6029: Collective
6031: Input Parameters:
6032: + ts - the `TS` context
6033: . J - Jacobian matrix (not altered in this routine)
6034: - B - newly computed Jacobian matrix to use with preconditioner
6036: Level: intermediate
6038: Notes:
6039: This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
6040: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
6041: and multiple fields are involved.
6043: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
6044: structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
6045: usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
6046: `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.
6048: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
6049: @*/
6050: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
6051: {
6052: MatColoring mc = NULL;
6053: ISColoring iscoloring = NULL;
6054: MatFDColoring matfdcoloring = NULL;
6056: PetscFunctionBegin;
6057: /* Generate new coloring after eliminating zeros in the matrix */
6058: PetscCall(MatEliminateZeros(B, PETSC_TRUE));
6059: PetscCall(MatColoringCreate(B, &mc));
6060: PetscCall(MatColoringSetDistance(mc, 2));
6061: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
6062: PetscCall(MatColoringSetFromOptions(mc));
6063: PetscCall(MatColoringApply(mc, &iscoloring));
6064: PetscCall(MatColoringDestroy(&mc));
6065: /* Replace the old coloring with the new one */
6066: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
6067: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode (*)(void))SNESTSFormFunction, (void *)ts));
6068: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
6069: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
6070: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
6071: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
6072: PetscCall(ISColoringDestroy(&iscoloring));
6073: PetscFunctionReturn(PETSC_SUCCESS);
6074: }