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
38: . -ts_max_steps <steps> - maximum number of time-steps to take
39: . -ts_init_time <time> - initial time to start computation
40: . -ts_final_time <time> - final time to compute to (deprecated: use `-ts_max_time`)
41: . -ts_dt <dt> - initial time step
42: . -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
43: . -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
44: . -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
45: . -ts_error_if_step_fails <true,false> - Error if no step succeeds
46: . -ts_rtol <rtol> - relative tolerance for local truncation error
47: . -ts_atol <atol> - Absolute tolerance for local truncation error
48: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function
49: . -ts_rhs_jacobian_test_mult_transpose - test the Jacobian at each iteration against finite difference with RHS function
50: . -ts_adjoint_solve <yes,no> - After solving the ODE/DAE solve the adjoint problem (requires `-ts_save_trajectory`)
51: . -ts_fd_color - Use finite differences with coloring to compute IJacobian
52: . -ts_monitor - print information at each timestep
53: . -ts_monitor_cancel - Cancel all monitors
54: . -ts_monitor_lg_solution - Monitor solution graphically
55: . -ts_monitor_lg_error - Monitor error graphically
56: . -ts_monitor_error - Monitors norm of error
57: . -ts_monitor_lg_timestep - Monitor timestep size graphically
58: . -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
59: . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
60: . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
61: . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
62: . -ts_monitor_draw_solution - Monitor solution graphically
63: . -ts_monitor_draw_solution_phase <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
64: . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
65: . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
66: . -ts_monitor_solution_interval <interval> - output once every interval (default=1) time steps. Use -1 to only output at the end of the simulation
67: . -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)
68: . -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
69: - -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
71: Level: beginner
73: Notes:
74: See `SNESSetFromOptions()` and `KSPSetFromOptions()` for how to control the nonlinear and linear solves used by the time-stepper.
76: Certain `SNES` options get reset for each new nonlinear solver, for example `-snes_lag_jacobian its` and `-snes_lag_preconditioner its`, in order
77: to retain them over the multiple nonlinear solves that `TS` uses you must also provide `-snes_lag_jacobian_persists true` and
78: `-snes_lag_preconditioner_persists true`
80: Developer Notes:
81: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
83: .seealso: [](ch_ts), `TS`, `TSGetType()`
84: @*/
85: PetscErrorCode TSSetFromOptions(TS ts)
86: {
87: PetscBool opt, flg, tflg;
88: char monfilename[PETSC_MAX_PATH_LEN];
89: PetscReal time_step, tspan[100];
90: PetscInt nt = PETSC_STATIC_ARRAY_LENGTH(tspan);
91: TSExactFinalTimeOption eftopt;
92: char dir[16];
93: TSIFunctionFn *ifun;
94: const char *defaultType;
95: char typeName[256];
97: PetscFunctionBegin;
100: PetscCall(TSRegisterAll());
101: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
103: PetscObjectOptionsBegin((PetscObject)ts);
104: if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
105: else defaultType = ifun ? TSBEULER : TSEULER;
106: PetscCall(PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, 256, &opt));
107: if (opt) PetscCall(TSSetType(ts, typeName));
108: else PetscCall(TSSetType(ts, defaultType));
110: /* Handle generic TS options */
111: PetscCall(PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL));
112: PetscCall(PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL));
113: PetscCall(PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", tspan, &nt, &flg));
114: if (flg) PetscCall(TSSetTimeSpan(ts, nt, tspan));
115: PetscCall(PetscOptionsInt("-ts_max_steps", "Maximum number of time steps", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL));
116: PetscCall(PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL));
117: PetscCall(PetscOptionsReal("-ts_dt", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg));
118: if (flg) PetscCall(TSSetTimeStep(ts, time_step));
119: PetscCall(PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg));
120: if (flg) PetscCall(TSSetExactFinalTime(ts, eftopt));
121: PetscCall(PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, &flg));
122: if (flg) PetscCall(TSSetMaxSNESFailures(ts, ts->max_snes_failures));
123: PetscCall(PetscOptionsInt("-ts_max_reject", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, &flg));
124: if (flg) PetscCall(TSSetMaxStepRejections(ts, ts->max_reject));
125: PetscCall(PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL));
126: PetscCall(PetscOptionsBoundedReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL, 0));
127: PetscCall(PetscOptionsBoundedReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL, 0));
129: PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL));
130: 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));
131: PetscCall(PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL));
132: #if defined(PETSC_HAVE_SAWS)
133: {
134: PetscBool set;
135: flg = PETSC_FALSE;
136: PetscCall(PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set));
137: if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)ts, flg));
138: }
139: #endif
141: /* Monitor options */
142: PetscCall(PetscOptionsInt("-ts_monitor_frequency", "Number of time steps between monitor output", "TSMonitorSetFrequency", ts->monitorFrequency, &ts->monitorFrequency, NULL));
143: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL));
144: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL));
145: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, NULL));
146: PetscCall(TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL));
148: PetscCall(PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg));
149: if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)ts, monfilename));
151: PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt));
152: if (opt) {
153: PetscInt howoften = 1;
154: DM dm;
155: PetscBool net;
157: PetscCall(PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL));
158: PetscCall(TSGetDM(ts, &dm));
159: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net));
160: if (net) {
161: TSMonitorLGCtxNetwork ctx;
162: PetscCall(TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx));
163: PetscCall(TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxNetworkDestroy));
164: PetscCall(PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL));
165: } else {
166: TSMonitorLGCtx ctx;
167: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
168: PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
169: }
170: }
172: PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt));
173: if (opt) {
174: TSMonitorLGCtx ctx;
175: PetscInt howoften = 1;
177: PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL));
178: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
179: PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
180: }
181: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL));
183: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt));
184: if (opt) {
185: TSMonitorLGCtx ctx;
186: PetscInt howoften = 1;
188: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
189: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
190: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
191: }
192: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt));
193: if (opt) {
194: TSMonitorLGCtx ctx;
195: PetscInt howoften = 1;
197: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
198: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
199: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
200: ctx->semilogy = PETSC_TRUE;
201: }
203: PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt));
204: if (opt) {
205: TSMonitorLGCtx ctx;
206: PetscInt howoften = 1;
208: PetscCall(PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL));
209: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
210: PetscCall(TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
211: }
212: PetscCall(PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt));
213: if (opt) {
214: TSMonitorLGCtx ctx;
215: PetscInt howoften = 1;
217: PetscCall(PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL));
218: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
219: PetscCall(TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
220: }
221: PetscCall(PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt));
222: if (opt) {
223: TSMonitorSPEigCtx ctx;
224: PetscInt howoften = 1;
226: PetscCall(PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL));
227: PetscCall(TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
228: PetscCall(TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscCtxDestroyFn *)TSMonitorSPEigCtxDestroy));
229: }
230: PetscCall(PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase space from the DMSwarm", "TSMonitorSPSwarm", &opt));
231: if (opt) {
232: TSMonitorSPCtx ctx;
233: PetscInt howoften = 1, retain = 0;
234: PetscBool phase = PETSC_TRUE, create = PETSC_TRUE, multispecies = PETSC_FALSE;
236: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
237: if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
238: create = PETSC_FALSE;
239: break;
240: }
241: if (create) {
242: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase space from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL));
243: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL));
244: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL));
245: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_multi_species", "Color particles by particle species", "TSMonitorSPSwarm", multispecies, &multispecies, NULL));
246: PetscCall(TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, multispecies, &ctx));
247: PetscCall(TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscCtxDestroyFn *)TSMonitorSPCtxDestroy));
248: }
249: }
250: PetscCall(PetscOptionsName("-ts_monitor_hg_swarm", "Display particle histogram from the DMSwarm", "TSMonitorHGSwarm", &opt));
251: if (opt) {
252: TSMonitorHGCtx ctx;
253: PetscInt howoften = 1, Ns = 1;
254: PetscBool velocity = PETSC_FALSE, create = PETSC_TRUE;
256: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
257: if (ts->monitor[i] == TSMonitorHGSwarmSolution) {
258: create = PETSC_FALSE;
259: break;
260: }
261: if (create) {
262: DM sw, dm;
263: PetscInt Nc, Nb;
265: PetscCall(TSGetDM(ts, &sw));
266: PetscCall(DMSwarmGetCellDM(sw, &dm));
267: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc));
268: Nb = PetscMin(20, PetscMax(10, Nc));
269: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL));
270: PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL));
271: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL));
272: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL));
273: PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx));
274: PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscCtxDestroyFn *)TSMonitorHGCtxDestroy));
275: }
276: }
277: opt = PETSC_FALSE;
278: PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt));
279: if (opt) {
280: TSMonitorDrawCtx ctx;
281: PetscInt howoften = 1;
283: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL));
284: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
285: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
286: }
287: opt = PETSC_FALSE;
288: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt));
289: if (opt) {
290: TSMonitorDrawCtx ctx;
291: PetscReal bounds[4];
292: PetscInt n = 4;
293: PetscDraw draw;
294: PetscDrawAxis axis;
296: PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL));
297: PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field");
298: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx));
299: PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw));
300: PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis));
301: PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]));
302: PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2"));
303: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
304: }
305: opt = PETSC_FALSE;
306: PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt));
307: if (opt) {
308: TSMonitorDrawCtx ctx;
309: PetscInt howoften = 1;
311: PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL));
312: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
313: PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
314: }
315: opt = PETSC_FALSE;
316: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt));
317: if (opt) {
318: TSMonitorDrawCtx ctx;
319: PetscInt howoften = 1;
321: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL));
322: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
323: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
324: }
326: opt = PETSC_FALSE;
327: 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));
328: if (flg) {
329: TSMonitorVTKCtx ctx;
331: PetscCall(TSMonitorSolutionVTKCtxCreate(monfilename, &ctx));
332: PetscCall(PetscOptionsInt("-ts_monitor_solution_vtk_interval", "Save every interval time step (-1 for last step only)", NULL, ctx->interval, &ctx->interval, NULL));
333: PetscCall(TSMonitorSet(ts, (PetscErrorCode (*)(TS, PetscInt, PetscReal, Vec, void *))TSMonitorSolutionVTK, ctx, (PetscCtxDestroyFn *)TSMonitorSolutionVTKDestroy));
334: }
336: PetscCall(PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg));
337: if (flg) {
338: TSMonitorDMDARayCtx *rayctx;
339: int ray = 0;
340: DMDirection ddir;
341: DM da;
342: PetscMPIInt rank;
344: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
345: if (dir[0] == 'x') ddir = DM_X;
346: else if (dir[0] == 'y') ddir = DM_Y;
347: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
348: sscanf(dir + 2, "%d", &ray);
350: PetscCall(PetscInfo(ts, "Displaying DMDA ray %c = %d\n", dir[0], ray));
351: PetscCall(PetscNew(&rayctx));
352: PetscCall(TSGetDM(ts, &da));
353: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
354: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank));
355: if (rank == 0) PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer));
356: rayctx->lgctx = NULL;
357: PetscCall(TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy));
358: }
359: PetscCall(PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg));
360: if (flg) {
361: TSMonitorDMDARayCtx *rayctx;
362: int ray = 0;
363: DMDirection ddir;
364: DM da;
365: PetscInt howoften = 1;
367: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
368: if (dir[0] == 'x') ddir = DM_X;
369: else if (dir[0] == 'y') ddir = DM_Y;
370: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
371: sscanf(dir + 2, "%d", &ray);
373: PetscCall(PetscInfo(ts, "Displaying LG DMDA ray %c = %d\n", dir[0], ray));
374: PetscCall(PetscNew(&rayctx));
375: PetscCall(TSGetDM(ts, &da));
376: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
377: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx));
378: PetscCall(TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy));
379: }
381: PetscCall(PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt));
382: if (opt) {
383: TSMonitorEnvelopeCtx ctx;
385: PetscCall(TSMonitorEnvelopeCtxCreate(ts, &ctx));
386: PetscCall(TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscCtxDestroyFn *)TSMonitorEnvelopeCtxDestroy));
387: }
388: flg = PETSC_FALSE;
389: PetscCall(PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt));
390: if (opt && flg) PetscCall(TSMonitorCancel(ts));
392: flg = PETSC_FALSE;
393: PetscCall(PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeIJacobianDefaultColor", flg, &flg, NULL));
394: if (flg) {
395: DM dm;
397: PetscCall(TSGetDM(ts, &dm));
398: PetscCall(DMTSUnsetIJacobianContext_Internal(dm));
399: PetscCall(TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL));
400: PetscCall(PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n"));
401: }
403: /* Handle specific TS options */
404: PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);
406: /* Handle TSAdapt options */
407: PetscCall(TSGetAdapt(ts, &ts->adapt));
408: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
409: PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject));
411: /* TS trajectory must be set after TS, since it may use some TS options above */
412: tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
413: PetscCall(PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL));
414: if (tflg) PetscCall(TSSetSaveTrajectory(ts));
416: PetscCall(TSAdjointSetFromOptions(ts, PetscOptionsObject));
418: /* process any options handlers added with PetscObjectAddOptionsHandler() */
419: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject));
420: PetscOptionsEnd();
422: if (ts->trajectory) PetscCall(TSTrajectorySetFromOptions(ts->trajectory, ts));
424: /* why do we have to do this here and not during TSSetUp? */
425: PetscCall(TSGetSNES(ts, &ts->snes));
426: if (ts->problem_type == TS_LINEAR) {
427: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
428: if (!flg) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
429: }
430: PetscCall(SNESSetFromOptions(ts->snes));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*@
435: TSGetTrajectory - Gets the trajectory from a `TS` if it exists
437: Collective
439: Input Parameter:
440: . ts - the `TS` context obtained from `TSCreate()`
442: Output Parameter:
443: . tr - the `TSTrajectory` object, if it exists
445: Level: advanced
447: Note:
448: This routine should be called after all `TS` options have been set
450: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectoryCreate()`
451: @*/
452: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
453: {
454: PetscFunctionBegin;
456: *tr = ts->trajectory;
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /*@
461: TSSetSaveTrajectory - Causes the `TS` to save its solutions as it iterates forward in time in a `TSTrajectory` object
463: Collective
465: Input Parameter:
466: . ts - the `TS` context obtained from `TSCreate()`
468: Options Database Keys:
469: + -ts_save_trajectory - saves the trajectory to a file
470: - -ts_trajectory_type type - set trajectory type
472: Level: intermediate
474: Notes:
475: This routine should be called after all `TS` options have been set
477: The `TSTRAJECTORYVISUALIZATION` files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
478: MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m
480: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`
481: @*/
482: PetscErrorCode TSSetSaveTrajectory(TS ts)
483: {
484: PetscFunctionBegin;
486: if (!ts->trajectory) PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: /*@
491: TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object
493: Collective
495: Input Parameter:
496: . ts - the `TS` context obtained from `TSCreate()`
498: Level: intermediate
500: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
501: @*/
502: PetscErrorCode TSResetTrajectory(TS ts)
503: {
504: PetscFunctionBegin;
506: if (ts->trajectory) {
507: PetscCall(TSTrajectoryDestroy(&ts->trajectory));
508: PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
509: }
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: /*@
514: TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from a `TS`
516: Collective
518: Input Parameter:
519: . ts - the `TS` context obtained from `TSCreate()`
521: Level: intermediate
523: .seealso: [](ch_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
524: @*/
525: PetscErrorCode TSRemoveTrajectory(TS ts)
526: {
527: PetscFunctionBegin;
529: if (ts->trajectory) PetscCall(TSTrajectoryDestroy(&ts->trajectory));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /*@
534: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
535: set with `TSSetRHSJacobian()`.
537: Collective
539: Input Parameters:
540: + ts - the `TS` context
541: . t - current timestep
542: - U - input vector
544: Output Parameters:
545: + A - Jacobian matrix
546: - B - optional preconditioning matrix
548: Level: developer
550: Note:
551: Most users should not need to explicitly call this routine, as it
552: is used internally within the nonlinear solvers.
554: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()`
555: @*/
556: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
557: {
558: PetscObjectState Ustate;
559: PetscObjectId Uid;
560: DM dm;
561: DMTS tsdm;
562: TSRHSJacobianFn *rhsjacobianfunc;
563: void *ctx;
564: TSRHSFunctionFn *rhsfunction;
566: PetscFunctionBegin;
569: PetscCheckSameComm(ts, 1, U, 3);
570: PetscCall(TSGetDM(ts, &dm));
571: PetscCall(DMGetDMTS(dm, &tsdm));
572: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
573: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx));
574: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
575: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
577: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) PetscFunctionReturn(PETSC_SUCCESS);
579: 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);
580: if (rhsjacobianfunc) {
581: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
582: PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
583: ts->rhsjacs++;
584: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
585: } else {
586: PetscCall(MatZeroEntries(A));
587: if (B && A != B) PetscCall(MatZeroEntries(B));
588: }
589: ts->rhsjacobian.time = t;
590: ts->rhsjacobian.shift = 0;
591: ts->rhsjacobian.scale = 1.;
592: PetscCall(PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid));
593: PetscCall(PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: /*@
598: TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS`
600: Collective
602: Input Parameters:
603: + ts - the `TS` context
604: . t - current time
605: - U - state vector
607: Output Parameter:
608: . y - right-hand side
610: Level: developer
612: Note:
613: Most users should not need to explicitly call this routine, as it
614: is used internally within the nonlinear solvers.
616: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
617: @*/
618: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
619: {
620: TSRHSFunctionFn *rhsfunction;
621: TSIFunctionFn *ifunction;
622: void *ctx;
623: DM dm;
625: PetscFunctionBegin;
629: PetscCall(TSGetDM(ts, &dm));
630: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx));
631: PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
633: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
635: if (rhsfunction) {
636: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, y, 0));
637: PetscCall(VecLockReadPush(U));
638: PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
639: PetscCall(VecLockReadPop(U));
640: ts->rhsfuncs++;
641: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, y, 0));
642: } else PetscCall(VecZeroEntries(y));
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: /*@
647: TSComputeSolutionFunction - Evaluates the solution function.
649: Collective
651: Input Parameters:
652: + ts - the `TS` context
653: - t - current time
655: Output Parameter:
656: . U - the solution
658: Level: developer
660: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
661: @*/
662: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
663: {
664: TSSolutionFn *solutionfunction;
665: void *ctx;
666: DM dm;
668: PetscFunctionBegin;
671: PetscCall(TSGetDM(ts, &dm));
672: PetscCall(DMTSGetSolutionFunction(dm, &solutionfunction, &ctx));
673: if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
676: /*@
677: TSComputeForcingFunction - Evaluates the forcing function.
679: Collective
681: Input Parameters:
682: + ts - the `TS` context
683: - t - current time
685: Output Parameter:
686: . U - the function value
688: Level: developer
690: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
691: @*/
692: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
693: {
694: void *ctx;
695: DM dm;
696: TSForcingFn *forcing;
698: PetscFunctionBegin;
701: PetscCall(TSGetDM(ts, &dm));
702: PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));
704: if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
709: {
710: Mat A, B;
711: TSIJacobianFn *ijacobian;
713: PetscFunctionBegin;
714: if (Arhs) *Arhs = NULL;
715: if (Brhs) *Brhs = NULL;
716: PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL));
717: if (Arhs) {
718: if (!ts->Arhs) {
719: if (ijacobian) {
720: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs));
721: PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN));
722: } else {
723: ts->Arhs = A;
724: PetscCall(PetscObjectReference((PetscObject)A));
725: }
726: } else {
727: PetscBool flg;
728: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
729: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
730: if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
731: PetscCall(PetscObjectDereference((PetscObject)ts->Arhs));
732: ts->Arhs = A;
733: PetscCall(PetscObjectReference((PetscObject)A));
734: }
735: }
736: *Arhs = ts->Arhs;
737: }
738: if (Brhs) {
739: if (!ts->Brhs) {
740: if (A != B) {
741: if (ijacobian) {
742: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs));
743: } else {
744: ts->Brhs = B;
745: PetscCall(PetscObjectReference((PetscObject)B));
746: }
747: } else {
748: PetscCall(PetscObjectReference((PetscObject)ts->Arhs));
749: ts->Brhs = ts->Arhs;
750: }
751: }
752: *Brhs = ts->Brhs;
753: }
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: /*@
758: TSComputeIFunction - Evaluates the DAE residual written in the implicit form F(t,U,Udot)=0
760: Collective
762: Input Parameters:
763: + ts - the `TS` context
764: . t - current time
765: . U - state vector
766: . Udot - time derivative of state vector
767: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSFunction should be kept separate
769: Output Parameter:
770: . Y - right-hand side
772: Level: developer
774: Note:
775: Most users should not need to explicitly call this routine, as it
776: is used internally within the nonlinear solvers.
778: If the user did not write their equations in implicit form, this
779: function recasts them in implicit form.
781: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
782: @*/
783: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
784: {
785: TSIFunctionFn *ifunction;
786: TSRHSFunctionFn *rhsfunction;
787: void *ctx;
788: DM dm;
790: PetscFunctionBegin;
796: PetscCall(TSGetDM(ts, &dm));
797: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
798: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
800: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
802: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, Udot, Y));
803: if (ifunction) {
804: PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
805: ts->ifuncs++;
806: }
807: if (imex) {
808: if (!ifunction) PetscCall(VecCopy(Udot, Y));
809: } else if (rhsfunction) {
810: if (ifunction) {
811: Vec Frhs;
813: PetscCall(DMGetGlobalVector(dm, &Frhs));
814: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
815: PetscCall(VecAXPY(Y, -1, Frhs));
816: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
817: } else {
818: PetscCall(TSComputeRHSFunction(ts, t, U, Y));
819: PetscCall(VecAYPX(Y, -1, Udot));
820: }
821: }
822: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, Udot, Y));
823: PetscFunctionReturn(PETSC_SUCCESS);
824: }
826: /*
827: TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call `TSComputeRHSJacobian()` on it.
829: Note:
830: This routine is needed when one switches from `TSComputeIJacobian()` to `TSComputeRHSJacobian()` because the Jacobian matrix may be shifted or scaled in `TSComputeIJacobian()`.
832: */
833: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
834: {
835: PetscFunctionBegin;
837: PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
838: PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");
840: if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
841: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
842: if (B && B == ts->Brhs && A != B) {
843: if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
844: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
845: }
846: ts->rhsjacobian.shift = 0;
847: ts->rhsjacobian.scale = 1.;
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: /*@
852: TSComputeIJacobian - Evaluates the Jacobian of the DAE
854: Collective
856: Input Parameters:
857: + ts - the `TS` context
858: . t - current timestep
859: . U - state vector
860: . Udot - time derivative of state vector
861: . shift - shift to apply, see note below
862: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSJacobian should be kept separate
864: Output Parameters:
865: + A - Jacobian matrix
866: - B - matrix from which the preconditioner is constructed; often the same as `A`
868: Level: developer
870: Notes:
871: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
872: .vb
873: dF/dU + shift*dF/dUdot
874: .ve
875: Most users should not need to explicitly call this routine, as it
876: is used internally within the nonlinear solvers.
878: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
879: @*/
880: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
881: {
882: TSIJacobianFn *ijacobian;
883: TSRHSJacobianFn *rhsjacobian;
884: DM dm;
885: void *ctx;
887: PetscFunctionBegin;
894: PetscCall(TSGetDM(ts, &dm));
895: PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
896: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
898: PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
900: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
901: if (ijacobian) {
902: PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
903: ts->ijacs++;
904: }
905: if (imex) {
906: if (!ijacobian) { /* system was written as Udot = G(t,U) */
907: PetscBool assembled;
908: if (rhsjacobian) {
909: Mat Arhs = NULL;
910: PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
911: if (A == Arhs) {
912: 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 */
913: ts->rhsjacobian.time = PETSC_MIN_REAL;
914: }
915: }
916: PetscCall(MatZeroEntries(A));
917: PetscCall(MatAssembled(A, &assembled));
918: if (!assembled) {
919: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
920: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
921: }
922: PetscCall(MatShift(A, shift));
923: if (A != B) {
924: PetscCall(MatZeroEntries(B));
925: PetscCall(MatAssembled(B, &assembled));
926: if (!assembled) {
927: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
928: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
929: }
930: PetscCall(MatShift(B, shift));
931: }
932: }
933: } else {
934: Mat Arhs = NULL, Brhs = NULL;
936: /* RHSJacobian needs to be converted to part of IJacobian if exists */
937: if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
938: if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
939: PetscObjectState Ustate;
940: PetscObjectId Uid;
941: TSRHSFunctionFn *rhsfunction;
943: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
944: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
945: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
946: if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
947: ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */
948: PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
949: if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
950: } else {
951: PetscBool flg;
953: if (ts->rhsjacobian.reuse) { /* Undo the damage */
954: /* MatScale has a short path for this case.
955: However, this code path is taken the first time TSComputeRHSJacobian is called
956: and the matrices have not been assembled yet */
957: PetscCall(TSRecoverRHSJacobian(ts, A, B));
958: }
959: PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
960: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
961: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
962: if (!flg) {
963: PetscCall(MatScale(A, -1));
964: PetscCall(MatShift(A, shift));
965: }
966: if (A != B) {
967: PetscCall(MatScale(B, -1));
968: PetscCall(MatShift(B, shift));
969: }
970: }
971: ts->rhsjacobian.scale = -1;
972: ts->rhsjacobian.shift = shift;
973: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
974: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
975: PetscCall(MatZeroEntries(A));
976: PetscCall(MatShift(A, shift));
977: if (A != B) {
978: PetscCall(MatZeroEntries(B));
979: PetscCall(MatShift(B, shift));
980: }
981: }
982: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
983: PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
984: if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
985: }
986: }
987: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
988: PetscFunctionReturn(PETSC_SUCCESS);
989: }
991: /*@C
992: TSSetRHSFunction - Sets the routine for evaluating the function,
993: where U_t = G(t,u).
995: Logically Collective
997: Input Parameters:
998: + ts - the `TS` context obtained from `TSCreate()`
999: . r - vector to put the computed right-hand side (or `NULL` to have it created)
1000: . f - routine for evaluating the right-hand-side function
1001: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
1003: Level: beginner
1005: Note:
1006: You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.
1008: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1009: @*/
1010: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, void *ctx)
1011: {
1012: SNES snes;
1013: Vec ralloc = NULL;
1014: DM dm;
1016: PetscFunctionBegin;
1020: PetscCall(TSGetDM(ts, &dm));
1021: PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1022: PetscCall(TSGetSNES(ts, &snes));
1023: if (!r && !ts->dm && ts->vec_sol) {
1024: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1025: r = ralloc;
1026: }
1027: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1028: PetscCall(VecDestroy(&ralloc));
1029: PetscFunctionReturn(PETSC_SUCCESS);
1030: }
1032: /*@C
1033: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1035: Logically Collective
1037: Input Parameters:
1038: + ts - the `TS` context obtained from `TSCreate()`
1039: . f - routine for evaluating the solution
1040: - ctx - [optional] user-defined context for private data for the
1041: function evaluation routine (may be `NULL`)
1043: Options Database Keys:
1044: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1045: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
1047: Level: intermediate
1049: Notes:
1050: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1051: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1052: create closed-form solutions with non-physical forcing terms.
1054: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1056: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1057: @*/
1058: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, void *ctx)
1059: {
1060: DM dm;
1062: PetscFunctionBegin;
1064: PetscCall(TSGetDM(ts, &dm));
1065: PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }
1069: /*@C
1070: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1072: Logically Collective
1074: Input Parameters:
1075: + ts - the `TS` context obtained from `TSCreate()`
1076: . func - routine for evaluating the forcing function
1077: - ctx - [optional] user-defined context for private data for the function evaluation routine
1078: (may be `NULL`)
1080: Level: intermediate
1082: Notes:
1083: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1084: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1085: definition of the problem you are solving and hence possibly introducing bugs.
1087: This replaces the ODE F(u,u_t,t) = 0 the `TS` is solving with F(u,u_t,t) - func(t) = 0
1089: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1090: parameters can be passed in the ctx variable.
1092: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1094: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
1095: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1096: @*/
1097: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, void *ctx)
1098: {
1099: DM dm;
1101: PetscFunctionBegin;
1103: PetscCall(TSGetDM(ts, &dm));
1104: PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1108: /*@C
1109: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1110: where U_t = G(U,t), as well as the location to store the matrix.
1112: Logically Collective
1114: Input Parameters:
1115: + ts - the `TS` context obtained from `TSCreate()`
1116: . Amat - (approximate) location to store Jacobian matrix entries computed by `f`
1117: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1118: . f - the Jacobian evaluation routine
1119: - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1121: Level: beginner
1123: Notes:
1124: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1126: The `TS` solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f()`
1127: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1129: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
1130: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
1131: @*/
1132: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, void *ctx)
1133: {
1134: SNES snes;
1135: DM dm;
1136: TSIJacobianFn *ijacobian;
1138: PetscFunctionBegin;
1142: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1143: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1145: PetscCall(TSGetDM(ts, &dm));
1146: PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1147: PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1148: PetscCall(TSGetSNES(ts, &snes));
1149: if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1150: if (Amat) {
1151: PetscCall(PetscObjectReference((PetscObject)Amat));
1152: PetscCall(MatDestroy(&ts->Arhs));
1153: ts->Arhs = Amat;
1154: }
1155: if (Pmat) {
1156: PetscCall(PetscObjectReference((PetscObject)Pmat));
1157: PetscCall(MatDestroy(&ts->Brhs));
1158: ts->Brhs = Pmat;
1159: }
1160: PetscFunctionReturn(PETSC_SUCCESS);
1161: }
1163: /*@C
1164: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1166: Logically Collective
1168: Input Parameters:
1169: + ts - the `TS` context obtained from `TSCreate()`
1170: . r - vector to hold the residual (or `NULL` to have it created internally)
1171: . f - the function evaluation routine
1172: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1174: Level: beginner
1176: Note:
1177: The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function.
1179: .seealso: [](ch_ts), `TS`, `TSIFunctionFn`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`,
1180: `TSSetIJacobian()`
1181: @*/
1182: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunctionFn *f, void *ctx)
1183: {
1184: SNES snes;
1185: Vec ralloc = NULL;
1186: DM dm;
1188: PetscFunctionBegin;
1192: PetscCall(TSGetDM(ts, &dm));
1193: PetscCall(DMTSSetIFunction(dm, f, ctx));
1195: PetscCall(TSGetSNES(ts, &snes));
1196: if (!r && !ts->dm && ts->vec_sol) {
1197: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1198: r = ralloc;
1199: }
1200: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1201: PetscCall(VecDestroy(&ralloc));
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: /*@C
1206: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.
1208: Not Collective
1210: Input Parameter:
1211: . ts - the `TS` context
1213: Output Parameters:
1214: + r - vector to hold residual (or `NULL`)
1215: . func - the function to compute residual (or `NULL`)
1216: - ctx - the function context (or `NULL`)
1218: Level: advanced
1220: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1221: @*/
1222: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunctionFn **func, void **ctx)
1223: {
1224: SNES snes;
1225: DM dm;
1227: PetscFunctionBegin;
1229: PetscCall(TSGetSNES(ts, &snes));
1230: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1231: PetscCall(TSGetDM(ts, &dm));
1232: PetscCall(DMTSGetIFunction(dm, func, ctx));
1233: PetscFunctionReturn(PETSC_SUCCESS);
1234: }
1236: /*@C
1237: TSGetRHSFunction - Returns the vector where the right-hand side is stored and the function/context to compute it.
1239: Not Collective
1241: Input Parameter:
1242: . ts - the `TS` context
1244: Output Parameters:
1245: + r - vector to hold computed right-hand side (or `NULL`)
1246: . func - the function to compute right-hand side (or `NULL`)
1247: - ctx - the function context (or `NULL`)
1249: Level: advanced
1251: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1252: @*/
1253: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunctionFn **func, void **ctx)
1254: {
1255: SNES snes;
1256: DM dm;
1258: PetscFunctionBegin;
1260: PetscCall(TSGetSNES(ts, &snes));
1261: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1262: PetscCall(TSGetDM(ts, &dm));
1263: PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1264: PetscFunctionReturn(PETSC_SUCCESS);
1265: }
1267: /*@C
1268: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1269: provided with `TSSetIFunction()`.
1271: Logically Collective
1273: Input Parameters:
1274: + ts - the `TS` context obtained from `TSCreate()`
1275: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
1276: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1277: . f - the Jacobian evaluation routine
1278: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1280: Level: beginner
1282: Notes:
1283: The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1285: If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
1286: space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.
1288: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1289: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1290: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1291: a and vector W depend on the integration method, step size, and past states. For example with
1292: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1293: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1295: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1297: The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f`
1298: You should not assume the values are the same in the next call to `f` as you set them in the previous call.
1300: In case `TSSetRHSJacobian()` is also used in conjunction with a fully-implicit solver,
1301: multilevel linear solvers, e.g. `PCMG`, will likely not work due to the way `TS` handles rhs matrices.
1303: .seealso: [](ch_ts), `TS`, `TSIJacobianFn`, `TSSetIFunction()`, `TSSetRHSJacobian()`,
1304: `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1305: @*/
1306: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobianFn *f, void *ctx)
1307: {
1308: SNES snes;
1309: DM dm;
1311: PetscFunctionBegin;
1315: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1316: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1318: PetscCall(TSGetDM(ts, &dm));
1319: PetscCall(DMTSSetIJacobian(dm, f, ctx));
1321: PetscCall(TSGetSNES(ts, &snes));
1322: PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1323: PetscFunctionReturn(PETSC_SUCCESS);
1324: }
1326: /*@
1327: TSRHSJacobianSetReuse - restore the RHS Jacobian before calling the user-provided `TSRHSJacobianFn` function again
1329: Logically Collective
1331: Input Parameters:
1332: + ts - `TS` context obtained from `TSCreate()`
1333: - reuse - `PETSC_TRUE` if the RHS Jacobian
1335: Level: intermediate
1337: Notes:
1338: Without this flag, `TS` will change the sign and shift the RHS Jacobian for a
1339: finite-time-step implicit solve, in which case the user function will need to recompute the
1340: entire Jacobian. The `reuse `flag must be set if the evaluation function assumes that the
1341: matrix entries have not been changed by the `TS`.
1343: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1344: @*/
1345: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1346: {
1347: PetscFunctionBegin;
1348: ts->rhsjacobian.reuse = reuse;
1349: PetscFunctionReturn(PETSC_SUCCESS);
1350: }
1352: /*@C
1353: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1355: Logically Collective
1357: Input Parameters:
1358: + ts - the `TS` context obtained from `TSCreate()`
1359: . F - vector to hold the residual (or `NULL` to have it created internally)
1360: . fun - the function evaluation routine
1361: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1363: Level: beginner
1365: .seealso: [](ch_ts), `TS`, `TSI2FunctionFn`, `TSSetI2Jacobian()`, `TSSetIFunction()`,
1366: `TSCreate()`, `TSSetRHSFunction()`
1367: @*/
1368: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2FunctionFn *fun, void *ctx)
1369: {
1370: DM dm;
1372: PetscFunctionBegin;
1375: PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1376: PetscCall(TSGetDM(ts, &dm));
1377: PetscCall(DMTSSetI2Function(dm, fun, ctx));
1378: PetscFunctionReturn(PETSC_SUCCESS);
1379: }
1381: /*@C
1382: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.
1384: Not Collective
1386: Input Parameter:
1387: . ts - the `TS` context
1389: Output Parameters:
1390: + r - vector to hold residual (or `NULL`)
1391: . fun - the function to compute residual (or `NULL`)
1392: - ctx - the function context (or `NULL`)
1394: Level: advanced
1396: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1397: @*/
1398: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2FunctionFn **fun, void **ctx)
1399: {
1400: SNES snes;
1401: DM dm;
1403: PetscFunctionBegin;
1405: PetscCall(TSGetSNES(ts, &snes));
1406: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1407: PetscCall(TSGetDM(ts, &dm));
1408: PetscCall(DMTSGetI2Function(dm, fun, ctx));
1409: PetscFunctionReturn(PETSC_SUCCESS);
1410: }
1412: /*@C
1413: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1414: where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.
1416: Logically Collective
1418: Input Parameters:
1419: + ts - the `TS` context obtained from `TSCreate()`
1420: . J - matrix to hold the Jacobian values
1421: . P - matrix for constructing the preconditioner (may be same as `J`)
1422: . jac - the Jacobian evaluation routine, see `TSI2JacobianFn` for the calling sequence
1423: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1425: Level: beginner
1427: Notes:
1428: The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1430: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1431: 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.
1432: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1433: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1435: .seealso: [](ch_ts), `TS`, `TSI2JacobianFn`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1436: @*/
1437: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2JacobianFn *jac, void *ctx)
1438: {
1439: DM dm;
1441: PetscFunctionBegin;
1445: PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1446: PetscCall(TSGetDM(ts, &dm));
1447: PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1448: PetscFunctionReturn(PETSC_SUCCESS);
1449: }
1451: /*@C
1452: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1454: Not Collective, but parallel objects are returned if `TS` is parallel
1456: Input Parameter:
1457: . ts - The `TS` context obtained from `TSCreate()`
1459: Output Parameters:
1460: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1461: . P - The matrix from which the preconditioner is constructed, often the same as `J`
1462: . jac - The function to compute the Jacobian matrices
1463: - ctx - User-defined context for Jacobian evaluation routine
1465: Level: advanced
1467: Note:
1468: You can pass in `NULL` for any return argument you do not need.
1470: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1471: @*/
1472: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2JacobianFn **jac, void **ctx)
1473: {
1474: SNES snes;
1475: DM dm;
1477: PetscFunctionBegin;
1478: PetscCall(TSGetSNES(ts, &snes));
1479: PetscCall(SNESSetUpMatrices(snes));
1480: PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1481: PetscCall(TSGetDM(ts, &dm));
1482: PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1483: PetscFunctionReturn(PETSC_SUCCESS);
1484: }
1486: /*@
1487: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1489: Collective
1491: Input Parameters:
1492: + ts - the `TS` context
1493: . t - current time
1494: . U - state vector
1495: . V - time derivative of state vector (U_t)
1496: - A - second time derivative of state vector (U_tt)
1498: Output Parameter:
1499: . F - the residual vector
1501: Level: developer
1503: Note:
1504: Most users should not need to explicitly call this routine, as it
1505: is used internally within the nonlinear solvers.
1507: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1508: @*/
1509: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1510: {
1511: DM dm;
1512: TSI2FunctionFn *I2Function;
1513: void *ctx;
1514: TSRHSFunctionFn *rhsfunction;
1516: PetscFunctionBegin;
1523: PetscCall(TSGetDM(ts, &dm));
1524: PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1525: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
1527: if (!I2Function) {
1528: PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1529: PetscFunctionReturn(PETSC_SUCCESS);
1530: }
1532: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, V, F));
1534: PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));
1536: if (rhsfunction) {
1537: Vec Frhs;
1539: PetscCall(DMGetGlobalVector(dm, &Frhs));
1540: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1541: PetscCall(VecAXPY(F, -1, Frhs));
1542: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
1543: }
1545: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, V, F));
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1549: /*@
1550: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1552: Collective
1554: Input Parameters:
1555: + ts - the `TS` context
1556: . t - current timestep
1557: . U - state vector
1558: . V - time derivative of state vector
1559: . A - second time derivative of state vector
1560: . shiftV - shift to apply, see note below
1561: - shiftA - shift to apply, see note below
1563: Output Parameters:
1564: + J - Jacobian matrix
1565: - P - optional preconditioning matrix
1567: Level: developer
1569: Notes:
1570: If F(t,U,V,A)=0 is the DAE, the required Jacobian is
1572: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1574: Most users should not need to explicitly call this routine, as it
1575: is used internally within the nonlinear solvers.
1577: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1578: @*/
1579: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1580: {
1581: DM dm;
1582: TSI2JacobianFn *I2Jacobian;
1583: void *ctx;
1584: TSRHSJacobianFn *rhsjacobian;
1586: PetscFunctionBegin;
1594: PetscCall(TSGetDM(ts, &dm));
1595: PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1596: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1598: if (!I2Jacobian) {
1599: PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1600: PetscFunctionReturn(PETSC_SUCCESS);
1601: }
1603: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, J, P));
1604: PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1605: if (rhsjacobian) {
1606: Mat Jrhs, Prhs;
1607: PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1608: PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1609: PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1610: if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1611: }
1613: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, J, P));
1614: PetscFunctionReturn(PETSC_SUCCESS);
1615: }
1617: /*@C
1618: TSSetTransientVariable - sets function to transform from state to transient variables
1620: Logically Collective
1622: Input Parameters:
1623: + ts - time stepping context on which to change the transient variable
1624: . tvar - a function that transforms to transient variables, see `TSTransientVariableFn` for the calling sequence
1625: - ctx - a context for tvar
1627: Level: advanced
1629: Notes:
1630: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1631: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1632: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1633: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1634: evaluated via the chain rule, as in
1635: .vb
1636: dF/dP + shift * dF/dCdot dC/dP.
1637: .ve
1639: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSTransientVariableFn`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1640: @*/
1641: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariableFn *tvar, void *ctx)
1642: {
1643: DM dm;
1645: PetscFunctionBegin;
1647: PetscCall(TSGetDM(ts, &dm));
1648: PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1649: PetscFunctionReturn(PETSC_SUCCESS);
1650: }
1652: /*@
1653: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1655: Logically Collective
1657: Input Parameters:
1658: + ts - TS on which to compute
1659: - U - state vector to be transformed to transient variables
1661: Output Parameter:
1662: . C - transient (conservative) variable
1664: Level: developer
1666: Developer Notes:
1667: If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed.
1668: This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are
1669: being used.
1671: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1672: @*/
1673: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1674: {
1675: DM dm;
1676: DMTS dmts;
1678: PetscFunctionBegin;
1681: PetscCall(TSGetDM(ts, &dm));
1682: PetscCall(DMGetDMTS(dm, &dmts));
1683: if (dmts->ops->transientvar) {
1685: PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1686: }
1687: PetscFunctionReturn(PETSC_SUCCESS);
1688: }
1690: /*@
1691: TSHasTransientVariable - determine whether transient variables have been set
1693: Logically Collective
1695: Input Parameter:
1696: . ts - `TS` on which to compute
1698: Output Parameter:
1699: . has - `PETSC_TRUE` if transient variables have been set
1701: Level: developer
1703: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1704: @*/
1705: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1706: {
1707: DM dm;
1708: DMTS dmts;
1710: PetscFunctionBegin;
1712: PetscCall(TSGetDM(ts, &dm));
1713: PetscCall(DMGetDMTS(dm, &dmts));
1714: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1715: PetscFunctionReturn(PETSC_SUCCESS);
1716: }
1718: /*@
1719: TS2SetSolution - Sets the initial solution and time derivative vectors
1720: for use by the `TS` routines handling second order equations.
1722: Logically Collective
1724: Input Parameters:
1725: + ts - the `TS` context obtained from `TSCreate()`
1726: . u - the solution vector
1727: - v - the time derivative vector
1729: Level: beginner
1731: .seealso: [](ch_ts), `TS`
1732: @*/
1733: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1734: {
1735: PetscFunctionBegin;
1739: PetscCall(TSSetSolution(ts, u));
1740: PetscCall(PetscObjectReference((PetscObject)v));
1741: PetscCall(VecDestroy(&ts->vec_dot));
1742: ts->vec_dot = v;
1743: PetscFunctionReturn(PETSC_SUCCESS);
1744: }
1746: /*@
1747: TS2GetSolution - Returns the solution and time derivative at the present timestep
1748: for second order equations.
1750: Not Collective
1752: Input Parameter:
1753: . ts - the `TS` context obtained from `TSCreate()`
1755: Output Parameters:
1756: + u - the vector containing the solution
1757: - v - the vector containing the time derivative
1759: Level: intermediate
1761: Notes:
1762: It is valid to call this routine inside the function
1763: that you are evaluating in order to move to the new timestep. This vector not
1764: changed until the solution at the next timestep has been calculated.
1766: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1767: @*/
1768: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1769: {
1770: PetscFunctionBegin;
1772: if (u) PetscAssertPointer(u, 2);
1773: if (v) PetscAssertPointer(v, 3);
1774: if (u) *u = ts->vec_sol;
1775: if (v) *v = ts->vec_dot;
1776: PetscFunctionReturn(PETSC_SUCCESS);
1777: }
1779: /*@
1780: TSLoad - Loads a `TS` that has been stored in binary with `TSView()`.
1782: Collective
1784: Input Parameters:
1785: + ts - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1786: some related function before a call to `TSLoad()`.
1787: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1789: Level: intermediate
1791: Note:
1792: The type is determined by the data in the file, any type set into the `TS` before this call is ignored.
1794: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1795: @*/
1796: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1797: {
1798: PetscBool isbinary;
1799: PetscInt classid;
1800: char type[256];
1801: DMTS sdm;
1802: DM dm;
1804: PetscFunctionBegin;
1807: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1808: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1810: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1811: PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1812: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1813: PetscCall(TSSetType(ts, type));
1814: PetscTryTypeMethod(ts, load, viewer);
1815: PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1816: PetscCall(DMLoad(dm, viewer));
1817: PetscCall(TSSetDM(ts, dm));
1818: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1819: PetscCall(VecLoad(ts->vec_sol, viewer));
1820: PetscCall(DMGetDMTS(ts->dm, &sdm));
1821: PetscCall(DMTSLoad(sdm, viewer));
1822: PetscFunctionReturn(PETSC_SUCCESS);
1823: }
1825: #include <petscdraw.h>
1826: #if defined(PETSC_HAVE_SAWS)
1827: #include <petscviewersaws.h>
1828: #endif
1830: /*@
1831: TSViewFromOptions - View a `TS` based on values in the options database
1833: Collective
1835: Input Parameters:
1836: + ts - the `TS` context
1837: . obj - Optional object that provides the prefix for the options database keys
1838: - name - command line option string to be passed by user
1840: Level: intermediate
1842: .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1843: @*/
1844: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1845: {
1846: PetscFunctionBegin;
1848: PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1849: PetscFunctionReturn(PETSC_SUCCESS);
1850: }
1852: /*@
1853: TSView - Prints the `TS` data structure.
1855: Collective
1857: Input Parameters:
1858: + ts - the `TS` context obtained from `TSCreate()`
1859: - viewer - visualization context
1861: Options Database Key:
1862: . -ts_view - calls `TSView()` at end of `TSStep()`
1864: Level: beginner
1866: Notes:
1867: The available visualization contexts include
1868: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1869: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1870: output where only the first processor opens
1871: the file. All other processors send their
1872: data to the first processor to print.
1874: The user can open an alternative visualization context with
1875: `PetscViewerASCIIOpen()` - output to a specified file.
1877: In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).
1879: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1880: @*/
1881: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1882: {
1883: TSType type;
1884: PetscBool iascii, isstring, isundials, isbinary, isdraw;
1885: DMTS sdm;
1886: #if defined(PETSC_HAVE_SAWS)
1887: PetscBool issaws;
1888: #endif
1890: PetscFunctionBegin;
1892: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1894: PetscCheckSameComm(ts, 1, viewer, 2);
1896: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1897: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1898: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1899: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1900: #if defined(PETSC_HAVE_SAWS)
1901: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1902: #endif
1903: if (iascii) {
1904: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer));
1905: if (ts->ops->view) {
1906: PetscCall(PetscViewerASCIIPushTab(viewer));
1907: PetscUseTypeMethod(ts, view, viewer);
1908: PetscCall(PetscViewerASCIIPopTab(viewer));
1909: }
1910: if (ts->max_steps < PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps));
1911: if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time));
1912: if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs));
1913: if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs));
1914: if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs));
1915: if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs));
1916: if (ts->usessnes) {
1917: PetscBool lin;
1918: if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its));
1919: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its));
1920: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
1921: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures));
1922: }
1923: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject));
1924: if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, "));
1925: else PetscCall(PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol));
1926: if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of absolute error tolerances\n"));
1927: else PetscCall(PetscViewerASCIIPrintf(viewer, " using absolute error tolerance of %g\n", (double)ts->atol));
1928: PetscCall(PetscViewerASCIIPushTab(viewer));
1929: PetscCall(TSAdaptView(ts->adapt, viewer));
1930: PetscCall(PetscViewerASCIIPopTab(viewer));
1931: } else if (isstring) {
1932: PetscCall(TSGetType(ts, &type));
1933: PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type));
1934: PetscTryTypeMethod(ts, view, viewer);
1935: } else if (isbinary) {
1936: PetscInt classid = TS_FILE_CLASSID;
1937: MPI_Comm comm;
1938: PetscMPIInt rank;
1939: char type[256];
1941: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
1942: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1943: if (rank == 0) {
1944: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1945: PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
1946: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1947: }
1948: PetscTryTypeMethod(ts, view, viewer);
1949: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1950: PetscCall(DMView(ts->dm, viewer));
1951: PetscCall(VecView(ts->vec_sol, viewer));
1952: PetscCall(DMGetDMTS(ts->dm, &sdm));
1953: PetscCall(DMTSView(sdm, viewer));
1954: } else if (isdraw) {
1955: PetscDraw draw;
1956: char str[36];
1957: PetscReal x, y, bottom, h;
1959: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1960: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1961: PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
1962: PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
1963: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
1964: bottom = y - h;
1965: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1966: PetscTryTypeMethod(ts, view, viewer);
1967: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1968: if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
1969: PetscCall(PetscDrawPopCurrentPoint(draw));
1970: #if defined(PETSC_HAVE_SAWS)
1971: } else if (issaws) {
1972: PetscMPIInt rank;
1973: const char *name;
1975: PetscCall(PetscObjectGetName((PetscObject)ts, &name));
1976: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1977: if (!((PetscObject)ts)->amsmem && rank == 0) {
1978: char dir[1024];
1980: PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
1981: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
1982: PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
1983: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
1984: PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
1985: }
1986: PetscTryTypeMethod(ts, view, viewer);
1987: #endif
1988: }
1989: if (ts->snes && ts->usessnes) {
1990: PetscCall(PetscViewerASCIIPushTab(viewer));
1991: PetscCall(SNESView(ts->snes, viewer));
1992: PetscCall(PetscViewerASCIIPopTab(viewer));
1993: }
1994: PetscCall(DMGetDMTS(ts->dm, &sdm));
1995: PetscCall(DMTSView(sdm, viewer));
1997: PetscCall(PetscViewerASCIIPushTab(viewer));
1998: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials));
1999: PetscCall(PetscViewerASCIIPopTab(viewer));
2000: PetscFunctionReturn(PETSC_SUCCESS);
2001: }
2003: /*@
2004: TSSetApplicationContext - Sets an optional user-defined context for
2005: the timesteppers.
2007: Logically Collective
2009: Input Parameters:
2010: + ts - the `TS` context obtained from `TSCreate()`
2011: - ctx - user context
2013: Level: intermediate
2015: Fortran Notes:
2016: You must write a Fortran interface definition for this
2017: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
2019: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2020: @*/
2021: PetscErrorCode TSSetApplicationContext(TS ts, void *ctx)
2022: {
2023: PetscFunctionBegin;
2025: ts->ctx = ctx;
2026: PetscFunctionReturn(PETSC_SUCCESS);
2027: }
2029: /*@
2030: TSGetApplicationContext - Gets the user-defined context for the
2031: timestepper that was set with `TSSetApplicationContext()`
2033: Not Collective
2035: Input Parameter:
2036: . ts - the `TS` context obtained from `TSCreate()`
2038: Output Parameter:
2039: . ctx - user context
2041: Level: intermediate
2043: Fortran Notes:
2044: You must write a Fortran interface definition for this
2045: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
2047: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2048: @*/
2049: PetscErrorCode TSGetApplicationContext(TS ts, void *ctx)
2050: {
2051: PetscFunctionBegin;
2053: *(void **)ctx = ts->ctx;
2054: PetscFunctionReturn(PETSC_SUCCESS);
2055: }
2057: /*@
2058: TSGetStepNumber - Gets the number of time steps completed.
2060: Not Collective
2062: Input Parameter:
2063: . ts - the `TS` context obtained from `TSCreate()`
2065: Output Parameter:
2066: . steps - number of steps completed so far
2068: Level: intermediate
2070: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2071: @*/
2072: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2073: {
2074: PetscFunctionBegin;
2076: PetscAssertPointer(steps, 2);
2077: *steps = ts->steps;
2078: PetscFunctionReturn(PETSC_SUCCESS);
2079: }
2081: /*@
2082: TSSetStepNumber - Sets the number of steps completed.
2084: Logically Collective
2086: Input Parameters:
2087: + ts - the `TS` context
2088: - steps - number of steps completed so far
2090: Level: developer
2092: Note:
2093: For most uses of the `TS` solvers the user need not explicitly call
2094: `TSSetStepNumber()`, as the step counter is appropriately updated in
2095: `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2096: reinitialize timestepping by setting the step counter to zero (and time
2097: to the initial time) to solve a similar problem with different initial
2098: conditions or parameters. Other possible use case is to continue
2099: timestepping from a previously interrupted run in such a way that `TS`
2100: monitors will be called with a initial nonzero step counter.
2102: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2103: @*/
2104: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2105: {
2106: PetscFunctionBegin;
2109: PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2110: ts->steps = steps;
2111: PetscFunctionReturn(PETSC_SUCCESS);
2112: }
2114: /*@
2115: TSSetTimeStep - Allows one to reset the timestep at any time,
2116: useful for simple pseudo-timestepping codes.
2118: Logically Collective
2120: Input Parameters:
2121: + ts - the `TS` context obtained from `TSCreate()`
2122: - time_step - the size of the timestep
2124: Level: intermediate
2126: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2127: @*/
2128: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2129: {
2130: PetscFunctionBegin;
2133: ts->time_step = time_step;
2134: PetscFunctionReturn(PETSC_SUCCESS);
2135: }
2137: /*@
2138: TSSetExactFinalTime - Determines whether to adapt the final time step to
2139: match the exact final time, interpolate solution to the exact final time,
2140: or just return at the final time `TS` computed.
2142: Logically Collective
2144: Input Parameters:
2145: + ts - the time-step context
2146: - eftopt - exact final time option
2147: .vb
2148: TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded
2149: TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2150: TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2151: .ve
2153: Options Database Key:
2154: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step at runtime
2156: Level: beginner
2158: Note:
2159: If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2160: then the final time you selected.
2162: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2163: @*/
2164: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2165: {
2166: PetscFunctionBegin;
2169: ts->exact_final_time = eftopt;
2170: PetscFunctionReturn(PETSC_SUCCESS);
2171: }
2173: /*@
2174: TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`
2176: Not Collective
2178: Input Parameter:
2179: . ts - the `TS` context
2181: Output Parameter:
2182: . eftopt - exact final time option
2184: Level: beginner
2186: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2187: @*/
2188: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2189: {
2190: PetscFunctionBegin;
2192: PetscAssertPointer(eftopt, 2);
2193: *eftopt = ts->exact_final_time;
2194: PetscFunctionReturn(PETSC_SUCCESS);
2195: }
2197: /*@
2198: TSGetTimeStep - Gets the current timestep size.
2200: Not Collective
2202: Input Parameter:
2203: . ts - the `TS` context obtained from `TSCreate()`
2205: Output Parameter:
2206: . dt - the current timestep size
2208: Level: intermediate
2210: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2211: @*/
2212: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2213: {
2214: PetscFunctionBegin;
2216: PetscAssertPointer(dt, 2);
2217: *dt = ts->time_step;
2218: PetscFunctionReturn(PETSC_SUCCESS);
2219: }
2221: /*@
2222: TSGetSolution - Returns the solution at the present timestep. It
2223: is valid to call this routine inside the function that you are evaluating
2224: in order to move to the new timestep. This vector not changed until
2225: the solution at the next timestep has been calculated.
2227: Not Collective, but v returned is parallel if ts is parallel
2229: Input Parameter:
2230: . ts - the `TS` context obtained from `TSCreate()`
2232: Output Parameter:
2233: . v - the vector containing the solution
2235: Level: intermediate
2237: Note:
2238: If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
2239: final time. It returns the solution at the next timestep.
2241: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2242: @*/
2243: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2244: {
2245: PetscFunctionBegin;
2247: PetscAssertPointer(v, 2);
2248: *v = ts->vec_sol;
2249: PetscFunctionReturn(PETSC_SUCCESS);
2250: }
2252: /*@
2253: TSGetSolutionComponents - Returns any solution components at the present
2254: timestep, if available for the time integration method being used.
2255: Solution components are quantities that share the same size and
2256: structure as the solution vector.
2258: Not Collective, but v returned is parallel if ts is parallel
2260: Input Parameters:
2261: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2262: . n - If v is `NULL`, then the number of solution components is
2263: returned through n, else the n-th solution component is
2264: returned in v.
2265: - v - the vector containing the n-th solution component
2266: (may be `NULL` to use this function to find out
2267: the number of solutions components).
2269: Level: advanced
2271: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2272: @*/
2273: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2274: {
2275: PetscFunctionBegin;
2277: if (!ts->ops->getsolutioncomponents) *n = 0;
2278: else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2279: PetscFunctionReturn(PETSC_SUCCESS);
2280: }
2282: /*@
2283: TSGetAuxSolution - Returns an auxiliary solution at the present
2284: timestep, if available for the time integration method being used.
2286: Not Collective, but v returned is parallel if ts is parallel
2288: Input Parameters:
2289: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2290: - v - the vector containing the auxiliary solution
2292: Level: intermediate
2294: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2295: @*/
2296: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2297: {
2298: PetscFunctionBegin;
2300: if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2301: else PetscCall(VecZeroEntries(*v));
2302: PetscFunctionReturn(PETSC_SUCCESS);
2303: }
2305: /*@
2306: TSGetTimeError - Returns the estimated error vector, if the chosen
2307: `TSType` has an error estimation functionality and `TSSetTimeError()` was called
2309: Not Collective, but v returned is parallel if ts is parallel
2311: Input Parameters:
2312: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2313: . n - current estimate (n=0) or previous one (n=-1)
2314: - v - the vector containing the error (same size as the solution).
2316: Level: intermediate
2318: Note:
2319: MUST call after `TSSetUp()`
2321: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2322: @*/
2323: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2324: {
2325: PetscFunctionBegin;
2327: if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2328: else PetscCall(VecZeroEntries(*v));
2329: PetscFunctionReturn(PETSC_SUCCESS);
2330: }
2332: /*@
2333: TSSetTimeError - Sets the estimated error vector, if the chosen
2334: `TSType` has an error estimation functionality. This can be used
2335: to restart such a time integrator with a given error vector.
2337: Not Collective, but v returned is parallel if ts is parallel
2339: Input Parameters:
2340: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2341: - v - the vector containing the error (same size as the solution).
2343: Level: intermediate
2345: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
2346: @*/
2347: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2348: {
2349: PetscFunctionBegin;
2351: PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2352: PetscTryTypeMethod(ts, settimeerror, v);
2353: PetscFunctionReturn(PETSC_SUCCESS);
2354: }
2356: /* ----- Routines to initialize and destroy a timestepper ---- */
2357: /*@
2358: TSSetProblemType - Sets the type of problem to be solved.
2360: Not collective
2362: Input Parameters:
2363: + ts - The `TS`
2364: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2365: .vb
2366: U_t - A U = 0 (linear)
2367: U_t - A(t) U = 0 (linear)
2368: F(t,U,U_t) = 0 (nonlinear)
2369: .ve
2371: Level: beginner
2373: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2374: @*/
2375: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2376: {
2377: PetscFunctionBegin;
2379: ts->problem_type = type;
2380: if (type == TS_LINEAR) {
2381: SNES snes;
2382: PetscCall(TSGetSNES(ts, &snes));
2383: PetscCall(SNESSetType(snes, SNESKSPONLY));
2384: }
2385: PetscFunctionReturn(PETSC_SUCCESS);
2386: }
2388: /*@
2389: TSGetProblemType - Gets the type of problem to be solved.
2391: Not collective
2393: Input Parameter:
2394: . ts - The `TS`
2396: Output Parameter:
2397: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2398: .vb
2399: M U_t = A U
2400: M(t) U_t = A(t) U
2401: F(t,U,U_t)
2402: .ve
2404: Level: beginner
2406: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2407: @*/
2408: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2409: {
2410: PetscFunctionBegin;
2412: PetscAssertPointer(type, 2);
2413: *type = ts->problem_type;
2414: PetscFunctionReturn(PETSC_SUCCESS);
2415: }
2417: /*
2418: Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2419: */
2420: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2421: {
2422: PetscBool isnone;
2424: PetscFunctionBegin;
2425: PetscCall(TSGetAdapt(ts, &ts->adapt));
2426: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2428: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2429: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2430: else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2431: PetscFunctionReturn(PETSC_SUCCESS);
2432: }
2434: /*@
2435: TSSetUp - Sets up the internal data structures for the later use of a timestepper.
2437: Collective
2439: Input Parameter:
2440: . ts - the `TS` context obtained from `TSCreate()`
2442: Level: advanced
2444: Note:
2445: For basic use of the `TS` solvers the user need not explicitly call
2446: `TSSetUp()`, since these actions will automatically occur during
2447: the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this
2448: phase separately, `TSSetUp()` should be called after `TSCreate()`
2449: and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.
2451: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2452: @*/
2453: PetscErrorCode TSSetUp(TS ts)
2454: {
2455: DM dm;
2456: PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2457: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2458: TSIFunctionFn *ifun;
2459: TSIJacobianFn *ijac;
2460: TSI2JacobianFn *i2jac;
2461: TSRHSJacobianFn *rhsjac;
2463: PetscFunctionBegin;
2465: if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
2467: if (!((PetscObject)ts)->type_name) {
2468: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2469: PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2470: }
2472: if (!ts->vec_sol) {
2473: PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2474: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2475: }
2477: if (ts->tspan) {
2478: if (!ts->tspan->vecs_sol) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2479: }
2480: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2481: PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2482: ts->Jacp = ts->Jacprhs;
2483: }
2485: if (ts->quadraturets) {
2486: PetscCall(TSSetUp(ts->quadraturets));
2487: PetscCall(VecDestroy(&ts->vec_costintegrand));
2488: PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2489: }
2491: PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
2492: if (rhsjac == TSComputeRHSJacobianConstant) {
2493: Mat Amat, Pmat;
2494: SNES snes;
2495: PetscCall(TSGetSNES(ts, &snes));
2496: PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
2497: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2498: * have displaced the RHS matrix */
2499: if (Amat && Amat == ts->Arhs) {
2500: /* 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 */
2501: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
2502: PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
2503: PetscCall(MatDestroy(&Amat));
2504: }
2505: if (Pmat && Pmat == ts->Brhs) {
2506: PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
2507: PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
2508: PetscCall(MatDestroy(&Pmat));
2509: }
2510: }
2512: PetscCall(TSGetAdapt(ts, &ts->adapt));
2513: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2515: PetscTryTypeMethod(ts, setup);
2517: PetscCall(TSSetExactFinalTimeDefault(ts));
2519: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2520: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2521: */
2522: PetscCall(TSGetDM(ts, &dm));
2523: PetscCall(DMSNESGetFunction(dm, &func, NULL));
2524: if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));
2526: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2527: Otherwise, the SNES will use coloring internally to form the Jacobian.
2528: */
2529: PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2530: PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2531: PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2532: PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2533: if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));
2535: /* if time integration scheme has a starting method, call it */
2536: PetscTryTypeMethod(ts, startingmethod);
2538: ts->setupcalled = PETSC_TRUE;
2539: PetscFunctionReturn(PETSC_SUCCESS);
2540: }
2542: /*@
2543: TSReset - Resets a `TS` context and removes any allocated `Vec`s and `Mat`s.
2545: Collective
2547: Input Parameter:
2548: . ts - the `TS` context obtained from `TSCreate()`
2550: Level: beginner
2552: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()`
2553: @*/
2554: PetscErrorCode TSReset(TS ts)
2555: {
2556: TS_RHSSplitLink ilink = ts->tsrhssplit, next;
2558: PetscFunctionBegin;
2561: PetscTryTypeMethod(ts, reset);
2562: if (ts->snes) PetscCall(SNESReset(ts->snes));
2563: if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));
2565: PetscCall(MatDestroy(&ts->Arhs));
2566: PetscCall(MatDestroy(&ts->Brhs));
2567: PetscCall(VecDestroy(&ts->Frhs));
2568: PetscCall(VecDestroy(&ts->vec_sol));
2569: PetscCall(VecDestroy(&ts->vec_sol0));
2570: PetscCall(VecDestroy(&ts->vec_dot));
2571: PetscCall(VecDestroy(&ts->vatol));
2572: PetscCall(VecDestroy(&ts->vrtol));
2573: PetscCall(VecDestroyVecs(ts->nwork, &ts->work));
2575: PetscCall(MatDestroy(&ts->Jacprhs));
2576: PetscCall(MatDestroy(&ts->Jacp));
2577: if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2578: if (ts->quadraturets) {
2579: PetscCall(TSReset(ts->quadraturets));
2580: PetscCall(VecDestroy(&ts->vec_costintegrand));
2581: }
2582: while (ilink) {
2583: next = ilink->next;
2584: PetscCall(TSDestroy(&ilink->ts));
2585: PetscCall(PetscFree(ilink->splitname));
2586: PetscCall(ISDestroy(&ilink->is));
2587: PetscCall(PetscFree(ilink));
2588: ilink = next;
2589: }
2590: ts->tsrhssplit = NULL;
2591: ts->num_rhs_splits = 0;
2592: if (ts->tspan) {
2593: PetscCall(PetscFree(ts->tspan->span_times));
2594: PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2595: PetscCall(PetscFree(ts->tspan));
2596: }
2597: ts->rhsjacobian.time = PETSC_MIN_REAL;
2598: ts->rhsjacobian.scale = 1.0;
2599: ts->ijacobian.shift = 1.0;
2600: ts->setupcalled = PETSC_FALSE;
2601: PetscFunctionReturn(PETSC_SUCCESS);
2602: }
2604: static PetscErrorCode TSResizeReset(TS);
2606: /*@
2607: TSDestroy - Destroys the timestepper context that was created
2608: with `TSCreate()`.
2610: Collective
2612: Input Parameter:
2613: . ts - the `TS` context obtained from `TSCreate()`
2615: Level: beginner
2617: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2618: @*/
2619: PetscErrorCode TSDestroy(TS *ts)
2620: {
2621: PetscFunctionBegin;
2622: if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2624: if (--((PetscObject)*ts)->refct > 0) {
2625: *ts = NULL;
2626: PetscFunctionReturn(PETSC_SUCCESS);
2627: }
2629: PetscCall(TSReset(*ts));
2630: PetscCall(TSAdjointReset(*ts));
2631: if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2632: PetscCall(TSResizeReset(*ts));
2634: /* if memory was published with SAWs then destroy it */
2635: PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2636: PetscTryTypeMethod(*ts, destroy);
2638: PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));
2640: PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2641: PetscCall(TSEventDestroy(&(*ts)->event));
2643: PetscCall(SNESDestroy(&(*ts)->snes));
2644: PetscCall(SNESDestroy(&(*ts)->snesrhssplit));
2645: PetscCall(DMDestroy(&(*ts)->dm));
2646: PetscCall(TSMonitorCancel(*ts));
2647: PetscCall(TSAdjointMonitorCancel(*ts));
2649: PetscCall(TSDestroy(&(*ts)->quadraturets));
2650: PetscCall(PetscHeaderDestroy(ts));
2651: PetscFunctionReturn(PETSC_SUCCESS);
2652: }
2654: /*@
2655: TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2656: a `TS` (timestepper) context. Valid only for nonlinear problems.
2658: Not Collective, but snes is parallel if ts is parallel
2660: Input Parameter:
2661: . ts - the `TS` context obtained from `TSCreate()`
2663: Output Parameter:
2664: . snes - the nonlinear solver context
2666: Level: beginner
2668: Notes:
2669: The user can then directly manipulate the `SNES` context to set various
2670: options, etc. Likewise, the user can then extract and manipulate the
2671: `KSP`, and `PC` contexts as well.
2673: `TSGetSNES()` does not work for integrators that do not use `SNES`; in
2674: this case `TSGetSNES()` returns `NULL` in `snes`.
2676: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2677: @*/
2678: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2679: {
2680: PetscFunctionBegin;
2682: PetscAssertPointer(snes, 2);
2683: if (!ts->snes) {
2684: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2685: PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2686: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2687: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2688: if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2689: if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2690: }
2691: *snes = ts->snes;
2692: PetscFunctionReturn(PETSC_SUCCESS);
2693: }
2695: /*@
2696: TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the timestepping context
2698: Collective
2700: Input Parameters:
2701: + ts - the `TS` context obtained from `TSCreate()`
2702: - snes - the nonlinear solver context
2704: Level: developer
2706: Note:
2707: Most users should have the `TS` created by calling `TSGetSNES()`
2709: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2710: @*/
2711: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2712: {
2713: PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);
2715: PetscFunctionBegin;
2718: PetscCall(PetscObjectReference((PetscObject)snes));
2719: PetscCall(SNESDestroy(&ts->snes));
2721: ts->snes = snes;
2723: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2724: PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2725: if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2726: PetscFunctionReturn(PETSC_SUCCESS);
2727: }
2729: /*@
2730: TSGetKSP - Returns the `KSP` (linear solver) associated with
2731: a `TS` (timestepper) context.
2733: Not Collective, but `ksp` is parallel if `ts` is parallel
2735: Input Parameter:
2736: . ts - the `TS` context obtained from `TSCreate()`
2738: Output Parameter:
2739: . ksp - the nonlinear solver context
2741: Level: beginner
2743: Notes:
2744: The user can then directly manipulate the `KSP` context to set various
2745: options, etc. Likewise, the user can then extract and manipulate the
2746: `PC` context as well.
2748: `TSGetKSP()` does not work for integrators that do not use `KSP`;
2749: in this case `TSGetKSP()` returns `NULL` in `ksp`.
2751: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2752: @*/
2753: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2754: {
2755: SNES snes;
2757: PetscFunctionBegin;
2759: PetscAssertPointer(ksp, 2);
2760: PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2761: PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
2762: PetscCall(TSGetSNES(ts, &snes));
2763: PetscCall(SNESGetKSP(snes, ksp));
2764: PetscFunctionReturn(PETSC_SUCCESS);
2765: }
2767: /* ----------- Routines to set solver parameters ---------- */
2769: /*@
2770: TSSetMaxSteps - Sets the maximum number of steps to use.
2772: Logically Collective
2774: Input Parameters:
2775: + ts - the `TS` context obtained from `TSCreate()`
2776: - maxsteps - maximum number of steps to use
2778: Options Database Key:
2779: . -ts_max_steps <maxsteps> - Sets maxsteps
2781: Level: intermediate
2783: Note:
2784: Use `PETSC_DETERMINE` to reset the maximum number of steps to the default from when the object's type was set
2786: The default maximum number of steps is 5,000
2788: Fortran Note:
2789: Use `PETSC_DETERMINE_INTEGER`
2791: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2792: @*/
2793: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2794: {
2795: PetscFunctionBegin;
2798: if (maxsteps == PETSC_DETERMINE) {
2799: ts->max_steps = ts->default_max_steps;
2800: } else {
2801: PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2802: ts->max_steps = maxsteps;
2803: }
2804: PetscFunctionReturn(PETSC_SUCCESS);
2805: }
2807: /*@
2808: TSGetMaxSteps - Gets the maximum number of steps to use.
2810: Not Collective
2812: Input Parameter:
2813: . ts - the `TS` context obtained from `TSCreate()`
2815: Output Parameter:
2816: . maxsteps - maximum number of steps to use
2818: Level: advanced
2820: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2821: @*/
2822: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2823: {
2824: PetscFunctionBegin;
2826: PetscAssertPointer(maxsteps, 2);
2827: *maxsteps = ts->max_steps;
2828: PetscFunctionReturn(PETSC_SUCCESS);
2829: }
2831: /*@
2832: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
2834: Logically Collective
2836: Input Parameters:
2837: + ts - the `TS` context obtained from `TSCreate()`
2838: - maxtime - final time to step to
2840: Options Database Key:
2841: . -ts_max_time <maxtime> - Sets maxtime
2843: Level: intermediate
2845: Notes:
2846: Use `PETSC_DETERMINE` to reset the maximum time to the default from when the object's type was set
2848: The default maximum time is 5.0
2850: Fortran Note:
2851: Use `PETSC_DETERMINE_REAL`
2853: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2854: @*/
2855: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2856: {
2857: PetscFunctionBegin;
2860: if (maxtime == PETSC_DETERMINE) {
2861: ts->max_time = ts->default_max_time;
2862: } else {
2863: ts->max_time = maxtime;
2864: }
2865: PetscFunctionReturn(PETSC_SUCCESS);
2866: }
2868: /*@
2869: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
2871: Not Collective
2873: Input Parameter:
2874: . ts - the `TS` context obtained from `TSCreate()`
2876: Output Parameter:
2877: . maxtime - final time to step to
2879: Level: advanced
2881: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2882: @*/
2883: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2884: {
2885: PetscFunctionBegin;
2887: PetscAssertPointer(maxtime, 2);
2888: *maxtime = ts->max_time;
2889: PetscFunctionReturn(PETSC_SUCCESS);
2890: }
2892: // PetscClangLinter pragma disable: -fdoc-*
2893: /*@
2894: TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.
2896: Level: deprecated
2898: @*/
2899: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2900: {
2901: PetscFunctionBegin;
2903: PetscCall(TSSetTime(ts, initial_time));
2904: PetscCall(TSSetTimeStep(ts, time_step));
2905: PetscFunctionReturn(PETSC_SUCCESS);
2906: }
2908: // PetscClangLinter pragma disable: -fdoc-*
2909: /*@
2910: TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.
2912: Level: deprecated
2914: @*/
2915: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2916: {
2917: PetscFunctionBegin;
2919: if (maxsteps) {
2920: PetscAssertPointer(maxsteps, 2);
2921: *maxsteps = ts->max_steps;
2922: }
2923: if (maxtime) {
2924: PetscAssertPointer(maxtime, 3);
2925: *maxtime = ts->max_time;
2926: }
2927: PetscFunctionReturn(PETSC_SUCCESS);
2928: }
2930: // PetscClangLinter pragma disable: -fdoc-*
2931: /*@
2932: TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.
2934: Level: deprecated
2936: @*/
2937: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2938: {
2939: PetscFunctionBegin;
2940: if (maxsteps != PETSC_CURRENT) PetscCall(TSSetMaxSteps(ts, maxsteps));
2941: if (maxtime != (PetscReal)PETSC_CURRENT) PetscCall(TSSetMaxTime(ts, maxtime));
2942: PetscFunctionReturn(PETSC_SUCCESS);
2943: }
2945: // PetscClangLinter pragma disable: -fdoc-*
2946: /*@
2947: TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.
2949: Level: deprecated
2951: @*/
2952: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
2953: {
2954: return TSGetStepNumber(ts, steps);
2955: }
2957: // PetscClangLinter pragma disable: -fdoc-*
2958: /*@
2959: TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.
2961: Level: deprecated
2963: @*/
2964: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
2965: {
2966: return TSGetStepNumber(ts, steps);
2967: }
2969: /*@
2970: TSSetSolution - Sets the initial solution vector
2971: for use by the `TS` routines.
2973: Logically Collective
2975: Input Parameters:
2976: + ts - the `TS` context obtained from `TSCreate()`
2977: - u - the solution vector
2979: Level: beginner
2981: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
2982: @*/
2983: PetscErrorCode TSSetSolution(TS ts, Vec u)
2984: {
2985: DM dm;
2987: PetscFunctionBegin;
2990: PetscCall(PetscObjectReference((PetscObject)u));
2991: PetscCall(VecDestroy(&ts->vec_sol));
2992: ts->vec_sol = u;
2994: PetscCall(TSGetDM(ts, &dm));
2995: PetscCall(DMShellSetGlobalVector(dm, u));
2996: PetscFunctionReturn(PETSC_SUCCESS);
2997: }
2999: /*@C
3000: TSSetPreStep - Sets the general-purpose function
3001: called once at the beginning of each time step.
3003: Logically Collective
3005: Input Parameters:
3006: + ts - The `TS` context obtained from `TSCreate()`
3007: - func - The function
3009: Calling sequence of `func`:
3010: . ts - the `TS` context
3012: Level: intermediate
3014: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3015: @*/
3016: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
3017: {
3018: PetscFunctionBegin;
3020: ts->prestep = func;
3021: PetscFunctionReturn(PETSC_SUCCESS);
3022: }
3024: /*@
3025: TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`
3027: Collective
3029: Input Parameter:
3030: . ts - The `TS` context obtained from `TSCreate()`
3032: Level: developer
3034: Note:
3035: `TSPreStep()` is typically used within time stepping implementations,
3036: so most users would not generally call this routine themselves.
3038: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3039: @*/
3040: PetscErrorCode TSPreStep(TS ts)
3041: {
3042: PetscFunctionBegin;
3044: if (ts->prestep) {
3045: Vec U;
3046: PetscObjectId idprev;
3047: PetscBool sameObject;
3048: PetscObjectState sprev, spost;
3050: PetscCall(TSGetSolution(ts, &U));
3051: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3052: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3053: PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3054: PetscCall(TSGetSolution(ts, &U));
3055: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3056: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3057: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3058: }
3059: PetscFunctionReturn(PETSC_SUCCESS);
3060: }
3062: /*@C
3063: TSSetPreStage - Sets the general-purpose function
3064: called once at the beginning of each stage.
3066: Logically Collective
3068: Input Parameters:
3069: + ts - The `TS` context obtained from `TSCreate()`
3070: - func - The function
3072: Calling sequence of `func`:
3073: + ts - the `TS` context
3074: - stagetime - the stage time
3076: Level: intermediate
3078: Note:
3079: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3080: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3081: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3083: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3084: @*/
3085: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime))
3086: {
3087: PetscFunctionBegin;
3089: ts->prestage = func;
3090: PetscFunctionReturn(PETSC_SUCCESS);
3091: }
3093: /*@C
3094: TSSetPostStage - Sets the general-purpose function
3095: called once at the end of each stage.
3097: Logically Collective
3099: Input Parameters:
3100: + ts - The `TS` context obtained from `TSCreate()`
3101: - func - The function
3103: Calling sequence of `func`:
3104: + ts - the `TS` context
3105: . stagetime - the stage time
3106: . stageindex - the stage index
3107: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3109: Level: intermediate
3111: Note:
3112: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3113: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3114: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3116: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3117: @*/
3118: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
3119: {
3120: PetscFunctionBegin;
3122: ts->poststage = func;
3123: PetscFunctionReturn(PETSC_SUCCESS);
3124: }
3126: /*@C
3127: TSSetPostEvaluate - Sets the general-purpose function
3128: called at the end of each step evaluation.
3130: Logically Collective
3132: Input Parameters:
3133: + ts - The `TS` context obtained from `TSCreate()`
3134: - func - The function
3136: Calling sequence of `func`:
3137: . ts - the `TS` context
3139: Level: intermediate
3141: Note:
3142: The function set by `TSSetPostEvaluate()` is called after the solution is evaluated, or after the step rollback.
3143: Inside the `func` callback, the solution vector can be obtained with `TSGetSolution()`, and modified, if need be.
3144: The time step can be obtained with `TSGetTimeStep()`, and the time at the start of the step - via `TSGetTime()`.
3145: The potential changes to the solution vector introduced by event handling (`postevent()`) are not relevant for `TSSetPostEvaluate()`,
3146: but are relevant for `TSSetPostStep()`, according to the function call scheme in `TSSolve()`, as shown below
3147: .vb
3148: ...
3149: Step()
3150: PostEvaluate()
3151: EventHandling()
3152: step_rollback ? PostEvaluate() : PostStep()
3153: ...
3154: .ve
3155: where EventHandling() may result in one of the following three outcomes
3156: .vb
3157: (1) | successful step | solution intact
3158: (2) | successful step | solution modified by `postevent()`
3159: (3) | step_rollback | solution rolled back
3160: .ve
3162: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3163: @*/
3164: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
3165: {
3166: PetscFunctionBegin;
3168: ts->postevaluate = func;
3169: PetscFunctionReturn(PETSC_SUCCESS);
3170: }
3172: /*@
3173: TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`
3175: Collective
3177: Input Parameters:
3178: + ts - The `TS` context obtained from `TSCreate()`
3179: - stagetime - The absolute time of the current stage
3181: Level: developer
3183: Note:
3184: `TSPreStage()` is typically used within time stepping implementations,
3185: most users would not generally call this routine themselves.
3187: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3188: @*/
3189: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3190: {
3191: PetscFunctionBegin;
3193: if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3194: PetscFunctionReturn(PETSC_SUCCESS);
3195: }
3197: /*@
3198: TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`
3200: Collective
3202: Input Parameters:
3203: + ts - The `TS` context obtained from `TSCreate()`
3204: . stagetime - The absolute time of the current stage
3205: . stageindex - Stage number
3206: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3208: Level: developer
3210: Note:
3211: `TSPostStage()` is typically used within time stepping implementations,
3212: most users would not generally call this routine themselves.
3214: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3215: @*/
3216: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3217: {
3218: PetscFunctionBegin;
3220: if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3221: PetscFunctionReturn(PETSC_SUCCESS);
3222: }
3224: /*@
3225: TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`
3227: Collective
3229: Input Parameter:
3230: . ts - The `TS` context obtained from `TSCreate()`
3232: Level: developer
3234: Note:
3235: `TSPostEvaluate()` is typically used within time stepping implementations,
3236: most users would not generally call this routine themselves.
3238: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3239: @*/
3240: PetscErrorCode TSPostEvaluate(TS ts)
3241: {
3242: PetscFunctionBegin;
3244: if (ts->postevaluate) {
3245: Vec U;
3246: PetscObjectState sprev, spost;
3248: PetscCall(TSGetSolution(ts, &U));
3249: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3250: PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3251: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3252: if (sprev != spost) PetscCall(TSRestartStep(ts));
3253: }
3254: PetscFunctionReturn(PETSC_SUCCESS);
3255: }
3257: /*@C
3258: TSSetPostStep - Sets the general-purpose function
3259: called once at the end of each successful time step.
3261: Logically Collective
3263: Input Parameters:
3264: + ts - The `TS` context obtained from `TSCreate()`
3265: - func - The function
3267: Calling sequence of `func`:
3268: . ts - the `TS` context
3270: Level: intermediate
3272: Note:
3273: The function set by `TSSetPostStep()` is called after each successful step. If the event handler locates an event at the
3274: given step, and `postevent()` modifies the solution vector, the solution vector obtained by `TSGetSolution()` inside `func` will
3275: contain the changes. To get the solution without these changes, use `TSSetPostEvaluate()` to set the appropriate callback.
3276: The scheme of the relevant function calls in `TSSolve()` is shown below
3277: .vb
3278: ...
3279: Step()
3280: PostEvaluate()
3281: EventHandling()
3282: step_rollback ? PostEvaluate() : PostStep()
3283: ...
3284: .ve
3285: where EventHandling() may result in one of the following three outcomes
3286: .vb
3287: (1) | successful step | solution intact
3288: (2) | successful step | solution modified by `postevent()`
3289: (3) | step_rollback | solution rolled back
3290: .ve
3292: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3293: @*/
3294: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
3295: {
3296: PetscFunctionBegin;
3298: ts->poststep = func;
3299: PetscFunctionReturn(PETSC_SUCCESS);
3300: }
3302: /*@
3303: TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`
3305: Collective
3307: Input Parameter:
3308: . ts - The `TS` context obtained from `TSCreate()`
3310: Note:
3311: `TSPostStep()` is typically used within time stepping implementations,
3312: so most users would not generally call this routine themselves.
3314: Level: developer
3316: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()`
3317: @*/
3318: PetscErrorCode TSPostStep(TS ts)
3319: {
3320: PetscFunctionBegin;
3322: if (ts->poststep) {
3323: Vec U;
3324: PetscObjectId idprev;
3325: PetscBool sameObject;
3326: PetscObjectState sprev, spost;
3328: PetscCall(TSGetSolution(ts, &U));
3329: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3330: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3331: PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3332: PetscCall(TSGetSolution(ts, &U));
3333: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3334: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3335: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3336: }
3337: PetscFunctionReturn(PETSC_SUCCESS);
3338: }
3340: /*@
3341: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3343: Collective
3345: Input Parameters:
3346: + ts - time stepping context
3347: - t - time to interpolate to
3349: Output Parameter:
3350: . U - state at given time
3352: Level: intermediate
3354: Developer Notes:
3355: `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3357: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3358: @*/
3359: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3360: {
3361: PetscFunctionBegin;
3364: 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);
3365: PetscUseTypeMethod(ts, interpolate, t, U);
3366: PetscFunctionReturn(PETSC_SUCCESS);
3367: }
3369: /*@
3370: TSStep - Steps one time step
3372: Collective
3374: Input Parameter:
3375: . ts - the `TS` context obtained from `TSCreate()`
3377: Level: developer
3379: Notes:
3380: The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.
3382: The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
3383: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3385: This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
3386: time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.
3388: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3389: @*/
3390: PetscErrorCode TSStep(TS ts)
3391: {
3392: static PetscBool cite = PETSC_FALSE;
3393: PetscReal ptime;
3395: PetscFunctionBegin;
3397: PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3398: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3399: " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3400: " journal = {arXiv e-preprints},\n"
3401: " eprint = {1806.01437},\n"
3402: " archivePrefix = {arXiv},\n"
3403: " year = {2018}\n}\n",
3404: &cite));
3405: PetscCall(TSSetUp(ts));
3406: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3407: if (ts->tspan)
3408: ts->tspan->worktol = 0; /* In each step of TSSolve() 'tspan->worktol' will be meaningfully defined (later) only once:
3409: in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */
3411: 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>");
3412: 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()");
3413: 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");
3415: if (!ts->vec_sol0) PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sol0));
3416: PetscCall(VecCopy(ts->vec_sol, ts->vec_sol0));
3417: ts->time_step0 = ts->time_step;
3419: if (!ts->steps) ts->ptime_prev = ts->ptime;
3420: ptime = ts->ptime;
3422: ts->ptime_prev_rollback = ts->ptime_prev;
3423: ts->reason = TS_CONVERGED_ITERATING;
3425: PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3426: PetscUseTypeMethod(ts, step);
3427: PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));
3429: if (ts->reason >= 0) {
3430: ts->ptime_prev = ptime;
3431: ts->steps++;
3432: ts->steprollback = PETSC_FALSE;
3433: ts->steprestart = PETSC_FALSE;
3434: ts->stepresize = PETSC_FALSE;
3435: }
3437: if (ts->reason < 0 && ts->errorifstepfailed) {
3438: PetscCall(TSMonitorCancel(ts));
3439: 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]);
3440: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3441: }
3442: PetscFunctionReturn(PETSC_SUCCESS);
3443: }
3445: /*@
3446: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3447: at the end of a time step with a given order of accuracy.
3449: Collective
3451: Input Parameters:
3452: + ts - time stepping context
3453: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3455: Input/Output Parameter:
3456: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3457: on output, the actual order of the error evaluation
3459: Output Parameter:
3460: . wlte - the weighted local truncation error norm
3462: Level: advanced
3464: Note:
3465: If the timestepper cannot evaluate the error in a particular step
3466: (eg. in the first step or restart steps after event handling),
3467: this routine returns wlte=-1.0 .
3469: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3470: @*/
3471: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3472: {
3473: PetscFunctionBegin;
3477: if (order) PetscAssertPointer(order, 3);
3479: PetscAssertPointer(wlte, 4);
3480: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3481: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3482: PetscFunctionReturn(PETSC_SUCCESS);
3483: }
3485: /*@
3486: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3488: Collective
3490: Input Parameters:
3491: + ts - time stepping context
3492: . order - desired order of accuracy
3493: - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)
3495: Output Parameter:
3496: . U - state at the end of the current step
3498: Level: advanced
3500: Notes:
3501: This function cannot be called until all stages have been evaluated.
3503: 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.
3505: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3506: @*/
3507: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3508: {
3509: PetscFunctionBegin;
3513: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3514: PetscFunctionReturn(PETSC_SUCCESS);
3515: }
3517: /*@C
3518: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3520: Not collective
3522: Input Parameter:
3523: . ts - time stepping context
3525: Output Parameter:
3526: . initCondition - The function which computes an initial condition
3528: Calling sequence of `initCondition`:
3529: + ts - The timestepping context
3530: - u - The input vector in which the initial condition is stored
3532: Level: advanced
3534: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3535: @*/
3536: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3537: {
3538: PetscFunctionBegin;
3540: PetscAssertPointer(initCondition, 2);
3541: *initCondition = ts->ops->initcondition;
3542: PetscFunctionReturn(PETSC_SUCCESS);
3543: }
3545: /*@C
3546: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3548: Logically collective
3550: Input Parameters:
3551: + ts - time stepping context
3552: - initCondition - The function which computes an initial condition
3554: Calling sequence of `initCondition`:
3555: + ts - The timestepping context
3556: - e - The input vector in which the initial condition is to be stored
3558: Level: advanced
3560: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3561: @*/
3562: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3563: {
3564: PetscFunctionBegin;
3567: ts->ops->initcondition = initCondition;
3568: PetscFunctionReturn(PETSC_SUCCESS);
3569: }
3571: /*@
3572: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3574: Collective
3576: Input Parameters:
3577: + ts - time stepping context
3578: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3580: Level: advanced
3582: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3583: @*/
3584: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3585: {
3586: PetscFunctionBegin;
3589: PetscTryTypeMethod(ts, initcondition, u);
3590: PetscFunctionReturn(PETSC_SUCCESS);
3591: }
3593: /*@C
3594: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3596: Not collective
3598: Input Parameter:
3599: . ts - time stepping context
3601: Output Parameter:
3602: . exactError - The function which computes the solution error
3604: Calling sequence of `exactError`:
3605: + ts - The timestepping context
3606: . u - The approximate solution vector
3607: - e - The vector in which the error is stored
3609: Level: advanced
3611: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3612: @*/
3613: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3614: {
3615: PetscFunctionBegin;
3617: PetscAssertPointer(exactError, 2);
3618: *exactError = ts->ops->exacterror;
3619: PetscFunctionReturn(PETSC_SUCCESS);
3620: }
3622: /*@C
3623: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3625: Logically collective
3627: Input Parameters:
3628: + ts - time stepping context
3629: - exactError - The function which computes the solution error
3631: Calling sequence of `exactError`:
3632: + ts - The timestepping context
3633: . u - The approximate solution vector
3634: - e - The vector in which the error is stored
3636: Level: advanced
3638: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3639: @*/
3640: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3641: {
3642: PetscFunctionBegin;
3645: ts->ops->exacterror = exactError;
3646: PetscFunctionReturn(PETSC_SUCCESS);
3647: }
3649: /*@
3650: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3652: Collective
3654: Input Parameters:
3655: + ts - time stepping context
3656: . u - The approximate solution
3657: - e - The `Vec` used to store the error
3659: Level: advanced
3661: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3662: @*/
3663: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3664: {
3665: PetscFunctionBegin;
3669: PetscTryTypeMethod(ts, exacterror, u, e);
3670: PetscFunctionReturn(PETSC_SUCCESS);
3671: }
3673: /*@C
3674: TSSetResize - Sets the resize callbacks.
3676: Logically Collective
3678: Input Parameters:
3679: + ts - The `TS` context obtained from `TSCreate()`
3680: . rollback - Whether a resize will restart the step
3681: . setup - The setup function
3682: . transfer - The transfer function
3683: - ctx - [optional] The user-defined context
3685: Calling sequence of `setup`:
3686: + ts - the `TS` context
3687: . step - the current step
3688: . time - the current time
3689: . state - the current vector of state
3690: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3691: - ctx - user defined context
3693: Calling sequence of `transfer`:
3694: + ts - the `TS` context
3695: . nv - the number of vectors to be transferred
3696: . vecsin - array of vectors to be transferred
3697: . vecsout - array of transferred vectors
3698: - ctx - user defined context
3700: Notes:
3701: The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
3702: depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
3703: and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
3704: that the size of the discrete problem has changed.
3705: In both cases, the solver will collect the needed vectors that will be
3706: transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
3707: current solution vector, and other vectors needed by the specific solver used.
3708: For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3709: Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3710: will be automatically reset if the sizes are changed and they must be specified again by the user
3711: inside the `transfer` function.
3712: The input and output arrays passed to `transfer` are allocated by PETSc.
3713: Vectors in `vecsout` must be created by the user.
3714: Ownership of vectors in `vecsout` is transferred to PETSc.
3716: Level: advanced
3718: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3719: @*/
3720: 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)
3721: {
3722: PetscFunctionBegin;
3724: ts->resizerollback = rollback;
3725: ts->resizesetup = setup;
3726: ts->resizetransfer = transfer;
3727: ts->resizectx = ctx;
3728: PetscFunctionReturn(PETSC_SUCCESS);
3729: }
3731: /*
3732: TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.
3734: Collective
3736: Input Parameters:
3737: + ts - The `TS` context obtained from `TSCreate()`
3738: - 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.
3740: Level: developer
3742: Note:
3743: `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3744: used within time stepping implementations,
3745: so most users would not generally call this routine themselves.
3747: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3748: @*/
3749: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3750: {
3751: PetscFunctionBegin;
3753: PetscTryTypeMethod(ts, resizeregister, flg);
3754: /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3755: PetscFunctionReturn(PETSC_SUCCESS);
3756: }
3758: static PetscErrorCode TSResizeReset(TS ts)
3759: {
3760: PetscFunctionBegin;
3762: PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3763: PetscFunctionReturn(PETSC_SUCCESS);
3764: }
3766: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3767: {
3768: PetscFunctionBegin;
3771: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3772: if (ts->resizetransfer) {
3773: PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3774: PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3775: }
3776: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3777: PetscFunctionReturn(PETSC_SUCCESS);
3778: }
3780: /*@C
3781: TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.
3783: Collective
3785: Input Parameters:
3786: + ts - The `TS` context obtained from `TSCreate()`
3787: . name - A string identifying the vector
3788: - vec - The vector
3790: Level: developer
3792: Note:
3793: `TSResizeRegisterVec()` is typically used within time stepping implementations,
3794: so most users would not generally call this routine themselves.
3796: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3797: @*/
3798: PetscErrorCode TSResizeRegisterVec(TS ts, const char name[], Vec vec)
3799: {
3800: PetscFunctionBegin;
3802: PetscAssertPointer(name, 2);
3804: PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3805: PetscFunctionReturn(PETSC_SUCCESS);
3806: }
3808: /*@C
3809: TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.
3811: Collective
3813: Input Parameters:
3814: + ts - The `TS` context obtained from `TSCreate()`
3815: . name - A string identifying the vector
3816: - vec - The vector
3818: Level: developer
3820: Note:
3821: `TSResizeRetrieveVec()` is typically used within time stepping implementations,
3822: so most users would not generally call this routine themselves.
3824: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
3825: @*/
3826: PetscErrorCode TSResizeRetrieveVec(TS ts, const char name[], Vec *vec)
3827: {
3828: PetscFunctionBegin;
3830: PetscAssertPointer(name, 2);
3831: PetscAssertPointer(vec, 3);
3832: PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
3833: PetscFunctionReturn(PETSC_SUCCESS);
3834: }
3836: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
3837: {
3838: PetscInt cnt;
3839: PetscObjectList tmp;
3840: Vec *vecsin = NULL;
3841: const char **namesin = NULL;
3843: PetscFunctionBegin;
3844: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
3845: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
3846: if (names) PetscCall(PetscMalloc1(cnt, &namesin));
3847: if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
3848: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
3849: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
3850: if (vecs) vecsin[cnt] = (Vec)tmp->obj;
3851: if (names) namesin[cnt] = tmp->name;
3852: cnt++;
3853: }
3854: }
3855: if (nv) *nv = cnt;
3856: if (names) *names = namesin;
3857: if (vecs) *vecs = vecsin;
3858: PetscFunctionReturn(PETSC_SUCCESS);
3859: }
3861: /*@
3862: TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`
3864: Collective
3866: Input Parameter:
3867: . ts - The `TS` context obtained from `TSCreate()`
3869: Level: developer
3871: Note:
3872: `TSResize()` is typically used within time stepping implementations,
3873: so most users would not generally call this routine themselves.
3875: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3876: @*/
3877: PetscErrorCode TSResize(TS ts)
3878: {
3879: PetscInt nv = 0;
3880: const char **names = NULL;
3881: Vec *vecsin = NULL;
3882: const char *solname = "ts:vec_sol";
3884: PetscFunctionBegin;
3886: if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
3887: if (ts->resizesetup) {
3888: PetscCall(VecLockReadPush(ts->vec_sol));
3889: PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &ts->stepresize, ts->resizectx));
3890: PetscCall(VecLockReadPop(ts->vec_sol));
3891: if (ts->stepresize) {
3892: if (ts->resizerollback) {
3893: PetscCall(TSRollBack(ts));
3894: ts->time_step = ts->time_step0;
3895: }
3896: PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
3897: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
3898: }
3899: }
3901: PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
3902: if (nv) {
3903: Vec *vecsout, vecsol;
3905: /* Reset internal objects */
3906: PetscCall(TSReset(ts));
3908: /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
3909: PetscCall(PetscCalloc1(nv, &vecsout));
3910: PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
3911: for (PetscInt i = 0; i < nv; i++) {
3912: const char *name;
3913: char *oname;
3915: PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
3916: PetscCall(PetscStrallocpy(name, &oname));
3917: PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
3918: if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
3919: PetscCall(PetscFree(oname));
3920: PetscCall(VecDestroy(&vecsout[i]));
3921: }
3922: PetscCall(PetscFree(vecsout));
3923: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */
3925: PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
3926: if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
3927: PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
3928: }
3930: PetscCall(PetscFree(names));
3931: PetscCall(PetscFree(vecsin));
3932: PetscCall(TSResizeReset(ts));
3933: PetscFunctionReturn(PETSC_SUCCESS);
3934: }
3936: /*@
3937: TSSolve - Steps the requested number of timesteps.
3939: Collective
3941: Input Parameters:
3942: + ts - the `TS` context obtained from `TSCreate()`
3943: - u - the solution vector (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
3944: otherwise must contain the initial conditions and will contain the solution at the final requested time
3946: Level: beginner
3948: Notes:
3949: The final time returned by this function may be different from the time of the internally
3950: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
3951: stepped over the final time.
3953: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3954: @*/
3955: PetscErrorCode TSSolve(TS ts, Vec u)
3956: {
3957: Vec solution;
3959: PetscFunctionBegin;
3963: PetscCall(TSSetExactFinalTimeDefault(ts));
3964: 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 */
3965: if (!ts->vec_sol || u == ts->vec_sol) {
3966: PetscCall(VecDuplicate(u, &solution));
3967: PetscCall(TSSetSolution(ts, solution));
3968: PetscCall(VecDestroy(&solution)); /* grant ownership */
3969: }
3970: PetscCall(VecCopy(u, ts->vec_sol));
3971: PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
3972: } else if (u) PetscCall(TSSetSolution(ts, u));
3973: PetscCall(TSSetUp(ts));
3974: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3976: 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>");
3977: 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()");
3978: 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");
3979: PetscCheck(!(ts->tspan && ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP), PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "You must use TS_EXACTFINALTIME_MATCHSTEP when using time span");
3981: if (ts->tspan && PetscIsCloseAtTol(ts->ptime, ts->tspan->span_times[0], ts->tspan->reltol * ts->time_step + ts->tspan->abstol, 0)) { /* starting point in time span */
3982: PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]));
3983: ts->tspan->spanctr = 1;
3984: }
3986: if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));
3988: /* reset number of steps only when the step is not restarted. ARKIMEX
3989: restarts the step after an event. Resetting these counters in such case causes
3990: TSTrajectory to incorrectly save the output files
3991: */
3992: /* reset time step and iteration counters */
3993: if (!ts->steps) {
3994: ts->ksp_its = 0;
3995: ts->snes_its = 0;
3996: ts->num_snes_failures = 0;
3997: ts->reject = 0;
3998: ts->steprestart = PETSC_TRUE;
3999: ts->steprollback = PETSC_FALSE;
4000: ts->stepresize = PETSC_FALSE;
4001: ts->rhsjacobian.time = PETSC_MIN_REAL;
4002: }
4004: /* make sure initial time step does not overshoot final time or the next point in tspan */
4005: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
4006: PetscReal maxdt;
4007: PetscReal dt = ts->time_step;
4009: if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
4010: else maxdt = ts->max_time - ts->ptime;
4011: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
4012: }
4013: ts->reason = TS_CONVERGED_ITERATING;
4015: {
4016: PetscViewer viewer;
4017: PetscViewerFormat format;
4018: PetscBool flg;
4019: static PetscBool incall = PETSC_FALSE;
4021: if (!incall) {
4022: /* Estimate the convergence rate of the time discretization */
4023: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
4024: if (flg) {
4025: PetscConvEst conv;
4026: DM dm;
4027: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4028: PetscInt Nf;
4029: PetscBool checkTemporal = PETSC_TRUE;
4031: incall = PETSC_TRUE;
4032: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
4033: PetscCall(TSGetDM(ts, &dm));
4034: PetscCall(DMGetNumFields(dm, &Nf));
4035: PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
4036: PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
4037: PetscCall(PetscConvEstUseTS(conv, checkTemporal));
4038: PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
4039: PetscCall(PetscConvEstSetFromOptions(conv));
4040: PetscCall(PetscConvEstSetUp(conv));
4041: PetscCall(PetscConvEstGetConvRate(conv, alpha));
4042: PetscCall(PetscViewerPushFormat(viewer, format));
4043: PetscCall(PetscConvEstRateView(conv, alpha, viewer));
4044: PetscCall(PetscViewerPopFormat(viewer));
4045: PetscCall(PetscViewerDestroy(&viewer));
4046: PetscCall(PetscConvEstDestroy(&conv));
4047: PetscCall(PetscFree(alpha));
4048: incall = PETSC_FALSE;
4049: }
4050: }
4051: }
4053: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));
4055: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4056: PetscUseTypeMethod(ts, solve);
4057: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4058: ts->solvetime = ts->ptime;
4059: solution = ts->vec_sol;
4060: } else { /* Step the requested number of timesteps. */
4061: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4062: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4064: if (!ts->steps) {
4065: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4066: PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4067: }
4069: while (!ts->reason) {
4070: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4071: if (!ts->steprollback || (ts->stepresize && ts->resizerollback)) PetscCall(TSPreStep(ts));
4072: PetscCall(TSStep(ts));
4073: if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4074: if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4075: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4076: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4077: PetscCall(TSForwardCostIntegral(ts));
4078: if (ts->reason >= 0) ts->steps++;
4079: }
4080: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4081: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4082: PetscCall(TSForwardStep(ts));
4083: if (ts->reason >= 0) ts->steps++;
4084: }
4085: PetscCall(TSPostEvaluate(ts));
4086: 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. */
4087: if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4088: if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
4089: /* check convergence */
4090: if (!ts->reason) {
4091: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4092: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4093: }
4094: if (!ts->steprollback) {
4095: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4096: PetscCall(TSPostStep(ts));
4097: if (!ts->resizerollback) PetscCall(TSResize(ts));
4099: if (ts->tspan && ts->tspan->spanctr < ts->tspan->num_span_times) {
4100: PetscCheck(ts->tspan->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(tspan->worktol > 0) in TSSolve()");
4101: if (PetscIsCloseAtTol(ts->ptime, ts->tspan->span_times[ts->tspan->spanctr], ts->tspan->worktol, 0)) PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[ts->tspan->spanctr++]));
4102: }
4103: }
4104: }
4105: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4107: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4108: if (!u) u = ts->vec_sol;
4109: PetscCall(TSInterpolate(ts, ts->max_time, u));
4110: ts->solvetime = ts->max_time;
4111: solution = u;
4112: PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4113: } else {
4114: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4115: ts->solvetime = ts->ptime;
4116: solution = ts->vec_sol;
4117: }
4118: }
4120: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4121: PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4122: PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4123: if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4124: PetscFunctionReturn(PETSC_SUCCESS);
4125: }
4127: /*@
4128: TSGetTime - Gets the time of the most recently completed step.
4130: Not Collective
4132: Input Parameter:
4133: . ts - the `TS` context obtained from `TSCreate()`
4135: Output Parameter:
4136: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
4138: Level: beginner
4140: Note:
4141: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
4142: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
4144: .seealso: [](ch_ts), `TS`, `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4145: @*/
4146: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4147: {
4148: PetscFunctionBegin;
4150: PetscAssertPointer(t, 2);
4151: *t = ts->ptime;
4152: PetscFunctionReturn(PETSC_SUCCESS);
4153: }
4155: /*@
4156: TSGetPrevTime - Gets the starting time of the previously completed step.
4158: Not Collective
4160: Input Parameter:
4161: . ts - the `TS` context obtained from `TSCreate()`
4163: Output Parameter:
4164: . t - the previous time
4166: Level: beginner
4168: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4169: @*/
4170: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4171: {
4172: PetscFunctionBegin;
4174: PetscAssertPointer(t, 2);
4175: *t = ts->ptime_prev;
4176: PetscFunctionReturn(PETSC_SUCCESS);
4177: }
4179: /*@
4180: TSSetTime - Allows one to reset the time.
4182: Logically Collective
4184: Input Parameters:
4185: + ts - the `TS` context obtained from `TSCreate()`
4186: - t - the time
4188: Level: intermediate
4190: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4191: @*/
4192: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4193: {
4194: PetscFunctionBegin;
4197: ts->ptime = t;
4198: PetscFunctionReturn(PETSC_SUCCESS);
4199: }
4201: /*@
4202: TSSetOptionsPrefix - Sets the prefix used for searching for all
4203: TS options in the database.
4205: Logically Collective
4207: Input Parameters:
4208: + ts - The `TS` context
4209: - prefix - The prefix to prepend to all option names
4211: Level: advanced
4213: Note:
4214: A hyphen (-) must NOT be given at the beginning of the prefix name.
4215: The first character of all runtime options is AUTOMATICALLY the
4216: hyphen.
4218: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4219: @*/
4220: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4221: {
4222: SNES snes;
4224: PetscFunctionBegin;
4226: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4227: PetscCall(TSGetSNES(ts, &snes));
4228: PetscCall(SNESSetOptionsPrefix(snes, prefix));
4229: PetscFunctionReturn(PETSC_SUCCESS);
4230: }
4232: /*@
4233: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4234: TS options in the database.
4236: Logically Collective
4238: Input Parameters:
4239: + ts - The `TS` context
4240: - prefix - The prefix to prepend to all option names
4242: Level: advanced
4244: Note:
4245: A hyphen (-) must NOT be given at the beginning of the prefix name.
4246: The first character of all runtime options is AUTOMATICALLY the
4247: hyphen.
4249: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4250: @*/
4251: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4252: {
4253: SNES snes;
4255: PetscFunctionBegin;
4257: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4258: PetscCall(TSGetSNES(ts, &snes));
4259: PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4260: PetscFunctionReturn(PETSC_SUCCESS);
4261: }
4263: /*@
4264: TSGetOptionsPrefix - Sets the prefix used for searching for all
4265: `TS` options in the database.
4267: Not Collective
4269: Input Parameter:
4270: . ts - The `TS` context
4272: Output Parameter:
4273: . prefix - A pointer to the prefix string used
4275: Level: intermediate
4277: Fortran Notes:
4278: The user should pass in a string 'prefix' of
4279: sufficient length to hold the prefix.
4281: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4282: @*/
4283: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4284: {
4285: PetscFunctionBegin;
4287: PetscAssertPointer(prefix, 2);
4288: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4289: PetscFunctionReturn(PETSC_SUCCESS);
4290: }
4292: /*@C
4293: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4295: Not Collective, but parallel objects are returned if ts is parallel
4297: Input Parameter:
4298: . ts - The `TS` context obtained from `TSCreate()`
4300: Output Parameters:
4301: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`)
4302: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`)
4303: . func - Function to compute the Jacobian of the RHS (or `NULL`)
4304: - ctx - User-defined context for Jacobian evaluation routine (or `NULL`)
4306: Level: intermediate
4308: Note:
4309: You can pass in `NULL` for any return argument you do not need.
4311: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4313: @*/
4314: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, void **ctx)
4315: {
4316: DM dm;
4318: PetscFunctionBegin;
4319: if (Amat || Pmat) {
4320: SNES snes;
4321: PetscCall(TSGetSNES(ts, &snes));
4322: PetscCall(SNESSetUpMatrices(snes));
4323: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4324: }
4325: PetscCall(TSGetDM(ts, &dm));
4326: PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4327: PetscFunctionReturn(PETSC_SUCCESS);
4328: }
4330: /*@C
4331: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4333: Not Collective, but parallel objects are returned if ts is parallel
4335: Input Parameter:
4336: . ts - The `TS` context obtained from `TSCreate()`
4338: Output Parameters:
4339: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4340: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4341: . f - The function to compute the matrices
4342: - ctx - User-defined context for Jacobian evaluation routine
4344: Level: advanced
4346: Note:
4347: You can pass in `NULL` for any return argument you do not need.
4349: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4350: @*/
4351: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, void **ctx)
4352: {
4353: DM dm;
4355: PetscFunctionBegin;
4356: if (Amat || Pmat) {
4357: SNES snes;
4358: PetscCall(TSGetSNES(ts, &snes));
4359: PetscCall(SNESSetUpMatrices(snes));
4360: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4361: }
4362: PetscCall(TSGetDM(ts, &dm));
4363: PetscCall(DMTSGetIJacobian(dm, f, ctx));
4364: PetscFunctionReturn(PETSC_SUCCESS);
4365: }
4367: #include <petsc/private/dmimpl.h>
4368: /*@
4369: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4371: Logically Collective
4373: Input Parameters:
4374: + ts - the `TS` integrator object
4375: - dm - the dm, cannot be `NULL`
4377: Level: intermediate
4379: Notes:
4380: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4381: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4382: different problems using the same function space.
4384: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4385: @*/
4386: PetscErrorCode TSSetDM(TS ts, DM dm)
4387: {
4388: SNES snes;
4389: DMTS tsdm;
4391: PetscFunctionBegin;
4394: PetscCall(PetscObjectReference((PetscObject)dm));
4395: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4396: if (ts->dm->dmts && !dm->dmts) {
4397: PetscCall(DMCopyDMTS(ts->dm, dm));
4398: PetscCall(DMGetDMTS(ts->dm, &tsdm));
4399: /* Grant write privileges to the replacement DM */
4400: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4401: }
4402: PetscCall(DMDestroy(&ts->dm));
4403: }
4404: ts->dm = dm;
4406: PetscCall(TSGetSNES(ts, &snes));
4407: PetscCall(SNESSetDM(snes, dm));
4408: PetscFunctionReturn(PETSC_SUCCESS);
4409: }
4411: /*@
4412: TSGetDM - Gets the `DM` that may be used by some preconditioners
4414: Not Collective
4416: Input Parameter:
4417: . ts - the `TS`
4419: Output Parameter:
4420: . dm - the `DM`
4422: Level: intermediate
4424: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4425: @*/
4426: PetscErrorCode TSGetDM(TS ts, DM *dm)
4427: {
4428: PetscFunctionBegin;
4430: if (!ts->dm) {
4431: PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4432: if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4433: }
4434: *dm = ts->dm;
4435: PetscFunctionReturn(PETSC_SUCCESS);
4436: }
4438: /*@
4439: SNESTSFormFunction - Function to evaluate nonlinear residual
4441: Logically Collective
4443: Input Parameters:
4444: + snes - nonlinear solver
4445: . U - the current state at which to evaluate the residual
4446: - ctx - user context, must be a TS
4448: Output Parameter:
4449: . F - the nonlinear residual
4451: Level: advanced
4453: Note:
4454: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4455: It is most frequently passed to `MatFDColoringSetFunction()`.
4457: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4458: @*/
4459: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4460: {
4461: TS ts = (TS)ctx;
4463: PetscFunctionBegin;
4468: PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4469: PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4470: PetscFunctionReturn(PETSC_SUCCESS);
4471: }
4473: /*@
4474: SNESTSFormJacobian - Function to evaluate the Jacobian
4476: Collective
4478: Input Parameters:
4479: + snes - nonlinear solver
4480: . U - the current state at which to evaluate the residual
4481: - ctx - user context, must be a `TS`
4483: Output Parameters:
4484: + A - the Jacobian
4485: - B - the preconditioning matrix (may be the same as A)
4487: Level: developer
4489: Note:
4490: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4492: .seealso: [](ch_ts), `SNESSetJacobian()`
4493: @*/
4494: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4495: {
4496: TS ts = (TS)ctx;
4498: PetscFunctionBegin;
4504: PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4505: PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4506: PetscFunctionReturn(PETSC_SUCCESS);
4507: }
4509: /*@C
4510: TSComputeRHSFunctionLinear - Evaluate the right-hand side via the user-provided Jacobian, for linear problems Udot = A U only
4512: Collective
4514: Input Parameters:
4515: + ts - time stepping context
4516: . t - time at which to evaluate
4517: . U - state at which to evaluate
4518: - ctx - context
4520: Output Parameter:
4521: . F - right-hand side
4523: Level: intermediate
4525: Note:
4526: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right-hand side for linear problems.
4527: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4529: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4530: @*/
4531: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4532: {
4533: Mat Arhs, Brhs;
4535: PetscFunctionBegin;
4536: PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4537: /* undo the damage caused by shifting */
4538: PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4539: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4540: PetscCall(MatMult(Arhs, U, F));
4541: PetscFunctionReturn(PETSC_SUCCESS);
4542: }
4544: /*@C
4545: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
4547: Collective
4549: Input Parameters:
4550: + ts - time stepping context
4551: . t - time at which to evaluate
4552: . U - state at which to evaluate
4553: - ctx - context
4555: Output Parameters:
4556: + A - pointer to operator
4557: - B - pointer to preconditioning matrix
4559: Level: intermediate
4561: Note:
4562: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4564: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4565: @*/
4566: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4567: {
4568: PetscFunctionBegin;
4569: PetscFunctionReturn(PETSC_SUCCESS);
4570: }
4572: /*@C
4573: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4575: Collective
4577: Input Parameters:
4578: + ts - time stepping context
4579: . t - time at which to evaluate
4580: . U - state at which to evaluate
4581: . Udot - time derivative of state vector
4582: - ctx - context
4584: Output Parameter:
4585: . F - left hand side
4587: Level: intermediate
4589: Notes:
4590: 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
4591: user is required to write their own `TSComputeIFunction()`.
4592: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4593: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4595: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4597: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4598: @*/
4599: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4600: {
4601: Mat A, B;
4603: PetscFunctionBegin;
4604: PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4605: PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4606: PetscCall(MatMult(A, Udot, F));
4607: PetscFunctionReturn(PETSC_SUCCESS);
4608: }
4610: /*@C
4611: TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE
4613: Collective
4615: Input Parameters:
4616: + ts - time stepping context
4617: . t - time at which to evaluate
4618: . U - state at which to evaluate
4619: . Udot - time derivative of state vector
4620: . shift - shift to apply
4621: - ctx - context
4623: Output Parameters:
4624: + A - pointer to operator
4625: - B - pointer to matrix from which the preconditioner is built (often `A`)
4627: Level: advanced
4629: Notes:
4630: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4632: It is only appropriate for problems of the form
4634: $$
4635: M \dot{U} = F(U,t)
4636: $$
4638: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4639: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4640: an implicit operator of the form
4642: $$
4643: shift*M + J
4644: $$
4646: 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
4647: a copy of M or reassemble it when requested.
4649: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4650: @*/
4651: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4652: {
4653: PetscFunctionBegin;
4654: PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4655: ts->ijacobian.shift = shift;
4656: PetscFunctionReturn(PETSC_SUCCESS);
4657: }
4659: /*@
4660: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4662: Not Collective
4664: Input Parameter:
4665: . ts - the `TS` context
4667: Output Parameter:
4668: . equation_type - see `TSEquationType`
4670: Level: beginner
4672: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4673: @*/
4674: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4675: {
4676: PetscFunctionBegin;
4678: PetscAssertPointer(equation_type, 2);
4679: *equation_type = ts->equation_type;
4680: PetscFunctionReturn(PETSC_SUCCESS);
4681: }
4683: /*@
4684: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4686: Not Collective
4688: Input Parameters:
4689: + ts - the `TS` context
4690: - equation_type - see `TSEquationType`
4692: Level: advanced
4694: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4695: @*/
4696: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4697: {
4698: PetscFunctionBegin;
4700: ts->equation_type = equation_type;
4701: PetscFunctionReturn(PETSC_SUCCESS);
4702: }
4704: /*@
4705: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4707: Not Collective
4709: Input Parameter:
4710: . ts - the `TS` context
4712: Output Parameter:
4713: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4714: manual pages for the individual convergence tests for complete lists
4716: Level: beginner
4718: Note:
4719: Can only be called after the call to `TSSolve()` is complete.
4721: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4722: @*/
4723: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4724: {
4725: PetscFunctionBegin;
4727: PetscAssertPointer(reason, 2);
4728: *reason = ts->reason;
4729: PetscFunctionReturn(PETSC_SUCCESS);
4730: }
4732: /*@
4733: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4735: Logically Collective; reason must contain common value
4737: Input Parameters:
4738: + ts - the `TS` context
4739: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4740: manual pages for the individual convergence tests for complete lists
4742: Level: advanced
4744: Note:
4745: Can only be called while `TSSolve()` is active.
4747: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4748: @*/
4749: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4750: {
4751: PetscFunctionBegin;
4753: ts->reason = reason;
4754: PetscFunctionReturn(PETSC_SUCCESS);
4755: }
4757: /*@
4758: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4760: Not Collective
4762: Input Parameter:
4763: . ts - the `TS` context
4765: Output Parameter:
4766: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4768: Level: beginner
4770: Note:
4771: Can only be called after the call to `TSSolve()` is complete.
4773: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4774: @*/
4775: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4776: {
4777: PetscFunctionBegin;
4779: PetscAssertPointer(ftime, 2);
4780: *ftime = ts->solvetime;
4781: PetscFunctionReturn(PETSC_SUCCESS);
4782: }
4784: /*@
4785: TSGetSNESIterations - Gets the total number of nonlinear iterations
4786: used by the time integrator.
4788: Not Collective
4790: Input Parameter:
4791: . ts - `TS` context
4793: Output Parameter:
4794: . nits - number of nonlinear iterations
4796: Level: intermediate
4798: Note:
4799: This counter is reset to zero for each successive call to `TSSolve()`.
4801: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4802: @*/
4803: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4804: {
4805: PetscFunctionBegin;
4807: PetscAssertPointer(nits, 2);
4808: *nits = ts->snes_its;
4809: PetscFunctionReturn(PETSC_SUCCESS);
4810: }
4812: /*@
4813: TSGetKSPIterations - Gets the total number of linear iterations
4814: used by the time integrator.
4816: Not Collective
4818: Input Parameter:
4819: . ts - `TS` context
4821: Output Parameter:
4822: . lits - number of linear iterations
4824: Level: intermediate
4826: Note:
4827: This counter is reset to zero for each successive call to `TSSolve()`.
4829: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()`
4830: @*/
4831: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4832: {
4833: PetscFunctionBegin;
4835: PetscAssertPointer(lits, 2);
4836: *lits = ts->ksp_its;
4837: PetscFunctionReturn(PETSC_SUCCESS);
4838: }
4840: /*@
4841: TSGetStepRejections - Gets the total number of rejected steps.
4843: Not Collective
4845: Input Parameter:
4846: . ts - `TS` context
4848: Output Parameter:
4849: . rejects - number of steps rejected
4851: Level: intermediate
4853: Note:
4854: This counter is reset to zero for each successive call to `TSSolve()`.
4856: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4857: @*/
4858: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4859: {
4860: PetscFunctionBegin;
4862: PetscAssertPointer(rejects, 2);
4863: *rejects = ts->reject;
4864: PetscFunctionReturn(PETSC_SUCCESS);
4865: }
4867: /*@
4868: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
4870: Not Collective
4872: Input Parameter:
4873: . ts - `TS` context
4875: Output Parameter:
4876: . fails - number of failed nonlinear solves
4878: Level: intermediate
4880: Note:
4881: This counter is reset to zero for each successive call to `TSSolve()`.
4883: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4884: @*/
4885: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4886: {
4887: PetscFunctionBegin;
4889: PetscAssertPointer(fails, 2);
4890: *fails = ts->num_snes_failures;
4891: PetscFunctionReturn(PETSC_SUCCESS);
4892: }
4894: /*@
4895: TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails
4897: Not Collective
4899: Input Parameters:
4900: + ts - `TS` context
4901: - rejects - maximum number of rejected steps, pass `PETSC_UNLIMITED` for unlimited
4903: Options Database Key:
4904: . -ts_max_reject - Maximum number of step rejections before a step fails
4906: Level: intermediate
4908: Developer Note:
4909: The options database name is incorrect.
4911: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4912: @*/
4913: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4914: {
4915: PetscFunctionBegin;
4917: if (rejects == PETSC_UNLIMITED || rejects == -1) {
4918: ts->max_reject = PETSC_UNLIMITED;
4919: } else {
4920: PetscCheck(rejects >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of rejections");
4921: ts->max_reject = rejects;
4922: }
4923: PetscFunctionReturn(PETSC_SUCCESS);
4924: }
4926: /*@
4927: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves
4929: Not Collective
4931: Input Parameters:
4932: + ts - `TS` context
4933: - fails - maximum number of failed nonlinear solves, pass `PETSC_UNLIMITED` to allow any number of failures.
4935: Options Database Key:
4936: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
4938: Level: intermediate
4940: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4941: @*/
4942: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4943: {
4944: PetscFunctionBegin;
4946: if (fails == PETSC_UNLIMITED || fails == -1) {
4947: ts->max_snes_failures = PETSC_UNLIMITED;
4948: } else {
4949: PetscCheck(fails >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures");
4950: ts->max_snes_failures = fails;
4951: }
4952: PetscFunctionReturn(PETSC_SUCCESS);
4953: }
4955: /*@
4956: TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`
4958: Not Collective
4960: Input Parameters:
4961: + ts - `TS` context
4962: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
4964: Options Database Key:
4965: . -ts_error_if_step_fails - Error if no step succeeds
4967: Level: intermediate
4969: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
4970: @*/
4971: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
4972: {
4973: PetscFunctionBegin;
4975: ts->errorifstepfailed = err;
4976: PetscFunctionReturn(PETSC_SUCCESS);
4977: }
4979: /*@
4980: TSGetAdapt - Get the adaptive controller context for the current method
4982: Collective if controller has not yet been created
4984: Input Parameter:
4985: . ts - time stepping context
4987: Output Parameter:
4988: . adapt - adaptive controller
4990: Level: intermediate
4992: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
4993: @*/
4994: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
4995: {
4996: PetscFunctionBegin;
4998: PetscAssertPointer(adapt, 2);
4999: if (!ts->adapt) {
5000: PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
5001: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
5002: }
5003: *adapt = ts->adapt;
5004: PetscFunctionReturn(PETSC_SUCCESS);
5005: }
5007: /*@
5008: TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller
5010: Logically Collective
5012: Input Parameters:
5013: + ts - time integration context
5014: . atol - scalar absolute tolerances
5015: . vatol - vector of absolute tolerances or `NULL`, used in preference to `atol` if present
5016: . rtol - scalar relative tolerances
5017: - vrtol - vector of relative tolerances or `NULL`, used in preference to `rtol` if present
5019: Options Database Keys:
5020: + -ts_rtol <rtol> - relative tolerance for local truncation error
5021: - -ts_atol <atol> - Absolute tolerance for local truncation error
5023: Level: beginner
5025: Notes:
5026: `PETSC_CURRENT` or `PETSC_DETERMINE` may be used for `atol` or `rtol` to indicate the current value
5027: or the default value from when the object's type was set.
5029: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5030: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5031: computed only for the differential or the algebraic part then this can be done using the vector of
5032: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
5033: differential part and infinity for the algebraic part, the LTE calculation will include only the
5034: differential variables.
5036: Fortran Note:
5037: Use `PETSC_CURRENT_INTEGER` or `PETSC_DETERMINE_INTEGER`.
5039: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
5040: @*/
5041: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
5042: {
5043: PetscFunctionBegin;
5044: if (atol == (PetscReal)PETSC_DETERMINE) {
5045: ts->atol = ts->default_atol;
5046: } else if (atol != (PetscReal)PETSC_CURRENT) {
5047: PetscCheck(atol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)atol);
5048: ts->atol = atol;
5049: }
5051: if (vatol) {
5052: PetscCall(PetscObjectReference((PetscObject)vatol));
5053: PetscCall(VecDestroy(&ts->vatol));
5054: ts->vatol = vatol;
5055: }
5057: if (rtol == (PetscReal)PETSC_DETERMINE) {
5058: ts->rtol = ts->default_rtol;
5059: } else if (rtol != (PetscReal)PETSC_CURRENT) {
5060: PetscCheck(rtol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative", (double)rtol);
5061: ts->rtol = rtol;
5062: }
5064: if (vrtol) {
5065: PetscCall(PetscObjectReference((PetscObject)vrtol));
5066: PetscCall(VecDestroy(&ts->vrtol));
5067: ts->vrtol = vrtol;
5068: }
5069: PetscFunctionReturn(PETSC_SUCCESS);
5070: }
5072: /*@
5073: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
5075: Logically Collective
5077: Input Parameter:
5078: . ts - time integration context
5080: Output Parameters:
5081: + atol - scalar absolute tolerances, `NULL` to ignore
5082: . vatol - vector of absolute tolerances, `NULL` to ignore
5083: . rtol - scalar relative tolerances, `NULL` to ignore
5084: - vrtol - vector of relative tolerances, `NULL` to ignore
5086: Level: beginner
5088: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
5089: @*/
5090: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
5091: {
5092: PetscFunctionBegin;
5093: if (atol) *atol = ts->atol;
5094: if (vatol) *vatol = ts->vatol;
5095: if (rtol) *rtol = ts->rtol;
5096: if (vrtol) *vrtol = ts->vrtol;
5097: PetscFunctionReturn(PETSC_SUCCESS);
5098: }
5100: /*@
5101: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
5103: Collective
5105: Input Parameters:
5106: + ts - time stepping context
5107: . U - state vector, usually ts->vec_sol
5108: . Y - state vector to be compared to U
5109: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5111: Output Parameters:
5112: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5113: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5114: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5116: Options Database Key:
5117: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5119: Level: developer
5121: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5122: @*/
5123: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5124: {
5125: PetscInt norma_loc, norm_loc, normr_loc;
5127: PetscFunctionBegin;
5129: 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));
5130: if (wnormtype == NORM_2) {
5131: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5132: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5133: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5134: }
5135: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5136: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5137: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5138: PetscFunctionReturn(PETSC_SUCCESS);
5139: }
5141: /*@
5142: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5144: Collective
5146: Input Parameters:
5147: + ts - time stepping context
5148: . E - error vector
5149: . U - state vector, usually ts->vec_sol
5150: . Y - state vector, previous time step
5151: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5153: Output Parameters:
5154: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5155: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5156: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5158: Options Database Key:
5159: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5161: Level: developer
5163: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5164: @*/
5165: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5166: {
5167: PetscInt norma_loc, norm_loc, normr_loc;
5169: PetscFunctionBegin;
5171: 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));
5172: if (wnormtype == NORM_2) {
5173: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5174: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5175: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5176: }
5177: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5178: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5179: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5180: PetscFunctionReturn(PETSC_SUCCESS);
5181: }
5183: /*@
5184: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5186: Logically Collective
5188: Input Parameters:
5189: + ts - time stepping context
5190: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5192: Note:
5193: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5195: Level: intermediate
5197: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5198: @*/
5199: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5200: {
5201: PetscFunctionBegin;
5203: ts->cfltime_local = cfltime;
5204: ts->cfltime = -1.;
5205: PetscFunctionReturn(PETSC_SUCCESS);
5206: }
5208: /*@
5209: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5211: Collective
5213: Input Parameter:
5214: . ts - time stepping context
5216: Output Parameter:
5217: . cfltime - maximum stable time step for forward Euler
5219: Level: advanced
5221: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5222: @*/
5223: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5224: {
5225: PetscFunctionBegin;
5226: if (ts->cfltime < 0) PetscCallMPI(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5227: *cfltime = ts->cfltime;
5228: PetscFunctionReturn(PETSC_SUCCESS);
5229: }
5231: /*@
5232: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5234: Input Parameters:
5235: + ts - the `TS` context.
5236: . xl - lower bound.
5237: - xu - upper bound.
5239: Level: advanced
5241: Note:
5242: If this routine is not called then the lower and upper bounds are set to
5243: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5245: .seealso: [](ch_ts), `TS`
5246: @*/
5247: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5248: {
5249: SNES snes;
5251: PetscFunctionBegin;
5252: PetscCall(TSGetSNES(ts, &snes));
5253: PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5254: PetscFunctionReturn(PETSC_SUCCESS);
5255: }
5257: /*@
5258: TSComputeLinearStability - computes the linear stability function at a point
5260: Collective
5262: Input Parameters:
5263: + ts - the `TS` context
5264: . xr - real part of input argument
5265: - xi - imaginary part of input argument
5267: Output Parameters:
5268: + yr - real part of function value
5269: - yi - imaginary part of function value
5271: Level: developer
5273: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5274: @*/
5275: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5276: {
5277: PetscFunctionBegin;
5279: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5280: PetscFunctionReturn(PETSC_SUCCESS);
5281: }
5283: /*@
5284: TSRestartStep - Flags the solver to restart the next step
5286: Collective
5288: Input Parameter:
5289: . ts - the `TS` context obtained from `TSCreate()`
5291: Level: advanced
5293: Notes:
5294: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5295: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5296: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5297: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5298: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5299: discontinuous source terms).
5301: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5302: @*/
5303: PetscErrorCode TSRestartStep(TS ts)
5304: {
5305: PetscFunctionBegin;
5307: ts->steprestart = PETSC_TRUE;
5308: PetscFunctionReturn(PETSC_SUCCESS);
5309: }
5311: /*@
5312: TSRollBack - Rolls back one time step
5314: Collective
5316: Input Parameter:
5317: . ts - the `TS` context obtained from `TSCreate()`
5319: Level: advanced
5321: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5322: @*/
5323: PetscErrorCode TSRollBack(TS ts)
5324: {
5325: PetscFunctionBegin;
5327: PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5328: PetscTryTypeMethod(ts, rollback);
5329: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
5330: ts->time_step = ts->ptime - ts->ptime_prev;
5331: ts->ptime = ts->ptime_prev;
5332: ts->ptime_prev = ts->ptime_prev_rollback;
5333: ts->steps--;
5334: ts->steprollback = PETSC_TRUE;
5335: PetscFunctionReturn(PETSC_SUCCESS);
5336: }
5338: /*@
5339: TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step
5341: Not collective
5343: Input Parameter:
5344: . ts - the `TS` context obtained from `TSCreate()`
5346: Output Parameter:
5347: . flg - the rollback flag
5349: Level: advanced
5351: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
5352: @*/
5353: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
5354: {
5355: PetscFunctionBegin;
5357: PetscAssertPointer(flg, 2);
5358: *flg = ts->steprollback;
5359: PetscFunctionReturn(PETSC_SUCCESS);
5360: }
5362: /*@
5363: TSGetStepResize - Get the internal flag indicating if the current step is after a resize.
5365: Not collective
5367: Input Parameter:
5368: . ts - the `TS` context obtained from `TSCreate()`
5370: Output Parameter:
5371: . flg - the resize flag
5373: Level: advanced
5375: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetResize()`
5376: @*/
5377: PetscErrorCode TSGetStepResize(TS ts, PetscBool *flg)
5378: {
5379: PetscFunctionBegin;
5381: PetscAssertPointer(flg, 2);
5382: *flg = ts->stepresize;
5383: PetscFunctionReturn(PETSC_SUCCESS);
5384: }
5386: /*@
5387: TSGetStages - Get the number of stages and stage values
5389: Input Parameter:
5390: . ts - the `TS` context obtained from `TSCreate()`
5392: Output Parameters:
5393: + ns - the number of stages
5394: - Y - the current stage vectors
5396: Level: advanced
5398: Note:
5399: Both `ns` and `Y` can be `NULL`.
5401: .seealso: [](ch_ts), `TS`, `TSCreate()`
5402: @*/
5403: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5404: {
5405: PetscFunctionBegin;
5407: if (ns) PetscAssertPointer(ns, 2);
5408: if (Y) PetscAssertPointer(Y, 3);
5409: if (!ts->ops->getstages) {
5410: if (ns) *ns = 0;
5411: if (Y) *Y = NULL;
5412: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5413: PetscFunctionReturn(PETSC_SUCCESS);
5414: }
5416: /*@C
5417: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5419: Collective
5421: Input Parameters:
5422: + ts - the `TS` context
5423: . t - current timestep
5424: . U - state vector
5425: . Udot - time derivative of state vector
5426: . shift - shift to apply, see note below
5427: - ctx - an optional user context
5429: Output Parameters:
5430: + J - Jacobian matrix (not altered in this routine)
5431: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
5433: Level: intermediate
5435: Notes:
5436: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5438: dF/dU + shift*dF/dUdot
5440: Most users should not need to explicitly call this routine, as it
5441: is used internally within the nonlinear solvers.
5443: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5444: routine, then it will try to get the coloring from the matrix. This requires that the
5445: matrix have nonzero entries precomputed.
5447: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5448: @*/
5449: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5450: {
5451: SNES snes;
5452: MatFDColoring color;
5453: PetscBool hascolor, matcolor = PETSC_FALSE;
5455: PetscFunctionBegin;
5456: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5457: PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5458: if (!color) {
5459: DM dm;
5460: ISColoring iscoloring;
5462: PetscCall(TSGetDM(ts, &dm));
5463: PetscCall(DMHasColoring(dm, &hascolor));
5464: if (hascolor && !matcolor) {
5465: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5466: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5467: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode (*)(void))SNESTSFormFunction, (void *)ts));
5468: PetscCall(MatFDColoringSetFromOptions(color));
5469: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5470: PetscCall(ISColoringDestroy(&iscoloring));
5471: } else {
5472: MatColoring mc;
5474: PetscCall(MatColoringCreate(B, &mc));
5475: PetscCall(MatColoringSetDistance(mc, 2));
5476: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5477: PetscCall(MatColoringSetFromOptions(mc));
5478: PetscCall(MatColoringApply(mc, &iscoloring));
5479: PetscCall(MatColoringDestroy(&mc));
5480: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5481: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode (*)(void))SNESTSFormFunction, (void *)ts));
5482: PetscCall(MatFDColoringSetFromOptions(color));
5483: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5484: PetscCall(ISColoringDestroy(&iscoloring));
5485: }
5486: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5487: PetscCall(PetscObjectDereference((PetscObject)color));
5488: }
5489: PetscCall(TSGetSNES(ts, &snes));
5490: PetscCall(MatFDColoringApply(B, color, U, snes));
5491: if (J != B) {
5492: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5493: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5494: }
5495: PetscFunctionReturn(PETSC_SUCCESS);
5496: }
5498: /*@C
5499: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5501: Input Parameters:
5502: + ts - the `TS` context
5503: - func - function called within `TSFunctionDomainError()`
5505: Calling sequence of `func`:
5506: + ts - the `TS` context
5507: . time - the current time (of the stage)
5508: . state - the state to check if it is valid
5509: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable
5511: Level: intermediate
5513: Notes:
5514: If an implicit ODE solver is being used then, in addition to providing this routine, the
5515: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5516: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5517: Use `TSGetSNES()` to obtain the `SNES` object
5519: Developer Notes:
5520: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5521: since one takes a function pointer and the other does not.
5523: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5524: @*/
5525: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5526: {
5527: PetscFunctionBegin;
5529: ts->functiondomainerror = func;
5530: PetscFunctionReturn(PETSC_SUCCESS);
5531: }
5533: /*@
5534: TSFunctionDomainError - Checks if the current state is valid
5536: Input Parameters:
5537: + ts - the `TS` context
5538: . stagetime - time of the simulation
5539: - Y - state vector to check.
5541: Output Parameter:
5542: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5544: Level: developer
5546: Note:
5547: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5548: to check if the current state is valid.
5550: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5551: @*/
5552: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5553: {
5554: PetscFunctionBegin;
5556: *accept = PETSC_TRUE;
5557: if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5558: PetscFunctionReturn(PETSC_SUCCESS);
5559: }
5561: /*@
5562: TSClone - This function clones a time step `TS` object.
5564: Collective
5566: Input Parameter:
5567: . tsin - The input `TS`
5569: Output Parameter:
5570: . tsout - The output `TS` (cloned)
5572: Level: developer
5574: Notes:
5575: 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.
5576: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5578: 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
5579: .vb
5580: SNES snes_dup = NULL;
5581: TSGetSNES(ts,&snes_dup);
5582: TSSetSNES(ts,snes_dup);
5583: .ve
5585: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5586: @*/
5587: PetscErrorCode TSClone(TS tsin, TS *tsout)
5588: {
5589: TS t;
5590: SNES snes_start;
5591: DM dm;
5592: TSType type;
5594: PetscFunctionBegin;
5595: PetscAssertPointer(tsin, 1);
5596: *tsout = NULL;
5598: PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));
5600: /* General TS description */
5601: t->numbermonitors = 0;
5602: t->monitorFrequency = 1;
5603: t->setupcalled = 0;
5604: t->ksp_its = 0;
5605: t->snes_its = 0;
5606: t->nwork = 0;
5607: t->rhsjacobian.time = PETSC_MIN_REAL;
5608: t->rhsjacobian.scale = 1.;
5609: t->ijacobian.shift = 1.;
5611: PetscCall(TSGetSNES(tsin, &snes_start));
5612: PetscCall(TSSetSNES(t, snes_start));
5614: PetscCall(TSGetDM(tsin, &dm));
5615: PetscCall(TSSetDM(t, dm));
5617: t->adapt = tsin->adapt;
5618: PetscCall(PetscObjectReference((PetscObject)t->adapt));
5620: t->trajectory = tsin->trajectory;
5621: PetscCall(PetscObjectReference((PetscObject)t->trajectory));
5623: t->event = tsin->event;
5624: if (t->event) t->event->refct++;
5626: t->problem_type = tsin->problem_type;
5627: t->ptime = tsin->ptime;
5628: t->ptime_prev = tsin->ptime_prev;
5629: t->time_step = tsin->time_step;
5630: t->max_time = tsin->max_time;
5631: t->steps = tsin->steps;
5632: t->max_steps = tsin->max_steps;
5633: t->equation_type = tsin->equation_type;
5634: t->atol = tsin->atol;
5635: t->rtol = tsin->rtol;
5636: t->max_snes_failures = tsin->max_snes_failures;
5637: t->max_reject = tsin->max_reject;
5638: t->errorifstepfailed = tsin->errorifstepfailed;
5640: PetscCall(TSGetType(tsin, &type));
5641: PetscCall(TSSetType(t, type));
5643: t->vec_sol = NULL;
5645: t->cfltime = tsin->cfltime;
5646: t->cfltime_local = tsin->cfltime_local;
5647: t->exact_final_time = tsin->exact_final_time;
5649: t->ops[0] = tsin->ops[0];
5651: if (((PetscObject)tsin)->fortran_func_pointers) {
5652: PetscInt i;
5653: PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers));
5654: for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5655: }
5656: *tsout = t;
5657: PetscFunctionReturn(PETSC_SUCCESS);
5658: }
5660: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5661: {
5662: TS ts = (TS)ctx;
5664: PetscFunctionBegin;
5665: PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5666: PetscFunctionReturn(PETSC_SUCCESS);
5667: }
5669: /*@
5670: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5672: Logically Collective
5674: Input Parameter:
5675: . ts - the time stepping routine
5677: Output Parameter:
5678: . flg - `PETSC_TRUE` if the multiply is likely correct
5680: Options Database Key:
5681: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5683: Level: advanced
5685: Note:
5686: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5688: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5689: @*/
5690: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5691: {
5692: Mat J, B;
5693: TSRHSJacobianFn *func;
5694: void *ctx;
5696: PetscFunctionBegin;
5697: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5698: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5699: PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5700: PetscFunctionReturn(PETSC_SUCCESS);
5701: }
5703: /*@
5704: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5706: Logically Collective
5708: Input Parameter:
5709: . ts - the time stepping routine
5711: Output Parameter:
5712: . flg - `PETSC_TRUE` if the multiply is likely correct
5714: Options Database Key:
5715: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5717: Level: advanced
5719: Notes:
5720: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5722: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5723: @*/
5724: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5725: {
5726: Mat J, B;
5727: void *ctx;
5728: TSRHSJacobianFn *func;
5730: PetscFunctionBegin;
5731: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5732: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5733: PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5734: PetscFunctionReturn(PETSC_SUCCESS);
5735: }
5737: /*@
5738: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
5740: Logically Collective
5742: Input Parameters:
5743: + ts - timestepping context
5744: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5746: Options Database Key:
5747: . -ts_use_splitrhsfunction - <true,false>
5749: Level: intermediate
5751: Note:
5752: This is only for multirate methods
5754: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5755: @*/
5756: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5757: {
5758: PetscFunctionBegin;
5760: ts->use_splitrhsfunction = use_splitrhsfunction;
5761: PetscFunctionReturn(PETSC_SUCCESS);
5762: }
5764: /*@
5765: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
5767: Not Collective
5769: Input Parameter:
5770: . ts - timestepping context
5772: Output Parameter:
5773: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5775: Level: intermediate
5777: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5778: @*/
5779: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5780: {
5781: PetscFunctionBegin;
5783: *use_splitrhsfunction = ts->use_splitrhsfunction;
5784: PetscFunctionReturn(PETSC_SUCCESS);
5785: }
5787: /*@
5788: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
5790: Logically Collective
5792: Input Parameters:
5793: + ts - the time-stepper
5794: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
5796: Level: intermediate
5798: Note:
5799: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
5801: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
5802: @*/
5803: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5804: {
5805: PetscFunctionBegin;
5807: ts->axpy_pattern = str;
5808: PetscFunctionReturn(PETSC_SUCCESS);
5809: }
5811: /*@
5812: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
5814: Collective
5816: Input Parameters:
5817: + ts - the time-stepper
5818: . n - number of the time points (>=2)
5819: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5821: Options Database Key:
5822: . -ts_time_span <t0,...tf> - Sets the time span
5824: Level: intermediate
5826: Notes:
5827: The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
5828: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
5829: The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may
5830: pressure the memory system when using a large number of span points.
5832: .seealso: [](ch_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()`
5833: @*/
5834: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
5835: {
5836: PetscFunctionBegin;
5838: PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
5839: if (ts->tspan && n != ts->tspan->num_span_times) {
5840: PetscCall(PetscFree(ts->tspan->span_times));
5841: PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
5842: PetscCall(PetscMalloc1(n, &ts->tspan->span_times));
5843: }
5844: if (!ts->tspan) {
5845: TSTimeSpan tspan;
5846: PetscCall(PetscNew(&tspan));
5847: PetscCall(PetscMalloc1(n, &tspan->span_times));
5848: tspan->reltol = 1e-6;
5849: tspan->abstol = 10 * PETSC_MACHINE_EPSILON;
5850: tspan->worktol = 0;
5851: ts->tspan = tspan;
5852: }
5853: ts->tspan->num_span_times = n;
5854: PetscCall(PetscArraycpy(ts->tspan->span_times, span_times, n));
5855: PetscCall(TSSetTime(ts, ts->tspan->span_times[0]));
5856: PetscCall(TSSetMaxTime(ts, ts->tspan->span_times[n - 1]));
5857: PetscFunctionReturn(PETSC_SUCCESS);
5858: }
5860: /*@C
5861: TSGetTimeSpan - gets the time span set with `TSSetTimeSpan()`
5863: Not Collective
5865: Input Parameter:
5866: . ts - the time-stepper
5868: Output Parameters:
5869: + n - number of the time points (>=2)
5870: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5872: Level: beginner
5874: Note:
5875: The values obtained are valid until the `TS` object is destroyed.
5877: Both `n` and `span_times` can be `NULL`.
5879: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()`
5880: @*/
5881: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal *span_times[])
5882: {
5883: PetscFunctionBegin;
5885: if (n) PetscAssertPointer(n, 2);
5886: if (span_times) PetscAssertPointer(span_times, 3);
5887: if (!ts->tspan) {
5888: if (n) *n = 0;
5889: if (span_times) *span_times = NULL;
5890: } else {
5891: if (n) *n = ts->tspan->num_span_times;
5892: if (span_times) *span_times = ts->tspan->span_times;
5893: }
5894: PetscFunctionReturn(PETSC_SUCCESS);
5895: }
5897: /*@
5898: TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.
5900: Input Parameter:
5901: . ts - the `TS` context obtained from `TSCreate()`
5903: Output Parameters:
5904: + nsol - the number of solutions
5905: - Sols - the solution vectors
5907: Level: intermediate
5909: Notes:
5910: Both `nsol` and `Sols` can be `NULL`.
5912: Some time points in the time span may be skipped by `TS` so that `nsol` is less than the number of points specified by `TSSetTimeSpan()`.
5913: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span.
5915: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`
5916: @*/
5917: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
5918: {
5919: PetscFunctionBegin;
5921: if (nsol) PetscAssertPointer(nsol, 2);
5922: if (Sols) PetscAssertPointer(Sols, 3);
5923: if (!ts->tspan) {
5924: if (nsol) *nsol = 0;
5925: if (Sols) *Sols = NULL;
5926: } else {
5927: if (nsol) *nsol = ts->tspan->spanctr;
5928: if (Sols) *Sols = ts->tspan->vecs_sol;
5929: }
5930: PetscFunctionReturn(PETSC_SUCCESS);
5931: }
5933: /*@
5934: TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
5936: Collective
5938: Input Parameters:
5939: + ts - the `TS` context
5940: . J - Jacobian matrix (not altered in this routine)
5941: - B - newly computed Jacobian matrix to use with preconditioner
5943: Level: intermediate
5945: Notes:
5946: This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
5947: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
5948: and multiple fields are involved.
5950: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
5951: structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
5952: usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
5953: `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.
5955: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5956: @*/
5957: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
5958: {
5959: MatColoring mc = NULL;
5960: ISColoring iscoloring = NULL;
5961: MatFDColoring matfdcoloring = NULL;
5963: PetscFunctionBegin;
5964: /* Generate new coloring after eliminating zeros in the matrix */
5965: PetscCall(MatEliminateZeros(B, PETSC_TRUE));
5966: PetscCall(MatColoringCreate(B, &mc));
5967: PetscCall(MatColoringSetDistance(mc, 2));
5968: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5969: PetscCall(MatColoringSetFromOptions(mc));
5970: PetscCall(MatColoringApply(mc, &iscoloring));
5971: PetscCall(MatColoringDestroy(&mc));
5972: /* Replace the old coloring with the new one */
5973: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
5974: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode (*)(void))SNESTSFormFunction, (void *)ts));
5975: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
5976: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
5977: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
5978: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
5979: PetscCall(ISColoringDestroy(&iscoloring));
5980: PetscFunctionReturn(PETSC_SUCCESS);
5981: }