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
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_envelope - determine maximum and minimum value of each component of the solution over the solution time
70: Level: beginner
72: Notes:
73: See `SNESSetFromOptions()` and `KSPSetFromOptions()` for how to control the nonlinear and linear solves used by the time-stepper.
75: Certain `SNES` options get reset for each new nonlinear solver, for example `-snes_lag_jacobian its` and `-snes_lag_preconditioner its`, in order
76: to retain them over the multiple nonlinear solves that `TS` uses you mush also provide `-snes_lag_jacobian_persists true` and
77: `-snes_lag_preconditioner_persists true`
79: Developer Notes:
80: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
82: .seealso: [](ch_ts), `TS`, `TSGetType()`
83: @*/
84: PetscErrorCode TSSetFromOptions(TS ts)
85: {
86: PetscBool opt, flg, tflg;
87: char monfilename[PETSC_MAX_PATH_LEN];
88: PetscReal time_step, tspan[100];
89: PetscInt nt = PETSC_STATIC_ARRAY_LENGTH(tspan);
90: TSExactFinalTimeOption eftopt;
91: char dir[16];
92: TSIFunctionFn *ifun;
93: const char *defaultType;
94: char typeName[256];
96: PetscFunctionBegin;
99: PetscCall(TSRegisterAll());
100: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
102: PetscObjectOptionsBegin((PetscObject)ts);
103: if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
104: else defaultType = ifun ? TSBEULER : TSEULER;
105: PetscCall(PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, 256, &opt));
106: if (opt) PetscCall(TSSetType(ts, typeName));
107: else PetscCall(TSSetType(ts, defaultType));
109: /* Handle generic TS options */
110: PetscCall(PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL));
111: PetscCall(PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL));
112: PetscCall(PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", tspan, &nt, &flg));
113: if (flg) PetscCall(TSSetTimeSpan(ts, nt, tspan));
114: PetscCall(PetscOptionsInt("-ts_max_steps", "Maximum number of time steps", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL));
115: PetscCall(PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL));
116: PetscCall(PetscOptionsReal("-ts_dt", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg));
117: if (flg) PetscCall(TSSetTimeStep(ts, time_step));
118: PetscCall(PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg));
119: if (flg) PetscCall(TSSetExactFinalTime(ts, eftopt));
120: PetscCall(PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, NULL));
121: PetscCall(PetscOptionsInt("-ts_max_reject", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, NULL));
122: PetscCall(PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL));
123: PetscCall(PetscOptionsReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL));
124: PetscCall(PetscOptionsReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL));
126: PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL));
127: 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));
128: PetscCall(PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL));
129: #if defined(PETSC_HAVE_SAWS)
130: {
131: PetscBool set;
132: flg = PETSC_FALSE;
133: PetscCall(PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set));
134: if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)ts, flg));
135: }
136: #endif
138: /* Monitor options */
139: PetscCall(PetscOptionsInt("-ts_monitor_frequency", "Number of time steps between monitor output", "TSMonitorSetFrequency", ts->monitorFrequency, &ts->monitorFrequency, NULL));
140: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL));
141: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL));
142: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, NULL));
143: PetscCall(TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL));
145: PetscCall(PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg));
146: if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)ts, monfilename));
148: PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt));
149: if (opt) {
150: PetscInt howoften = 1;
151: DM dm;
152: PetscBool net;
154: PetscCall(PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL));
155: PetscCall(TSGetDM(ts, &dm));
156: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net));
157: if (net) {
158: TSMonitorLGCtxNetwork ctx;
159: PetscCall(TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx));
160: PetscCall(TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxNetworkDestroy));
161: PetscCall(PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL));
162: } else {
163: TSMonitorLGCtx ctx;
164: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
165: PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
166: }
167: }
169: PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt));
170: if (opt) {
171: TSMonitorLGCtx ctx;
172: PetscInt howoften = 1;
174: PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL));
175: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
176: PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
177: }
178: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL));
180: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt));
181: if (opt) {
182: TSMonitorLGCtx ctx;
183: PetscInt howoften = 1;
185: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
186: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
187: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
188: }
189: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt));
190: if (opt) {
191: TSMonitorLGCtx ctx;
192: PetscInt howoften = 1;
194: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
195: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
196: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
197: ctx->semilogy = PETSC_TRUE;
198: }
200: PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt));
201: if (opt) {
202: TSMonitorLGCtx ctx;
203: PetscInt howoften = 1;
205: PetscCall(PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL));
206: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
207: PetscCall(TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
208: }
209: PetscCall(PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt));
210: if (opt) {
211: TSMonitorLGCtx ctx;
212: PetscInt howoften = 1;
214: PetscCall(PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL));
215: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
216: PetscCall(TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
217: }
218: PetscCall(PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt));
219: if (opt) {
220: TSMonitorSPEigCtx ctx;
221: PetscInt howoften = 1;
223: PetscCall(PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL));
224: PetscCall(TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
225: PetscCall(TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscErrorCode(*)(void **))TSMonitorSPEigCtxDestroy));
226: }
227: PetscCall(PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase space from the DMSwarm", "TSMonitorSPSwarm", &opt));
228: if (opt) {
229: TSMonitorSPCtx ctx;
230: PetscInt howoften = 1, retain = 0;
231: PetscBool phase = PETSC_TRUE, create = PETSC_TRUE, multispecies = PETSC_FALSE;
233: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
234: if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
235: create = PETSC_FALSE;
236: break;
237: }
238: if (create) {
239: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase space from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL));
240: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL));
241: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL));
242: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_multi_species", "Color particles by particle species", "TSMonitorSPSwarm", multispecies, &multispecies, NULL));
243: PetscCall(TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, multispecies, &ctx));
244: PetscCall(TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorSPCtxDestroy));
245: }
246: }
247: PetscCall(PetscOptionsName("-ts_monitor_hg_swarm", "Display particle histogram from the DMSwarm", "TSMonitorHGSwarm", &opt));
248: if (opt) {
249: TSMonitorHGCtx ctx;
250: PetscInt howoften = 1, Ns = 1;
251: PetscBool velocity = PETSC_FALSE, create = PETSC_TRUE;
253: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
254: if (ts->monitor[i] == TSMonitorHGSwarmSolution) {
255: create = PETSC_FALSE;
256: break;
257: }
258: if (create) {
259: DM sw, dm;
260: PetscInt Nc, Nb;
262: PetscCall(TSGetDM(ts, &sw));
263: PetscCall(DMSwarmGetCellDM(sw, &dm));
264: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc));
265: Nb = PetscMin(20, PetscMax(10, Nc));
266: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL));
267: PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL));
268: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL));
269: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL));
270: PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx));
271: PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorHGCtxDestroy));
272: }
273: }
274: opt = PETSC_FALSE;
275: PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt));
276: if (opt) {
277: TSMonitorDrawCtx ctx;
278: PetscInt howoften = 1;
280: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL));
281: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
282: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
283: }
284: opt = PETSC_FALSE;
285: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt));
286: if (opt) {
287: TSMonitorDrawCtx ctx;
288: PetscReal bounds[4];
289: PetscInt n = 4;
290: PetscDraw draw;
291: PetscDrawAxis axis;
293: PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL));
294: PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field");
295: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx));
296: PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw));
297: PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis));
298: PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]));
299: PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2"));
300: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
301: }
302: opt = PETSC_FALSE;
303: PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt));
304: if (opt) {
305: TSMonitorDrawCtx ctx;
306: PetscInt howoften = 1;
308: PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL));
309: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
310: PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
311: }
312: opt = PETSC_FALSE;
313: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt));
314: if (opt) {
315: TSMonitorDrawCtx ctx;
316: PetscInt howoften = 1;
318: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL));
319: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
320: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
321: }
323: opt = PETSC_FALSE;
324: 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));
325: if (flg) {
326: const char *ptr = NULL, *ptr2 = NULL;
327: char *filetemplate;
328: PetscCheck(monfilename[0], PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03" PetscInt_FMT ".vts");
329: /* Do some cursory validation of the input. */
330: PetscCall(PetscStrstr(monfilename, "%", (char **)&ptr));
331: PetscCheck(ptr, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03" PetscInt_FMT ".vts");
332: for (ptr++; ptr && *ptr; ptr++) {
333: PetscCall(PetscStrchr("DdiouxX", *ptr, (char **)&ptr2));
334: PetscCheck(ptr2 || (*ptr >= '0' && *ptr <= '9'), PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Invalid file template argument to -ts_monitor_solution_vtk, should look like filename-%%03" PetscInt_FMT ".vts");
335: if (ptr2) break;
336: }
337: PetscCall(PetscStrallocpy(monfilename, &filetemplate));
338: PetscCall(TSMonitorSet(ts, TSMonitorSolutionVTK, filetemplate, (PetscErrorCode(*)(void **))TSMonitorSolutionVTKDestroy));
339: }
341: PetscCall(PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg));
342: if (flg) {
343: TSMonitorDMDARayCtx *rayctx;
344: int ray = 0;
345: DMDirection ddir;
346: DM da;
347: PetscMPIInt rank;
349: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
350: if (dir[0] == 'x') ddir = DM_X;
351: else if (dir[0] == 'y') ddir = DM_Y;
352: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
353: sscanf(dir + 2, "%d", &ray);
355: PetscCall(PetscInfo(((PetscObject)ts), "Displaying DMDA ray %c = %d\n", dir[0], ray));
356: PetscCall(PetscNew(&rayctx));
357: PetscCall(TSGetDM(ts, &da));
358: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
359: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank));
360: if (rank == 0) PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer));
361: rayctx->lgctx = NULL;
362: PetscCall(TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy));
363: }
364: PetscCall(PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg));
365: if (flg) {
366: TSMonitorDMDARayCtx *rayctx;
367: int ray = 0;
368: DMDirection ddir;
369: DM da;
370: PetscInt howoften = 1;
372: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
373: if (dir[0] == 'x') ddir = DM_X;
374: else if (dir[0] == 'y') ddir = DM_Y;
375: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
376: sscanf(dir + 2, "%d", &ray);
378: PetscCall(PetscInfo(((PetscObject)ts), "Displaying LG DMDA ray %c = %d\n", dir[0], ray));
379: PetscCall(PetscNew(&rayctx));
380: PetscCall(TSGetDM(ts, &da));
381: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
382: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx));
383: PetscCall(TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy));
384: }
386: PetscCall(PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt));
387: if (opt) {
388: TSMonitorEnvelopeCtx ctx;
390: PetscCall(TSMonitorEnvelopeCtxCreate(ts, &ctx));
391: PetscCall(TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscErrorCode(*)(void **))TSMonitorEnvelopeCtxDestroy));
392: }
393: flg = PETSC_FALSE;
394: PetscCall(PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt));
395: if (opt && flg) PetscCall(TSMonitorCancel(ts));
397: flg = PETSC_FALSE;
398: PetscCall(PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeIJacobianDefaultColor", flg, &flg, NULL));
399: if (flg) {
400: DM dm;
402: PetscCall(TSGetDM(ts, &dm));
403: PetscCall(DMTSUnsetIJacobianContext_Internal(dm));
404: PetscCall(TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL));
405: PetscCall(PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n"));
406: }
408: /* Handle specific TS options */
409: PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);
411: /* Handle TSAdapt options */
412: PetscCall(TSGetAdapt(ts, &ts->adapt));
413: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
414: PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject));
416: /* TS trajectory must be set after TS, since it may use some TS options above */
417: tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
418: PetscCall(PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL));
419: if (tflg) PetscCall(TSSetSaveTrajectory(ts));
421: PetscCall(TSAdjointSetFromOptions(ts, PetscOptionsObject));
423: /* process any options handlers added with PetscObjectAddOptionsHandler() */
424: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject));
425: PetscOptionsEnd();
427: if (ts->trajectory) PetscCall(TSTrajectorySetFromOptions(ts->trajectory, ts));
429: /* why do we have to do this here and not during TSSetUp? */
430: PetscCall(TSGetSNES(ts, &ts->snes));
431: if (ts->problem_type == TS_LINEAR) {
432: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
433: if (!flg) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
434: }
435: PetscCall(SNESSetFromOptions(ts->snes));
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: /*@
440: TSGetTrajectory - Gets the trajectory from a `TS` if it exists
442: Collective
444: Input Parameter:
445: . ts - the `TS` context obtained from `TSCreate()`
447: Output Parameter:
448: . tr - the `TSTrajectory` object, if it exists
450: Level: advanced
452: Note:
453: This routine should be called after all `TS` options have been set
455: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectoryCreate()`
456: @*/
457: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
458: {
459: PetscFunctionBegin;
461: *tr = ts->trajectory;
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: /*@
466: TSSetSaveTrajectory - Causes the `TS` to save its solutions as it iterates forward in time in a `TSTrajectory` object
468: Collective
470: Input Parameter:
471: . ts - the `TS` context obtained from `TSCreate()`
473: Options Database Keys:
474: + -ts_save_trajectory - saves the trajectory to a file
475: - -ts_trajectory_type type - set trajectory type
477: Level: intermediate
479: Notes:
480: This routine should be called after all `TS` options have been set
482: The `TSTRAJECTORYVISUALIZATION` files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
483: MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m
485: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`
486: @*/
487: PetscErrorCode TSSetSaveTrajectory(TS ts)
488: {
489: PetscFunctionBegin;
491: if (!ts->trajectory) PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@
496: TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object
498: Collective
500: Input Parameter:
501: . ts - the `TS` context obtained from `TSCreate()`
503: Level: intermediate
505: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
506: @*/
507: PetscErrorCode TSResetTrajectory(TS ts)
508: {
509: PetscFunctionBegin;
511: if (ts->trajectory) {
512: PetscCall(TSTrajectoryDestroy(&ts->trajectory));
513: PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
514: }
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: /*@
519: TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from a `TS`
521: Collective
523: Input Parameter:
524: . ts - the `TS` context obtained from `TSCreate()`
526: Level: intermediate
528: .seealso: [](ch_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
529: @*/
530: PetscErrorCode TSRemoveTrajectory(TS ts)
531: {
532: PetscFunctionBegin;
534: if (ts->trajectory) PetscCall(TSTrajectoryDestroy(&ts->trajectory));
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: /*@
539: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
540: set with `TSSetRHSJacobian()`.
542: Collective
544: Input Parameters:
545: + ts - the `TS` context
546: . t - current timestep
547: - U - input vector
549: Output Parameters:
550: + A - Jacobian matrix
551: - B - optional preconditioning matrix
553: Level: developer
555: Note:
556: Most users should not need to explicitly call this routine, as it
557: is used internally within the nonlinear solvers.
559: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()`
560: @*/
561: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
562: {
563: PetscObjectState Ustate;
564: PetscObjectId Uid;
565: DM dm;
566: DMTS tsdm;
567: TSRHSJacobianFn *rhsjacobianfunc;
568: void *ctx;
569: TSRHSFunctionFn *rhsfunction;
571: PetscFunctionBegin;
574: PetscCheckSameComm(ts, 1, U, 3);
575: PetscCall(TSGetDM(ts, &dm));
576: PetscCall(DMGetDMTS(dm, &tsdm));
577: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
578: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx));
579: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
580: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
582: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) PetscFunctionReturn(PETSC_SUCCESS);
584: PetscCheck(ts->rhsjacobian.shift == 0.0 || !ts->rhsjacobian.reuse, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Should not call TSComputeRHSJacobian() on a shifted matrix (shift=%lf) when RHSJacobian is reusable.", (double)ts->rhsjacobian.shift);
585: if (rhsjacobianfunc) {
586: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
587: PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
588: ts->rhsjacs++;
589: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
590: } else {
591: PetscCall(MatZeroEntries(A));
592: if (B && A != B) PetscCall(MatZeroEntries(B));
593: }
594: ts->rhsjacobian.time = t;
595: ts->rhsjacobian.shift = 0;
596: ts->rhsjacobian.scale = 1.;
597: PetscCall(PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid));
598: PetscCall(PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate));
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: /*@
603: TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS`
605: Collective
607: Input Parameters:
608: + ts - the `TS` context
609: . t - current time
610: - U - state vector
612: Output Parameter:
613: . y - right-hand side
615: Level: developer
617: Note:
618: Most users should not need to explicitly call this routine, as it
619: is used internally within the nonlinear solvers.
621: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
622: @*/
623: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
624: {
625: TSRHSFunctionFn *rhsfunction;
626: TSIFunctionFn *ifunction;
627: void *ctx;
628: DM dm;
630: PetscFunctionBegin;
634: PetscCall(TSGetDM(ts, &dm));
635: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx));
636: PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
638: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
640: if (rhsfunction) {
641: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, y, 0));
642: PetscCall(VecLockReadPush(U));
643: PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
644: PetscCall(VecLockReadPop(U));
645: ts->rhsfuncs++;
646: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, y, 0));
647: } else PetscCall(VecZeroEntries(y));
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
651: /*@
652: TSComputeSolutionFunction - Evaluates the solution function.
654: Collective
656: Input Parameters:
657: + ts - the `TS` context
658: - t - current time
660: Output Parameter:
661: . U - the solution
663: Level: developer
665: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
666: @*/
667: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
668: {
669: TSSolutionFn *solutionfunction;
670: void *ctx;
671: DM dm;
673: PetscFunctionBegin;
676: PetscCall(TSGetDM(ts, &dm));
677: PetscCall(DMTSGetSolutionFunction(dm, &solutionfunction, &ctx));
678: if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
681: /*@
682: TSComputeForcingFunction - Evaluates the forcing function.
684: Collective
686: Input Parameters:
687: + ts - the `TS` context
688: - t - current time
690: Output Parameter:
691: . U - the function value
693: Level: developer
695: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
696: @*/
697: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
698: {
699: void *ctx;
700: DM dm;
701: TSForcingFn *forcing;
703: PetscFunctionBegin;
706: PetscCall(TSGetDM(ts, &dm));
707: PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));
709: if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: static PetscErrorCode TSGetRHSVec_Private(TS ts, Vec *Frhs)
714: {
715: Vec F;
717: PetscFunctionBegin;
718: *Frhs = NULL;
719: PetscCall(TSGetIFunction(ts, &F, NULL, NULL));
720: if (!ts->Frhs) PetscCall(VecDuplicate(F, &ts->Frhs));
721: *Frhs = ts->Frhs;
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
726: {
727: Mat A, B;
728: TSIJacobianFn *ijacobian;
730: PetscFunctionBegin;
731: if (Arhs) *Arhs = NULL;
732: if (Brhs) *Brhs = NULL;
733: PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL));
734: if (Arhs) {
735: if (!ts->Arhs) {
736: if (ijacobian) {
737: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs));
738: PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN));
739: } else {
740: ts->Arhs = A;
741: PetscCall(PetscObjectReference((PetscObject)A));
742: }
743: } else {
744: PetscBool flg;
745: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
746: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
747: if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
748: PetscCall(PetscObjectDereference((PetscObject)ts->Arhs));
749: ts->Arhs = A;
750: PetscCall(PetscObjectReference((PetscObject)A));
751: }
752: }
753: *Arhs = ts->Arhs;
754: }
755: if (Brhs) {
756: if (!ts->Brhs) {
757: if (A != B) {
758: if (ijacobian) {
759: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs));
760: } else {
761: ts->Brhs = B;
762: PetscCall(PetscObjectReference((PetscObject)B));
763: }
764: } else {
765: PetscCall(PetscObjectReference((PetscObject)ts->Arhs));
766: ts->Brhs = ts->Arhs;
767: }
768: }
769: *Brhs = ts->Brhs;
770: }
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: /*@
775: TSComputeIFunction - Evaluates the DAE residual written in the implicit form F(t,U,Udot)=0
777: Collective
779: Input Parameters:
780: + ts - the `TS` context
781: . t - current time
782: . U - state vector
783: . Udot - time derivative of state vector
784: - imex - flag indicates if the method is `TSIMEX` so that the RHSFunction should be kept separate
786: Output Parameter:
787: . Y - right-hand side
789: Level: developer
791: Note:
792: Most users should not need to explicitly call this routine, as it
793: is used internally within the nonlinear solvers.
795: If the user did not write their equations in implicit form, this
796: function recasts them in implicit form.
798: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
799: @*/
800: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
801: {
802: TSIFunctionFn *ifunction;
803: TSRHSFunctionFn *rhsfunction;
804: void *ctx;
805: DM dm;
807: PetscFunctionBegin;
813: PetscCall(TSGetDM(ts, &dm));
814: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
815: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
817: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
819: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, Udot, Y));
820: if (ifunction) {
821: PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
822: ts->ifuncs++;
823: }
824: if (imex) {
825: if (!ifunction) PetscCall(VecCopy(Udot, Y));
826: } else if (rhsfunction) {
827: if (ifunction) {
828: Vec Frhs;
829: PetscCall(TSGetRHSVec_Private(ts, &Frhs));
830: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
831: PetscCall(VecAXPY(Y, -1, Frhs));
832: } else {
833: PetscCall(TSComputeRHSFunction(ts, t, U, Y));
834: PetscCall(VecAYPX(Y, -1, Udot));
835: }
836: }
837: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, Udot, Y));
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: /*
842: TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call `TSComputeRHSJacobian()` on it.
844: Note:
845: This routine is needed when one switches from `TSComputeIJacobian()` to `TSComputeRHSJacobian()` because the Jacobian matrix may be shifted or scaled in `TSComputeIJacobian()`.
847: */
848: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
849: {
850: PetscFunctionBegin;
852: PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
853: PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");
855: if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
856: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
857: if (B && B == ts->Brhs && A != B) {
858: if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
859: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
860: }
861: ts->rhsjacobian.shift = 0;
862: ts->rhsjacobian.scale = 1.;
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: /*@
867: TSComputeIJacobian - Evaluates the Jacobian of the DAE
869: Collective
871: Input Parameters:
872: + ts - the `TS` context
873: . t - current timestep
874: . U - state vector
875: . Udot - time derivative of state vector
876: . shift - shift to apply, see note below
877: - imex - flag indicates if the method is `TSIMEX` so that the RHSJacobian should be kept separate
879: Output Parameters:
880: + A - Jacobian matrix
881: - B - matrix from which the preconditioner is constructed; often the same as `A`
883: Level: developer
885: Notes:
886: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
887: .vb
888: dF/dU + shift*dF/dUdot
889: .ve
890: Most users should not need to explicitly call this routine, as it
891: is used internally within the nonlinear solvers.
893: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
894: @*/
895: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
896: {
897: TSIJacobianFn *ijacobian;
898: TSRHSJacobianFn *rhsjacobian;
899: DM dm;
900: void *ctx;
902: PetscFunctionBegin;
909: PetscCall(TSGetDM(ts, &dm));
910: PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
911: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
913: PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
915: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
916: if (ijacobian) {
917: PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
918: ts->ijacs++;
919: }
920: if (imex) {
921: if (!ijacobian) { /* system was written as Udot = G(t,U) */
922: PetscBool assembled;
923: if (rhsjacobian) {
924: Mat Arhs = NULL;
925: PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
926: if (A == Arhs) {
927: 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 */
928: ts->rhsjacobian.time = PETSC_MIN_REAL;
929: }
930: }
931: PetscCall(MatZeroEntries(A));
932: PetscCall(MatAssembled(A, &assembled));
933: if (!assembled) {
934: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
935: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
936: }
937: PetscCall(MatShift(A, shift));
938: if (A != B) {
939: PetscCall(MatZeroEntries(B));
940: PetscCall(MatAssembled(B, &assembled));
941: if (!assembled) {
942: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
943: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
944: }
945: PetscCall(MatShift(B, shift));
946: }
947: }
948: } else {
949: Mat Arhs = NULL, Brhs = NULL;
951: /* RHSJacobian needs to be converted to part of IJacobian if exists */
952: if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
953: if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
954: PetscObjectState Ustate;
955: PetscObjectId Uid;
956: TSRHSFunctionFn *rhsfunction;
958: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
959: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
960: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
961: if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
962: ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */
963: PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
964: if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
965: } else {
966: PetscBool flg;
968: if (ts->rhsjacobian.reuse) { /* Undo the damage */
969: /* MatScale has a short path for this case.
970: However, this code path is taken the first time TSComputeRHSJacobian is called
971: and the matrices have not been assembled yet */
972: PetscCall(TSRecoverRHSJacobian(ts, A, B));
973: }
974: PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
975: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
976: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
977: if (!flg) {
978: PetscCall(MatScale(A, -1));
979: PetscCall(MatShift(A, shift));
980: }
981: if (A != B) {
982: PetscCall(MatScale(B, -1));
983: PetscCall(MatShift(B, shift));
984: }
985: }
986: ts->rhsjacobian.scale = -1;
987: ts->rhsjacobian.shift = shift;
988: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
989: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
990: PetscCall(MatZeroEntries(A));
991: PetscCall(MatShift(A, shift));
992: if (A != B) {
993: PetscCall(MatZeroEntries(B));
994: PetscCall(MatShift(B, shift));
995: }
996: }
997: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
998: PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
999: if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
1000: }
1001: }
1002: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
1003: PetscFunctionReturn(PETSC_SUCCESS);
1004: }
1006: /*@C
1007: TSSetRHSFunction - Sets the routine for evaluating the function,
1008: where U_t = G(t,u).
1010: Logically Collective
1012: Input Parameters:
1013: + ts - the `TS` context obtained from `TSCreate()`
1014: . r - vector to put the computed right-hand side (or `NULL` to have it created)
1015: . f - routine for evaluating the right-hand-side function
1016: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
1018: Level: beginner
1020: Note:
1021: You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.
1023: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1024: @*/
1025: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, void *ctx)
1026: {
1027: SNES snes;
1028: Vec ralloc = NULL;
1029: DM dm;
1031: PetscFunctionBegin;
1035: PetscCall(TSGetDM(ts, &dm));
1036: PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1037: PetscCall(TSGetSNES(ts, &snes));
1038: if (!r && !ts->dm && ts->vec_sol) {
1039: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1040: r = ralloc;
1041: }
1042: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1043: PetscCall(VecDestroy(&ralloc));
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: /*@C
1048: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1050: Logically Collective
1052: Input Parameters:
1053: + ts - the `TS` context obtained from `TSCreate()`
1054: . f - routine for evaluating the solution
1055: - ctx - [optional] user-defined context for private data for the
1056: function evaluation routine (may be `NULL`)
1058: Options Database Keys:
1059: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1060: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
1062: Level: intermediate
1064: Notes:
1065: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1066: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1067: create closed-form solutions with non-physical forcing terms.
1069: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1071: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1072: @*/
1073: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, void *ctx)
1074: {
1075: DM dm;
1077: PetscFunctionBegin;
1079: PetscCall(TSGetDM(ts, &dm));
1080: PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1081: PetscFunctionReturn(PETSC_SUCCESS);
1082: }
1084: /*@C
1085: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1087: Logically Collective
1089: Input Parameters:
1090: + ts - the `TS` context obtained from `TSCreate()`
1091: . func - routine for evaluating the forcing function
1092: - ctx - [optional] user-defined context for private data for the function evaluation routine
1093: (may be `NULL`)
1095: Level: intermediate
1097: Notes:
1098: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1099: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1100: definition of the problem you are solving and hence possibly introducing bugs.
1102: This replaces the ODE F(u,u_t,t) = 0 the `TS` is solving with F(u,u_t,t) - func(t) = 0
1104: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1105: parameters can be passed in the ctx variable.
1107: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1109: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
1110: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1111: @*/
1112: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, void *ctx)
1113: {
1114: DM dm;
1116: PetscFunctionBegin;
1118: PetscCall(TSGetDM(ts, &dm));
1119: PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1120: PetscFunctionReturn(PETSC_SUCCESS);
1121: }
1123: /*@C
1124: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1125: where U_t = G(U,t), as well as the location to store the matrix.
1127: Logically Collective
1129: Input Parameters:
1130: + ts - the `TS` context obtained from `TSCreate()`
1131: . Amat - (approximate) location to store Jacobian matrix entries computed by `f`
1132: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1133: . f - the Jacobian evaluation routine
1134: - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1136: Level: beginner
1138: Notes:
1139: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1141: The `TS` solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f()`
1142: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1144: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
1145: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
1146: @*/
1147: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, void *ctx)
1148: {
1149: SNES snes;
1150: DM dm;
1151: TSIJacobianFn *ijacobian;
1153: PetscFunctionBegin;
1157: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1158: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1160: PetscCall(TSGetDM(ts, &dm));
1161: PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1162: PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1163: PetscCall(TSGetSNES(ts, &snes));
1164: if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1165: if (Amat) {
1166: PetscCall(PetscObjectReference((PetscObject)Amat));
1167: PetscCall(MatDestroy(&ts->Arhs));
1168: ts->Arhs = Amat;
1169: }
1170: if (Pmat) {
1171: PetscCall(PetscObjectReference((PetscObject)Pmat));
1172: PetscCall(MatDestroy(&ts->Brhs));
1173: ts->Brhs = Pmat;
1174: }
1175: PetscFunctionReturn(PETSC_SUCCESS);
1176: }
1178: /*@C
1179: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1181: Logically Collective
1183: Input Parameters:
1184: + ts - the `TS` context obtained from `TSCreate()`
1185: . r - vector to hold the residual (or `NULL` to have it created internally)
1186: . f - the function evaluation routine
1187: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1189: Level: beginner
1191: Note:
1192: The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function.
1194: .seealso: [](ch_ts), `TS`, `TSIFunctionFn`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`,
1195: `TSSetIJacobian()`
1196: @*/
1197: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunctionFn *f, void *ctx)
1198: {
1199: SNES snes;
1200: Vec ralloc = NULL;
1201: DM dm;
1203: PetscFunctionBegin;
1207: PetscCall(TSGetDM(ts, &dm));
1208: PetscCall(DMTSSetIFunction(dm, f, ctx));
1210: PetscCall(TSGetSNES(ts, &snes));
1211: if (!r && !ts->dm && ts->vec_sol) {
1212: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1213: r = ralloc;
1214: }
1215: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1216: PetscCall(VecDestroy(&ralloc));
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: /*@C
1221: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.
1223: Not Collective
1225: Input Parameter:
1226: . ts - the `TS` context
1228: Output Parameters:
1229: + r - vector to hold residual (or `NULL`)
1230: . func - the function to compute residual (or `NULL`)
1231: - ctx - the function context (or `NULL`)
1233: Level: advanced
1235: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1236: @*/
1237: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunctionFn **func, void **ctx)
1238: {
1239: SNES snes;
1240: DM dm;
1242: PetscFunctionBegin;
1244: PetscCall(TSGetSNES(ts, &snes));
1245: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1246: PetscCall(TSGetDM(ts, &dm));
1247: PetscCall(DMTSGetIFunction(dm, func, ctx));
1248: PetscFunctionReturn(PETSC_SUCCESS);
1249: }
1251: /*@C
1252: TSGetRHSFunction - Returns the vector where the right-hand side is stored and the function/context to compute it.
1254: Not Collective
1256: Input Parameter:
1257: . ts - the `TS` context
1259: Output Parameters:
1260: + r - vector to hold computed right-hand side (or `NULL`)
1261: . func - the function to compute right-hand side (or `NULL`)
1262: - ctx - the function context (or `NULL`)
1264: Level: advanced
1266: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1267: @*/
1268: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunctionFn **func, void **ctx)
1269: {
1270: SNES snes;
1271: DM dm;
1273: PetscFunctionBegin;
1275: PetscCall(TSGetSNES(ts, &snes));
1276: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1277: PetscCall(TSGetDM(ts, &dm));
1278: PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1279: PetscFunctionReturn(PETSC_SUCCESS);
1280: }
1282: /*@C
1283: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1284: provided with `TSSetIFunction()`.
1286: Logically Collective
1288: Input Parameters:
1289: + ts - the `TS` context obtained from `TSCreate()`
1290: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
1291: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1292: . f - the Jacobian evaluation routine
1293: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1295: Level: beginner
1297: Notes:
1298: The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1300: If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
1301: space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.
1303: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1304: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1305: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1306: a and vector W depend on the integration method, step size, and past states. For example with
1307: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1308: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1310: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1312: The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f`
1313: You should not assume the values are the same in the next call to `f` as you set them in the previous call.
1315: .seealso: [](ch_ts), `TS`, `TSIJacobianFn`, `TSSetIFunction()`, `TSSetRHSJacobian()`,
1316: `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1317: @*/
1318: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobianFn *f, void *ctx)
1319: {
1320: SNES snes;
1321: DM dm;
1323: PetscFunctionBegin;
1327: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1328: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1330: PetscCall(TSGetDM(ts, &dm));
1331: PetscCall(DMTSSetIJacobian(dm, f, ctx));
1333: PetscCall(TSGetSNES(ts, &snes));
1334: PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1335: PetscFunctionReturn(PETSC_SUCCESS);
1336: }
1338: /*@
1339: TSRHSJacobianSetReuse - restore the RHS Jacobian before calling the user-provided `TSRHSJacobianFn` function again
1341: Logically Collective
1343: Input Parameters:
1344: + ts - `TS` context obtained from `TSCreate()`
1345: - reuse - `PETSC_TRUE` if the RHS Jacobian
1347: Level: intermediate
1349: Notes:
1350: Without this flag, `TS` will change the sign and shift the RHS Jacobian for a
1351: finite-time-step implicit solve, in which case the user function will need to recompute the
1352: entire Jacobian. The `reuse `flag must be set if the evaluation function assumes that the
1353: matrix entries have not been changed by the `TS`.
1355: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1356: @*/
1357: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1358: {
1359: PetscFunctionBegin;
1360: ts->rhsjacobian.reuse = reuse;
1361: PetscFunctionReturn(PETSC_SUCCESS);
1362: }
1364: /*@C
1365: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1367: Logically Collective
1369: Input Parameters:
1370: + ts - the `TS` context obtained from `TSCreate()`
1371: . F - vector to hold the residual (or `NULL` to have it created internally)
1372: . fun - the function evaluation routine
1373: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1375: Level: beginner
1377: .seealso: [](ch_ts), `TS`, `TSI2FunctionFn`, `TSSetI2Jacobian()`, `TSSetIFunction()`,
1378: `TSCreate()`, `TSSetRHSFunction()`
1379: @*/
1380: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2FunctionFn *fun, void *ctx)
1381: {
1382: DM dm;
1384: PetscFunctionBegin;
1387: PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1388: PetscCall(TSGetDM(ts, &dm));
1389: PetscCall(DMTSSetI2Function(dm, fun, ctx));
1390: PetscFunctionReturn(PETSC_SUCCESS);
1391: }
1393: /*@C
1394: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.
1396: Not Collective
1398: Input Parameter:
1399: . ts - the `TS` context
1401: Output Parameters:
1402: + r - vector to hold residual (or `NULL`)
1403: . fun - the function to compute residual (or `NULL`)
1404: - ctx - the function context (or `NULL`)
1406: Level: advanced
1408: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1409: @*/
1410: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2FunctionFn **fun, void **ctx)
1411: {
1412: SNES snes;
1413: DM dm;
1415: PetscFunctionBegin;
1417: PetscCall(TSGetSNES(ts, &snes));
1418: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1419: PetscCall(TSGetDM(ts, &dm));
1420: PetscCall(DMTSGetI2Function(dm, fun, ctx));
1421: PetscFunctionReturn(PETSC_SUCCESS);
1422: }
1424: /*@C
1425: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1426: where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.
1428: Logically Collective
1430: Input Parameters:
1431: + ts - the `TS` context obtained from `TSCreate()`
1432: . J - matrix to hold the Jacobian values
1433: . P - matrix for constructing the preconditioner (may be same as `J`)
1434: . jac - the Jacobian evaluation routine, see `TSI2JacobianFn` for the calling sequence
1435: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1437: Level: beginner
1439: Notes:
1440: The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1442: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1443: 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.
1444: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1445: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1447: .seealso: [](ch_ts), `TS`, `TSI2JacobianFn`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1448: @*/
1449: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2JacobianFn *jac, void *ctx)
1450: {
1451: DM dm;
1453: PetscFunctionBegin;
1457: PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1458: PetscCall(TSGetDM(ts, &dm));
1459: PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1460: PetscFunctionReturn(PETSC_SUCCESS);
1461: }
1463: /*@C
1464: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1466: Not Collective, but parallel objects are returned if `TS` is parallel
1468: Input Parameter:
1469: . ts - The `TS` context obtained from `TSCreate()`
1471: Output Parameters:
1472: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1473: . P - The matrix from which the preconditioner is constructed, often the same as `J`
1474: . jac - The function to compute the Jacobian matrices
1475: - ctx - User-defined context for Jacobian evaluation routine
1477: Level: advanced
1479: Note:
1480: You can pass in `NULL` for any return argument you do not need.
1482: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1483: @*/
1484: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2JacobianFn **jac, void **ctx)
1485: {
1486: SNES snes;
1487: DM dm;
1489: PetscFunctionBegin;
1490: PetscCall(TSGetSNES(ts, &snes));
1491: PetscCall(SNESSetUpMatrices(snes));
1492: PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1493: PetscCall(TSGetDM(ts, &dm));
1494: PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1495: PetscFunctionReturn(PETSC_SUCCESS);
1496: }
1498: /*@
1499: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1501: Collective
1503: Input Parameters:
1504: + ts - the `TS` context
1505: . t - current time
1506: . U - state vector
1507: . V - time derivative of state vector (U_t)
1508: - A - second time derivative of state vector (U_tt)
1510: Output Parameter:
1511: . F - the residual vector
1513: Level: developer
1515: Note:
1516: Most users should not need to explicitly call this routine, as it
1517: is used internally within the nonlinear solvers.
1519: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1520: @*/
1521: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1522: {
1523: DM dm;
1524: TSI2FunctionFn *I2Function;
1525: void *ctx;
1526: TSRHSFunctionFn *rhsfunction;
1528: PetscFunctionBegin;
1535: PetscCall(TSGetDM(ts, &dm));
1536: PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1537: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
1539: if (!I2Function) {
1540: PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1541: PetscFunctionReturn(PETSC_SUCCESS);
1542: }
1544: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, V, F));
1546: PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));
1548: if (rhsfunction) {
1549: Vec Frhs;
1550: PetscCall(TSGetRHSVec_Private(ts, &Frhs));
1551: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1552: PetscCall(VecAXPY(F, -1, Frhs));
1553: }
1555: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, V, F));
1556: PetscFunctionReturn(PETSC_SUCCESS);
1557: }
1559: /*@
1560: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1562: Collective
1564: Input Parameters:
1565: + ts - the `TS` context
1566: . t - current timestep
1567: . U - state vector
1568: . V - time derivative of state vector
1569: . A - second time derivative of state vector
1570: . shiftV - shift to apply, see note below
1571: - shiftA - shift to apply, see note below
1573: Output Parameters:
1574: + J - Jacobian matrix
1575: - P - optional preconditioning matrix
1577: Level: developer
1579: Notes:
1580: If F(t,U,V,A)=0 is the DAE, the required Jacobian is
1582: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1584: Most users should not need to explicitly call this routine, as it
1585: is used internally within the nonlinear solvers.
1587: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1588: @*/
1589: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1590: {
1591: DM dm;
1592: TSI2JacobianFn *I2Jacobian;
1593: void *ctx;
1594: TSRHSJacobianFn *rhsjacobian;
1596: PetscFunctionBegin;
1604: PetscCall(TSGetDM(ts, &dm));
1605: PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1606: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1608: if (!I2Jacobian) {
1609: PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1610: PetscFunctionReturn(PETSC_SUCCESS);
1611: }
1613: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, J, P));
1614: PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1615: if (rhsjacobian) {
1616: Mat Jrhs, Prhs;
1617: PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1618: PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1619: PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1620: if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1621: }
1623: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, J, P));
1624: PetscFunctionReturn(PETSC_SUCCESS);
1625: }
1627: /*@C
1628: TSSetTransientVariable - sets function to transform from state to transient variables
1630: Logically Collective
1632: Input Parameters:
1633: + ts - time stepping context on which to change the transient variable
1634: . tvar - a function that transforms to transient variables, see `TSTransientVariableFn` for the calling sequence
1635: - ctx - a context for tvar
1637: Level: advanced
1639: Notes:
1640: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1641: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1642: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1643: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1644: evaluated via the chain rule, as in
1645: .vb
1646: dF/dP + shift * dF/dCdot dC/dP.
1647: .ve
1649: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSTransientVariableFn`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1650: @*/
1651: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariableFn *tvar, void *ctx)
1652: {
1653: DM dm;
1655: PetscFunctionBegin;
1657: PetscCall(TSGetDM(ts, &dm));
1658: PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1659: PetscFunctionReturn(PETSC_SUCCESS);
1660: }
1662: /*@
1663: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1665: Logically Collective
1667: Input Parameters:
1668: + ts - TS on which to compute
1669: - U - state vector to be transformed to transient variables
1671: Output Parameter:
1672: . C - transient (conservative) variable
1674: Level: developer
1676: Developer Notes:
1677: If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed.
1678: This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are
1679: being used.
1681: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1682: @*/
1683: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1684: {
1685: DM dm;
1686: DMTS dmts;
1688: PetscFunctionBegin;
1691: PetscCall(TSGetDM(ts, &dm));
1692: PetscCall(DMGetDMTS(dm, &dmts));
1693: if (dmts->ops->transientvar) {
1695: PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1696: }
1697: PetscFunctionReturn(PETSC_SUCCESS);
1698: }
1700: /*@
1701: TSHasTransientVariable - determine whether transient variables have been set
1703: Logically Collective
1705: Input Parameter:
1706: . ts - `TS` on which to compute
1708: Output Parameter:
1709: . has - `PETSC_TRUE` if transient variables have been set
1711: Level: developer
1713: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1714: @*/
1715: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1716: {
1717: DM dm;
1718: DMTS dmts;
1720: PetscFunctionBegin;
1722: PetscCall(TSGetDM(ts, &dm));
1723: PetscCall(DMGetDMTS(dm, &dmts));
1724: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1725: PetscFunctionReturn(PETSC_SUCCESS);
1726: }
1728: /*@
1729: TS2SetSolution - Sets the initial solution and time derivative vectors
1730: for use by the `TS` routines handling second order equations.
1732: Logically Collective
1734: Input Parameters:
1735: + ts - the `TS` context obtained from `TSCreate()`
1736: . u - the solution vector
1737: - v - the time derivative vector
1739: Level: beginner
1741: .seealso: [](ch_ts), `TS`
1742: @*/
1743: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1744: {
1745: PetscFunctionBegin;
1749: PetscCall(TSSetSolution(ts, u));
1750: PetscCall(PetscObjectReference((PetscObject)v));
1751: PetscCall(VecDestroy(&ts->vec_dot));
1752: ts->vec_dot = v;
1753: PetscFunctionReturn(PETSC_SUCCESS);
1754: }
1756: /*@
1757: TS2GetSolution - Returns the solution and time derivative at the present timestep
1758: for second order equations.
1760: Not Collective
1762: Input Parameter:
1763: . ts - the `TS` context obtained from `TSCreate()`
1765: Output Parameters:
1766: + u - the vector containing the solution
1767: - v - the vector containing the time derivative
1769: Level: intermediate
1771: Notes:
1772: It is valid to call this routine inside the function
1773: that you are evaluating in order to move to the new timestep. This vector not
1774: changed until the solution at the next timestep has been calculated.
1776: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1777: @*/
1778: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1779: {
1780: PetscFunctionBegin;
1782: if (u) PetscAssertPointer(u, 2);
1783: if (v) PetscAssertPointer(v, 3);
1784: if (u) *u = ts->vec_sol;
1785: if (v) *v = ts->vec_dot;
1786: PetscFunctionReturn(PETSC_SUCCESS);
1787: }
1789: /*@C
1790: TSLoad - Loads a `TS` that has been stored in binary with `TSView()`.
1792: Collective
1794: Input Parameters:
1795: + ts - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1796: some related function before a call to `TSLoad()`.
1797: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1799: Level: intermediate
1801: Note:
1802: The type is determined by the data in the file, any type set into the `TS` before this call is ignored.
1804: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1805: @*/
1806: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1807: {
1808: PetscBool isbinary;
1809: PetscInt classid;
1810: char type[256];
1811: DMTS sdm;
1812: DM dm;
1814: PetscFunctionBegin;
1817: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1818: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1820: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1821: PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1822: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1823: PetscCall(TSSetType(ts, type));
1824: PetscTryTypeMethod(ts, load, viewer);
1825: PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1826: PetscCall(DMLoad(dm, viewer));
1827: PetscCall(TSSetDM(ts, dm));
1828: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1829: PetscCall(VecLoad(ts->vec_sol, viewer));
1830: PetscCall(DMGetDMTS(ts->dm, &sdm));
1831: PetscCall(DMTSLoad(sdm, viewer));
1832: PetscFunctionReturn(PETSC_SUCCESS);
1833: }
1835: #include <petscdraw.h>
1836: #if defined(PETSC_HAVE_SAWS)
1837: #include <petscviewersaws.h>
1838: #endif
1840: /*@C
1841: TSViewFromOptions - View a `TS` based on values in the options database
1843: Collective
1845: Input Parameters:
1846: + ts - the `TS` context
1847: . obj - Optional object that provides the prefix for the options database keys
1848: - name - command line option string to be passed by user
1850: Level: intermediate
1852: .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1853: @*/
1854: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1855: {
1856: PetscFunctionBegin;
1858: PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1859: PetscFunctionReturn(PETSC_SUCCESS);
1860: }
1862: /*@C
1863: TSView - Prints the `TS` data structure.
1865: Collective
1867: Input Parameters:
1868: + ts - the `TS` context obtained from `TSCreate()`
1869: - viewer - visualization context
1871: Options Database Key:
1872: . -ts_view - calls `TSView()` at end of `TSStep()`
1874: Level: beginner
1876: Notes:
1877: The available visualization contexts include
1878: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1879: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1880: output where only the first processor opens
1881: the file. All other processors send their
1882: data to the first processor to print.
1884: The user can open an alternative visualization context with
1885: `PetscViewerASCIIOpen()` - output to a specified file.
1887: In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).
1889: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1890: @*/
1891: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1892: {
1893: TSType type;
1894: PetscBool iascii, isstring, isundials, isbinary, isdraw;
1895: DMTS sdm;
1896: #if defined(PETSC_HAVE_SAWS)
1897: PetscBool issaws;
1898: #endif
1900: PetscFunctionBegin;
1902: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1904: PetscCheckSameComm(ts, 1, viewer, 2);
1906: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1907: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1908: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1909: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1910: #if defined(PETSC_HAVE_SAWS)
1911: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1912: #endif
1913: if (iascii) {
1914: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer));
1915: if (ts->ops->view) {
1916: PetscCall(PetscViewerASCIIPushTab(viewer));
1917: PetscUseTypeMethod(ts, view, viewer);
1918: PetscCall(PetscViewerASCIIPopTab(viewer));
1919: }
1920: if (ts->max_steps < PETSC_MAX_INT) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps));
1921: if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time));
1922: if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs));
1923: if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs));
1924: if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs));
1925: if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs));
1926: if (ts->usessnes) {
1927: PetscBool lin;
1928: if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its));
1929: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its));
1930: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
1931: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures));
1932: }
1933: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject));
1934: if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, "));
1935: else PetscCall(PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol));
1936: if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of absolute error tolerances\n"));
1937: else PetscCall(PetscViewerASCIIPrintf(viewer, " using absolute error tolerance of %g\n", (double)ts->atol));
1938: PetscCall(PetscViewerASCIIPushTab(viewer));
1939: PetscCall(TSAdaptView(ts->adapt, viewer));
1940: PetscCall(PetscViewerASCIIPopTab(viewer));
1941: } else if (isstring) {
1942: PetscCall(TSGetType(ts, &type));
1943: PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type));
1944: PetscTryTypeMethod(ts, view, viewer);
1945: } else if (isbinary) {
1946: PetscInt classid = TS_FILE_CLASSID;
1947: MPI_Comm comm;
1948: PetscMPIInt rank;
1949: char type[256];
1951: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
1952: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1953: if (rank == 0) {
1954: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1955: PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
1956: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1957: }
1958: PetscTryTypeMethod(ts, view, viewer);
1959: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1960: PetscCall(DMView(ts->dm, viewer));
1961: PetscCall(VecView(ts->vec_sol, viewer));
1962: PetscCall(DMGetDMTS(ts->dm, &sdm));
1963: PetscCall(DMTSView(sdm, viewer));
1964: } else if (isdraw) {
1965: PetscDraw draw;
1966: char str[36];
1967: PetscReal x, y, bottom, h;
1969: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1970: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1971: PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
1972: PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
1973: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
1974: bottom = y - h;
1975: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1976: PetscTryTypeMethod(ts, view, viewer);
1977: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1978: if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
1979: PetscCall(PetscDrawPopCurrentPoint(draw));
1980: #if defined(PETSC_HAVE_SAWS)
1981: } else if (issaws) {
1982: PetscMPIInt rank;
1983: const char *name;
1985: PetscCall(PetscObjectGetName((PetscObject)ts, &name));
1986: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1987: if (!((PetscObject)ts)->amsmem && rank == 0) {
1988: char dir[1024];
1990: PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
1991: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
1992: PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
1993: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
1994: PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
1995: }
1996: PetscTryTypeMethod(ts, view, viewer);
1997: #endif
1998: }
1999: if (ts->snes && ts->usessnes) {
2000: PetscCall(PetscViewerASCIIPushTab(viewer));
2001: PetscCall(SNESView(ts->snes, viewer));
2002: PetscCall(PetscViewerASCIIPopTab(viewer));
2003: }
2004: PetscCall(DMGetDMTS(ts->dm, &sdm));
2005: PetscCall(DMTSView(sdm, viewer));
2007: PetscCall(PetscViewerASCIIPushTab(viewer));
2008: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials));
2009: PetscCall(PetscViewerASCIIPopTab(viewer));
2010: PetscFunctionReturn(PETSC_SUCCESS);
2011: }
2013: /*@
2014: TSSetApplicationContext - Sets an optional user-defined context for
2015: the timesteppers.
2017: Logically Collective
2019: Input Parameters:
2020: + ts - the `TS` context obtained from `TSCreate()`
2021: - usrP - user context
2023: Level: intermediate
2025: Fortran Notes:
2026: You must write a Fortran interface definition for this
2027: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
2029: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2030: @*/
2031: PetscErrorCode TSSetApplicationContext(TS ts, void *usrP)
2032: {
2033: PetscFunctionBegin;
2035: ts->user = usrP;
2036: PetscFunctionReturn(PETSC_SUCCESS);
2037: }
2039: /*@
2040: TSGetApplicationContext - Gets the user-defined context for the
2041: timestepper that was set with `TSSetApplicationContext()`
2043: Not Collective
2045: Input Parameter:
2046: . ts - the `TS` context obtained from `TSCreate()`
2048: Output Parameter:
2049: . usrP - user context
2051: Level: intermediate
2053: Fortran Notes:
2054: You must write a Fortran interface definition for this
2055: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
2057: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2058: @*/
2059: PetscErrorCode TSGetApplicationContext(TS ts, void *usrP)
2060: {
2061: PetscFunctionBegin;
2063: *(void **)usrP = ts->user;
2064: PetscFunctionReturn(PETSC_SUCCESS);
2065: }
2067: /*@
2068: TSGetStepNumber - Gets the number of time steps completed.
2070: Not Collective
2072: Input Parameter:
2073: . ts - the `TS` context obtained from `TSCreate()`
2075: Output Parameter:
2076: . steps - number of steps completed so far
2078: Level: intermediate
2080: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2081: @*/
2082: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2083: {
2084: PetscFunctionBegin;
2086: PetscAssertPointer(steps, 2);
2087: *steps = ts->steps;
2088: PetscFunctionReturn(PETSC_SUCCESS);
2089: }
2091: /*@
2092: TSSetStepNumber - Sets the number of steps completed.
2094: Logically Collective
2096: Input Parameters:
2097: + ts - the `TS` context
2098: - steps - number of steps completed so far
2100: Level: developer
2102: Note:
2103: For most uses of the `TS` solvers the user need not explicitly call
2104: `TSSetStepNumber()`, as the step counter is appropriately updated in
2105: `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2106: reinitialize timestepping by setting the step counter to zero (and time
2107: to the initial time) to solve a similar problem with different initial
2108: conditions or parameters. Other possible use case is to continue
2109: timestepping from a previously interrupted run in such a way that `TS`
2110: monitors will be called with a initial nonzero step counter.
2112: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2113: @*/
2114: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2115: {
2116: PetscFunctionBegin;
2119: PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2120: ts->steps = steps;
2121: PetscFunctionReturn(PETSC_SUCCESS);
2122: }
2124: /*@
2125: TSSetTimeStep - Allows one to reset the timestep at any time,
2126: useful for simple pseudo-timestepping codes.
2128: Logically Collective
2130: Input Parameters:
2131: + ts - the `TS` context obtained from `TSCreate()`
2132: - time_step - the size of the timestep
2134: Level: intermediate
2136: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2137: @*/
2138: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2139: {
2140: PetscFunctionBegin;
2143: ts->time_step = time_step;
2144: PetscFunctionReturn(PETSC_SUCCESS);
2145: }
2147: /*@
2148: TSSetExactFinalTime - Determines whether to adapt the final time step to
2149: match the exact final time, interpolate solution to the exact final time,
2150: or just return at the final time `TS` computed.
2152: Logically Collective
2154: Input Parameters:
2155: + ts - the time-step context
2156: - eftopt - exact final time option
2157: .vb
2158: TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded
2159: TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2160: TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2161: .ve
2163: Options Database Key:
2164: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step at runtime
2166: Level: beginner
2168: Note:
2169: If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2170: then the final time you selected.
2172: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2173: @*/
2174: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2175: {
2176: PetscFunctionBegin;
2179: ts->exact_final_time = eftopt;
2180: PetscFunctionReturn(PETSC_SUCCESS);
2181: }
2183: /*@
2184: TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`
2186: Not Collective
2188: Input Parameter:
2189: . ts - the `TS` context
2191: Output Parameter:
2192: . eftopt - exact final time option
2194: Level: beginner
2196: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2197: @*/
2198: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2199: {
2200: PetscFunctionBegin;
2202: PetscAssertPointer(eftopt, 2);
2203: *eftopt = ts->exact_final_time;
2204: PetscFunctionReturn(PETSC_SUCCESS);
2205: }
2207: /*@
2208: TSGetTimeStep - Gets the current timestep size.
2210: Not Collective
2212: Input Parameter:
2213: . ts - the `TS` context obtained from `TSCreate()`
2215: Output Parameter:
2216: . dt - the current timestep size
2218: Level: intermediate
2220: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2221: @*/
2222: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2223: {
2224: PetscFunctionBegin;
2226: PetscAssertPointer(dt, 2);
2227: *dt = ts->time_step;
2228: PetscFunctionReturn(PETSC_SUCCESS);
2229: }
2231: /*@
2232: TSGetSolution - Returns the solution at the present timestep. It
2233: is valid to call this routine inside the function that you are evaluating
2234: in order to move to the new timestep. This vector not changed until
2235: the solution at the next timestep has been calculated.
2237: Not Collective, but v returned is parallel if ts is parallel
2239: Input Parameter:
2240: . ts - the `TS` context obtained from `TSCreate()`
2242: Output Parameter:
2243: . v - the vector containing the solution
2245: Level: intermediate
2247: Note:
2248: If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
2249: final time. It returns the solution at the next timestep.
2251: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2252: @*/
2253: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2254: {
2255: PetscFunctionBegin;
2257: PetscAssertPointer(v, 2);
2258: *v = ts->vec_sol;
2259: PetscFunctionReturn(PETSC_SUCCESS);
2260: }
2262: /*@
2263: TSGetSolutionComponents - Returns any solution components at the present
2264: timestep, if available for the time integration method being used.
2265: Solution components are quantities that share the same size and
2266: structure as the solution vector.
2268: Not Collective, but v returned is parallel if ts is parallel
2270: Input Parameters:
2271: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2272: . n - If v is `NULL`, then the number of solution components is
2273: returned through n, else the n-th solution component is
2274: returned in v.
2275: - v - the vector containing the n-th solution component
2276: (may be `NULL` to use this function to find out
2277: the number of solutions components).
2279: Level: advanced
2281: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2282: @*/
2283: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2284: {
2285: PetscFunctionBegin;
2287: if (!ts->ops->getsolutioncomponents) *n = 0;
2288: else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2289: PetscFunctionReturn(PETSC_SUCCESS);
2290: }
2292: /*@
2293: TSGetAuxSolution - Returns an auxiliary solution at the present
2294: timestep, if available for the time integration method being used.
2296: Not Collective, but v returned is parallel if ts is parallel
2298: Input Parameters:
2299: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2300: - v - the vector containing the auxiliary solution
2302: Level: intermediate
2304: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2305: @*/
2306: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2307: {
2308: PetscFunctionBegin;
2310: if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2311: else PetscCall(VecZeroEntries(*v));
2312: PetscFunctionReturn(PETSC_SUCCESS);
2313: }
2315: /*@
2316: TSGetTimeError - Returns the estimated error vector, if the chosen
2317: `TSType` has an error estimation functionality and `TSSetTimeError()` was called
2319: Not Collective, but v returned is parallel if ts is parallel
2321: Input Parameters:
2322: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2323: . n - current estimate (n=0) or previous one (n=-1)
2324: - v - the vector containing the error (same size as the solution).
2326: Level: intermediate
2328: Note:
2329: MUST call after `TSSetUp()`
2331: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2332: @*/
2333: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2334: {
2335: PetscFunctionBegin;
2337: if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2338: else PetscCall(VecZeroEntries(*v));
2339: PetscFunctionReturn(PETSC_SUCCESS);
2340: }
2342: /*@
2343: TSSetTimeError - Sets the estimated error vector, if the chosen
2344: `TSType` has an error estimation functionality. This can be used
2345: to restart such a time integrator with a given error vector.
2347: Not Collective, but v returned is parallel if ts is parallel
2349: Input Parameters:
2350: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2351: - v - the vector containing the error (same size as the solution).
2353: Level: intermediate
2355: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
2356: @*/
2357: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2358: {
2359: PetscFunctionBegin;
2361: PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2362: PetscTryTypeMethod(ts, settimeerror, v);
2363: PetscFunctionReturn(PETSC_SUCCESS);
2364: }
2366: /* ----- Routines to initialize and destroy a timestepper ---- */
2367: /*@
2368: TSSetProblemType - Sets the type of problem to be solved.
2370: Not collective
2372: Input Parameters:
2373: + ts - The `TS`
2374: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2375: .vb
2376: U_t - A U = 0 (linear)
2377: U_t - A(t) U = 0 (linear)
2378: F(t,U,U_t) = 0 (nonlinear)
2379: .ve
2381: Level: beginner
2383: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2384: @*/
2385: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2386: {
2387: PetscFunctionBegin;
2389: ts->problem_type = type;
2390: if (type == TS_LINEAR) {
2391: SNES snes;
2392: PetscCall(TSGetSNES(ts, &snes));
2393: PetscCall(SNESSetType(snes, SNESKSPONLY));
2394: }
2395: PetscFunctionReturn(PETSC_SUCCESS);
2396: }
2398: /*@C
2399: TSGetProblemType - Gets the type of problem to be solved.
2401: Not collective
2403: Input Parameter:
2404: . ts - The `TS`
2406: Output Parameter:
2407: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2408: .vb
2409: M U_t = A U
2410: M(t) U_t = A(t) U
2411: F(t,U,U_t)
2412: .ve
2414: Level: beginner
2416: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2417: @*/
2418: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2419: {
2420: PetscFunctionBegin;
2422: PetscAssertPointer(type, 2);
2423: *type = ts->problem_type;
2424: PetscFunctionReturn(PETSC_SUCCESS);
2425: }
2427: /*
2428: Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2429: */
2430: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2431: {
2432: PetscBool isnone;
2434: PetscFunctionBegin;
2435: PetscCall(TSGetAdapt(ts, &ts->adapt));
2436: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2438: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2439: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2440: else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2441: PetscFunctionReturn(PETSC_SUCCESS);
2442: }
2444: /*@
2445: TSSetUp - Sets up the internal data structures for the later use of a timestepper.
2447: Collective
2449: Input Parameter:
2450: . ts - the `TS` context obtained from `TSCreate()`
2452: Level: advanced
2454: Note:
2455: For basic use of the `TS` solvers the user need not explicitly call
2456: `TSSetUp()`, since these actions will automatically occur during
2457: the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this
2458: phase separately, `TSSetUp()` should be called after `TSCreate()`
2459: and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.
2461: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2462: @*/
2463: PetscErrorCode TSSetUp(TS ts)
2464: {
2465: DM dm;
2466: PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2467: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2468: TSIFunctionFn *ifun;
2469: TSIJacobianFn *ijac;
2470: TSI2JacobianFn *i2jac;
2471: TSRHSJacobianFn *rhsjac;
2473: PetscFunctionBegin;
2475: if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
2477: if (!((PetscObject)ts)->type_name) {
2478: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2479: PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2480: }
2482: if (!ts->vec_sol) {
2483: PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2484: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2485: }
2487: if (ts->tspan) {
2488: if (!ts->tspan->vecs_sol) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2489: }
2490: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2491: PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2492: ts->Jacp = ts->Jacprhs;
2493: }
2495: if (ts->quadraturets) {
2496: PetscCall(TSSetUp(ts->quadraturets));
2497: PetscCall(VecDestroy(&ts->vec_costintegrand));
2498: PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2499: }
2501: PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
2502: if (rhsjac == TSComputeRHSJacobianConstant) {
2503: Mat Amat, Pmat;
2504: SNES snes;
2505: PetscCall(TSGetSNES(ts, &snes));
2506: PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
2507: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2508: * have displaced the RHS matrix */
2509: if (Amat && Amat == ts->Arhs) {
2510: /* 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 */
2511: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
2512: PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
2513: PetscCall(MatDestroy(&Amat));
2514: }
2515: if (Pmat && Pmat == ts->Brhs) {
2516: PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
2517: PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
2518: PetscCall(MatDestroy(&Pmat));
2519: }
2520: }
2522: PetscCall(TSGetAdapt(ts, &ts->adapt));
2523: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2525: PetscTryTypeMethod(ts, setup);
2527: PetscCall(TSSetExactFinalTimeDefault(ts));
2529: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2530: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2531: */
2532: PetscCall(TSGetDM(ts, &dm));
2533: PetscCall(DMSNESGetFunction(dm, &func, NULL));
2534: if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));
2536: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2537: Otherwise, the SNES will use coloring internally to form the Jacobian.
2538: */
2539: PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2540: PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2541: PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2542: PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2543: if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));
2545: /* if time integration scheme has a starting method, call it */
2546: PetscTryTypeMethod(ts, startingmethod);
2548: ts->setupcalled = PETSC_TRUE;
2549: PetscFunctionReturn(PETSC_SUCCESS);
2550: }
2552: /*@
2553: TSReset - Resets a `TS` context and removes any allocated `Vec`s and `Mat`s.
2555: Collective
2557: Input Parameter:
2558: . ts - the `TS` context obtained from `TSCreate()`
2560: Level: beginner
2562: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()`
2563: @*/
2564: PetscErrorCode TSReset(TS ts)
2565: {
2566: TS_RHSSplitLink ilink = ts->tsrhssplit, next;
2568: PetscFunctionBegin;
2571: PetscTryTypeMethod(ts, reset);
2572: if (ts->snes) PetscCall(SNESReset(ts->snes));
2573: if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));
2575: PetscCall(MatDestroy(&ts->Arhs));
2576: PetscCall(MatDestroy(&ts->Brhs));
2577: PetscCall(VecDestroy(&ts->Frhs));
2578: PetscCall(VecDestroy(&ts->vec_sol));
2579: PetscCall(VecDestroy(&ts->vec_sol0));
2580: PetscCall(VecDestroy(&ts->vec_dot));
2581: PetscCall(VecDestroy(&ts->vatol));
2582: PetscCall(VecDestroy(&ts->vrtol));
2583: PetscCall(VecDestroyVecs(ts->nwork, &ts->work));
2585: PetscCall(MatDestroy(&ts->Jacprhs));
2586: PetscCall(MatDestroy(&ts->Jacp));
2587: if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2588: if (ts->quadraturets) {
2589: PetscCall(TSReset(ts->quadraturets));
2590: PetscCall(VecDestroy(&ts->vec_costintegrand));
2591: }
2592: while (ilink) {
2593: next = ilink->next;
2594: PetscCall(TSDestroy(&ilink->ts));
2595: PetscCall(PetscFree(ilink->splitname));
2596: PetscCall(ISDestroy(&ilink->is));
2597: PetscCall(PetscFree(ilink));
2598: ilink = next;
2599: }
2600: ts->tsrhssplit = NULL;
2601: ts->num_rhs_splits = 0;
2602: if (ts->tspan) {
2603: PetscCall(PetscFree(ts->tspan->span_times));
2604: PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2605: PetscCall(PetscFree(ts->tspan));
2606: }
2607: ts->rhsjacobian.time = PETSC_MIN_REAL;
2608: ts->rhsjacobian.scale = 1.0;
2609: ts->ijacobian.shift = 1.0;
2610: ts->setupcalled = PETSC_FALSE;
2611: PetscFunctionReturn(PETSC_SUCCESS);
2612: }
2614: static PetscErrorCode TSResizeReset(TS);
2616: /*@C
2617: TSDestroy - Destroys the timestepper context that was created
2618: with `TSCreate()`.
2620: Collective
2622: Input Parameter:
2623: . ts - the `TS` context obtained from `TSCreate()`
2625: Level: beginner
2627: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2628: @*/
2629: PetscErrorCode TSDestroy(TS *ts)
2630: {
2631: PetscFunctionBegin;
2632: if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2634: if (--((PetscObject)*ts)->refct > 0) {
2635: *ts = NULL;
2636: PetscFunctionReturn(PETSC_SUCCESS);
2637: }
2639: PetscCall(TSReset(*ts));
2640: PetscCall(TSAdjointReset(*ts));
2641: if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2642: PetscCall(TSResizeReset(*ts));
2644: /* if memory was published with SAWs then destroy it */
2645: PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2646: PetscTryTypeMethod(*ts, destroy);
2648: PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));
2650: PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2651: PetscCall(TSEventDestroy(&(*ts)->event));
2653: PetscCall(SNESDestroy(&(*ts)->snes));
2654: PetscCall(DMDestroy(&(*ts)->dm));
2655: PetscCall(TSMonitorCancel(*ts));
2656: PetscCall(TSAdjointMonitorCancel(*ts));
2658: PetscCall(TSDestroy(&(*ts)->quadraturets));
2659: PetscCall(PetscHeaderDestroy(ts));
2660: PetscFunctionReturn(PETSC_SUCCESS);
2661: }
2663: /*@
2664: TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2665: a `TS` (timestepper) context. Valid only for nonlinear problems.
2667: Not Collective, but snes is parallel if ts is parallel
2669: Input Parameter:
2670: . ts - the `TS` context obtained from `TSCreate()`
2672: Output Parameter:
2673: . snes - the nonlinear solver context
2675: Level: beginner
2677: Notes:
2678: The user can then directly manipulate the `SNES` context to set various
2679: options, etc. Likewise, the user can then extract and manipulate the
2680: `KSP`, and `PC` contexts as well.
2682: `TSGetSNES()` does not work for integrators that do not use `SNES`; in
2683: this case `TSGetSNES()` returns `NULL` in `snes`.
2685: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2686: @*/
2687: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2688: {
2689: PetscFunctionBegin;
2691: PetscAssertPointer(snes, 2);
2692: if (!ts->snes) {
2693: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2694: PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2695: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2696: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2697: if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2698: if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2699: }
2700: *snes = ts->snes;
2701: PetscFunctionReturn(PETSC_SUCCESS);
2702: }
2704: /*@
2705: TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the timestepping context
2707: Collective
2709: Input Parameters:
2710: + ts - the `TS` context obtained from `TSCreate()`
2711: - snes - the nonlinear solver context
2713: Level: developer
2715: Note:
2716: Most users should have the `TS` created by calling `TSGetSNES()`
2718: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2719: @*/
2720: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2721: {
2722: PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);
2724: PetscFunctionBegin;
2727: PetscCall(PetscObjectReference((PetscObject)snes));
2728: PetscCall(SNESDestroy(&ts->snes));
2730: ts->snes = snes;
2732: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2733: PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2734: if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2735: PetscFunctionReturn(PETSC_SUCCESS);
2736: }
2738: /*@
2739: TSGetKSP - Returns the `KSP` (linear solver) associated with
2740: a `TS` (timestepper) context.
2742: Not Collective, but `ksp` is parallel if `ts` is parallel
2744: Input Parameter:
2745: . ts - the `TS` context obtained from `TSCreate()`
2747: Output Parameter:
2748: . ksp - the nonlinear solver context
2750: Level: beginner
2752: Notes:
2753: The user can then directly manipulate the `KSP` context to set various
2754: options, etc. Likewise, the user can then extract and manipulate the
2755: `PC` context as well.
2757: `TSGetKSP()` does not work for integrators that do not use `KSP`;
2758: in this case `TSGetKSP()` returns `NULL` in `ksp`.
2760: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2761: @*/
2762: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2763: {
2764: SNES snes;
2766: PetscFunctionBegin;
2768: PetscAssertPointer(ksp, 2);
2769: PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2770: PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
2771: PetscCall(TSGetSNES(ts, &snes));
2772: PetscCall(SNESGetKSP(snes, ksp));
2773: PetscFunctionReturn(PETSC_SUCCESS);
2774: }
2776: /* ----------- Routines to set solver parameters ---------- */
2778: /*@
2779: TSSetMaxSteps - Sets the maximum number of steps to use.
2781: Logically Collective
2783: Input Parameters:
2784: + ts - the `TS` context obtained from `TSCreate()`
2785: - maxsteps - maximum number of steps to use
2787: Options Database Key:
2788: . -ts_max_steps <maxsteps> - Sets maxsteps
2790: Level: intermediate
2792: Note:
2793: The default maximum number of steps is 5000
2795: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2796: @*/
2797: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2798: {
2799: PetscFunctionBegin;
2802: PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2803: ts->max_steps = maxsteps;
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: The default maximum time is 5.0
2848: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2849: @*/
2850: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2851: {
2852: PetscFunctionBegin;
2855: ts->max_time = maxtime;
2856: PetscFunctionReturn(PETSC_SUCCESS);
2857: }
2859: /*@
2860: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
2862: Not Collective
2864: Input Parameter:
2865: . ts - the `TS` context obtained from `TSCreate()`
2867: Output Parameter:
2868: . maxtime - final time to step to
2870: Level: advanced
2872: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2873: @*/
2874: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2875: {
2876: PetscFunctionBegin;
2878: PetscAssertPointer(maxtime, 2);
2879: *maxtime = ts->max_time;
2880: PetscFunctionReturn(PETSC_SUCCESS);
2881: }
2883: // PetscClangLinter pragma disable: -fdoc-*
2884: /*@
2885: TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.
2887: Level: deprecated
2889: @*/
2890: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2891: {
2892: PetscFunctionBegin;
2894: PetscCall(TSSetTime(ts, initial_time));
2895: PetscCall(TSSetTimeStep(ts, time_step));
2896: PetscFunctionReturn(PETSC_SUCCESS);
2897: }
2899: // PetscClangLinter pragma disable: -fdoc-*
2900: /*@
2901: TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.
2903: Level: deprecated
2905: @*/
2906: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2907: {
2908: PetscFunctionBegin;
2910: if (maxsteps) {
2911: PetscAssertPointer(maxsteps, 2);
2912: *maxsteps = ts->max_steps;
2913: }
2914: if (maxtime) {
2915: PetscAssertPointer(maxtime, 3);
2916: *maxtime = ts->max_time;
2917: }
2918: PetscFunctionReturn(PETSC_SUCCESS);
2919: }
2921: // PetscClangLinter pragma disable: -fdoc-*
2922: /*@
2923: TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.
2925: Level: deprecated
2927: @*/
2928: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2929: {
2930: PetscFunctionBegin;
2934: if (maxsteps >= 0) ts->max_steps = maxsteps;
2935: if (maxtime != (PetscReal)PETSC_DEFAULT) ts->max_time = maxtime;
2936: PetscFunctionReturn(PETSC_SUCCESS);
2937: }
2939: // PetscClangLinter pragma disable: -fdoc-*
2940: /*@
2941: TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.
2943: Level: deprecated
2945: @*/
2946: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
2947: {
2948: return TSGetStepNumber(ts, steps);
2949: }
2951: // PetscClangLinter pragma disable: -fdoc-*
2952: /*@
2953: TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.
2955: Level: deprecated
2957: @*/
2958: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
2959: {
2960: return TSGetStepNumber(ts, steps);
2961: }
2963: /*@
2964: TSSetSolution - Sets the initial solution vector
2965: for use by the `TS` routines.
2967: Logically Collective
2969: Input Parameters:
2970: + ts - the `TS` context obtained from `TSCreate()`
2971: - u - the solution vector
2973: Level: beginner
2975: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
2976: @*/
2977: PetscErrorCode TSSetSolution(TS ts, Vec u)
2978: {
2979: DM dm;
2981: PetscFunctionBegin;
2984: PetscCall(PetscObjectReference((PetscObject)u));
2985: PetscCall(VecDestroy(&ts->vec_sol));
2986: ts->vec_sol = u;
2988: PetscCall(TSGetDM(ts, &dm));
2989: PetscCall(DMShellSetGlobalVector(dm, u));
2990: PetscFunctionReturn(PETSC_SUCCESS);
2991: }
2993: /*@C
2994: TSSetPreStep - Sets the general-purpose function
2995: called once at the beginning of each time step.
2997: Logically Collective
2999: Input Parameters:
3000: + ts - The `TS` context obtained from `TSCreate()`
3001: - func - The function
3003: Calling sequence of `func`:
3004: . ts - the `TS` context
3006: Level: intermediate
3008: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3009: @*/
3010: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
3011: {
3012: PetscFunctionBegin;
3014: ts->prestep = func;
3015: PetscFunctionReturn(PETSC_SUCCESS);
3016: }
3018: /*@
3019: TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`
3021: Collective
3023: Input Parameter:
3024: . ts - The `TS` context obtained from `TSCreate()`
3026: Level: developer
3028: Note:
3029: `TSPreStep()` is typically used within time stepping implementations,
3030: so most users would not generally call this routine themselves.
3032: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3033: @*/
3034: PetscErrorCode TSPreStep(TS ts)
3035: {
3036: PetscFunctionBegin;
3038: if (ts->prestep) {
3039: Vec U;
3040: PetscObjectId idprev;
3041: PetscBool sameObject;
3042: PetscObjectState sprev, spost;
3044: PetscCall(TSGetSolution(ts, &U));
3045: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3046: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3047: PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3048: PetscCall(TSGetSolution(ts, &U));
3049: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3050: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3051: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3052: }
3053: PetscFunctionReturn(PETSC_SUCCESS);
3054: }
3056: /*@C
3057: TSSetPreStage - Sets the general-purpose function
3058: called once at the beginning of each stage.
3060: Logically Collective
3062: Input Parameters:
3063: + ts - The `TS` context obtained from `TSCreate()`
3064: - func - The function
3066: Calling sequence of `func`:
3067: + ts - the `TS` context
3068: - stagetime - the stage time
3070: Level: intermediate
3072: Note:
3073: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3074: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3075: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3077: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3078: @*/
3079: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime))
3080: {
3081: PetscFunctionBegin;
3083: ts->prestage = func;
3084: PetscFunctionReturn(PETSC_SUCCESS);
3085: }
3087: /*@C
3088: TSSetPostStage - Sets the general-purpose function, provided with `TSSetPostStep()`,
3089: called once at the end of each stage.
3091: Logically Collective
3093: Input Parameters:
3094: + ts - The `TS` context obtained from `TSCreate()`
3095: - func - The function
3097: Calling sequence of `func`:
3098: + ts - the `TS` context
3099: . stagetime - the stage time
3100: . stageindex - the stage index
3101: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3103: Level: intermediate
3105: Note:
3106: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3107: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3108: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3110: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3111: @*/
3112: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
3113: {
3114: PetscFunctionBegin;
3116: ts->poststage = func;
3117: PetscFunctionReturn(PETSC_SUCCESS);
3118: }
3120: /*@C
3121: TSSetPostEvaluate - Sets the general-purpose function
3122: called once at the end of each step evaluation.
3124: Logically Collective
3126: Input Parameters:
3127: + ts - The `TS` context obtained from `TSCreate()`
3128: - func - The function
3130: Calling sequence of `func`:
3131: . ts - the `TS` context
3133: Level: intermediate
3135: Note:
3136: Semantically, `TSSetPostEvaluate()` differs from `TSSetPostStep()` since the function it sets is called before event-handling
3137: thus guaranteeing the same solution (computed by the time-stepper) will be passed to it. On the other hand, `TSPostStep()`
3138: may be passed a different solution, possibly changed by the event handler. `TSPostEvaluate()` is called after the next step
3139: solution is evaluated allowing to modify it, if need be. The solution can be obtained with `TSGetSolution()`, the time step
3140: with `TSGetTimeStep()`, and the time at the start of the step is available via `TSGetTime()`
3142: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3143: @*/
3144: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
3145: {
3146: PetscFunctionBegin;
3148: ts->postevaluate = func;
3149: PetscFunctionReturn(PETSC_SUCCESS);
3150: }
3152: /*@
3153: TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`
3155: Collective
3157: Input Parameters:
3158: + ts - The `TS` context obtained from `TSCreate()`
3159: - stagetime - The absolute time of the current stage
3161: Level: developer
3163: Note:
3164: `TSPreStage()` is typically used within time stepping implementations,
3165: most users would not generally call this routine themselves.
3167: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3168: @*/
3169: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3170: {
3171: PetscFunctionBegin;
3173: if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3174: PetscFunctionReturn(PETSC_SUCCESS);
3175: }
3177: /*@
3178: TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`
3180: Collective
3182: Input Parameters:
3183: + ts - The `TS` context obtained from `TSCreate()`
3184: . stagetime - The absolute time of the current stage
3185: . stageindex - Stage number
3186: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3188: Level: developer
3190: Note:
3191: `TSPostStage()` is typically used within time stepping implementations,
3192: most users would not generally call this routine themselves.
3194: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3195: @*/
3196: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3197: {
3198: PetscFunctionBegin;
3200: if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3201: PetscFunctionReturn(PETSC_SUCCESS);
3202: }
3204: /*@
3205: TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`
3207: Collective
3209: Input Parameter:
3210: . ts - The `TS` context obtained from `TSCreate()`
3212: Level: developer
3214: Note:
3215: `TSPostEvaluate()` is typically used within time stepping implementations,
3216: most users would not generally call this routine themselves.
3218: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3219: @*/
3220: PetscErrorCode TSPostEvaluate(TS ts)
3221: {
3222: PetscFunctionBegin;
3224: if (ts->postevaluate) {
3225: Vec U;
3226: PetscObjectState sprev, spost;
3228: PetscCall(TSGetSolution(ts, &U));
3229: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3230: PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3231: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3232: if (sprev != spost) PetscCall(TSRestartStep(ts));
3233: }
3234: PetscFunctionReturn(PETSC_SUCCESS);
3235: }
3237: /*@C
3238: TSSetPostStep - Sets the general-purpose function
3239: called once at the end of each time step.
3241: Logically Collective
3243: Input Parameters:
3244: + ts - The `TS` context obtained from `TSCreate()`
3245: - func - The function
3247: Calling sequence of `func`:
3248: . ts - the `TS` context
3250: Level: intermediate
3252: Note:
3253: The function set by `TSSetPostStep()` is called after each successful step. The solution vector
3254: obtained by `TSGetSolution()` may be different than that computed at the step end if the event handler
3255: locates an event and `TSPostEvent()` modifies it. Use `TSSetPostEvaluate()` if an unmodified solution is needed instead.
3257: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3258: @*/
3259: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
3260: {
3261: PetscFunctionBegin;
3263: ts->poststep = func;
3264: PetscFunctionReturn(PETSC_SUCCESS);
3265: }
3267: /*@
3268: TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`
3270: Collective
3272: Input Parameter:
3273: . ts - The `TS` context obtained from `TSCreate()`
3275: Note:
3276: `TSPostStep()` is typically used within time stepping implementations,
3277: so most users would not generally call this routine themselves.
3279: Level: developer
3281: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()`
3282: @*/
3283: PetscErrorCode TSPostStep(TS ts)
3284: {
3285: PetscFunctionBegin;
3287: if (ts->poststep) {
3288: Vec U;
3289: PetscObjectId idprev;
3290: PetscBool sameObject;
3291: PetscObjectState sprev, spost;
3293: PetscCall(TSGetSolution(ts, &U));
3294: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3295: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3296: PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3297: PetscCall(TSGetSolution(ts, &U));
3298: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3299: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3300: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3301: }
3302: PetscFunctionReturn(PETSC_SUCCESS);
3303: }
3305: /*@
3306: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3308: Collective
3310: Input Parameters:
3311: + ts - time stepping context
3312: - t - time to interpolate to
3314: Output Parameter:
3315: . U - state at given time
3317: Level: intermediate
3319: Developer Notes:
3320: `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3322: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3323: @*/
3324: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3325: {
3326: PetscFunctionBegin;
3329: 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);
3330: PetscUseTypeMethod(ts, interpolate, t, U);
3331: PetscFunctionReturn(PETSC_SUCCESS);
3332: }
3334: /*@
3335: TSStep - Steps one time step
3337: Collective
3339: Input Parameter:
3340: . ts - the `TS` context obtained from `TSCreate()`
3342: Level: developer
3344: Notes:
3345: The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.
3347: The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
3348: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3350: This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
3351: time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.
3353: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3354: @*/
3355: PetscErrorCode TSStep(TS ts)
3356: {
3357: static PetscBool cite = PETSC_FALSE;
3358: PetscReal ptime;
3360: PetscFunctionBegin;
3362: PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3363: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3364: " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3365: " journal = {arXiv e-preprints},\n"
3366: " eprint = {1806.01437},\n"
3367: " archivePrefix = {arXiv},\n"
3368: " year = {2018}\n}\n",
3369: &cite));
3370: PetscCall(TSSetUp(ts));
3371: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3372: if (ts->tspan)
3373: ts->tspan->worktol = 0; /* In each step of TSSolve() 'tspan->worktol' will be meaningfully defined (later) only once:
3374: in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */
3376: PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_MAX_INT, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
3377: 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()");
3378: 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");
3380: if (!ts->vec_sol0) PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sol0));
3381: PetscCall(VecCopy(ts->vec_sol, ts->vec_sol0));
3382: ts->time_step0 = ts->time_step;
3384: if (!ts->steps) ts->ptime_prev = ts->ptime;
3385: ptime = ts->ptime;
3387: ts->ptime_prev_rollback = ts->ptime_prev;
3388: ts->reason = TS_CONVERGED_ITERATING;
3390: PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3391: PetscUseTypeMethod(ts, step);
3392: PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));
3394: if (ts->reason >= 0) {
3395: ts->ptime_prev = ptime;
3396: ts->steps++;
3397: ts->steprollback = PETSC_FALSE;
3398: ts->steprestart = PETSC_FALSE;
3399: }
3401: if (ts->reason < 0 && ts->errorifstepfailed) {
3402: PetscCall(TSMonitorCancel(ts));
3403: 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 make negative to attempt recovery", TSConvergedReasons[ts->reason]);
3404: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3405: }
3406: PetscFunctionReturn(PETSC_SUCCESS);
3407: }
3409: /*@
3410: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3411: at the end of a time step with a given order of accuracy.
3413: Collective
3415: Input Parameters:
3416: + ts - time stepping context
3417: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3419: Input/Output Parameter:
3420: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3421: on output, the actual order of the error evaluation
3423: Output Parameter:
3424: . wlte - the weighted local truncation error norm
3426: Level: advanced
3428: Note:
3429: If the timestepper cannot evaluate the error in a particular step
3430: (eg. in the first step or restart steps after event handling),
3431: this routine returns wlte=-1.0 .
3433: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3434: @*/
3435: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3436: {
3437: PetscFunctionBegin;
3441: if (order) PetscAssertPointer(order, 3);
3443: PetscAssertPointer(wlte, 4);
3444: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3445: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3446: PetscFunctionReturn(PETSC_SUCCESS);
3447: }
3449: /*@
3450: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3452: Collective
3454: Input Parameters:
3455: + ts - time stepping context
3456: . order - desired order of accuracy
3457: - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)
3459: Output Parameter:
3460: . U - state at the end of the current step
3462: Level: advanced
3464: Notes:
3465: This function cannot be called until all stages have been evaluated.
3467: 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.
3469: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3470: @*/
3471: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3472: {
3473: PetscFunctionBegin;
3477: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3478: PetscFunctionReturn(PETSC_SUCCESS);
3479: }
3481: /*@C
3482: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3484: Not collective
3486: Input Parameter:
3487: . ts - time stepping context
3489: Output Parameter:
3490: . initCondition - The function which computes an initial condition
3492: Calling sequence of `initCondition`:
3493: + ts - The timestepping context
3494: - u - The input vector in which the initial condition is stored
3496: Level: advanced
3498: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3499: @*/
3500: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3501: {
3502: PetscFunctionBegin;
3504: PetscAssertPointer(initCondition, 2);
3505: *initCondition = ts->ops->initcondition;
3506: PetscFunctionReturn(PETSC_SUCCESS);
3507: }
3509: /*@C
3510: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3512: Logically collective
3514: Input Parameters:
3515: + ts - time stepping context
3516: - initCondition - The function which computes an initial condition
3518: Calling sequence of `initCondition`:
3519: + ts - The timestepping context
3520: - e - The input vector in which the initial condition is to be stored
3522: Level: advanced
3524: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3525: @*/
3526: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3527: {
3528: PetscFunctionBegin;
3531: ts->ops->initcondition = initCondition;
3532: PetscFunctionReturn(PETSC_SUCCESS);
3533: }
3535: /*@
3536: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3538: Collective
3540: Input Parameters:
3541: + ts - time stepping context
3542: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3544: Level: advanced
3546: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3547: @*/
3548: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3549: {
3550: PetscFunctionBegin;
3553: PetscTryTypeMethod(ts, initcondition, u);
3554: PetscFunctionReturn(PETSC_SUCCESS);
3555: }
3557: /*@C
3558: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3560: Not collective
3562: Input Parameter:
3563: . ts - time stepping context
3565: Output Parameter:
3566: . exactError - The function which computes the solution error
3568: Calling sequence of `exactError`:
3569: + ts - The timestepping context
3570: . u - The approximate solution vector
3571: - e - The vector in which the error is stored
3573: Level: advanced
3575: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3576: @*/
3577: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3578: {
3579: PetscFunctionBegin;
3581: PetscAssertPointer(exactError, 2);
3582: *exactError = ts->ops->exacterror;
3583: PetscFunctionReturn(PETSC_SUCCESS);
3584: }
3586: /*@C
3587: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3589: Logically collective
3591: Input Parameters:
3592: + ts - time stepping context
3593: - exactError - The function which computes the solution error
3595: Calling sequence of `exactError`:
3596: + ts - The timestepping context
3597: . u - The approximate solution vector
3598: - e - The vector in which the error is stored
3600: Level: advanced
3602: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3603: @*/
3604: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3605: {
3606: PetscFunctionBegin;
3609: ts->ops->exacterror = exactError;
3610: PetscFunctionReturn(PETSC_SUCCESS);
3611: }
3613: /*@
3614: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3616: Collective
3618: Input Parameters:
3619: + ts - time stepping context
3620: . u - The approximate solution
3621: - e - The `Vec` used to store the error
3623: Level: advanced
3625: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3626: @*/
3627: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3628: {
3629: PetscFunctionBegin;
3633: PetscTryTypeMethod(ts, exacterror, u, e);
3634: PetscFunctionReturn(PETSC_SUCCESS);
3635: }
3637: /*@C
3638: TSSetResize - Sets the resize callbacks.
3640: Logically Collective
3642: Input Parameters:
3643: + ts - The `TS` context obtained from `TSCreate()`
3644: . rollback - Whether a resize will restart the step
3645: . setup - The setup function
3646: . transfer - The transfer function
3647: - ctx - [optional] The user-defined context
3649: Calling sequence of `setup`:
3650: + ts - the `TS` context
3651: . step - the current step
3652: . time - the current time
3653: . state - the current vector of state
3654: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3655: - ctx - user defined context
3657: Calling sequence of `transfer`:
3658: + ts - the `TS` context
3659: . nv - the number of vectors to be transferred
3660: . vecsin - array of vectors to be transferred
3661: . vecsout - array of transferred vectors
3662: - ctx - user defined context
3664: Notes:
3665: The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
3666: depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
3667: and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
3668: that the size of the discrete problem has changed.
3669: In both cases, the solver will collect the needed vectors that will be
3670: transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
3671: current solution vector, and other vectors needed by the specific solver used.
3672: For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3673: Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3674: will be automatically reset if the sizes are changed and they must be specified again by the user
3675: inside the `transfer` function.
3676: The input and output arrays passed to `transfer` are allocated by PETSc.
3677: Vectors in `vecsout` must be created by the user.
3678: Ownership of vectors in `vecsout` is transferred to PETSc.
3680: Level: advanced
3682: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3683: @*/
3684: 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)
3685: {
3686: PetscFunctionBegin;
3688: ts->resizerollback = rollback;
3689: ts->resizesetup = setup;
3690: ts->resizetransfer = transfer;
3691: ts->resizectx = ctx;
3692: PetscFunctionReturn(PETSC_SUCCESS);
3693: }
3695: /*
3696: TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.
3698: Collective
3700: Input Parameters:
3701: + ts - The `TS` context obtained from `TSCreate()`
3702: - 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.
3704: Level: developer
3706: Note:
3707: `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3708: used within time stepping implementations,
3709: so most users would not generally call this routine themselves.
3711: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3712: @*/
3713: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3714: {
3715: PetscFunctionBegin;
3717: PetscTryTypeMethod(ts, resizeregister, flg);
3718: /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3719: PetscFunctionReturn(PETSC_SUCCESS);
3720: }
3722: static PetscErrorCode TSResizeReset(TS ts)
3723: {
3724: PetscFunctionBegin;
3726: PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3727: PetscFunctionReturn(PETSC_SUCCESS);
3728: }
3730: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3731: {
3732: PetscFunctionBegin;
3735: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3736: if (ts->resizetransfer) {
3737: PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3738: PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3739: }
3740: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3741: PetscFunctionReturn(PETSC_SUCCESS);
3742: }
3744: /*@C
3745: TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.
3747: Collective
3749: Input Parameters:
3750: + ts - The `TS` context obtained from `TSCreate()`
3751: . name - A string identifying the vector
3752: - vec - The vector
3754: Level: developer
3756: Note:
3757: `TSResizeRegisterVec()` is typically used within time stepping implementations,
3758: so most users would not generally call this routine themselves.
3760: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3761: @*/
3762: PetscErrorCode TSResizeRegisterVec(TS ts, const char *name, Vec vec)
3763: {
3764: PetscFunctionBegin;
3766: PetscAssertPointer(name, 2);
3768: PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3769: PetscFunctionReturn(PETSC_SUCCESS);
3770: }
3772: /*@C
3773: TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.
3775: Collective
3777: Input Parameters:
3778: + ts - The `TS` context obtained from `TSCreate()`
3779: . name - A string identifying the vector
3780: - vec - The vector
3782: Level: developer
3784: Note:
3785: `TSResizeRetrieveVec()` is typically used within time stepping implementations,
3786: so most users would not generally call this routine themselves.
3788: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
3789: @*/
3790: PetscErrorCode TSResizeRetrieveVec(TS ts, const char *name, Vec *vec)
3791: {
3792: PetscFunctionBegin;
3794: PetscAssertPointer(name, 2);
3795: PetscAssertPointer(vec, 3);
3796: PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
3797: PetscFunctionReturn(PETSC_SUCCESS);
3798: }
3800: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
3801: {
3802: PetscInt cnt;
3803: PetscObjectList tmp;
3804: Vec *vecsin = NULL;
3805: const char **namesin = NULL;
3807: PetscFunctionBegin;
3808: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
3809: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
3810: if (names) PetscCall(PetscMalloc1(cnt, &namesin));
3811: if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
3812: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
3813: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
3814: if (vecs) vecsin[cnt] = (Vec)tmp->obj;
3815: if (names) namesin[cnt] = tmp->name;
3816: cnt++;
3817: }
3818: }
3819: if (nv) *nv = cnt;
3820: if (names) *names = namesin;
3821: if (vecs) *vecs = vecsin;
3822: PetscFunctionReturn(PETSC_SUCCESS);
3823: }
3825: /*@
3826: TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`
3828: Collective
3830: Input Parameter:
3831: . ts - The `TS` context obtained from `TSCreate()`
3833: Level: developer
3835: Note:
3836: `TSResize()` is typically used within time stepping implementations,
3837: so most users would not generally call this routine themselves.
3839: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3840: @*/
3841: PetscErrorCode TSResize(TS ts)
3842: {
3843: PetscInt nv = 0;
3844: const char **names = NULL;
3845: Vec *vecsin = NULL;
3846: const char *solname = "ts:vec_sol";
3848: PetscFunctionBegin;
3850: if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
3851: if (ts->resizesetup) {
3852: PetscBool flg = PETSC_FALSE;
3854: PetscCall(VecLockReadPush(ts->vec_sol));
3855: PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &flg, ts->resizectx));
3856: PetscCall(VecLockReadPop(ts->vec_sol));
3857: if (flg) {
3858: if (ts->resizerollback) {
3859: PetscCall(TSRollBack(ts));
3860: ts->time_step = ts->time_step0;
3861: }
3862: PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
3863: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
3864: }
3865: }
3867: PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
3868: if (nv) {
3869: Vec *vecsout, vecsol;
3871: /* Reset internal objects */
3872: PetscCall(TSReset(ts));
3874: /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
3875: PetscCall(PetscCalloc1(nv, &vecsout));
3876: PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
3877: for (PetscInt i = 0; i < nv; i++) {
3878: const char *name;
3879: char *oname;
3881: PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
3882: PetscCall(PetscStrallocpy(name, &oname));
3883: PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
3884: if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
3885: PetscCall(PetscFree(oname));
3886: PetscCall(VecDestroy(&vecsout[i]));
3887: }
3888: PetscCall(PetscFree(vecsout));
3889: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */
3891: PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
3892: if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
3893: PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
3894: }
3896: PetscCall(PetscFree(names));
3897: PetscCall(PetscFree(vecsin));
3898: PetscCall(TSResizeReset(ts));
3899: PetscFunctionReturn(PETSC_SUCCESS);
3900: }
3902: /*@
3903: TSSolve - Steps the requested number of timesteps.
3905: Collective
3907: Input Parameters:
3908: + ts - the `TS` context obtained from `TSCreate()`
3909: - u - the solution vector (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
3910: otherwise must contain the initial conditions and will contain the solution at the final requested time
3912: Level: beginner
3914: Notes:
3915: The final time returned by this function may be different from the time of the internally
3916: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
3917: stepped over the final time.
3919: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3920: @*/
3921: PetscErrorCode TSSolve(TS ts, Vec u)
3922: {
3923: Vec solution;
3925: PetscFunctionBegin;
3929: PetscCall(TSSetExactFinalTimeDefault(ts));
3930: 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 */
3931: if (!ts->vec_sol || u == ts->vec_sol) {
3932: PetscCall(VecDuplicate(u, &solution));
3933: PetscCall(TSSetSolution(ts, solution));
3934: PetscCall(VecDestroy(&solution)); /* grant ownership */
3935: }
3936: PetscCall(VecCopy(u, ts->vec_sol));
3937: PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
3938: } else if (u) PetscCall(TSSetSolution(ts, u));
3939: PetscCall(TSSetUp(ts));
3940: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3942: PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_MAX_INT, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
3943: 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()");
3944: 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");
3945: 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");
3947: 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 */
3948: PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]));
3949: ts->tspan->spanctr = 1;
3950: }
3952: if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));
3954: /* reset number of steps only when the step is not restarted. ARKIMEX
3955: restarts the step after an event. Resetting these counters in such case causes
3956: TSTrajectory to incorrectly save the output files
3957: */
3958: /* reset time step and iteration counters */
3959: if (!ts->steps) {
3960: ts->ksp_its = 0;
3961: ts->snes_its = 0;
3962: ts->num_snes_failures = 0;
3963: ts->reject = 0;
3964: ts->steprestart = PETSC_TRUE;
3965: ts->steprollback = PETSC_FALSE;
3966: ts->rhsjacobian.time = PETSC_MIN_REAL;
3967: }
3969: /* make sure initial time step does not overshoot final time or the next point in tspan */
3970: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
3971: PetscReal maxdt;
3972: PetscReal dt = ts->time_step;
3974: if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
3975: else maxdt = ts->max_time - ts->ptime;
3976: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
3977: }
3978: ts->reason = TS_CONVERGED_ITERATING;
3980: {
3981: PetscViewer viewer;
3982: PetscViewerFormat format;
3983: PetscBool flg;
3984: static PetscBool incall = PETSC_FALSE;
3986: if (!incall) {
3987: /* Estimate the convergence rate of the time discretization */
3988: PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
3989: if (flg) {
3990: PetscConvEst conv;
3991: DM dm;
3992: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
3993: PetscInt Nf;
3994: PetscBool checkTemporal = PETSC_TRUE;
3996: incall = PETSC_TRUE;
3997: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
3998: PetscCall(TSGetDM(ts, &dm));
3999: PetscCall(DMGetNumFields(dm, &Nf));
4000: PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
4001: PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
4002: PetscCall(PetscConvEstUseTS(conv, checkTemporal));
4003: PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
4004: PetscCall(PetscConvEstSetFromOptions(conv));
4005: PetscCall(PetscConvEstSetUp(conv));
4006: PetscCall(PetscConvEstGetConvRate(conv, alpha));
4007: PetscCall(PetscViewerPushFormat(viewer, format));
4008: PetscCall(PetscConvEstRateView(conv, alpha, viewer));
4009: PetscCall(PetscViewerPopFormat(viewer));
4010: PetscCall(PetscOptionsRestoreViewer(&viewer));
4011: PetscCall(PetscConvEstDestroy(&conv));
4012: PetscCall(PetscFree(alpha));
4013: incall = PETSC_FALSE;
4014: }
4015: }
4016: }
4018: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));
4020: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4021: PetscUseTypeMethod(ts, solve);
4022: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4023: ts->solvetime = ts->ptime;
4024: solution = ts->vec_sol;
4025: } else { /* Step the requested number of timesteps. */
4026: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4027: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4029: if (!ts->steps) {
4030: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4031: PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4032: }
4034: while (!ts->reason) {
4035: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4036: if (!ts->steprollback) PetscCall(TSPreStep(ts));
4037: PetscCall(TSStep(ts));
4038: if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4039: if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4040: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4041: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4042: PetscCall(TSForwardCostIntegral(ts));
4043: if (ts->reason >= 0) ts->steps++;
4044: }
4045: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4046: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4047: PetscCall(TSForwardStep(ts));
4048: if (ts->reason >= 0) ts->steps++;
4049: }
4050: PetscCall(TSPostEvaluate(ts));
4051: 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. */
4052: if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4053: if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
4054: /* check convergence */
4055: if (!ts->reason) {
4056: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4057: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4058: }
4059: if (!ts->steprollback) {
4060: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4061: PetscCall(TSPostStep(ts));
4062: if (!ts->resizerollback) PetscCall(TSResize(ts));
4064: if (ts->tspan && ts->tspan->spanctr < ts->tspan->num_span_times) {
4065: PetscCheck(ts->tspan->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(tspan->worktol > 0) in TSSolve()");
4066: 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++]));
4067: }
4068: }
4069: }
4070: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4072: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4073: if (!u) u = ts->vec_sol;
4074: PetscCall(TSInterpolate(ts, ts->max_time, u));
4075: ts->solvetime = ts->max_time;
4076: solution = u;
4077: PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4078: } else {
4079: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4080: ts->solvetime = ts->ptime;
4081: solution = ts->vec_sol;
4082: }
4083: }
4085: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4086: PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4087: PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4088: if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4089: PetscFunctionReturn(PETSC_SUCCESS);
4090: }
4092: /*@
4093: TSGetTime - Gets the time of the most recently completed step.
4095: Not Collective
4097: Input Parameter:
4098: . ts - the `TS` context obtained from `TSCreate()`
4100: Output Parameter:
4101: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
4103: Level: beginner
4105: Note:
4106: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
4107: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
4109: .seealso: [](ch_ts), `TS`, ``TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4110: @*/
4111: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4112: {
4113: PetscFunctionBegin;
4115: PetscAssertPointer(t, 2);
4116: *t = ts->ptime;
4117: PetscFunctionReturn(PETSC_SUCCESS);
4118: }
4120: /*@
4121: TSGetPrevTime - Gets the starting time of the previously completed step.
4123: Not Collective
4125: Input Parameter:
4126: . ts - the `TS` context obtained from `TSCreate()`
4128: Output Parameter:
4129: . t - the previous time
4131: Level: beginner
4133: .seealso: [](ch_ts), `TS`, ``TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4134: @*/
4135: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4136: {
4137: PetscFunctionBegin;
4139: PetscAssertPointer(t, 2);
4140: *t = ts->ptime_prev;
4141: PetscFunctionReturn(PETSC_SUCCESS);
4142: }
4144: /*@
4145: TSSetTime - Allows one to reset the time.
4147: Logically Collective
4149: Input Parameters:
4150: + ts - the `TS` context obtained from `TSCreate()`
4151: - t - the time
4153: Level: intermediate
4155: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4156: @*/
4157: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4158: {
4159: PetscFunctionBegin;
4162: ts->ptime = t;
4163: PetscFunctionReturn(PETSC_SUCCESS);
4164: }
4166: /*@C
4167: TSSetOptionsPrefix - Sets the prefix used for searching for all
4168: TS options in the database.
4170: Logically Collective
4172: Input Parameters:
4173: + ts - The `TS` context
4174: - prefix - The prefix to prepend to all option names
4176: Level: advanced
4178: Note:
4179: A hyphen (-) must NOT be given at the beginning of the prefix name.
4180: The first character of all runtime options is AUTOMATICALLY the
4181: hyphen.
4183: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4184: @*/
4185: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4186: {
4187: SNES snes;
4189: PetscFunctionBegin;
4191: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4192: PetscCall(TSGetSNES(ts, &snes));
4193: PetscCall(SNESSetOptionsPrefix(snes, prefix));
4194: PetscFunctionReturn(PETSC_SUCCESS);
4195: }
4197: /*@C
4198: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4199: TS options in the database.
4201: Logically Collective
4203: Input Parameters:
4204: + ts - The `TS` context
4205: - prefix - The prefix to prepend to all option names
4207: Level: advanced
4209: Note:
4210: A hyphen (-) must NOT be given at the beginning of the prefix name.
4211: The first character of all runtime options is AUTOMATICALLY the
4212: hyphen.
4214: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4215: @*/
4216: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4217: {
4218: SNES snes;
4220: PetscFunctionBegin;
4222: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4223: PetscCall(TSGetSNES(ts, &snes));
4224: PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4225: PetscFunctionReturn(PETSC_SUCCESS);
4226: }
4228: /*@C
4229: TSGetOptionsPrefix - Sets the prefix used for searching for all
4230: `TS` options in the database.
4232: Not Collective
4234: Input Parameter:
4235: . ts - The `TS` context
4237: Output Parameter:
4238: . prefix - A pointer to the prefix string used
4240: Level: intermediate
4242: Fortran Notes:
4243: The user should pass in a string 'prefix' of
4244: sufficient length to hold the prefix.
4246: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4247: @*/
4248: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4249: {
4250: PetscFunctionBegin;
4252: PetscAssertPointer(prefix, 2);
4253: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4254: PetscFunctionReturn(PETSC_SUCCESS);
4255: }
4257: /*@C
4258: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4260: Not Collective, but parallel objects are returned if ts is parallel
4262: Input Parameter:
4263: . ts - The `TS` context obtained from `TSCreate()`
4265: Output Parameters:
4266: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`)
4267: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`)
4268: . func - Function to compute the Jacobian of the RHS (or `NULL`)
4269: - ctx - User-defined context for Jacobian evaluation routine (or `NULL`)
4271: Level: intermediate
4273: Note:
4274: You can pass in `NULL` for any return argument you do not need.
4276: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4278: @*/
4279: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, void **ctx)
4280: {
4281: DM dm;
4283: PetscFunctionBegin;
4284: if (Amat || Pmat) {
4285: SNES snes;
4286: PetscCall(TSGetSNES(ts, &snes));
4287: PetscCall(SNESSetUpMatrices(snes));
4288: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4289: }
4290: PetscCall(TSGetDM(ts, &dm));
4291: PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4292: PetscFunctionReturn(PETSC_SUCCESS);
4293: }
4295: /*@C
4296: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4298: Not Collective, but parallel objects are returned if ts is parallel
4300: Input Parameter:
4301: . ts - The `TS` context obtained from `TSCreate()`
4303: Output Parameters:
4304: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4305: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4306: . f - The function to compute the matrices
4307: - ctx - User-defined context for Jacobian evaluation routine
4309: Level: advanced
4311: Note:
4312: You can pass in `NULL` for any return argument you do not need.
4314: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4315: @*/
4316: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, void **ctx)
4317: {
4318: DM dm;
4320: PetscFunctionBegin;
4321: if (Amat || Pmat) {
4322: SNES snes;
4323: PetscCall(TSGetSNES(ts, &snes));
4324: PetscCall(SNESSetUpMatrices(snes));
4325: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4326: }
4327: PetscCall(TSGetDM(ts, &dm));
4328: PetscCall(DMTSGetIJacobian(dm, f, ctx));
4329: PetscFunctionReturn(PETSC_SUCCESS);
4330: }
4332: #include <petsc/private/dmimpl.h>
4333: /*@
4334: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4336: Logically Collective
4338: Input Parameters:
4339: + ts - the `TS` integrator object
4340: - dm - the dm, cannot be `NULL`
4342: Level: intermediate
4344: Notes:
4345: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4346: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4347: different problems using the same function space.
4349: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4350: @*/
4351: PetscErrorCode TSSetDM(TS ts, DM dm)
4352: {
4353: SNES snes;
4354: DMTS tsdm;
4356: PetscFunctionBegin;
4359: PetscCall(PetscObjectReference((PetscObject)dm));
4360: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4361: if (ts->dm->dmts && !dm->dmts) {
4362: PetscCall(DMCopyDMTS(ts->dm, dm));
4363: PetscCall(DMGetDMTS(ts->dm, &tsdm));
4364: /* Grant write privileges to the replacement DM */
4365: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4366: }
4367: PetscCall(DMDestroy(&ts->dm));
4368: }
4369: ts->dm = dm;
4371: PetscCall(TSGetSNES(ts, &snes));
4372: PetscCall(SNESSetDM(snes, dm));
4373: PetscFunctionReturn(PETSC_SUCCESS);
4374: }
4376: /*@
4377: TSGetDM - Gets the `DM` that may be used by some preconditioners
4379: Not Collective
4381: Input Parameter:
4382: . ts - the `TS`
4384: Output Parameter:
4385: . dm - the `DM`
4387: Level: intermediate
4389: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4390: @*/
4391: PetscErrorCode TSGetDM(TS ts, DM *dm)
4392: {
4393: PetscFunctionBegin;
4395: if (!ts->dm) {
4396: PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4397: if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4398: }
4399: *dm = ts->dm;
4400: PetscFunctionReturn(PETSC_SUCCESS);
4401: }
4403: /*@
4404: SNESTSFormFunction - Function to evaluate nonlinear residual
4406: Logically Collective
4408: Input Parameters:
4409: + snes - nonlinear solver
4410: . U - the current state at which to evaluate the residual
4411: - ctx - user context, must be a TS
4413: Output Parameter:
4414: . F - the nonlinear residual
4416: Level: advanced
4418: Note:
4419: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4420: It is most frequently passed to `MatFDColoringSetFunction()`.
4422: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4423: @*/
4424: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4425: {
4426: TS ts = (TS)ctx;
4428: PetscFunctionBegin;
4433: PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4434: PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4435: PetscFunctionReturn(PETSC_SUCCESS);
4436: }
4438: /*@
4439: SNESTSFormJacobian - Function to evaluate the Jacobian
4441: 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 Parameters:
4449: + A - the Jacobian
4450: - B - the preconditioning matrix (may be the same as A)
4452: Level: developer
4454: Note:
4455: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4457: .seealso: [](ch_ts), `SNESSetJacobian()`
4458: @*/
4459: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4460: {
4461: TS ts = (TS)ctx;
4463: PetscFunctionBegin;
4469: PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4470: PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4471: PetscFunctionReturn(PETSC_SUCCESS);
4472: }
4474: /*@C
4475: TSComputeRHSFunctionLinear - Evaluate the right-hand side via the user-provided Jacobian, for linear problems Udot = A U only
4477: Collective
4479: Input Parameters:
4480: + ts - time stepping context
4481: . t - time at which to evaluate
4482: . U - state at which to evaluate
4483: - ctx - context
4485: Output Parameter:
4486: . F - right-hand side
4488: Level: intermediate
4490: Note:
4491: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right-hand side for linear problems.
4492: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4494: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4495: @*/
4496: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4497: {
4498: Mat Arhs, Brhs;
4500: PetscFunctionBegin;
4501: PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4502: /* undo the damage caused by shifting */
4503: PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4504: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4505: PetscCall(MatMult(Arhs, U, F));
4506: PetscFunctionReturn(PETSC_SUCCESS);
4507: }
4509: /*@C
4510: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
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 Parameters:
4521: + A - pointer to operator
4522: - B - pointer to preconditioning matrix
4524: Level: intermediate
4526: Note:
4527: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4529: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4530: @*/
4531: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4532: {
4533: PetscFunctionBegin;
4534: PetscFunctionReturn(PETSC_SUCCESS);
4535: }
4537: /*@C
4538: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4540: Collective
4542: Input Parameters:
4543: + ts - time stepping context
4544: . t - time at which to evaluate
4545: . U - state at which to evaluate
4546: . Udot - time derivative of state vector
4547: - ctx - context
4549: Output Parameter:
4550: . F - left hand side
4552: Level: intermediate
4554: Notes:
4555: 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
4556: user is required to write their own `TSComputeIFunction()`.
4557: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4558: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4560: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4562: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4563: @*/
4564: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4565: {
4566: Mat A, B;
4568: PetscFunctionBegin;
4569: PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4570: PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4571: PetscCall(MatMult(A, Udot, F));
4572: PetscFunctionReturn(PETSC_SUCCESS);
4573: }
4575: /*@C
4576: TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE
4578: Collective
4580: Input Parameters:
4581: + ts - time stepping context
4582: . t - time at which to evaluate
4583: . U - state at which to evaluate
4584: . Udot - time derivative of state vector
4585: . shift - shift to apply
4586: - ctx - context
4588: Output Parameters:
4589: + A - pointer to operator
4590: - B - pointer to matrix from which the preconditioner is built (often `A`)
4592: Level: advanced
4594: Notes:
4595: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4597: It is only appropriate for problems of the form
4599: $$
4600: M \dot{U} = F(U,t)
4601: $$
4603: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4604: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4605: an implicit operator of the form
4607: $$
4608: shift*M + J
4609: $$
4611: 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
4612: a copy of M or reassemble it when requested.
4614: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4615: @*/
4616: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4617: {
4618: PetscFunctionBegin;
4619: PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4620: ts->ijacobian.shift = shift;
4621: PetscFunctionReturn(PETSC_SUCCESS);
4622: }
4624: /*@
4625: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4627: Not Collective
4629: Input Parameter:
4630: . ts - the `TS` context
4632: Output Parameter:
4633: . equation_type - see `TSEquationType`
4635: Level: beginner
4637: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4638: @*/
4639: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4640: {
4641: PetscFunctionBegin;
4643: PetscAssertPointer(equation_type, 2);
4644: *equation_type = ts->equation_type;
4645: PetscFunctionReturn(PETSC_SUCCESS);
4646: }
4648: /*@
4649: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4651: Not Collective
4653: Input Parameters:
4654: + ts - the `TS` context
4655: - equation_type - see `TSEquationType`
4657: Level: advanced
4659: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4660: @*/
4661: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4662: {
4663: PetscFunctionBegin;
4665: ts->equation_type = equation_type;
4666: PetscFunctionReturn(PETSC_SUCCESS);
4667: }
4669: /*@
4670: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4672: Not Collective
4674: Input Parameter:
4675: . ts - the `TS` context
4677: Output Parameter:
4678: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4679: manual pages for the individual convergence tests for complete lists
4681: Level: beginner
4683: Note:
4684: Can only be called after the call to `TSSolve()` is complete.
4686: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4687: @*/
4688: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4689: {
4690: PetscFunctionBegin;
4692: PetscAssertPointer(reason, 2);
4693: *reason = ts->reason;
4694: PetscFunctionReturn(PETSC_SUCCESS);
4695: }
4697: /*@
4698: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4700: Logically Collective; reason must contain common value
4702: Input Parameters:
4703: + ts - the `TS` context
4704: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4705: manual pages for the individual convergence tests for complete lists
4707: Level: advanced
4709: Note:
4710: Can only be called while `TSSolve()` is active.
4712: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4713: @*/
4714: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4715: {
4716: PetscFunctionBegin;
4718: ts->reason = reason;
4719: PetscFunctionReturn(PETSC_SUCCESS);
4720: }
4722: /*@
4723: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4725: Not Collective
4727: Input Parameter:
4728: . ts - the `TS` context
4730: Output Parameter:
4731: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4733: Level: beginner
4735: Note:
4736: Can only be called after the call to `TSSolve()` is complete.
4738: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4739: @*/
4740: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4741: {
4742: PetscFunctionBegin;
4744: PetscAssertPointer(ftime, 2);
4745: *ftime = ts->solvetime;
4746: PetscFunctionReturn(PETSC_SUCCESS);
4747: }
4749: /*@
4750: TSGetSNESIterations - Gets the total number of nonlinear iterations
4751: used by the time integrator.
4753: Not Collective
4755: Input Parameter:
4756: . ts - `TS` context
4758: Output Parameter:
4759: . nits - number of nonlinear iterations
4761: Level: intermediate
4763: Note:
4764: This counter is reset to zero for each successive call to `TSSolve()`.
4766: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4767: @*/
4768: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4769: {
4770: PetscFunctionBegin;
4772: PetscAssertPointer(nits, 2);
4773: *nits = ts->snes_its;
4774: PetscFunctionReturn(PETSC_SUCCESS);
4775: }
4777: /*@
4778: TSGetKSPIterations - Gets the total number of linear iterations
4779: used by the time integrator.
4781: Not Collective
4783: Input Parameter:
4784: . ts - `TS` context
4786: Output Parameter:
4787: . lits - number of linear iterations
4789: Level: intermediate
4791: Note:
4792: This counter is reset to zero for each successive call to `TSSolve()`.
4794: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()`
4795: @*/
4796: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4797: {
4798: PetscFunctionBegin;
4800: PetscAssertPointer(lits, 2);
4801: *lits = ts->ksp_its;
4802: PetscFunctionReturn(PETSC_SUCCESS);
4803: }
4805: /*@
4806: TSGetStepRejections - Gets the total number of rejected steps.
4808: Not Collective
4810: Input Parameter:
4811: . ts - `TS` context
4813: Output Parameter:
4814: . rejects - number of steps rejected
4816: Level: intermediate
4818: Note:
4819: This counter is reset to zero for each successive call to `TSSolve()`.
4821: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4822: @*/
4823: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4824: {
4825: PetscFunctionBegin;
4827: PetscAssertPointer(rejects, 2);
4828: *rejects = ts->reject;
4829: PetscFunctionReturn(PETSC_SUCCESS);
4830: }
4832: /*@
4833: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
4835: Not Collective
4837: Input Parameter:
4838: . ts - `TS` context
4840: Output Parameter:
4841: . fails - number of failed nonlinear solves
4843: Level: intermediate
4845: Note:
4846: This counter is reset to zero for each successive call to `TSSolve()`.
4848: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4849: @*/
4850: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4851: {
4852: PetscFunctionBegin;
4854: PetscAssertPointer(fails, 2);
4855: *fails = ts->num_snes_failures;
4856: PetscFunctionReturn(PETSC_SUCCESS);
4857: }
4859: /*@
4860: TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails
4862: Not Collective
4864: Input Parameters:
4865: + ts - `TS` context
4866: - rejects - maximum number of rejected steps, pass -1 for unlimited
4868: Options Database Key:
4869: . -ts_max_reject - Maximum number of step rejections before a step fails
4871: Level: intermediate
4873: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4874: @*/
4875: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4876: {
4877: PetscFunctionBegin;
4879: ts->max_reject = rejects;
4880: PetscFunctionReturn(PETSC_SUCCESS);
4881: }
4883: /*@
4884: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves
4886: Not Collective
4888: Input Parameters:
4889: + ts - `TS` context
4890: - fails - maximum number of failed nonlinear solves, pass -1 for unlimited
4892: Options Database Key:
4893: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
4895: Level: intermediate
4897: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4898: @*/
4899: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4900: {
4901: PetscFunctionBegin;
4903: ts->max_snes_failures = fails;
4904: PetscFunctionReturn(PETSC_SUCCESS);
4905: }
4907: /*@
4908: TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`
4910: Not Collective
4912: Input Parameters:
4913: + ts - `TS` context
4914: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
4916: Options Database Key:
4917: . -ts_error_if_step_fails - Error if no step succeeds
4919: Level: intermediate
4921: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
4922: @*/
4923: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
4924: {
4925: PetscFunctionBegin;
4927: ts->errorifstepfailed = err;
4928: PetscFunctionReturn(PETSC_SUCCESS);
4929: }
4931: /*@
4932: TSGetAdapt - Get the adaptive controller context for the current method
4934: Collective if controller has not yet been created
4936: Input Parameter:
4937: . ts - time stepping context
4939: Output Parameter:
4940: . adapt - adaptive controller
4942: Level: intermediate
4944: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
4945: @*/
4946: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
4947: {
4948: PetscFunctionBegin;
4950: PetscAssertPointer(adapt, 2);
4951: if (!ts->adapt) {
4952: PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
4953: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
4954: }
4955: *adapt = ts->adapt;
4956: PetscFunctionReturn(PETSC_SUCCESS);
4957: }
4959: /*@
4960: TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller
4962: Logically Collective
4964: Input Parameters:
4965: + ts - time integration context
4966: . atol - scalar absolute tolerances, `PETSC_DECIDE` to leave current value
4967: . vatol - vector of absolute tolerances or `NULL`, used in preference to atol if present
4968: . rtol - scalar relative tolerances, `PETSC_DECIDE` to leave current value
4969: - vrtol - vector of relative tolerances or `NULL`, used in preference to atol if present
4971: Options Database Keys:
4972: + -ts_rtol <rtol> - relative tolerance for local truncation error
4973: - -ts_atol <atol> - Absolute tolerance for local truncation error
4975: Level: beginner
4977: Notes:
4978: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
4979: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
4980: computed only for the differential or the algebraic part then this can be done using the vector of
4981: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
4982: differential part and infinity for the algebraic part, the LTE calculation will include only the
4983: differential variables.
4985: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
4986: @*/
4987: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
4988: {
4989: PetscFunctionBegin;
4990: if (atol != (PetscReal)PETSC_DECIDE && atol != (PetscReal)PETSC_DEFAULT) ts->atol = atol;
4991: if (vatol) {
4992: PetscCall(PetscObjectReference((PetscObject)vatol));
4993: PetscCall(VecDestroy(&ts->vatol));
4994: ts->vatol = vatol;
4995: }
4996: if (rtol != (PetscReal)PETSC_DECIDE && rtol != (PetscReal)PETSC_DEFAULT) ts->rtol = rtol;
4997: if (vrtol) {
4998: PetscCall(PetscObjectReference((PetscObject)vrtol));
4999: PetscCall(VecDestroy(&ts->vrtol));
5000: ts->vrtol = vrtol;
5001: }
5002: PetscFunctionReturn(PETSC_SUCCESS);
5003: }
5005: /*@
5006: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
5008: Logically Collective
5010: Input Parameter:
5011: . ts - time integration context
5013: Output Parameters:
5014: + atol - scalar absolute tolerances, `NULL` to ignore
5015: . vatol - vector of absolute tolerances, `NULL` to ignore
5016: . rtol - scalar relative tolerances, `NULL` to ignore
5017: - vrtol - vector of relative tolerances, `NULL` to ignore
5019: Level: beginner
5021: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
5022: @*/
5023: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
5024: {
5025: PetscFunctionBegin;
5026: if (atol) *atol = ts->atol;
5027: if (vatol) *vatol = ts->vatol;
5028: if (rtol) *rtol = ts->rtol;
5029: if (vrtol) *vrtol = ts->vrtol;
5030: PetscFunctionReturn(PETSC_SUCCESS);
5031: }
5033: /*@
5034: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
5036: Collective
5038: Input Parameters:
5039: + ts - time stepping context
5040: . U - state vector, usually ts->vec_sol
5041: . Y - state vector to be compared to U
5042: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5044: Output Parameters:
5045: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5046: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5047: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5049: Options Database Key:
5050: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5052: Level: developer
5054: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5055: @*/
5056: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5057: {
5058: PetscInt norma_loc, norm_loc, normr_loc;
5060: PetscFunctionBegin;
5062: 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));
5063: if (wnormtype == NORM_2) {
5064: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5065: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5066: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5067: }
5068: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5069: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5070: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5071: PetscFunctionReturn(PETSC_SUCCESS);
5072: }
5074: /*@
5075: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5077: Collective
5079: Input Parameters:
5080: + ts - time stepping context
5081: . E - error vector
5082: . U - state vector, usually ts->vec_sol
5083: . Y - state vector, previous time step
5084: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5086: Output Parameters:
5087: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5088: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5089: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5091: Options Database Key:
5092: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5094: Level: developer
5096: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5097: @*/
5098: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5099: {
5100: PetscInt norma_loc, norm_loc, normr_loc;
5102: PetscFunctionBegin;
5104: 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));
5105: if (wnormtype == NORM_2) {
5106: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5107: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5108: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5109: }
5110: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5111: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5112: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5113: PetscFunctionReturn(PETSC_SUCCESS);
5114: }
5116: /*@
5117: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5119: Logically Collective
5121: Input Parameters:
5122: + ts - time stepping context
5123: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5125: Note:
5126: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5128: Level: intermediate
5130: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5131: @*/
5132: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5133: {
5134: PetscFunctionBegin;
5136: ts->cfltime_local = cfltime;
5137: ts->cfltime = -1.;
5138: PetscFunctionReturn(PETSC_SUCCESS);
5139: }
5141: /*@
5142: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5144: Collective
5146: Input Parameter:
5147: . ts - time stepping context
5149: Output Parameter:
5150: . cfltime - maximum stable time step for forward Euler
5152: Level: advanced
5154: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5155: @*/
5156: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5157: {
5158: PetscFunctionBegin;
5159: if (ts->cfltime < 0) PetscCall(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5160: *cfltime = ts->cfltime;
5161: PetscFunctionReturn(PETSC_SUCCESS);
5162: }
5164: /*@
5165: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5167: Input Parameters:
5168: + ts - the `TS` context.
5169: . xl - lower bound.
5170: - xu - upper bound.
5172: Level: advanced
5174: Note:
5175: If this routine is not called then the lower and upper bounds are set to
5176: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5178: .seealso: [](ch_ts), `TS`
5179: @*/
5180: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5181: {
5182: SNES snes;
5184: PetscFunctionBegin;
5185: PetscCall(TSGetSNES(ts, &snes));
5186: PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5187: PetscFunctionReturn(PETSC_SUCCESS);
5188: }
5190: /*@
5191: TSComputeLinearStability - computes the linear stability function at a point
5193: Collective
5195: Input Parameters:
5196: + ts - the `TS` context
5197: . xr - real part of input argument
5198: - xi - imaginary part of input argument
5200: Output Parameters:
5201: + yr - real part of function value
5202: - yi - imaginary part of function value
5204: Level: developer
5206: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5207: @*/
5208: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5209: {
5210: PetscFunctionBegin;
5212: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5213: PetscFunctionReturn(PETSC_SUCCESS);
5214: }
5216: /*@
5217: TSRestartStep - Flags the solver to restart the next step
5219: Collective
5221: Input Parameter:
5222: . ts - the `TS` context obtained from `TSCreate()`
5224: Level: advanced
5226: Notes:
5227: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5228: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5229: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5230: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5231: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5232: discontinuous source terms).
5234: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5235: @*/
5236: PetscErrorCode TSRestartStep(TS ts)
5237: {
5238: PetscFunctionBegin;
5240: ts->steprestart = PETSC_TRUE;
5241: PetscFunctionReturn(PETSC_SUCCESS);
5242: }
5244: /*@
5245: TSRollBack - Rolls back one time step
5247: Collective
5249: Input Parameter:
5250: . ts - the `TS` context obtained from `TSCreate()`
5252: Level: advanced
5254: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5255: @*/
5256: PetscErrorCode TSRollBack(TS ts)
5257: {
5258: PetscFunctionBegin;
5260: PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5261: PetscTryTypeMethod(ts, rollback);
5262: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
5263: ts->time_step = ts->ptime - ts->ptime_prev;
5264: ts->ptime = ts->ptime_prev;
5265: ts->ptime_prev = ts->ptime_prev_rollback;
5266: ts->steps--;
5267: ts->steprollback = PETSC_TRUE;
5268: PetscFunctionReturn(PETSC_SUCCESS);
5269: }
5271: /*@
5272: TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step
5274: Not collective
5276: Input Parameter:
5277: . ts - the `TS` context obtained from `TSCreate()`
5279: Output Parameter:
5280: . flg - the rollback flag
5282: Level: advanced
5284: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
5285: @*/
5286: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
5287: {
5288: PetscFunctionBegin;
5290: PetscAssertPointer(flg, 2);
5291: *flg = ts->steprollback;
5292: PetscFunctionReturn(PETSC_SUCCESS);
5293: }
5295: /*@
5296: TSGetStages - Get the number of stages and stage values
5298: Input Parameter:
5299: . ts - the `TS` context obtained from `TSCreate()`
5301: Output Parameters:
5302: + ns - the number of stages
5303: - Y - the current stage vectors
5305: Level: advanced
5307: Note:
5308: Both `ns` and `Y` can be `NULL`.
5310: .seealso: [](ch_ts), `TS`, `TSCreate()`
5311: @*/
5312: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5313: {
5314: PetscFunctionBegin;
5316: if (ns) PetscAssertPointer(ns, 2);
5317: if (Y) PetscAssertPointer(Y, 3);
5318: if (!ts->ops->getstages) {
5319: if (ns) *ns = 0;
5320: if (Y) *Y = NULL;
5321: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5322: PetscFunctionReturn(PETSC_SUCCESS);
5323: }
5325: /*@C
5326: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5328: Collective
5330: Input Parameters:
5331: + ts - the `TS` context
5332: . t - current timestep
5333: . U - state vector
5334: . Udot - time derivative of state vector
5335: . shift - shift to apply, see note below
5336: - ctx - an optional user context
5338: Output Parameters:
5339: + J - Jacobian matrix (not altered in this routine)
5340: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
5342: Level: intermediate
5344: Notes:
5345: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5347: dF/dU + shift*dF/dUdot
5349: Most users should not need to explicitly call this routine, as it
5350: is used internally within the nonlinear solvers.
5352: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5353: routine, then it will try to get the coloring from the matrix. This requires that the
5354: matrix have nonzero entries precomputed.
5356: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5357: @*/
5358: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5359: {
5360: SNES snes;
5361: MatFDColoring color;
5362: PetscBool hascolor, matcolor = PETSC_FALSE;
5364: PetscFunctionBegin;
5365: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5366: PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5367: if (!color) {
5368: DM dm;
5369: ISColoring iscoloring;
5371: PetscCall(TSGetDM(ts, &dm));
5372: PetscCall(DMHasColoring(dm, &hascolor));
5373: if (hascolor && !matcolor) {
5374: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5375: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5376: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5377: PetscCall(MatFDColoringSetFromOptions(color));
5378: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5379: PetscCall(ISColoringDestroy(&iscoloring));
5380: } else {
5381: MatColoring mc;
5383: PetscCall(MatColoringCreate(B, &mc));
5384: PetscCall(MatColoringSetDistance(mc, 2));
5385: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5386: PetscCall(MatColoringSetFromOptions(mc));
5387: PetscCall(MatColoringApply(mc, &iscoloring));
5388: PetscCall(MatColoringDestroy(&mc));
5389: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5390: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5391: PetscCall(MatFDColoringSetFromOptions(color));
5392: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5393: PetscCall(ISColoringDestroy(&iscoloring));
5394: }
5395: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5396: PetscCall(PetscObjectDereference((PetscObject)color));
5397: }
5398: PetscCall(TSGetSNES(ts, &snes));
5399: PetscCall(MatFDColoringApply(B, color, U, snes));
5400: if (J != B) {
5401: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5402: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5403: }
5404: PetscFunctionReturn(PETSC_SUCCESS);
5405: }
5407: /*@C
5408: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5410: Input Parameters:
5411: + ts - the `TS` context
5412: - func - function called within `TSFunctionDomainError()`
5414: Calling sequence of `func`:
5415: + ts - the `TS` context
5416: . time - the current time (of the stage)
5417: . state - the state to check if it is valid
5418: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable
5420: Level: intermediate
5422: Notes:
5423: If an implicit ODE solver is being used then, in addition to providing this routine, the
5424: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5425: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5426: Use `TSGetSNES()` to obtain the `SNES` object
5428: Developer Notes:
5429: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5430: since one takes a function pointer and the other does not.
5432: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5433: @*/
5434: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5435: {
5436: PetscFunctionBegin;
5438: ts->functiondomainerror = func;
5439: PetscFunctionReturn(PETSC_SUCCESS);
5440: }
5442: /*@
5443: TSFunctionDomainError - Checks if the current state is valid
5445: Input Parameters:
5446: + ts - the `TS` context
5447: . stagetime - time of the simulation
5448: - Y - state vector to check.
5450: Output Parameter:
5451: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5453: Level: developer
5455: Note:
5456: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5457: to check if the current state is valid.
5459: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5460: @*/
5461: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5462: {
5463: PetscFunctionBegin;
5465: *accept = PETSC_TRUE;
5466: if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5467: PetscFunctionReturn(PETSC_SUCCESS);
5468: }
5470: /*@C
5471: TSClone - This function clones a time step `TS` object.
5473: Collective
5475: Input Parameter:
5476: . tsin - The input `TS`
5478: Output Parameter:
5479: . tsout - The output `TS` (cloned)
5481: Level: developer
5483: Notes:
5484: 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.
5485: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5487: 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
5488: .vb
5489: SNES snes_dup = NULL;
5490: TSGetSNES(ts,&snes_dup);
5491: TSSetSNES(ts,snes_dup);
5492: .ve
5494: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5495: @*/
5496: PetscErrorCode TSClone(TS tsin, TS *tsout)
5497: {
5498: TS t;
5499: SNES snes_start;
5500: DM dm;
5501: TSType type;
5503: PetscFunctionBegin;
5504: PetscAssertPointer(tsin, 1);
5505: *tsout = NULL;
5507: PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));
5509: /* General TS description */
5510: t->numbermonitors = 0;
5511: t->monitorFrequency = 1;
5512: t->setupcalled = 0;
5513: t->ksp_its = 0;
5514: t->snes_its = 0;
5515: t->nwork = 0;
5516: t->rhsjacobian.time = PETSC_MIN_REAL;
5517: t->rhsjacobian.scale = 1.;
5518: t->ijacobian.shift = 1.;
5520: PetscCall(TSGetSNES(tsin, &snes_start));
5521: PetscCall(TSSetSNES(t, snes_start));
5523: PetscCall(TSGetDM(tsin, &dm));
5524: PetscCall(TSSetDM(t, dm));
5526: t->adapt = tsin->adapt;
5527: PetscCall(PetscObjectReference((PetscObject)t->adapt));
5529: t->trajectory = tsin->trajectory;
5530: PetscCall(PetscObjectReference((PetscObject)t->trajectory));
5532: t->event = tsin->event;
5533: if (t->event) t->event->refct++;
5535: t->problem_type = tsin->problem_type;
5536: t->ptime = tsin->ptime;
5537: t->ptime_prev = tsin->ptime_prev;
5538: t->time_step = tsin->time_step;
5539: t->max_time = tsin->max_time;
5540: t->steps = tsin->steps;
5541: t->max_steps = tsin->max_steps;
5542: t->equation_type = tsin->equation_type;
5543: t->atol = tsin->atol;
5544: t->rtol = tsin->rtol;
5545: t->max_snes_failures = tsin->max_snes_failures;
5546: t->max_reject = tsin->max_reject;
5547: t->errorifstepfailed = tsin->errorifstepfailed;
5549: PetscCall(TSGetType(tsin, &type));
5550: PetscCall(TSSetType(t, type));
5552: t->vec_sol = NULL;
5554: t->cfltime = tsin->cfltime;
5555: t->cfltime_local = tsin->cfltime_local;
5556: t->exact_final_time = tsin->exact_final_time;
5558: t->ops[0] = tsin->ops[0];
5560: if (((PetscObject)tsin)->fortran_func_pointers) {
5561: PetscInt i;
5562: PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers));
5563: for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5564: }
5565: *tsout = t;
5566: PetscFunctionReturn(PETSC_SUCCESS);
5567: }
5569: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5570: {
5571: TS ts = (TS)ctx;
5573: PetscFunctionBegin;
5574: PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5575: PetscFunctionReturn(PETSC_SUCCESS);
5576: }
5578: /*@
5579: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5581: Logically Collective
5583: Input Parameter:
5584: . ts - the time stepping routine
5586: Output Parameter:
5587: . flg - `PETSC_TRUE` if the multiply is likely correct
5589: Options Database Key:
5590: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5592: Level: advanced
5594: Note:
5595: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5597: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5598: @*/
5599: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5600: {
5601: Mat J, B;
5602: TSRHSJacobianFn *func;
5603: void *ctx;
5605: PetscFunctionBegin;
5606: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5607: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5608: PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5609: PetscFunctionReturn(PETSC_SUCCESS);
5610: }
5612: /*@C
5613: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5615: Logically Collective
5617: Input Parameter:
5618: . ts - the time stepping routine
5620: Output Parameter:
5621: . flg - `PETSC_TRUE` if the multiply is likely correct
5623: Options Database Key:
5624: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5626: Level: advanced
5628: Notes:
5629: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5631: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5632: @*/
5633: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5634: {
5635: Mat J, B;
5636: void *ctx;
5637: TSRHSJacobianFn *func;
5639: PetscFunctionBegin;
5640: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5641: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5642: PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5643: PetscFunctionReturn(PETSC_SUCCESS);
5644: }
5646: /*@
5647: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
5649: Logically Collective
5651: Input Parameters:
5652: + ts - timestepping context
5653: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5655: Options Database Key:
5656: . -ts_use_splitrhsfunction - <true,false>
5658: Level: intermediate
5660: Note:
5661: This is only for multirate methods
5663: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5664: @*/
5665: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5666: {
5667: PetscFunctionBegin;
5669: ts->use_splitrhsfunction = use_splitrhsfunction;
5670: PetscFunctionReturn(PETSC_SUCCESS);
5671: }
5673: /*@
5674: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
5676: Not Collective
5678: Input Parameter:
5679: . ts - timestepping context
5681: Output Parameter:
5682: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5684: Level: intermediate
5686: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5687: @*/
5688: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5689: {
5690: PetscFunctionBegin;
5692: *use_splitrhsfunction = ts->use_splitrhsfunction;
5693: PetscFunctionReturn(PETSC_SUCCESS);
5694: }
5696: /*@
5697: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
5699: Logically Collective
5701: Input Parameters:
5702: + ts - the time-stepper
5703: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
5705: Level: intermediate
5707: Note:
5708: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
5710: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
5711: @*/
5712: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5713: {
5714: PetscFunctionBegin;
5716: ts->axpy_pattern = str;
5717: PetscFunctionReturn(PETSC_SUCCESS);
5718: }
5720: /*@
5721: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
5723: Collective
5725: Input Parameters:
5726: + ts - the time-stepper
5727: . n - number of the time points (>=2)
5728: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5730: Options Database Key:
5731: . -ts_time_span <t0,...tf> - Sets the time span
5733: Level: intermediate
5735: Notes:
5736: The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
5737: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
5738: The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may
5739: pressure the memory system when using a large number of span points.
5741: .seealso: [](ch_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()`
5742: @*/
5743: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
5744: {
5745: PetscFunctionBegin;
5747: PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
5748: if (ts->tspan && n != ts->tspan->num_span_times) {
5749: PetscCall(PetscFree(ts->tspan->span_times));
5750: PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
5751: PetscCall(PetscMalloc1(n, &ts->tspan->span_times));
5752: }
5753: if (!ts->tspan) {
5754: TSTimeSpan tspan;
5755: PetscCall(PetscNew(&tspan));
5756: PetscCall(PetscMalloc1(n, &tspan->span_times));
5757: tspan->reltol = 1e-6;
5758: tspan->abstol = 10 * PETSC_MACHINE_EPSILON;
5759: tspan->worktol = 0;
5760: ts->tspan = tspan;
5761: }
5762: ts->tspan->num_span_times = n;
5763: PetscCall(PetscArraycpy(ts->tspan->span_times, span_times, n));
5764: PetscCall(TSSetTime(ts, ts->tspan->span_times[0]));
5765: PetscCall(TSSetMaxTime(ts, ts->tspan->span_times[n - 1]));
5766: PetscFunctionReturn(PETSC_SUCCESS);
5767: }
5769: /*@C
5770: TSGetTimeSpan - gets the time span set with `TSSetTimeSpan()`
5772: Not Collective
5774: Input Parameter:
5775: . ts - the time-stepper
5777: Output Parameters:
5778: + n - number of the time points (>=2)
5779: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5781: Level: beginner
5783: Note:
5784: The values obtained are valid until the `TS` object is destroyed.
5786: Both `n` and `span_times` can be `NULL`.
5788: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()`
5789: @*/
5790: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal **span_times)
5791: {
5792: PetscFunctionBegin;
5794: if (n) PetscAssertPointer(n, 2);
5795: if (span_times) PetscAssertPointer(span_times, 3);
5796: if (!ts->tspan) {
5797: if (n) *n = 0;
5798: if (span_times) *span_times = NULL;
5799: } else {
5800: if (n) *n = ts->tspan->num_span_times;
5801: if (span_times) *span_times = ts->tspan->span_times;
5802: }
5803: PetscFunctionReturn(PETSC_SUCCESS);
5804: }
5806: /*@
5807: TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.
5809: Input Parameter:
5810: . ts - the `TS` context obtained from `TSCreate()`
5812: Output Parameters:
5813: + nsol - the number of solutions
5814: - Sols - the solution vectors
5816: Level: intermediate
5818: Notes:
5819: Both `nsol` and `Sols` can be `NULL`.
5821: 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()`.
5822: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span.
5824: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`
5825: @*/
5826: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
5827: {
5828: PetscFunctionBegin;
5830: if (nsol) PetscAssertPointer(nsol, 2);
5831: if (Sols) PetscAssertPointer(Sols, 3);
5832: if (!ts->tspan) {
5833: if (nsol) *nsol = 0;
5834: if (Sols) *Sols = NULL;
5835: } else {
5836: if (nsol) *nsol = ts->tspan->spanctr;
5837: if (Sols) *Sols = ts->tspan->vecs_sol;
5838: }
5839: PetscFunctionReturn(PETSC_SUCCESS);
5840: }
5842: /*@C
5843: TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
5845: Collective
5847: Input Parameters:
5848: + ts - the `TS` context
5849: . J - Jacobian matrix (not altered in this routine)
5850: - B - newly computed Jacobian matrix to use with preconditioner
5852: Level: intermediate
5854: Notes:
5855: This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
5856: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
5857: and multiple fields are involved.
5859: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
5860: structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
5861: usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
5862: `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.
5864: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5865: @*/
5866: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
5867: {
5868: MatColoring mc = NULL;
5869: ISColoring iscoloring = NULL;
5870: MatFDColoring matfdcoloring = NULL;
5872: PetscFunctionBegin;
5873: /* Generate new coloring after eliminating zeros in the matrix */
5874: PetscCall(MatEliminateZeros(B, PETSC_TRUE));
5875: PetscCall(MatColoringCreate(B, &mc));
5876: PetscCall(MatColoringSetDistance(mc, 2));
5877: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5878: PetscCall(MatColoringSetFromOptions(mc));
5879: PetscCall(MatColoringApply(mc, &iscoloring));
5880: PetscCall(MatColoringDestroy(&mc));
5881: /* Replace the old coloring with the new one */
5882: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
5883: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5884: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
5885: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
5886: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
5887: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
5888: PetscCall(ISColoringDestroy(&iscoloring));
5889: PetscFunctionReturn(PETSC_SUCCESS);
5890: }