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: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
714: {
715: Mat A, B;
716: TSIJacobianFn *ijacobian;
718: PetscFunctionBegin;
719: if (Arhs) *Arhs = NULL;
720: if (Brhs) *Brhs = NULL;
721: PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL));
722: if (Arhs) {
723: if (!ts->Arhs) {
724: if (ijacobian) {
725: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs));
726: PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN));
727: } else {
728: ts->Arhs = A;
729: PetscCall(PetscObjectReference((PetscObject)A));
730: }
731: } else {
732: PetscBool flg;
733: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
734: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
735: if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
736: PetscCall(PetscObjectDereference((PetscObject)ts->Arhs));
737: ts->Arhs = A;
738: PetscCall(PetscObjectReference((PetscObject)A));
739: }
740: }
741: *Arhs = ts->Arhs;
742: }
743: if (Brhs) {
744: if (!ts->Brhs) {
745: if (A != B) {
746: if (ijacobian) {
747: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs));
748: } else {
749: ts->Brhs = B;
750: PetscCall(PetscObjectReference((PetscObject)B));
751: }
752: } else {
753: PetscCall(PetscObjectReference((PetscObject)ts->Arhs));
754: ts->Brhs = ts->Arhs;
755: }
756: }
757: *Brhs = ts->Brhs;
758: }
759: PetscFunctionReturn(PETSC_SUCCESS);
760: }
762: /*@
763: TSComputeIFunction - Evaluates the DAE residual written in the implicit form F(t,U,Udot)=0
765: Collective
767: Input Parameters:
768: + ts - the `TS` context
769: . t - current time
770: . U - state vector
771: . Udot - time derivative of state vector
772: - imex - flag indicates if the method is `TSIMEX` so that the RHSFunction should be kept separate
774: Output Parameter:
775: . Y - right-hand side
777: Level: developer
779: Note:
780: Most users should not need to explicitly call this routine, as it
781: is used internally within the nonlinear solvers.
783: If the user did not write their equations in implicit form, this
784: function recasts them in implicit form.
786: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
787: @*/
788: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
789: {
790: TSIFunctionFn *ifunction;
791: TSRHSFunctionFn *rhsfunction;
792: void *ctx;
793: DM dm;
795: PetscFunctionBegin;
801: PetscCall(TSGetDM(ts, &dm));
802: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
803: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
805: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
807: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, Udot, Y));
808: if (ifunction) {
809: PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
810: ts->ifuncs++;
811: }
812: if (imex) {
813: if (!ifunction) PetscCall(VecCopy(Udot, Y));
814: } else if (rhsfunction) {
815: if (ifunction) {
816: Vec Frhs;
818: PetscCall(DMGetGlobalVector(dm, &Frhs));
819: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
820: PetscCall(VecAXPY(Y, -1, Frhs));
821: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
822: } else {
823: PetscCall(TSComputeRHSFunction(ts, t, U, Y));
824: PetscCall(VecAYPX(Y, -1, Udot));
825: }
826: }
827: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, Udot, Y));
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }
831: /*
832: TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call `TSComputeRHSJacobian()` on it.
834: Note:
835: This routine is needed when one switches from `TSComputeIJacobian()` to `TSComputeRHSJacobian()` because the Jacobian matrix may be shifted or scaled in `TSComputeIJacobian()`.
837: */
838: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
839: {
840: PetscFunctionBegin;
842: PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
843: PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");
845: if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
846: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
847: if (B && B == ts->Brhs && A != B) {
848: if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
849: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
850: }
851: ts->rhsjacobian.shift = 0;
852: ts->rhsjacobian.scale = 1.;
853: PetscFunctionReturn(PETSC_SUCCESS);
854: }
856: /*@
857: TSComputeIJacobian - Evaluates the Jacobian of the DAE
859: Collective
861: Input Parameters:
862: + ts - the `TS` context
863: . t - current timestep
864: . U - state vector
865: . Udot - time derivative of state vector
866: . shift - shift to apply, see note below
867: - imex - flag indicates if the method is `TSIMEX` so that the RHSJacobian should be kept separate
869: Output Parameters:
870: + A - Jacobian matrix
871: - B - matrix from which the preconditioner is constructed; often the same as `A`
873: Level: developer
875: Notes:
876: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
877: .vb
878: dF/dU + shift*dF/dUdot
879: .ve
880: Most users should not need to explicitly call this routine, as it
881: is used internally within the nonlinear solvers.
883: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
884: @*/
885: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
886: {
887: TSIJacobianFn *ijacobian;
888: TSRHSJacobianFn *rhsjacobian;
889: DM dm;
890: void *ctx;
892: PetscFunctionBegin;
899: PetscCall(TSGetDM(ts, &dm));
900: PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
901: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
903: PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
905: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
906: if (ijacobian) {
907: PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
908: ts->ijacs++;
909: }
910: if (imex) {
911: if (!ijacobian) { /* system was written as Udot = G(t,U) */
912: PetscBool assembled;
913: if (rhsjacobian) {
914: Mat Arhs = NULL;
915: PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
916: if (A == Arhs) {
917: PetscCheck(rhsjacobian != TSComputeRHSJacobianConstant, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Unsupported operation! cannot use TSComputeRHSJacobianConstant"); /* there is no way to reconstruct shift*M-J since J cannot be reevaluated */
918: ts->rhsjacobian.time = PETSC_MIN_REAL;
919: }
920: }
921: PetscCall(MatZeroEntries(A));
922: PetscCall(MatAssembled(A, &assembled));
923: if (!assembled) {
924: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
925: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
926: }
927: PetscCall(MatShift(A, shift));
928: if (A != B) {
929: PetscCall(MatZeroEntries(B));
930: PetscCall(MatAssembled(B, &assembled));
931: if (!assembled) {
932: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
933: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
934: }
935: PetscCall(MatShift(B, shift));
936: }
937: }
938: } else {
939: Mat Arhs = NULL, Brhs = NULL;
941: /* RHSJacobian needs to be converted to part of IJacobian if exists */
942: if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
943: if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
944: PetscObjectState Ustate;
945: PetscObjectId Uid;
946: TSRHSFunctionFn *rhsfunction;
948: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
949: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
950: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
951: if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
952: ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */
953: PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
954: if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
955: } else {
956: PetscBool flg;
958: if (ts->rhsjacobian.reuse) { /* Undo the damage */
959: /* MatScale has a short path for this case.
960: However, this code path is taken the first time TSComputeRHSJacobian is called
961: and the matrices have not been assembled yet */
962: PetscCall(TSRecoverRHSJacobian(ts, A, B));
963: }
964: PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
965: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
966: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
967: if (!flg) {
968: PetscCall(MatScale(A, -1));
969: PetscCall(MatShift(A, shift));
970: }
971: if (A != B) {
972: PetscCall(MatScale(B, -1));
973: PetscCall(MatShift(B, shift));
974: }
975: }
976: ts->rhsjacobian.scale = -1;
977: ts->rhsjacobian.shift = shift;
978: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
979: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
980: PetscCall(MatZeroEntries(A));
981: PetscCall(MatShift(A, shift));
982: if (A != B) {
983: PetscCall(MatZeroEntries(B));
984: PetscCall(MatShift(B, shift));
985: }
986: }
987: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
988: PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
989: if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
990: }
991: }
992: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
993: PetscFunctionReturn(PETSC_SUCCESS);
994: }
996: /*@C
997: TSSetRHSFunction - Sets the routine for evaluating the function,
998: where U_t = G(t,u).
1000: Logically Collective
1002: Input Parameters:
1003: + ts - the `TS` context obtained from `TSCreate()`
1004: . r - vector to put the computed right-hand side (or `NULL` to have it created)
1005: . f - routine for evaluating the right-hand-side function
1006: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
1008: Level: beginner
1010: Note:
1011: You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.
1013: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1014: @*/
1015: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, void *ctx)
1016: {
1017: SNES snes;
1018: Vec ralloc = NULL;
1019: DM dm;
1021: PetscFunctionBegin;
1025: PetscCall(TSGetDM(ts, &dm));
1026: PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1027: PetscCall(TSGetSNES(ts, &snes));
1028: if (!r && !ts->dm && ts->vec_sol) {
1029: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1030: r = ralloc;
1031: }
1032: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1033: PetscCall(VecDestroy(&ralloc));
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: /*@C
1038: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1040: Logically Collective
1042: Input Parameters:
1043: + ts - the `TS` context obtained from `TSCreate()`
1044: . f - routine for evaluating the solution
1045: - ctx - [optional] user-defined context for private data for the
1046: function evaluation routine (may be `NULL`)
1048: Options Database Keys:
1049: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1050: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
1052: Level: intermediate
1054: Notes:
1055: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1056: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1057: create closed-form solutions with non-physical forcing terms.
1059: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1061: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1062: @*/
1063: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, void *ctx)
1064: {
1065: DM dm;
1067: PetscFunctionBegin;
1069: PetscCall(TSGetDM(ts, &dm));
1070: PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: /*@C
1075: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1077: Logically Collective
1079: Input Parameters:
1080: + ts - the `TS` context obtained from `TSCreate()`
1081: . func - routine for evaluating the forcing function
1082: - ctx - [optional] user-defined context for private data for the function evaluation routine
1083: (may be `NULL`)
1085: Level: intermediate
1087: Notes:
1088: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1089: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1090: definition of the problem you are solving and hence possibly introducing bugs.
1092: This replaces the ODE F(u,u_t,t) = 0 the `TS` is solving with F(u,u_t,t) - func(t) = 0
1094: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1095: parameters can be passed in the ctx variable.
1097: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1099: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
1100: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1101: @*/
1102: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, void *ctx)
1103: {
1104: DM dm;
1106: PetscFunctionBegin;
1108: PetscCall(TSGetDM(ts, &dm));
1109: PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1110: PetscFunctionReturn(PETSC_SUCCESS);
1111: }
1113: /*@C
1114: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1115: where U_t = G(U,t), as well as the location to store the matrix.
1117: Logically Collective
1119: Input Parameters:
1120: + ts - the `TS` context obtained from `TSCreate()`
1121: . Amat - (approximate) location to store Jacobian matrix entries computed by `f`
1122: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1123: . f - the Jacobian evaluation routine
1124: - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1126: Level: beginner
1128: Notes:
1129: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1131: The `TS` solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f()`
1132: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1134: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
1135: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
1136: @*/
1137: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, void *ctx)
1138: {
1139: SNES snes;
1140: DM dm;
1141: TSIJacobianFn *ijacobian;
1143: PetscFunctionBegin;
1147: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1148: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1150: PetscCall(TSGetDM(ts, &dm));
1151: PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1152: PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1153: PetscCall(TSGetSNES(ts, &snes));
1154: if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1155: if (Amat) {
1156: PetscCall(PetscObjectReference((PetscObject)Amat));
1157: PetscCall(MatDestroy(&ts->Arhs));
1158: ts->Arhs = Amat;
1159: }
1160: if (Pmat) {
1161: PetscCall(PetscObjectReference((PetscObject)Pmat));
1162: PetscCall(MatDestroy(&ts->Brhs));
1163: ts->Brhs = Pmat;
1164: }
1165: PetscFunctionReturn(PETSC_SUCCESS);
1166: }
1168: /*@C
1169: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1171: Logically Collective
1173: Input Parameters:
1174: + ts - the `TS` context obtained from `TSCreate()`
1175: . r - vector to hold the residual (or `NULL` to have it created internally)
1176: . f - the function evaluation routine
1177: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1179: Level: beginner
1181: Note:
1182: The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function.
1184: .seealso: [](ch_ts), `TS`, `TSIFunctionFn`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`,
1185: `TSSetIJacobian()`
1186: @*/
1187: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunctionFn *f, void *ctx)
1188: {
1189: SNES snes;
1190: Vec ralloc = NULL;
1191: DM dm;
1193: PetscFunctionBegin;
1197: PetscCall(TSGetDM(ts, &dm));
1198: PetscCall(DMTSSetIFunction(dm, f, ctx));
1200: PetscCall(TSGetSNES(ts, &snes));
1201: if (!r && !ts->dm && ts->vec_sol) {
1202: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1203: r = ralloc;
1204: }
1205: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1206: PetscCall(VecDestroy(&ralloc));
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: /*@C
1211: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.
1213: Not Collective
1215: Input Parameter:
1216: . ts - the `TS` context
1218: Output Parameters:
1219: + r - vector to hold residual (or `NULL`)
1220: . func - the function to compute residual (or `NULL`)
1221: - ctx - the function context (or `NULL`)
1223: Level: advanced
1225: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1226: @*/
1227: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunctionFn **func, void **ctx)
1228: {
1229: SNES snes;
1230: DM dm;
1232: PetscFunctionBegin;
1234: PetscCall(TSGetSNES(ts, &snes));
1235: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1236: PetscCall(TSGetDM(ts, &dm));
1237: PetscCall(DMTSGetIFunction(dm, func, ctx));
1238: PetscFunctionReturn(PETSC_SUCCESS);
1239: }
1241: /*@C
1242: TSGetRHSFunction - Returns the vector where the right-hand side is stored and the function/context to compute it.
1244: Not Collective
1246: Input Parameter:
1247: . ts - the `TS` context
1249: Output Parameters:
1250: + r - vector to hold computed right-hand side (or `NULL`)
1251: . func - the function to compute right-hand side (or `NULL`)
1252: - ctx - the function context (or `NULL`)
1254: Level: advanced
1256: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1257: @*/
1258: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunctionFn **func, void **ctx)
1259: {
1260: SNES snes;
1261: DM dm;
1263: PetscFunctionBegin;
1265: PetscCall(TSGetSNES(ts, &snes));
1266: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1267: PetscCall(TSGetDM(ts, &dm));
1268: PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1269: PetscFunctionReturn(PETSC_SUCCESS);
1270: }
1272: /*@C
1273: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1274: provided with `TSSetIFunction()`.
1276: Logically Collective
1278: Input Parameters:
1279: + ts - the `TS` context obtained from `TSCreate()`
1280: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
1281: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1282: . f - the Jacobian evaluation routine
1283: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1285: Level: beginner
1287: Notes:
1288: The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1290: If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
1291: space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.
1293: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1294: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1295: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1296: a and vector W depend on the integration method, step size, and past states. For example with
1297: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1298: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1300: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1302: The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f`
1303: You should not assume the values are the same in the next call to `f` as you set them in the previous call.
1305: In case `TSSetRHSJacobian()` is also used in conjuction with a fully-implicit solver,
1306: multilevel linear solvers, e.g. `PCMG`, will likely not work due to the way `TS` handles rhs matrices.
1308: .seealso: [](ch_ts), `TS`, `TSIJacobianFn`, `TSSetIFunction()`, `TSSetRHSJacobian()`,
1309: `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1310: @*/
1311: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobianFn *f, void *ctx)
1312: {
1313: SNES snes;
1314: DM dm;
1316: PetscFunctionBegin;
1320: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1321: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1323: PetscCall(TSGetDM(ts, &dm));
1324: PetscCall(DMTSSetIJacobian(dm, f, ctx));
1326: PetscCall(TSGetSNES(ts, &snes));
1327: PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1328: PetscFunctionReturn(PETSC_SUCCESS);
1329: }
1331: /*@
1332: TSRHSJacobianSetReuse - restore the RHS Jacobian before calling the user-provided `TSRHSJacobianFn` function again
1334: Logically Collective
1336: Input Parameters:
1337: + ts - `TS` context obtained from `TSCreate()`
1338: - reuse - `PETSC_TRUE` if the RHS Jacobian
1340: Level: intermediate
1342: Notes:
1343: Without this flag, `TS` will change the sign and shift the RHS Jacobian for a
1344: finite-time-step implicit solve, in which case the user function will need to recompute the
1345: entire Jacobian. The `reuse `flag must be set if the evaluation function assumes that the
1346: matrix entries have not been changed by the `TS`.
1348: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1349: @*/
1350: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1351: {
1352: PetscFunctionBegin;
1353: ts->rhsjacobian.reuse = reuse;
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: /*@C
1358: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1360: Logically Collective
1362: Input Parameters:
1363: + ts - the `TS` context obtained from `TSCreate()`
1364: . F - vector to hold the residual (or `NULL` to have it created internally)
1365: . fun - the function evaluation routine
1366: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1368: Level: beginner
1370: .seealso: [](ch_ts), `TS`, `TSI2FunctionFn`, `TSSetI2Jacobian()`, `TSSetIFunction()`,
1371: `TSCreate()`, `TSSetRHSFunction()`
1372: @*/
1373: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2FunctionFn *fun, void *ctx)
1374: {
1375: DM dm;
1377: PetscFunctionBegin;
1380: PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1381: PetscCall(TSGetDM(ts, &dm));
1382: PetscCall(DMTSSetI2Function(dm, fun, ctx));
1383: PetscFunctionReturn(PETSC_SUCCESS);
1384: }
1386: /*@C
1387: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.
1389: Not Collective
1391: Input Parameter:
1392: . ts - the `TS` context
1394: Output Parameters:
1395: + r - vector to hold residual (or `NULL`)
1396: . fun - the function to compute residual (or `NULL`)
1397: - ctx - the function context (or `NULL`)
1399: Level: advanced
1401: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1402: @*/
1403: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2FunctionFn **fun, void **ctx)
1404: {
1405: SNES snes;
1406: DM dm;
1408: PetscFunctionBegin;
1410: PetscCall(TSGetSNES(ts, &snes));
1411: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1412: PetscCall(TSGetDM(ts, &dm));
1413: PetscCall(DMTSGetI2Function(dm, fun, ctx));
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: /*@C
1418: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1419: where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.
1421: Logically Collective
1423: Input Parameters:
1424: + ts - the `TS` context obtained from `TSCreate()`
1425: . J - matrix to hold the Jacobian values
1426: . P - matrix for constructing the preconditioner (may be same as `J`)
1427: . jac - the Jacobian evaluation routine, see `TSI2JacobianFn` for the calling sequence
1428: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1430: Level: beginner
1432: Notes:
1433: The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1435: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1436: the Jacobian of G(U) = F(t,U,W+v*U,W'+a*U) where F(t,U,U_t,U_tt) = 0 is the DAE to be solved.
1437: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1438: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1440: .seealso: [](ch_ts), `TS`, `TSI2JacobianFn`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1441: @*/
1442: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2JacobianFn *jac, void *ctx)
1443: {
1444: DM dm;
1446: PetscFunctionBegin;
1450: PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1451: PetscCall(TSGetDM(ts, &dm));
1452: PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1453: PetscFunctionReturn(PETSC_SUCCESS);
1454: }
1456: /*@C
1457: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1459: Not Collective, but parallel objects are returned if `TS` is parallel
1461: Input Parameter:
1462: . ts - The `TS` context obtained from `TSCreate()`
1464: Output Parameters:
1465: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1466: . P - The matrix from which the preconditioner is constructed, often the same as `J`
1467: . jac - The function to compute the Jacobian matrices
1468: - ctx - User-defined context for Jacobian evaluation routine
1470: Level: advanced
1472: Note:
1473: You can pass in `NULL` for any return argument you do not need.
1475: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1476: @*/
1477: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2JacobianFn **jac, void **ctx)
1478: {
1479: SNES snes;
1480: DM dm;
1482: PetscFunctionBegin;
1483: PetscCall(TSGetSNES(ts, &snes));
1484: PetscCall(SNESSetUpMatrices(snes));
1485: PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1486: PetscCall(TSGetDM(ts, &dm));
1487: PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1488: PetscFunctionReturn(PETSC_SUCCESS);
1489: }
1491: /*@
1492: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1494: Collective
1496: Input Parameters:
1497: + ts - the `TS` context
1498: . t - current time
1499: . U - state vector
1500: . V - time derivative of state vector (U_t)
1501: - A - second time derivative of state vector (U_tt)
1503: Output Parameter:
1504: . F - the residual vector
1506: Level: developer
1508: Note:
1509: Most users should not need to explicitly call this routine, as it
1510: is used internally within the nonlinear solvers.
1512: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1513: @*/
1514: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1515: {
1516: DM dm;
1517: TSI2FunctionFn *I2Function;
1518: void *ctx;
1519: TSRHSFunctionFn *rhsfunction;
1521: PetscFunctionBegin;
1528: PetscCall(TSGetDM(ts, &dm));
1529: PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1530: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
1532: if (!I2Function) {
1533: PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1534: PetscFunctionReturn(PETSC_SUCCESS);
1535: }
1537: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, V, F));
1539: PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));
1541: if (rhsfunction) {
1542: Vec Frhs;
1544: PetscCall(DMGetGlobalVector(dm, &Frhs));
1545: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1546: PetscCall(VecAXPY(F, -1, Frhs));
1547: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
1548: }
1550: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, V, F));
1551: PetscFunctionReturn(PETSC_SUCCESS);
1552: }
1554: /*@
1555: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1557: Collective
1559: Input Parameters:
1560: + ts - the `TS` context
1561: . t - current timestep
1562: . U - state vector
1563: . V - time derivative of state vector
1564: . A - second time derivative of state vector
1565: . shiftV - shift to apply, see note below
1566: - shiftA - shift to apply, see note below
1568: Output Parameters:
1569: + J - Jacobian matrix
1570: - P - optional preconditioning matrix
1572: Level: developer
1574: Notes:
1575: If F(t,U,V,A)=0 is the DAE, the required Jacobian is
1577: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1579: Most users should not need to explicitly call this routine, as it
1580: is used internally within the nonlinear solvers.
1582: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1583: @*/
1584: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1585: {
1586: DM dm;
1587: TSI2JacobianFn *I2Jacobian;
1588: void *ctx;
1589: TSRHSJacobianFn *rhsjacobian;
1591: PetscFunctionBegin;
1599: PetscCall(TSGetDM(ts, &dm));
1600: PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1601: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1603: if (!I2Jacobian) {
1604: PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1605: PetscFunctionReturn(PETSC_SUCCESS);
1606: }
1608: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, J, P));
1609: PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1610: if (rhsjacobian) {
1611: Mat Jrhs, Prhs;
1612: PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1613: PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1614: PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1615: if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1616: }
1618: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, J, P));
1619: PetscFunctionReturn(PETSC_SUCCESS);
1620: }
1622: /*@C
1623: TSSetTransientVariable - sets function to transform from state to transient variables
1625: Logically Collective
1627: Input Parameters:
1628: + ts - time stepping context on which to change the transient variable
1629: . tvar - a function that transforms to transient variables, see `TSTransientVariableFn` for the calling sequence
1630: - ctx - a context for tvar
1632: Level: advanced
1634: Notes:
1635: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1636: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1637: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1638: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1639: evaluated via the chain rule, as in
1640: .vb
1641: dF/dP + shift * dF/dCdot dC/dP.
1642: .ve
1644: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSTransientVariableFn`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1645: @*/
1646: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariableFn *tvar, void *ctx)
1647: {
1648: DM dm;
1650: PetscFunctionBegin;
1652: PetscCall(TSGetDM(ts, &dm));
1653: PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1654: PetscFunctionReturn(PETSC_SUCCESS);
1655: }
1657: /*@
1658: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1660: Logically Collective
1662: Input Parameters:
1663: + ts - TS on which to compute
1664: - U - state vector to be transformed to transient variables
1666: Output Parameter:
1667: . C - transient (conservative) variable
1669: Level: developer
1671: Developer Notes:
1672: If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed.
1673: This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are
1674: being used.
1676: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1677: @*/
1678: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1679: {
1680: DM dm;
1681: DMTS dmts;
1683: PetscFunctionBegin;
1686: PetscCall(TSGetDM(ts, &dm));
1687: PetscCall(DMGetDMTS(dm, &dmts));
1688: if (dmts->ops->transientvar) {
1690: PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1691: }
1692: PetscFunctionReturn(PETSC_SUCCESS);
1693: }
1695: /*@
1696: TSHasTransientVariable - determine whether transient variables have been set
1698: Logically Collective
1700: Input Parameter:
1701: . ts - `TS` on which to compute
1703: Output Parameter:
1704: . has - `PETSC_TRUE` if transient variables have been set
1706: Level: developer
1708: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1709: @*/
1710: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1711: {
1712: DM dm;
1713: DMTS dmts;
1715: PetscFunctionBegin;
1717: PetscCall(TSGetDM(ts, &dm));
1718: PetscCall(DMGetDMTS(dm, &dmts));
1719: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1720: PetscFunctionReturn(PETSC_SUCCESS);
1721: }
1723: /*@
1724: TS2SetSolution - Sets the initial solution and time derivative vectors
1725: for use by the `TS` routines handling second order equations.
1727: Logically Collective
1729: Input Parameters:
1730: + ts - the `TS` context obtained from `TSCreate()`
1731: . u - the solution vector
1732: - v - the time derivative vector
1734: Level: beginner
1736: .seealso: [](ch_ts), `TS`
1737: @*/
1738: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1739: {
1740: PetscFunctionBegin;
1744: PetscCall(TSSetSolution(ts, u));
1745: PetscCall(PetscObjectReference((PetscObject)v));
1746: PetscCall(VecDestroy(&ts->vec_dot));
1747: ts->vec_dot = v;
1748: PetscFunctionReturn(PETSC_SUCCESS);
1749: }
1751: /*@
1752: TS2GetSolution - Returns the solution and time derivative at the present timestep
1753: for second order equations.
1755: Not Collective
1757: Input Parameter:
1758: . ts - the `TS` context obtained from `TSCreate()`
1760: Output Parameters:
1761: + u - the vector containing the solution
1762: - v - the vector containing the time derivative
1764: Level: intermediate
1766: Notes:
1767: It is valid to call this routine inside the function
1768: that you are evaluating in order to move to the new timestep. This vector not
1769: changed until the solution at the next timestep has been calculated.
1771: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1772: @*/
1773: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1774: {
1775: PetscFunctionBegin;
1777: if (u) PetscAssertPointer(u, 2);
1778: if (v) PetscAssertPointer(v, 3);
1779: if (u) *u = ts->vec_sol;
1780: if (v) *v = ts->vec_dot;
1781: PetscFunctionReturn(PETSC_SUCCESS);
1782: }
1784: /*@C
1785: TSLoad - Loads a `TS` that has been stored in binary with `TSView()`.
1787: Collective
1789: Input Parameters:
1790: + ts - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1791: some related function before a call to `TSLoad()`.
1792: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1794: Level: intermediate
1796: Note:
1797: The type is determined by the data in the file, any type set into the `TS` before this call is ignored.
1799: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1800: @*/
1801: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1802: {
1803: PetscBool isbinary;
1804: PetscInt classid;
1805: char type[256];
1806: DMTS sdm;
1807: DM dm;
1809: PetscFunctionBegin;
1812: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1813: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1815: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1816: PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1817: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1818: PetscCall(TSSetType(ts, type));
1819: PetscTryTypeMethod(ts, load, viewer);
1820: PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1821: PetscCall(DMLoad(dm, viewer));
1822: PetscCall(TSSetDM(ts, dm));
1823: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1824: PetscCall(VecLoad(ts->vec_sol, viewer));
1825: PetscCall(DMGetDMTS(ts->dm, &sdm));
1826: PetscCall(DMTSLoad(sdm, viewer));
1827: PetscFunctionReturn(PETSC_SUCCESS);
1828: }
1830: #include <petscdraw.h>
1831: #if defined(PETSC_HAVE_SAWS)
1832: #include <petscviewersaws.h>
1833: #endif
1835: /*@C
1836: TSViewFromOptions - View a `TS` based on values in the options database
1838: Collective
1840: Input Parameters:
1841: + ts - the `TS` context
1842: . obj - Optional object that provides the prefix for the options database keys
1843: - name - command line option string to be passed by user
1845: Level: intermediate
1847: .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1848: @*/
1849: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1850: {
1851: PetscFunctionBegin;
1853: PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1854: PetscFunctionReturn(PETSC_SUCCESS);
1855: }
1857: /*@C
1858: TSView - Prints the `TS` data structure.
1860: Collective
1862: Input Parameters:
1863: + ts - the `TS` context obtained from `TSCreate()`
1864: - viewer - visualization context
1866: Options Database Key:
1867: . -ts_view - calls `TSView()` at end of `TSStep()`
1869: Level: beginner
1871: Notes:
1872: The available visualization contexts include
1873: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1874: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1875: output where only the first processor opens
1876: the file. All other processors send their
1877: data to the first processor to print.
1879: The user can open an alternative visualization context with
1880: `PetscViewerASCIIOpen()` - output to a specified file.
1882: In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).
1884: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1885: @*/
1886: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1887: {
1888: TSType type;
1889: PetscBool iascii, isstring, isundials, isbinary, isdraw;
1890: DMTS sdm;
1891: #if defined(PETSC_HAVE_SAWS)
1892: PetscBool issaws;
1893: #endif
1895: PetscFunctionBegin;
1897: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1899: PetscCheckSameComm(ts, 1, viewer, 2);
1901: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1902: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1903: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1904: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1905: #if defined(PETSC_HAVE_SAWS)
1906: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1907: #endif
1908: if (iascii) {
1909: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer));
1910: if (ts->ops->view) {
1911: PetscCall(PetscViewerASCIIPushTab(viewer));
1912: PetscUseTypeMethod(ts, view, viewer);
1913: PetscCall(PetscViewerASCIIPopTab(viewer));
1914: }
1915: if (ts->max_steps < PETSC_MAX_INT) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps));
1916: if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time));
1917: if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs));
1918: if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs));
1919: if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs));
1920: if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs));
1921: if (ts->usessnes) {
1922: PetscBool lin;
1923: if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its));
1924: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its));
1925: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
1926: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures));
1927: }
1928: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject));
1929: if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, "));
1930: else PetscCall(PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol));
1931: if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of absolute error tolerances\n"));
1932: else PetscCall(PetscViewerASCIIPrintf(viewer, " using absolute error tolerance of %g\n", (double)ts->atol));
1933: PetscCall(PetscViewerASCIIPushTab(viewer));
1934: PetscCall(TSAdaptView(ts->adapt, viewer));
1935: PetscCall(PetscViewerASCIIPopTab(viewer));
1936: } else if (isstring) {
1937: PetscCall(TSGetType(ts, &type));
1938: PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type));
1939: PetscTryTypeMethod(ts, view, viewer);
1940: } else if (isbinary) {
1941: PetscInt classid = TS_FILE_CLASSID;
1942: MPI_Comm comm;
1943: PetscMPIInt rank;
1944: char type[256];
1946: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
1947: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1948: if (rank == 0) {
1949: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
1950: PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
1951: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
1952: }
1953: PetscTryTypeMethod(ts, view, viewer);
1954: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1955: PetscCall(DMView(ts->dm, viewer));
1956: PetscCall(VecView(ts->vec_sol, viewer));
1957: PetscCall(DMGetDMTS(ts->dm, &sdm));
1958: PetscCall(DMTSView(sdm, viewer));
1959: } else if (isdraw) {
1960: PetscDraw draw;
1961: char str[36];
1962: PetscReal x, y, bottom, h;
1964: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1965: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
1966: PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
1967: PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
1968: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
1969: bottom = y - h;
1970: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
1971: PetscTryTypeMethod(ts, view, viewer);
1972: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
1973: if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
1974: PetscCall(PetscDrawPopCurrentPoint(draw));
1975: #if defined(PETSC_HAVE_SAWS)
1976: } else if (issaws) {
1977: PetscMPIInt rank;
1978: const char *name;
1980: PetscCall(PetscObjectGetName((PetscObject)ts, &name));
1981: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1982: if (!((PetscObject)ts)->amsmem && rank == 0) {
1983: char dir[1024];
1985: PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
1986: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
1987: PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
1988: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
1989: PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
1990: }
1991: PetscTryTypeMethod(ts, view, viewer);
1992: #endif
1993: }
1994: if (ts->snes && ts->usessnes) {
1995: PetscCall(PetscViewerASCIIPushTab(viewer));
1996: PetscCall(SNESView(ts->snes, viewer));
1997: PetscCall(PetscViewerASCIIPopTab(viewer));
1998: }
1999: PetscCall(DMGetDMTS(ts->dm, &sdm));
2000: PetscCall(DMTSView(sdm, viewer));
2002: PetscCall(PetscViewerASCIIPushTab(viewer));
2003: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials));
2004: PetscCall(PetscViewerASCIIPopTab(viewer));
2005: PetscFunctionReturn(PETSC_SUCCESS);
2006: }
2008: /*@
2009: TSSetApplicationContext - Sets an optional user-defined context for
2010: the timesteppers.
2012: Logically Collective
2014: Input Parameters:
2015: + ts - the `TS` context obtained from `TSCreate()`
2016: - usrP - user context
2018: Level: intermediate
2020: Fortran Notes:
2021: You must write a Fortran interface definition for this
2022: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
2024: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2025: @*/
2026: PetscErrorCode TSSetApplicationContext(TS ts, void *usrP)
2027: {
2028: PetscFunctionBegin;
2030: ts->user = usrP;
2031: PetscFunctionReturn(PETSC_SUCCESS);
2032: }
2034: /*@
2035: TSGetApplicationContext - Gets the user-defined context for the
2036: timestepper that was set with `TSSetApplicationContext()`
2038: Not Collective
2040: Input Parameter:
2041: . ts - the `TS` context obtained from `TSCreate()`
2043: Output Parameter:
2044: . usrP - user context
2046: Level: intermediate
2048: Fortran Notes:
2049: You must write a Fortran interface definition for this
2050: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
2052: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2053: @*/
2054: PetscErrorCode TSGetApplicationContext(TS ts, void *usrP)
2055: {
2056: PetscFunctionBegin;
2058: *(void **)usrP = ts->user;
2059: PetscFunctionReturn(PETSC_SUCCESS);
2060: }
2062: /*@
2063: TSGetStepNumber - Gets the number of time steps completed.
2065: Not Collective
2067: Input Parameter:
2068: . ts - the `TS` context obtained from `TSCreate()`
2070: Output Parameter:
2071: . steps - number of steps completed so far
2073: Level: intermediate
2075: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2076: @*/
2077: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2078: {
2079: PetscFunctionBegin;
2081: PetscAssertPointer(steps, 2);
2082: *steps = ts->steps;
2083: PetscFunctionReturn(PETSC_SUCCESS);
2084: }
2086: /*@
2087: TSSetStepNumber - Sets the number of steps completed.
2089: Logically Collective
2091: Input Parameters:
2092: + ts - the `TS` context
2093: - steps - number of steps completed so far
2095: Level: developer
2097: Note:
2098: For most uses of the `TS` solvers the user need not explicitly call
2099: `TSSetStepNumber()`, as the step counter is appropriately updated in
2100: `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2101: reinitialize timestepping by setting the step counter to zero (and time
2102: to the initial time) to solve a similar problem with different initial
2103: conditions or parameters. Other possible use case is to continue
2104: timestepping from a previously interrupted run in such a way that `TS`
2105: monitors will be called with a initial nonzero step counter.
2107: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2108: @*/
2109: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2110: {
2111: PetscFunctionBegin;
2114: PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2115: ts->steps = steps;
2116: PetscFunctionReturn(PETSC_SUCCESS);
2117: }
2119: /*@
2120: TSSetTimeStep - Allows one to reset the timestep at any time,
2121: useful for simple pseudo-timestepping codes.
2123: Logically Collective
2125: Input Parameters:
2126: + ts - the `TS` context obtained from `TSCreate()`
2127: - time_step - the size of the timestep
2129: Level: intermediate
2131: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2132: @*/
2133: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2134: {
2135: PetscFunctionBegin;
2138: ts->time_step = time_step;
2139: PetscFunctionReturn(PETSC_SUCCESS);
2140: }
2142: /*@
2143: TSSetExactFinalTime - Determines whether to adapt the final time step to
2144: match the exact final time, interpolate solution to the exact final time,
2145: or just return at the final time `TS` computed.
2147: Logically Collective
2149: Input Parameters:
2150: + ts - the time-step context
2151: - eftopt - exact final time option
2152: .vb
2153: TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded
2154: TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2155: TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2156: .ve
2158: Options Database Key:
2159: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step at runtime
2161: Level: beginner
2163: Note:
2164: If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2165: then the final time you selected.
2167: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2168: @*/
2169: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2170: {
2171: PetscFunctionBegin;
2174: ts->exact_final_time = eftopt;
2175: PetscFunctionReturn(PETSC_SUCCESS);
2176: }
2178: /*@
2179: TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`
2181: Not Collective
2183: Input Parameter:
2184: . ts - the `TS` context
2186: Output Parameter:
2187: . eftopt - exact final time option
2189: Level: beginner
2191: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2192: @*/
2193: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2194: {
2195: PetscFunctionBegin;
2197: PetscAssertPointer(eftopt, 2);
2198: *eftopt = ts->exact_final_time;
2199: PetscFunctionReturn(PETSC_SUCCESS);
2200: }
2202: /*@
2203: TSGetTimeStep - Gets the current timestep size.
2205: Not Collective
2207: Input Parameter:
2208: . ts - the `TS` context obtained from `TSCreate()`
2210: Output Parameter:
2211: . dt - the current timestep size
2213: Level: intermediate
2215: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2216: @*/
2217: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2218: {
2219: PetscFunctionBegin;
2221: PetscAssertPointer(dt, 2);
2222: *dt = ts->time_step;
2223: PetscFunctionReturn(PETSC_SUCCESS);
2224: }
2226: /*@
2227: TSGetSolution - Returns the solution at the present timestep. It
2228: is valid to call this routine inside the function that you are evaluating
2229: in order to move to the new timestep. This vector not changed until
2230: the solution at the next timestep has been calculated.
2232: Not Collective, but v returned is parallel if ts is parallel
2234: Input Parameter:
2235: . ts - the `TS` context obtained from `TSCreate()`
2237: Output Parameter:
2238: . v - the vector containing the solution
2240: Level: intermediate
2242: Note:
2243: If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
2244: final time. It returns the solution at the next timestep.
2246: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2247: @*/
2248: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2249: {
2250: PetscFunctionBegin;
2252: PetscAssertPointer(v, 2);
2253: *v = ts->vec_sol;
2254: PetscFunctionReturn(PETSC_SUCCESS);
2255: }
2257: /*@
2258: TSGetSolutionComponents - Returns any solution components at the present
2259: timestep, if available for the time integration method being used.
2260: Solution components are quantities that share the same size and
2261: structure as the solution vector.
2263: Not Collective, but v returned is parallel if ts is parallel
2265: Input Parameters:
2266: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2267: . n - If v is `NULL`, then the number of solution components is
2268: returned through n, else the n-th solution component is
2269: returned in v.
2270: - v - the vector containing the n-th solution component
2271: (may be `NULL` to use this function to find out
2272: the number of solutions components).
2274: Level: advanced
2276: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2277: @*/
2278: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2279: {
2280: PetscFunctionBegin;
2282: if (!ts->ops->getsolutioncomponents) *n = 0;
2283: else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2284: PetscFunctionReturn(PETSC_SUCCESS);
2285: }
2287: /*@
2288: TSGetAuxSolution - Returns an auxiliary solution at the present
2289: timestep, if available for the time integration method being used.
2291: Not Collective, but v returned is parallel if ts is parallel
2293: Input Parameters:
2294: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2295: - v - the vector containing the auxiliary solution
2297: Level: intermediate
2299: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2300: @*/
2301: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2302: {
2303: PetscFunctionBegin;
2305: if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2306: else PetscCall(VecZeroEntries(*v));
2307: PetscFunctionReturn(PETSC_SUCCESS);
2308: }
2310: /*@
2311: TSGetTimeError - Returns the estimated error vector, if the chosen
2312: `TSType` has an error estimation functionality and `TSSetTimeError()` was called
2314: Not Collective, but v returned is parallel if ts is parallel
2316: Input Parameters:
2317: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2318: . n - current estimate (n=0) or previous one (n=-1)
2319: - v - the vector containing the error (same size as the solution).
2321: Level: intermediate
2323: Note:
2324: MUST call after `TSSetUp()`
2326: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2327: @*/
2328: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2329: {
2330: PetscFunctionBegin;
2332: if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2333: else PetscCall(VecZeroEntries(*v));
2334: PetscFunctionReturn(PETSC_SUCCESS);
2335: }
2337: /*@
2338: TSSetTimeError - Sets the estimated error vector, if the chosen
2339: `TSType` has an error estimation functionality. This can be used
2340: to restart such a time integrator with a given error vector.
2342: Not Collective, but v returned is parallel if ts is parallel
2344: Input Parameters:
2345: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2346: - v - the vector containing the error (same size as the solution).
2348: Level: intermediate
2350: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
2351: @*/
2352: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2353: {
2354: PetscFunctionBegin;
2356: PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2357: PetscTryTypeMethod(ts, settimeerror, v);
2358: PetscFunctionReturn(PETSC_SUCCESS);
2359: }
2361: /* ----- Routines to initialize and destroy a timestepper ---- */
2362: /*@
2363: TSSetProblemType - Sets the type of problem to be solved.
2365: Not collective
2367: Input Parameters:
2368: + ts - The `TS`
2369: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2370: .vb
2371: U_t - A U = 0 (linear)
2372: U_t - A(t) U = 0 (linear)
2373: F(t,U,U_t) = 0 (nonlinear)
2374: .ve
2376: Level: beginner
2378: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2379: @*/
2380: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2381: {
2382: PetscFunctionBegin;
2384: ts->problem_type = type;
2385: if (type == TS_LINEAR) {
2386: SNES snes;
2387: PetscCall(TSGetSNES(ts, &snes));
2388: PetscCall(SNESSetType(snes, SNESKSPONLY));
2389: }
2390: PetscFunctionReturn(PETSC_SUCCESS);
2391: }
2393: /*@C
2394: TSGetProblemType - Gets the type of problem to be solved.
2396: Not collective
2398: Input Parameter:
2399: . ts - The `TS`
2401: Output Parameter:
2402: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2403: .vb
2404: M U_t = A U
2405: M(t) U_t = A(t) U
2406: F(t,U,U_t)
2407: .ve
2409: Level: beginner
2411: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2412: @*/
2413: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2414: {
2415: PetscFunctionBegin;
2417: PetscAssertPointer(type, 2);
2418: *type = ts->problem_type;
2419: PetscFunctionReturn(PETSC_SUCCESS);
2420: }
2422: /*
2423: Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2424: */
2425: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2426: {
2427: PetscBool isnone;
2429: PetscFunctionBegin;
2430: PetscCall(TSGetAdapt(ts, &ts->adapt));
2431: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2433: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2434: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2435: else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2436: PetscFunctionReturn(PETSC_SUCCESS);
2437: }
2439: /*@
2440: TSSetUp - Sets up the internal data structures for the later use of a timestepper.
2442: Collective
2444: Input Parameter:
2445: . ts - the `TS` context obtained from `TSCreate()`
2447: Level: advanced
2449: Note:
2450: For basic use of the `TS` solvers the user need not explicitly call
2451: `TSSetUp()`, since these actions will automatically occur during
2452: the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this
2453: phase separately, `TSSetUp()` should be called after `TSCreate()`
2454: and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.
2456: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2457: @*/
2458: PetscErrorCode TSSetUp(TS ts)
2459: {
2460: DM dm;
2461: PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2462: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2463: TSIFunctionFn *ifun;
2464: TSIJacobianFn *ijac;
2465: TSI2JacobianFn *i2jac;
2466: TSRHSJacobianFn *rhsjac;
2468: PetscFunctionBegin;
2470: if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
2472: if (!((PetscObject)ts)->type_name) {
2473: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2474: PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2475: }
2477: if (!ts->vec_sol) {
2478: PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2479: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2480: }
2482: if (ts->tspan) {
2483: if (!ts->tspan->vecs_sol) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2484: }
2485: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2486: PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2487: ts->Jacp = ts->Jacprhs;
2488: }
2490: if (ts->quadraturets) {
2491: PetscCall(TSSetUp(ts->quadraturets));
2492: PetscCall(VecDestroy(&ts->vec_costintegrand));
2493: PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2494: }
2496: PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
2497: if (rhsjac == TSComputeRHSJacobianConstant) {
2498: Mat Amat, Pmat;
2499: SNES snes;
2500: PetscCall(TSGetSNES(ts, &snes));
2501: PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
2502: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2503: * have displaced the RHS matrix */
2504: if (Amat && Amat == ts->Arhs) {
2505: /* 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 */
2506: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
2507: PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
2508: PetscCall(MatDestroy(&Amat));
2509: }
2510: if (Pmat && Pmat == ts->Brhs) {
2511: PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
2512: PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
2513: PetscCall(MatDestroy(&Pmat));
2514: }
2515: }
2517: PetscCall(TSGetAdapt(ts, &ts->adapt));
2518: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2520: PetscTryTypeMethod(ts, setup);
2522: PetscCall(TSSetExactFinalTimeDefault(ts));
2524: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2525: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2526: */
2527: PetscCall(TSGetDM(ts, &dm));
2528: PetscCall(DMSNESGetFunction(dm, &func, NULL));
2529: if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));
2531: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2532: Otherwise, the SNES will use coloring internally to form the Jacobian.
2533: */
2534: PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2535: PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2536: PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2537: PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2538: if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));
2540: /* if time integration scheme has a starting method, call it */
2541: PetscTryTypeMethod(ts, startingmethod);
2543: ts->setupcalled = PETSC_TRUE;
2544: PetscFunctionReturn(PETSC_SUCCESS);
2545: }
2547: /*@
2548: TSReset - Resets a `TS` context and removes any allocated `Vec`s and `Mat`s.
2550: Collective
2552: Input Parameter:
2553: . ts - the `TS` context obtained from `TSCreate()`
2555: Level: beginner
2557: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()`
2558: @*/
2559: PetscErrorCode TSReset(TS ts)
2560: {
2561: TS_RHSSplitLink ilink = ts->tsrhssplit, next;
2563: PetscFunctionBegin;
2566: PetscTryTypeMethod(ts, reset);
2567: if (ts->snes) PetscCall(SNESReset(ts->snes));
2568: if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));
2570: PetscCall(MatDestroy(&ts->Arhs));
2571: PetscCall(MatDestroy(&ts->Brhs));
2572: PetscCall(VecDestroy(&ts->Frhs));
2573: PetscCall(VecDestroy(&ts->vec_sol));
2574: PetscCall(VecDestroy(&ts->vec_sol0));
2575: PetscCall(VecDestroy(&ts->vec_dot));
2576: PetscCall(VecDestroy(&ts->vatol));
2577: PetscCall(VecDestroy(&ts->vrtol));
2578: PetscCall(VecDestroyVecs(ts->nwork, &ts->work));
2580: PetscCall(MatDestroy(&ts->Jacprhs));
2581: PetscCall(MatDestroy(&ts->Jacp));
2582: if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2583: if (ts->quadraturets) {
2584: PetscCall(TSReset(ts->quadraturets));
2585: PetscCall(VecDestroy(&ts->vec_costintegrand));
2586: }
2587: while (ilink) {
2588: next = ilink->next;
2589: PetscCall(TSDestroy(&ilink->ts));
2590: PetscCall(PetscFree(ilink->splitname));
2591: PetscCall(ISDestroy(&ilink->is));
2592: PetscCall(PetscFree(ilink));
2593: ilink = next;
2594: }
2595: ts->tsrhssplit = NULL;
2596: ts->num_rhs_splits = 0;
2597: if (ts->tspan) {
2598: PetscCall(PetscFree(ts->tspan->span_times));
2599: PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2600: PetscCall(PetscFree(ts->tspan));
2601: }
2602: ts->rhsjacobian.time = PETSC_MIN_REAL;
2603: ts->rhsjacobian.scale = 1.0;
2604: ts->ijacobian.shift = 1.0;
2605: ts->setupcalled = PETSC_FALSE;
2606: PetscFunctionReturn(PETSC_SUCCESS);
2607: }
2609: static PetscErrorCode TSResizeReset(TS);
2611: /*@C
2612: TSDestroy - Destroys the timestepper context that was created
2613: with `TSCreate()`.
2615: Collective
2617: Input Parameter:
2618: . ts - the `TS` context obtained from `TSCreate()`
2620: Level: beginner
2622: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2623: @*/
2624: PetscErrorCode TSDestroy(TS *ts)
2625: {
2626: PetscFunctionBegin;
2627: if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2629: if (--((PetscObject)*ts)->refct > 0) {
2630: *ts = NULL;
2631: PetscFunctionReturn(PETSC_SUCCESS);
2632: }
2634: PetscCall(TSReset(*ts));
2635: PetscCall(TSAdjointReset(*ts));
2636: if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2637: PetscCall(TSResizeReset(*ts));
2639: /* if memory was published with SAWs then destroy it */
2640: PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2641: PetscTryTypeMethod(*ts, destroy);
2643: PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));
2645: PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2646: PetscCall(TSEventDestroy(&(*ts)->event));
2648: PetscCall(SNESDestroy(&(*ts)->snes));
2649: PetscCall(DMDestroy(&(*ts)->dm));
2650: PetscCall(TSMonitorCancel(*ts));
2651: PetscCall(TSAdjointMonitorCancel(*ts));
2653: PetscCall(TSDestroy(&(*ts)->quadraturets));
2654: PetscCall(PetscHeaderDestroy(ts));
2655: PetscFunctionReturn(PETSC_SUCCESS);
2656: }
2658: /*@
2659: TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2660: a `TS` (timestepper) context. Valid only for nonlinear problems.
2662: Not Collective, but snes is parallel if ts is parallel
2664: Input Parameter:
2665: . ts - the `TS` context obtained from `TSCreate()`
2667: Output Parameter:
2668: . snes - the nonlinear solver context
2670: Level: beginner
2672: Notes:
2673: The user can then directly manipulate the `SNES` context to set various
2674: options, etc. Likewise, the user can then extract and manipulate the
2675: `KSP`, and `PC` contexts as well.
2677: `TSGetSNES()` does not work for integrators that do not use `SNES`; in
2678: this case `TSGetSNES()` returns `NULL` in `snes`.
2680: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2681: @*/
2682: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2683: {
2684: PetscFunctionBegin;
2686: PetscAssertPointer(snes, 2);
2687: if (!ts->snes) {
2688: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2689: PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2690: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2691: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2692: if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2693: if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2694: }
2695: *snes = ts->snes;
2696: PetscFunctionReturn(PETSC_SUCCESS);
2697: }
2699: /*@
2700: TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the timestepping context
2702: Collective
2704: Input Parameters:
2705: + ts - the `TS` context obtained from `TSCreate()`
2706: - snes - the nonlinear solver context
2708: Level: developer
2710: Note:
2711: Most users should have the `TS` created by calling `TSGetSNES()`
2713: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2714: @*/
2715: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2716: {
2717: PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);
2719: PetscFunctionBegin;
2722: PetscCall(PetscObjectReference((PetscObject)snes));
2723: PetscCall(SNESDestroy(&ts->snes));
2725: ts->snes = snes;
2727: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2728: PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2729: if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2730: PetscFunctionReturn(PETSC_SUCCESS);
2731: }
2733: /*@
2734: TSGetKSP - Returns the `KSP` (linear solver) associated with
2735: a `TS` (timestepper) context.
2737: Not Collective, but `ksp` is parallel if `ts` is parallel
2739: Input Parameter:
2740: . ts - the `TS` context obtained from `TSCreate()`
2742: Output Parameter:
2743: . ksp - the nonlinear solver context
2745: Level: beginner
2747: Notes:
2748: The user can then directly manipulate the `KSP` context to set various
2749: options, etc. Likewise, the user can then extract and manipulate the
2750: `PC` context as well.
2752: `TSGetKSP()` does not work for integrators that do not use `KSP`;
2753: in this case `TSGetKSP()` returns `NULL` in `ksp`.
2755: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2756: @*/
2757: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2758: {
2759: SNES snes;
2761: PetscFunctionBegin;
2763: PetscAssertPointer(ksp, 2);
2764: PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2765: PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
2766: PetscCall(TSGetSNES(ts, &snes));
2767: PetscCall(SNESGetKSP(snes, ksp));
2768: PetscFunctionReturn(PETSC_SUCCESS);
2769: }
2771: /* ----------- Routines to set solver parameters ---------- */
2773: /*@
2774: TSSetMaxSteps - Sets the maximum number of steps to use.
2776: Logically Collective
2778: Input Parameters:
2779: + ts - the `TS` context obtained from `TSCreate()`
2780: - maxsteps - maximum number of steps to use
2782: Options Database Key:
2783: . -ts_max_steps <maxsteps> - Sets maxsteps
2785: Level: intermediate
2787: Note:
2788: The default maximum number of steps is 5000
2790: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2791: @*/
2792: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2793: {
2794: PetscFunctionBegin;
2797: PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2798: ts->max_steps = maxsteps;
2799: PetscFunctionReturn(PETSC_SUCCESS);
2800: }
2802: /*@
2803: TSGetMaxSteps - Gets the maximum number of steps to use.
2805: Not Collective
2807: Input Parameter:
2808: . ts - the `TS` context obtained from `TSCreate()`
2810: Output Parameter:
2811: . maxsteps - maximum number of steps to use
2813: Level: advanced
2815: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2816: @*/
2817: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2818: {
2819: PetscFunctionBegin;
2821: PetscAssertPointer(maxsteps, 2);
2822: *maxsteps = ts->max_steps;
2823: PetscFunctionReturn(PETSC_SUCCESS);
2824: }
2826: /*@
2827: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
2829: Logically Collective
2831: Input Parameters:
2832: + ts - the `TS` context obtained from `TSCreate()`
2833: - maxtime - final time to step to
2835: Options Database Key:
2836: . -ts_max_time <maxtime> - Sets maxtime
2838: Level: intermediate
2840: Notes:
2841: The default maximum time is 5.0
2843: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2844: @*/
2845: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2846: {
2847: PetscFunctionBegin;
2850: ts->max_time = maxtime;
2851: PetscFunctionReturn(PETSC_SUCCESS);
2852: }
2854: /*@
2855: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
2857: Not Collective
2859: Input Parameter:
2860: . ts - the `TS` context obtained from `TSCreate()`
2862: Output Parameter:
2863: . maxtime - final time to step to
2865: Level: advanced
2867: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2868: @*/
2869: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2870: {
2871: PetscFunctionBegin;
2873: PetscAssertPointer(maxtime, 2);
2874: *maxtime = ts->max_time;
2875: PetscFunctionReturn(PETSC_SUCCESS);
2876: }
2878: // PetscClangLinter pragma disable: -fdoc-*
2879: /*@
2880: TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.
2882: Level: deprecated
2884: @*/
2885: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2886: {
2887: PetscFunctionBegin;
2889: PetscCall(TSSetTime(ts, initial_time));
2890: PetscCall(TSSetTimeStep(ts, time_step));
2891: PetscFunctionReturn(PETSC_SUCCESS);
2892: }
2894: // PetscClangLinter pragma disable: -fdoc-*
2895: /*@
2896: TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.
2898: Level: deprecated
2900: @*/
2901: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2902: {
2903: PetscFunctionBegin;
2905: if (maxsteps) {
2906: PetscAssertPointer(maxsteps, 2);
2907: *maxsteps = ts->max_steps;
2908: }
2909: if (maxtime) {
2910: PetscAssertPointer(maxtime, 3);
2911: *maxtime = ts->max_time;
2912: }
2913: PetscFunctionReturn(PETSC_SUCCESS);
2914: }
2916: // PetscClangLinter pragma disable: -fdoc-*
2917: /*@
2918: TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.
2920: Level: deprecated
2922: @*/
2923: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2924: {
2925: PetscFunctionBegin;
2929: if (maxsteps >= 0) ts->max_steps = maxsteps;
2930: if (maxtime != (PetscReal)PETSC_DEFAULT) ts->max_time = maxtime;
2931: PetscFunctionReturn(PETSC_SUCCESS);
2932: }
2934: // PetscClangLinter pragma disable: -fdoc-*
2935: /*@
2936: TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.
2938: Level: deprecated
2940: @*/
2941: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
2942: {
2943: return TSGetStepNumber(ts, steps);
2944: }
2946: // PetscClangLinter pragma disable: -fdoc-*
2947: /*@
2948: TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.
2950: Level: deprecated
2952: @*/
2953: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
2954: {
2955: return TSGetStepNumber(ts, steps);
2956: }
2958: /*@
2959: TSSetSolution - Sets the initial solution vector
2960: for use by the `TS` routines.
2962: Logically Collective
2964: Input Parameters:
2965: + ts - the `TS` context obtained from `TSCreate()`
2966: - u - the solution vector
2968: Level: beginner
2970: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
2971: @*/
2972: PetscErrorCode TSSetSolution(TS ts, Vec u)
2973: {
2974: DM dm;
2976: PetscFunctionBegin;
2979: PetscCall(PetscObjectReference((PetscObject)u));
2980: PetscCall(VecDestroy(&ts->vec_sol));
2981: ts->vec_sol = u;
2983: PetscCall(TSGetDM(ts, &dm));
2984: PetscCall(DMShellSetGlobalVector(dm, u));
2985: PetscFunctionReturn(PETSC_SUCCESS);
2986: }
2988: /*@C
2989: TSSetPreStep - Sets the general-purpose function
2990: called once at the beginning of each time step.
2992: Logically Collective
2994: Input Parameters:
2995: + ts - The `TS` context obtained from `TSCreate()`
2996: - func - The function
2998: Calling sequence of `func`:
2999: . ts - the `TS` context
3001: Level: intermediate
3003: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3004: @*/
3005: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
3006: {
3007: PetscFunctionBegin;
3009: ts->prestep = func;
3010: PetscFunctionReturn(PETSC_SUCCESS);
3011: }
3013: /*@
3014: TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`
3016: Collective
3018: Input Parameter:
3019: . ts - The `TS` context obtained from `TSCreate()`
3021: Level: developer
3023: Note:
3024: `TSPreStep()` is typically used within time stepping implementations,
3025: so most users would not generally call this routine themselves.
3027: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3028: @*/
3029: PetscErrorCode TSPreStep(TS ts)
3030: {
3031: PetscFunctionBegin;
3033: if (ts->prestep) {
3034: Vec U;
3035: PetscObjectId idprev;
3036: PetscBool sameObject;
3037: PetscObjectState sprev, spost;
3039: PetscCall(TSGetSolution(ts, &U));
3040: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3041: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3042: PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3043: PetscCall(TSGetSolution(ts, &U));
3044: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3045: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3046: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3047: }
3048: PetscFunctionReturn(PETSC_SUCCESS);
3049: }
3051: /*@C
3052: TSSetPreStage - Sets the general-purpose function
3053: called once at the beginning of each stage.
3055: Logically Collective
3057: Input Parameters:
3058: + ts - The `TS` context obtained from `TSCreate()`
3059: - func - The function
3061: Calling sequence of `func`:
3062: + ts - the `TS` context
3063: - stagetime - the stage time
3065: Level: intermediate
3067: Note:
3068: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3069: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3070: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3072: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3073: @*/
3074: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime))
3075: {
3076: PetscFunctionBegin;
3078: ts->prestage = func;
3079: PetscFunctionReturn(PETSC_SUCCESS);
3080: }
3082: /*@C
3083: TSSetPostStage - Sets the general-purpose function, provided with `TSSetPostStep()`,
3084: called once at the end of each stage.
3086: Logically Collective
3088: Input Parameters:
3089: + ts - The `TS` context obtained from `TSCreate()`
3090: - func - The function
3092: Calling sequence of `func`:
3093: + ts - the `TS` context
3094: . stagetime - the stage time
3095: . stageindex - the stage index
3096: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3098: Level: intermediate
3100: Note:
3101: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3102: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3103: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3105: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3106: @*/
3107: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
3108: {
3109: PetscFunctionBegin;
3111: ts->poststage = func;
3112: PetscFunctionReturn(PETSC_SUCCESS);
3113: }
3115: /*@C
3116: TSSetPostEvaluate - Sets the general-purpose function
3117: called once at the end of each step evaluation.
3119: Logically Collective
3121: Input Parameters:
3122: + ts - The `TS` context obtained from `TSCreate()`
3123: - func - The function
3125: Calling sequence of `func`:
3126: . ts - the `TS` context
3128: Level: intermediate
3130: Note:
3131: Semantically, `TSSetPostEvaluate()` differs from `TSSetPostStep()` since the function it sets is called before event-handling
3132: thus guaranteeing the same solution (computed by the time-stepper) will be passed to it. On the other hand, `TSPostStep()`
3133: may be passed a different solution, possibly changed by the event handler. `TSPostEvaluate()` is called after the next step
3134: solution is evaluated allowing to modify it, if need be. The solution can be obtained with `TSGetSolution()`, the time step
3135: with `TSGetTimeStep()`, and the time at the start of the step is available via `TSGetTime()`
3137: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3138: @*/
3139: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
3140: {
3141: PetscFunctionBegin;
3143: ts->postevaluate = func;
3144: PetscFunctionReturn(PETSC_SUCCESS);
3145: }
3147: /*@
3148: TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`
3150: Collective
3152: Input Parameters:
3153: + ts - The `TS` context obtained from `TSCreate()`
3154: - stagetime - The absolute time of the current stage
3156: Level: developer
3158: Note:
3159: `TSPreStage()` is typically used within time stepping implementations,
3160: most users would not generally call this routine themselves.
3162: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3163: @*/
3164: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3165: {
3166: PetscFunctionBegin;
3168: if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3169: PetscFunctionReturn(PETSC_SUCCESS);
3170: }
3172: /*@
3173: TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`
3175: Collective
3177: Input Parameters:
3178: + ts - The `TS` context obtained from `TSCreate()`
3179: . stagetime - The absolute time of the current stage
3180: . stageindex - Stage number
3181: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3183: Level: developer
3185: Note:
3186: `TSPostStage()` is typically used within time stepping implementations,
3187: most users would not generally call this routine themselves.
3189: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3190: @*/
3191: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3192: {
3193: PetscFunctionBegin;
3195: if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3196: PetscFunctionReturn(PETSC_SUCCESS);
3197: }
3199: /*@
3200: TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`
3202: Collective
3204: Input Parameter:
3205: . ts - The `TS` context obtained from `TSCreate()`
3207: Level: developer
3209: Note:
3210: `TSPostEvaluate()` is typically used within time stepping implementations,
3211: most users would not generally call this routine themselves.
3213: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3214: @*/
3215: PetscErrorCode TSPostEvaluate(TS ts)
3216: {
3217: PetscFunctionBegin;
3219: if (ts->postevaluate) {
3220: Vec U;
3221: PetscObjectState sprev, spost;
3223: PetscCall(TSGetSolution(ts, &U));
3224: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3225: PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3226: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3227: if (sprev != spost) PetscCall(TSRestartStep(ts));
3228: }
3229: PetscFunctionReturn(PETSC_SUCCESS);
3230: }
3232: /*@C
3233: TSSetPostStep - Sets the general-purpose function
3234: called once at the end of each time step.
3236: Logically Collective
3238: Input Parameters:
3239: + ts - The `TS` context obtained from `TSCreate()`
3240: - func - The function
3242: Calling sequence of `func`:
3243: . ts - the `TS` context
3245: Level: intermediate
3247: Note:
3248: The function set by `TSSetPostStep()` is called after each successful step. The solution vector
3249: obtained by `TSGetSolution()` may be different than that computed at the step end if the event handler
3250: locates an event and `TSPostEvent()` modifies it. Use `TSSetPostEvaluate()` if an unmodified solution is needed instead.
3252: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3253: @*/
3254: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
3255: {
3256: PetscFunctionBegin;
3258: ts->poststep = func;
3259: PetscFunctionReturn(PETSC_SUCCESS);
3260: }
3262: /*@
3263: TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`
3265: Collective
3267: Input Parameter:
3268: . ts - The `TS` context obtained from `TSCreate()`
3270: Note:
3271: `TSPostStep()` is typically used within time stepping implementations,
3272: so most users would not generally call this routine themselves.
3274: Level: developer
3276: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()`
3277: @*/
3278: PetscErrorCode TSPostStep(TS ts)
3279: {
3280: PetscFunctionBegin;
3282: if (ts->poststep) {
3283: Vec U;
3284: PetscObjectId idprev;
3285: PetscBool sameObject;
3286: PetscObjectState sprev, spost;
3288: PetscCall(TSGetSolution(ts, &U));
3289: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3290: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3291: PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3292: PetscCall(TSGetSolution(ts, &U));
3293: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3294: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3295: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3296: }
3297: PetscFunctionReturn(PETSC_SUCCESS);
3298: }
3300: /*@
3301: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3303: Collective
3305: Input Parameters:
3306: + ts - time stepping context
3307: - t - time to interpolate to
3309: Output Parameter:
3310: . U - state at given time
3312: Level: intermediate
3314: Developer Notes:
3315: `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3317: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3318: @*/
3319: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3320: {
3321: PetscFunctionBegin;
3324: 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);
3325: PetscUseTypeMethod(ts, interpolate, t, U);
3326: PetscFunctionReturn(PETSC_SUCCESS);
3327: }
3329: /*@
3330: TSStep - Steps one time step
3332: Collective
3334: Input Parameter:
3335: . ts - the `TS` context obtained from `TSCreate()`
3337: Level: developer
3339: Notes:
3340: The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.
3342: The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
3343: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3345: This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
3346: time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.
3348: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3349: @*/
3350: PetscErrorCode TSStep(TS ts)
3351: {
3352: static PetscBool cite = PETSC_FALSE;
3353: PetscReal ptime;
3355: PetscFunctionBegin;
3357: PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3358: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3359: " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3360: " journal = {arXiv e-preprints},\n"
3361: " eprint = {1806.01437},\n"
3362: " archivePrefix = {arXiv},\n"
3363: " year = {2018}\n}\n",
3364: &cite));
3365: PetscCall(TSSetUp(ts));
3366: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3367: if (ts->tspan)
3368: ts->tspan->worktol = 0; /* In each step of TSSolve() 'tspan->worktol' will be meaningfully defined (later) only once:
3369: in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */
3371: 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>");
3372: 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()");
3373: 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");
3375: if (!ts->vec_sol0) PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sol0));
3376: PetscCall(VecCopy(ts->vec_sol, ts->vec_sol0));
3377: ts->time_step0 = ts->time_step;
3379: if (!ts->steps) ts->ptime_prev = ts->ptime;
3380: ptime = ts->ptime;
3382: ts->ptime_prev_rollback = ts->ptime_prev;
3383: ts->reason = TS_CONVERGED_ITERATING;
3385: PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3386: PetscUseTypeMethod(ts, step);
3387: PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));
3389: if (ts->reason >= 0) {
3390: ts->ptime_prev = ptime;
3391: ts->steps++;
3392: ts->steprollback = PETSC_FALSE;
3393: ts->steprestart = PETSC_FALSE;
3394: }
3396: if (ts->reason < 0 && ts->errorifstepfailed) {
3397: PetscCall(TSMonitorCancel(ts));
3398: 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]);
3399: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3400: }
3401: PetscFunctionReturn(PETSC_SUCCESS);
3402: }
3404: /*@
3405: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3406: at the end of a time step with a given order of accuracy.
3408: Collective
3410: Input Parameters:
3411: + ts - time stepping context
3412: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3414: Input/Output Parameter:
3415: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3416: on output, the actual order of the error evaluation
3418: Output Parameter:
3419: . wlte - the weighted local truncation error norm
3421: Level: advanced
3423: Note:
3424: If the timestepper cannot evaluate the error in a particular step
3425: (eg. in the first step or restart steps after event handling),
3426: this routine returns wlte=-1.0 .
3428: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3429: @*/
3430: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3431: {
3432: PetscFunctionBegin;
3436: if (order) PetscAssertPointer(order, 3);
3438: PetscAssertPointer(wlte, 4);
3439: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3440: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3441: PetscFunctionReturn(PETSC_SUCCESS);
3442: }
3444: /*@
3445: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3447: Collective
3449: Input Parameters:
3450: + ts - time stepping context
3451: . order - desired order of accuracy
3452: - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)
3454: Output Parameter:
3455: . U - state at the end of the current step
3457: Level: advanced
3459: Notes:
3460: This function cannot be called until all stages have been evaluated.
3462: 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.
3464: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3465: @*/
3466: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3467: {
3468: PetscFunctionBegin;
3472: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3473: PetscFunctionReturn(PETSC_SUCCESS);
3474: }
3476: /*@C
3477: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3479: Not collective
3481: Input Parameter:
3482: . ts - time stepping context
3484: Output Parameter:
3485: . initCondition - The function which computes an initial condition
3487: Calling sequence of `initCondition`:
3488: + ts - The timestepping context
3489: - u - The input vector in which the initial condition is stored
3491: Level: advanced
3493: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3494: @*/
3495: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3496: {
3497: PetscFunctionBegin;
3499: PetscAssertPointer(initCondition, 2);
3500: *initCondition = ts->ops->initcondition;
3501: PetscFunctionReturn(PETSC_SUCCESS);
3502: }
3504: /*@C
3505: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3507: Logically collective
3509: Input Parameters:
3510: + ts - time stepping context
3511: - initCondition - The function which computes an initial condition
3513: Calling sequence of `initCondition`:
3514: + ts - The timestepping context
3515: - e - The input vector in which the initial condition is to be stored
3517: Level: advanced
3519: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3520: @*/
3521: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3522: {
3523: PetscFunctionBegin;
3526: ts->ops->initcondition = initCondition;
3527: PetscFunctionReturn(PETSC_SUCCESS);
3528: }
3530: /*@
3531: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3533: Collective
3535: Input Parameters:
3536: + ts - time stepping context
3537: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3539: Level: advanced
3541: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3542: @*/
3543: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3544: {
3545: PetscFunctionBegin;
3548: PetscTryTypeMethod(ts, initcondition, u);
3549: PetscFunctionReturn(PETSC_SUCCESS);
3550: }
3552: /*@C
3553: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3555: Not collective
3557: Input Parameter:
3558: . ts - time stepping context
3560: Output Parameter:
3561: . exactError - The function which computes the solution error
3563: Calling sequence of `exactError`:
3564: + ts - The timestepping context
3565: . u - The approximate solution vector
3566: - e - The vector in which the error is stored
3568: Level: advanced
3570: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3571: @*/
3572: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3573: {
3574: PetscFunctionBegin;
3576: PetscAssertPointer(exactError, 2);
3577: *exactError = ts->ops->exacterror;
3578: PetscFunctionReturn(PETSC_SUCCESS);
3579: }
3581: /*@C
3582: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3584: Logically collective
3586: Input Parameters:
3587: + ts - time stepping context
3588: - exactError - The function which computes the solution error
3590: Calling sequence of `exactError`:
3591: + ts - The timestepping context
3592: . u - The approximate solution vector
3593: - e - The vector in which the error is stored
3595: Level: advanced
3597: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3598: @*/
3599: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3600: {
3601: PetscFunctionBegin;
3604: ts->ops->exacterror = exactError;
3605: PetscFunctionReturn(PETSC_SUCCESS);
3606: }
3608: /*@
3609: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3611: Collective
3613: Input Parameters:
3614: + ts - time stepping context
3615: . u - The approximate solution
3616: - e - The `Vec` used to store the error
3618: Level: advanced
3620: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3621: @*/
3622: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3623: {
3624: PetscFunctionBegin;
3628: PetscTryTypeMethod(ts, exacterror, u, e);
3629: PetscFunctionReturn(PETSC_SUCCESS);
3630: }
3632: /*@C
3633: TSSetResize - Sets the resize callbacks.
3635: Logically Collective
3637: Input Parameters:
3638: + ts - The `TS` context obtained from `TSCreate()`
3639: . rollback - Whether a resize will restart the step
3640: . setup - The setup function
3641: . transfer - The transfer function
3642: - ctx - [optional] The user-defined context
3644: Calling sequence of `setup`:
3645: + ts - the `TS` context
3646: . step - the current step
3647: . time - the current time
3648: . state - the current vector of state
3649: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3650: - ctx - user defined context
3652: Calling sequence of `transfer`:
3653: + ts - the `TS` context
3654: . nv - the number of vectors to be transferred
3655: . vecsin - array of vectors to be transferred
3656: . vecsout - array of transferred vectors
3657: - ctx - user defined context
3659: Notes:
3660: The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
3661: depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
3662: and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
3663: that the size of the discrete problem has changed.
3664: In both cases, the solver will collect the needed vectors that will be
3665: transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
3666: current solution vector, and other vectors needed by the specific solver used.
3667: For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3668: Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3669: will be automatically reset if the sizes are changed and they must be specified again by the user
3670: inside the `transfer` function.
3671: The input and output arrays passed to `transfer` are allocated by PETSc.
3672: Vectors in `vecsout` must be created by the user.
3673: Ownership of vectors in `vecsout` is transferred to PETSc.
3675: Level: advanced
3677: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3678: @*/
3679: 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)
3680: {
3681: PetscFunctionBegin;
3683: ts->resizerollback = rollback;
3684: ts->resizesetup = setup;
3685: ts->resizetransfer = transfer;
3686: ts->resizectx = ctx;
3687: PetscFunctionReturn(PETSC_SUCCESS);
3688: }
3690: /*
3691: TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.
3693: Collective
3695: Input Parameters:
3696: + ts - The `TS` context obtained from `TSCreate()`
3697: - 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.
3699: Level: developer
3701: Note:
3702: `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3703: used within time stepping implementations,
3704: so most users would not generally call this routine themselves.
3706: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3707: @*/
3708: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3709: {
3710: PetscFunctionBegin;
3712: PetscTryTypeMethod(ts, resizeregister, flg);
3713: /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3714: PetscFunctionReturn(PETSC_SUCCESS);
3715: }
3717: static PetscErrorCode TSResizeReset(TS ts)
3718: {
3719: PetscFunctionBegin;
3721: PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3722: PetscFunctionReturn(PETSC_SUCCESS);
3723: }
3725: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3726: {
3727: PetscFunctionBegin;
3730: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3731: if (ts->resizetransfer) {
3732: PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3733: PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3734: }
3735: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3736: PetscFunctionReturn(PETSC_SUCCESS);
3737: }
3739: /*@C
3740: TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.
3742: Collective
3744: Input Parameters:
3745: + ts - The `TS` context obtained from `TSCreate()`
3746: . name - A string identifying the vector
3747: - vec - The vector
3749: Level: developer
3751: Note:
3752: `TSResizeRegisterVec()` is typically used within time stepping implementations,
3753: so most users would not generally call this routine themselves.
3755: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3756: @*/
3757: PetscErrorCode TSResizeRegisterVec(TS ts, const char *name, Vec vec)
3758: {
3759: PetscFunctionBegin;
3761: PetscAssertPointer(name, 2);
3763: PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3764: PetscFunctionReturn(PETSC_SUCCESS);
3765: }
3767: /*@C
3768: TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.
3770: Collective
3772: Input Parameters:
3773: + ts - The `TS` context obtained from `TSCreate()`
3774: . name - A string identifying the vector
3775: - vec - The vector
3777: Level: developer
3779: Note:
3780: `TSResizeRetrieveVec()` is typically used within time stepping implementations,
3781: so most users would not generally call this routine themselves.
3783: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
3784: @*/
3785: PetscErrorCode TSResizeRetrieveVec(TS ts, const char *name, Vec *vec)
3786: {
3787: PetscFunctionBegin;
3789: PetscAssertPointer(name, 2);
3790: PetscAssertPointer(vec, 3);
3791: PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
3792: PetscFunctionReturn(PETSC_SUCCESS);
3793: }
3795: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
3796: {
3797: PetscInt cnt;
3798: PetscObjectList tmp;
3799: Vec *vecsin = NULL;
3800: const char **namesin = NULL;
3802: PetscFunctionBegin;
3803: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
3804: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
3805: if (names) PetscCall(PetscMalloc1(cnt, &namesin));
3806: if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
3807: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
3808: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
3809: if (vecs) vecsin[cnt] = (Vec)tmp->obj;
3810: if (names) namesin[cnt] = tmp->name;
3811: cnt++;
3812: }
3813: }
3814: if (nv) *nv = cnt;
3815: if (names) *names = namesin;
3816: if (vecs) *vecs = vecsin;
3817: PetscFunctionReturn(PETSC_SUCCESS);
3818: }
3820: /*@
3821: TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`
3823: Collective
3825: Input Parameter:
3826: . ts - The `TS` context obtained from `TSCreate()`
3828: Level: developer
3830: Note:
3831: `TSResize()` is typically used within time stepping implementations,
3832: so most users would not generally call this routine themselves.
3834: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3835: @*/
3836: PetscErrorCode TSResize(TS ts)
3837: {
3838: PetscInt nv = 0;
3839: const char **names = NULL;
3840: Vec *vecsin = NULL;
3841: const char *solname = "ts:vec_sol";
3843: PetscFunctionBegin;
3845: if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
3846: if (ts->resizesetup) {
3847: PetscBool flg = PETSC_FALSE;
3849: PetscCall(VecLockReadPush(ts->vec_sol));
3850: PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &flg, ts->resizectx));
3851: PetscCall(VecLockReadPop(ts->vec_sol));
3852: if (flg) {
3853: if (ts->resizerollback) {
3854: PetscCall(TSRollBack(ts));
3855: ts->time_step = ts->time_step0;
3856: }
3857: PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
3858: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
3859: }
3860: }
3862: PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
3863: if (nv) {
3864: Vec *vecsout, vecsol;
3866: /* Reset internal objects */
3867: PetscCall(TSReset(ts));
3869: /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
3870: PetscCall(PetscCalloc1(nv, &vecsout));
3871: PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
3872: for (PetscInt i = 0; i < nv; i++) {
3873: const char *name;
3874: char *oname;
3876: PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
3877: PetscCall(PetscStrallocpy(name, &oname));
3878: PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
3879: if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
3880: PetscCall(PetscFree(oname));
3881: PetscCall(VecDestroy(&vecsout[i]));
3882: }
3883: PetscCall(PetscFree(vecsout));
3884: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */
3886: PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
3887: if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
3888: PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
3889: }
3891: PetscCall(PetscFree(names));
3892: PetscCall(PetscFree(vecsin));
3893: PetscCall(TSResizeReset(ts));
3894: PetscFunctionReturn(PETSC_SUCCESS);
3895: }
3897: /*@
3898: TSSolve - Steps the requested number of timesteps.
3900: Collective
3902: Input Parameters:
3903: + ts - the `TS` context obtained from `TSCreate()`
3904: - u - the solution vector (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
3905: otherwise must contain the initial conditions and will contain the solution at the final requested time
3907: Level: beginner
3909: Notes:
3910: The final time returned by this function may be different from the time of the internally
3911: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
3912: stepped over the final time.
3914: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3915: @*/
3916: PetscErrorCode TSSolve(TS ts, Vec u)
3917: {
3918: Vec solution;
3920: PetscFunctionBegin;
3924: PetscCall(TSSetExactFinalTimeDefault(ts));
3925: 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 */
3926: if (!ts->vec_sol || u == ts->vec_sol) {
3927: PetscCall(VecDuplicate(u, &solution));
3928: PetscCall(TSSetSolution(ts, solution));
3929: PetscCall(VecDestroy(&solution)); /* grant ownership */
3930: }
3931: PetscCall(VecCopy(u, ts->vec_sol));
3932: PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
3933: } else if (u) PetscCall(TSSetSolution(ts, u));
3934: PetscCall(TSSetUp(ts));
3935: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3937: 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>");
3938: 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()");
3939: 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");
3940: 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");
3942: 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 */
3943: PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]));
3944: ts->tspan->spanctr = 1;
3945: }
3947: if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));
3949: /* reset number of steps only when the step is not restarted. ARKIMEX
3950: restarts the step after an event. Resetting these counters in such case causes
3951: TSTrajectory to incorrectly save the output files
3952: */
3953: /* reset time step and iteration counters */
3954: if (!ts->steps) {
3955: ts->ksp_its = 0;
3956: ts->snes_its = 0;
3957: ts->num_snes_failures = 0;
3958: ts->reject = 0;
3959: ts->steprestart = PETSC_TRUE;
3960: ts->steprollback = PETSC_FALSE;
3961: ts->rhsjacobian.time = PETSC_MIN_REAL;
3962: }
3964: /* make sure initial time step does not overshoot final time or the next point in tspan */
3965: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
3966: PetscReal maxdt;
3967: PetscReal dt = ts->time_step;
3969: if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
3970: else maxdt = ts->max_time - ts->ptime;
3971: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
3972: }
3973: ts->reason = TS_CONVERGED_ITERATING;
3975: {
3976: PetscViewer viewer;
3977: PetscViewerFormat format;
3978: PetscBool flg;
3979: static PetscBool incall = PETSC_FALSE;
3981: if (!incall) {
3982: /* Estimate the convergence rate of the time discretization */
3983: PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
3984: if (flg) {
3985: PetscConvEst conv;
3986: DM dm;
3987: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
3988: PetscInt Nf;
3989: PetscBool checkTemporal = PETSC_TRUE;
3991: incall = PETSC_TRUE;
3992: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
3993: PetscCall(TSGetDM(ts, &dm));
3994: PetscCall(DMGetNumFields(dm, &Nf));
3995: PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
3996: PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
3997: PetscCall(PetscConvEstUseTS(conv, checkTemporal));
3998: PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
3999: PetscCall(PetscConvEstSetFromOptions(conv));
4000: PetscCall(PetscConvEstSetUp(conv));
4001: PetscCall(PetscConvEstGetConvRate(conv, alpha));
4002: PetscCall(PetscViewerPushFormat(viewer, format));
4003: PetscCall(PetscConvEstRateView(conv, alpha, viewer));
4004: PetscCall(PetscViewerPopFormat(viewer));
4005: PetscCall(PetscOptionsRestoreViewer(&viewer));
4006: PetscCall(PetscConvEstDestroy(&conv));
4007: PetscCall(PetscFree(alpha));
4008: incall = PETSC_FALSE;
4009: }
4010: }
4011: }
4013: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));
4015: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4016: PetscUseTypeMethod(ts, solve);
4017: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4018: ts->solvetime = ts->ptime;
4019: solution = ts->vec_sol;
4020: } else { /* Step the requested number of timesteps. */
4021: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4022: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4024: if (!ts->steps) {
4025: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4026: PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4027: }
4029: while (!ts->reason) {
4030: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4031: if (!ts->steprollback) PetscCall(TSPreStep(ts));
4032: PetscCall(TSStep(ts));
4033: if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4034: if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4035: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4036: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4037: PetscCall(TSForwardCostIntegral(ts));
4038: if (ts->reason >= 0) ts->steps++;
4039: }
4040: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4041: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4042: PetscCall(TSForwardStep(ts));
4043: if (ts->reason >= 0) ts->steps++;
4044: }
4045: PetscCall(TSPostEvaluate(ts));
4046: 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. */
4047: if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4048: if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
4049: /* check convergence */
4050: if (!ts->reason) {
4051: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4052: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4053: }
4054: if (!ts->steprollback) {
4055: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4056: PetscCall(TSPostStep(ts));
4057: if (!ts->resizerollback) PetscCall(TSResize(ts));
4059: if (ts->tspan && ts->tspan->spanctr < ts->tspan->num_span_times) {
4060: PetscCheck(ts->tspan->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(tspan->worktol > 0) in TSSolve()");
4061: 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++]));
4062: }
4063: }
4064: }
4065: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4067: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4068: if (!u) u = ts->vec_sol;
4069: PetscCall(TSInterpolate(ts, ts->max_time, u));
4070: ts->solvetime = ts->max_time;
4071: solution = u;
4072: PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4073: } else {
4074: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4075: ts->solvetime = ts->ptime;
4076: solution = ts->vec_sol;
4077: }
4078: }
4080: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4081: PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4082: PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4083: if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4084: PetscFunctionReturn(PETSC_SUCCESS);
4085: }
4087: /*@
4088: TSGetTime - Gets the time of the most recently completed step.
4090: Not Collective
4092: Input Parameter:
4093: . ts - the `TS` context obtained from `TSCreate()`
4095: Output Parameter:
4096: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
4098: Level: beginner
4100: Note:
4101: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
4102: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
4104: .seealso: [](ch_ts), `TS`, `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4105: @*/
4106: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4107: {
4108: PetscFunctionBegin;
4110: PetscAssertPointer(t, 2);
4111: *t = ts->ptime;
4112: PetscFunctionReturn(PETSC_SUCCESS);
4113: }
4115: /*@
4116: TSGetPrevTime - Gets the starting time of the previously completed step.
4118: Not Collective
4120: Input Parameter:
4121: . ts - the `TS` context obtained from `TSCreate()`
4123: Output Parameter:
4124: . t - the previous time
4126: Level: beginner
4128: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4129: @*/
4130: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4131: {
4132: PetscFunctionBegin;
4134: PetscAssertPointer(t, 2);
4135: *t = ts->ptime_prev;
4136: PetscFunctionReturn(PETSC_SUCCESS);
4137: }
4139: /*@
4140: TSSetTime - Allows one to reset the time.
4142: Logically Collective
4144: Input Parameters:
4145: + ts - the `TS` context obtained from `TSCreate()`
4146: - t - the time
4148: Level: intermediate
4150: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4151: @*/
4152: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4153: {
4154: PetscFunctionBegin;
4157: ts->ptime = t;
4158: PetscFunctionReturn(PETSC_SUCCESS);
4159: }
4161: /*@C
4162: TSSetOptionsPrefix - Sets the prefix used for searching for all
4163: TS options in the database.
4165: Logically Collective
4167: Input Parameters:
4168: + ts - The `TS` context
4169: - prefix - The prefix to prepend to all option names
4171: Level: advanced
4173: Note:
4174: A hyphen (-) must NOT be given at the beginning of the prefix name.
4175: The first character of all runtime options is AUTOMATICALLY the
4176: hyphen.
4178: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4179: @*/
4180: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4181: {
4182: SNES snes;
4184: PetscFunctionBegin;
4186: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4187: PetscCall(TSGetSNES(ts, &snes));
4188: PetscCall(SNESSetOptionsPrefix(snes, prefix));
4189: PetscFunctionReturn(PETSC_SUCCESS);
4190: }
4192: /*@C
4193: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4194: TS options in the database.
4196: Logically Collective
4198: Input Parameters:
4199: + ts - The `TS` context
4200: - prefix - The prefix to prepend to all option names
4202: Level: advanced
4204: Note:
4205: A hyphen (-) must NOT be given at the beginning of the prefix name.
4206: The first character of all runtime options is AUTOMATICALLY the
4207: hyphen.
4209: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4210: @*/
4211: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4212: {
4213: SNES snes;
4215: PetscFunctionBegin;
4217: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4218: PetscCall(TSGetSNES(ts, &snes));
4219: PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4220: PetscFunctionReturn(PETSC_SUCCESS);
4221: }
4223: /*@C
4224: TSGetOptionsPrefix - Sets the prefix used for searching for all
4225: `TS` options in the database.
4227: Not Collective
4229: Input Parameter:
4230: . ts - The `TS` context
4232: Output Parameter:
4233: . prefix - A pointer to the prefix string used
4235: Level: intermediate
4237: Fortran Notes:
4238: The user should pass in a string 'prefix' of
4239: sufficient length to hold the prefix.
4241: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4242: @*/
4243: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4244: {
4245: PetscFunctionBegin;
4247: PetscAssertPointer(prefix, 2);
4248: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4249: PetscFunctionReturn(PETSC_SUCCESS);
4250: }
4252: /*@C
4253: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4255: Not Collective, but parallel objects are returned if ts is parallel
4257: Input Parameter:
4258: . ts - The `TS` context obtained from `TSCreate()`
4260: Output Parameters:
4261: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`)
4262: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`)
4263: . func - Function to compute the Jacobian of the RHS (or `NULL`)
4264: - ctx - User-defined context for Jacobian evaluation routine (or `NULL`)
4266: Level: intermediate
4268: Note:
4269: You can pass in `NULL` for any return argument you do not need.
4271: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4273: @*/
4274: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, void **ctx)
4275: {
4276: DM dm;
4278: PetscFunctionBegin;
4279: if (Amat || Pmat) {
4280: SNES snes;
4281: PetscCall(TSGetSNES(ts, &snes));
4282: PetscCall(SNESSetUpMatrices(snes));
4283: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4284: }
4285: PetscCall(TSGetDM(ts, &dm));
4286: PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4287: PetscFunctionReturn(PETSC_SUCCESS);
4288: }
4290: /*@C
4291: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4293: Not Collective, but parallel objects are returned if ts is parallel
4295: Input Parameter:
4296: . ts - The `TS` context obtained from `TSCreate()`
4298: Output Parameters:
4299: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4300: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4301: . f - The function to compute the matrices
4302: - ctx - User-defined context for Jacobian evaluation routine
4304: Level: advanced
4306: Note:
4307: You can pass in `NULL` for any return argument you do not need.
4309: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4310: @*/
4311: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, void **ctx)
4312: {
4313: DM dm;
4315: PetscFunctionBegin;
4316: if (Amat || Pmat) {
4317: SNES snes;
4318: PetscCall(TSGetSNES(ts, &snes));
4319: PetscCall(SNESSetUpMatrices(snes));
4320: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4321: }
4322: PetscCall(TSGetDM(ts, &dm));
4323: PetscCall(DMTSGetIJacobian(dm, f, ctx));
4324: PetscFunctionReturn(PETSC_SUCCESS);
4325: }
4327: #include <petsc/private/dmimpl.h>
4328: /*@
4329: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4331: Logically Collective
4333: Input Parameters:
4334: + ts - the `TS` integrator object
4335: - dm - the dm, cannot be `NULL`
4337: Level: intermediate
4339: Notes:
4340: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4341: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4342: different problems using the same function space.
4344: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4345: @*/
4346: PetscErrorCode TSSetDM(TS ts, DM dm)
4347: {
4348: SNES snes;
4349: DMTS tsdm;
4351: PetscFunctionBegin;
4354: PetscCall(PetscObjectReference((PetscObject)dm));
4355: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4356: if (ts->dm->dmts && !dm->dmts) {
4357: PetscCall(DMCopyDMTS(ts->dm, dm));
4358: PetscCall(DMGetDMTS(ts->dm, &tsdm));
4359: /* Grant write privileges to the replacement DM */
4360: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4361: }
4362: PetscCall(DMDestroy(&ts->dm));
4363: }
4364: ts->dm = dm;
4366: PetscCall(TSGetSNES(ts, &snes));
4367: PetscCall(SNESSetDM(snes, dm));
4368: PetscFunctionReturn(PETSC_SUCCESS);
4369: }
4371: /*@
4372: TSGetDM - Gets the `DM` that may be used by some preconditioners
4374: Not Collective
4376: Input Parameter:
4377: . ts - the `TS`
4379: Output Parameter:
4380: . dm - the `DM`
4382: Level: intermediate
4384: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4385: @*/
4386: PetscErrorCode TSGetDM(TS ts, DM *dm)
4387: {
4388: PetscFunctionBegin;
4390: if (!ts->dm) {
4391: PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4392: if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4393: }
4394: *dm = ts->dm;
4395: PetscFunctionReturn(PETSC_SUCCESS);
4396: }
4398: /*@
4399: SNESTSFormFunction - Function to evaluate nonlinear residual
4401: Logically Collective
4403: Input Parameters:
4404: + snes - nonlinear solver
4405: . U - the current state at which to evaluate the residual
4406: - ctx - user context, must be a TS
4408: Output Parameter:
4409: . F - the nonlinear residual
4411: Level: advanced
4413: Note:
4414: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4415: It is most frequently passed to `MatFDColoringSetFunction()`.
4417: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4418: @*/
4419: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4420: {
4421: TS ts = (TS)ctx;
4423: PetscFunctionBegin;
4428: PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4429: PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4430: PetscFunctionReturn(PETSC_SUCCESS);
4431: }
4433: /*@
4434: SNESTSFormJacobian - Function to evaluate the Jacobian
4436: Collective
4438: Input Parameters:
4439: + snes - nonlinear solver
4440: . U - the current state at which to evaluate the residual
4441: - ctx - user context, must be a `TS`
4443: Output Parameters:
4444: + A - the Jacobian
4445: - B - the preconditioning matrix (may be the same as A)
4447: Level: developer
4449: Note:
4450: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4452: .seealso: [](ch_ts), `SNESSetJacobian()`
4453: @*/
4454: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4455: {
4456: TS ts = (TS)ctx;
4458: PetscFunctionBegin;
4464: PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4465: PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4466: PetscFunctionReturn(PETSC_SUCCESS);
4467: }
4469: /*@C
4470: TSComputeRHSFunctionLinear - Evaluate the right-hand side via the user-provided Jacobian, for linear problems Udot = A U only
4472: Collective
4474: Input Parameters:
4475: + ts - time stepping context
4476: . t - time at which to evaluate
4477: . U - state at which to evaluate
4478: - ctx - context
4480: Output Parameter:
4481: . F - right-hand side
4483: Level: intermediate
4485: Note:
4486: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right-hand side for linear problems.
4487: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4489: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4490: @*/
4491: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4492: {
4493: Mat Arhs, Brhs;
4495: PetscFunctionBegin;
4496: PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4497: /* undo the damage caused by shifting */
4498: PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4499: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4500: PetscCall(MatMult(Arhs, U, F));
4501: PetscFunctionReturn(PETSC_SUCCESS);
4502: }
4504: /*@C
4505: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
4507: Collective
4509: Input Parameters:
4510: + ts - time stepping context
4511: . t - time at which to evaluate
4512: . U - state at which to evaluate
4513: - ctx - context
4515: Output Parameters:
4516: + A - pointer to operator
4517: - B - pointer to preconditioning matrix
4519: Level: intermediate
4521: Note:
4522: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4524: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4525: @*/
4526: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4527: {
4528: PetscFunctionBegin;
4529: PetscFunctionReturn(PETSC_SUCCESS);
4530: }
4532: /*@C
4533: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4535: Collective
4537: Input Parameters:
4538: + ts - time stepping context
4539: . t - time at which to evaluate
4540: . U - state at which to evaluate
4541: . Udot - time derivative of state vector
4542: - ctx - context
4544: Output Parameter:
4545: . F - left hand side
4547: Level: intermediate
4549: Notes:
4550: 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
4551: user is required to write their own `TSComputeIFunction()`.
4552: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4553: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4555: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4557: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4558: @*/
4559: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4560: {
4561: Mat A, B;
4563: PetscFunctionBegin;
4564: PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4565: PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4566: PetscCall(MatMult(A, Udot, F));
4567: PetscFunctionReturn(PETSC_SUCCESS);
4568: }
4570: /*@C
4571: TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE
4573: Collective
4575: Input Parameters:
4576: + ts - time stepping context
4577: . t - time at which to evaluate
4578: . U - state at which to evaluate
4579: . Udot - time derivative of state vector
4580: . shift - shift to apply
4581: - ctx - context
4583: Output Parameters:
4584: + A - pointer to operator
4585: - B - pointer to matrix from which the preconditioner is built (often `A`)
4587: Level: advanced
4589: Notes:
4590: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4592: It is only appropriate for problems of the form
4594: $$
4595: M \dot{U} = F(U,t)
4596: $$
4598: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4599: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4600: an implicit operator of the form
4602: $$
4603: shift*M + J
4604: $$
4606: 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
4607: a copy of M or reassemble it when requested.
4609: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4610: @*/
4611: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4612: {
4613: PetscFunctionBegin;
4614: PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4615: ts->ijacobian.shift = shift;
4616: PetscFunctionReturn(PETSC_SUCCESS);
4617: }
4619: /*@
4620: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4622: Not Collective
4624: Input Parameter:
4625: . ts - the `TS` context
4627: Output Parameter:
4628: . equation_type - see `TSEquationType`
4630: Level: beginner
4632: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4633: @*/
4634: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4635: {
4636: PetscFunctionBegin;
4638: PetscAssertPointer(equation_type, 2);
4639: *equation_type = ts->equation_type;
4640: PetscFunctionReturn(PETSC_SUCCESS);
4641: }
4643: /*@
4644: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4646: Not Collective
4648: Input Parameters:
4649: + ts - the `TS` context
4650: - equation_type - see `TSEquationType`
4652: Level: advanced
4654: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4655: @*/
4656: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4657: {
4658: PetscFunctionBegin;
4660: ts->equation_type = equation_type;
4661: PetscFunctionReturn(PETSC_SUCCESS);
4662: }
4664: /*@
4665: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4667: Not Collective
4669: Input Parameter:
4670: . ts - the `TS` context
4672: Output Parameter:
4673: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4674: manual pages for the individual convergence tests for complete lists
4676: Level: beginner
4678: Note:
4679: Can only be called after the call to `TSSolve()` is complete.
4681: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4682: @*/
4683: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4684: {
4685: PetscFunctionBegin;
4687: PetscAssertPointer(reason, 2);
4688: *reason = ts->reason;
4689: PetscFunctionReturn(PETSC_SUCCESS);
4690: }
4692: /*@
4693: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4695: Logically Collective; reason must contain common value
4697: Input Parameters:
4698: + ts - the `TS` context
4699: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4700: manual pages for the individual convergence tests for complete lists
4702: Level: advanced
4704: Note:
4705: Can only be called while `TSSolve()` is active.
4707: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4708: @*/
4709: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4710: {
4711: PetscFunctionBegin;
4713: ts->reason = reason;
4714: PetscFunctionReturn(PETSC_SUCCESS);
4715: }
4717: /*@
4718: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4720: Not Collective
4722: Input Parameter:
4723: . ts - the `TS` context
4725: Output Parameter:
4726: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4728: Level: beginner
4730: Note:
4731: Can only be called after the call to `TSSolve()` is complete.
4733: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4734: @*/
4735: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4736: {
4737: PetscFunctionBegin;
4739: PetscAssertPointer(ftime, 2);
4740: *ftime = ts->solvetime;
4741: PetscFunctionReturn(PETSC_SUCCESS);
4742: }
4744: /*@
4745: TSGetSNESIterations - Gets the total number of nonlinear iterations
4746: used by the time integrator.
4748: Not Collective
4750: Input Parameter:
4751: . ts - `TS` context
4753: Output Parameter:
4754: . nits - number of nonlinear iterations
4756: Level: intermediate
4758: Note:
4759: This counter is reset to zero for each successive call to `TSSolve()`.
4761: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4762: @*/
4763: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4764: {
4765: PetscFunctionBegin;
4767: PetscAssertPointer(nits, 2);
4768: *nits = ts->snes_its;
4769: PetscFunctionReturn(PETSC_SUCCESS);
4770: }
4772: /*@
4773: TSGetKSPIterations - Gets the total number of linear iterations
4774: used by the time integrator.
4776: Not Collective
4778: Input Parameter:
4779: . ts - `TS` context
4781: Output Parameter:
4782: . lits - number of linear iterations
4784: Level: intermediate
4786: Note:
4787: This counter is reset to zero for each successive call to `TSSolve()`.
4789: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()`
4790: @*/
4791: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4792: {
4793: PetscFunctionBegin;
4795: PetscAssertPointer(lits, 2);
4796: *lits = ts->ksp_its;
4797: PetscFunctionReturn(PETSC_SUCCESS);
4798: }
4800: /*@
4801: TSGetStepRejections - Gets the total number of rejected steps.
4803: Not Collective
4805: Input Parameter:
4806: . ts - `TS` context
4808: Output Parameter:
4809: . rejects - number of steps rejected
4811: Level: intermediate
4813: Note:
4814: This counter is reset to zero for each successive call to `TSSolve()`.
4816: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4817: @*/
4818: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4819: {
4820: PetscFunctionBegin;
4822: PetscAssertPointer(rejects, 2);
4823: *rejects = ts->reject;
4824: PetscFunctionReturn(PETSC_SUCCESS);
4825: }
4827: /*@
4828: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
4830: Not Collective
4832: Input Parameter:
4833: . ts - `TS` context
4835: Output Parameter:
4836: . fails - number of failed nonlinear solves
4838: Level: intermediate
4840: Note:
4841: This counter is reset to zero for each successive call to `TSSolve()`.
4843: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4844: @*/
4845: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4846: {
4847: PetscFunctionBegin;
4849: PetscAssertPointer(fails, 2);
4850: *fails = ts->num_snes_failures;
4851: PetscFunctionReturn(PETSC_SUCCESS);
4852: }
4854: /*@
4855: TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails
4857: Not Collective
4859: Input Parameters:
4860: + ts - `TS` context
4861: - rejects - maximum number of rejected steps, pass -1 for unlimited
4863: Options Database Key:
4864: . -ts_max_reject - Maximum number of step rejections before a step fails
4866: Level: intermediate
4868: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4869: @*/
4870: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4871: {
4872: PetscFunctionBegin;
4874: ts->max_reject = rejects;
4875: PetscFunctionReturn(PETSC_SUCCESS);
4876: }
4878: /*@
4879: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves
4881: Not Collective
4883: Input Parameters:
4884: + ts - `TS` context
4885: - fails - maximum number of failed nonlinear solves, pass -1 for unlimited
4887: Options Database Key:
4888: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
4890: Level: intermediate
4892: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4893: @*/
4894: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4895: {
4896: PetscFunctionBegin;
4898: ts->max_snes_failures = fails;
4899: PetscFunctionReturn(PETSC_SUCCESS);
4900: }
4902: /*@
4903: TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`
4905: Not Collective
4907: Input Parameters:
4908: + ts - `TS` context
4909: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
4911: Options Database Key:
4912: . -ts_error_if_step_fails - Error if no step succeeds
4914: Level: intermediate
4916: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
4917: @*/
4918: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
4919: {
4920: PetscFunctionBegin;
4922: ts->errorifstepfailed = err;
4923: PetscFunctionReturn(PETSC_SUCCESS);
4924: }
4926: /*@
4927: TSGetAdapt - Get the adaptive controller context for the current method
4929: Collective if controller has not yet been created
4931: Input Parameter:
4932: . ts - time stepping context
4934: Output Parameter:
4935: . adapt - adaptive controller
4937: Level: intermediate
4939: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
4940: @*/
4941: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
4942: {
4943: PetscFunctionBegin;
4945: PetscAssertPointer(adapt, 2);
4946: if (!ts->adapt) {
4947: PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
4948: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
4949: }
4950: *adapt = ts->adapt;
4951: PetscFunctionReturn(PETSC_SUCCESS);
4952: }
4954: /*@
4955: TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller
4957: Logically Collective
4959: Input Parameters:
4960: + ts - time integration context
4961: . atol - scalar absolute tolerances, `PETSC_DECIDE` to leave current value
4962: . vatol - vector of absolute tolerances or `NULL`, used in preference to atol if present
4963: . rtol - scalar relative tolerances, `PETSC_DECIDE` to leave current value
4964: - vrtol - vector of relative tolerances or `NULL`, used in preference to atol if present
4966: Options Database Keys:
4967: + -ts_rtol <rtol> - relative tolerance for local truncation error
4968: - -ts_atol <atol> - Absolute tolerance for local truncation error
4970: Level: beginner
4972: Notes:
4973: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
4974: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
4975: computed only for the differential or the algebraic part then this can be done using the vector of
4976: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
4977: differential part and infinity for the algebraic part, the LTE calculation will include only the
4978: differential variables.
4980: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
4981: @*/
4982: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
4983: {
4984: PetscFunctionBegin;
4985: if (atol != (PetscReal)PETSC_DECIDE && atol != (PetscReal)PETSC_DEFAULT) ts->atol = atol;
4986: if (vatol) {
4987: PetscCall(PetscObjectReference((PetscObject)vatol));
4988: PetscCall(VecDestroy(&ts->vatol));
4989: ts->vatol = vatol;
4990: }
4991: if (rtol != (PetscReal)PETSC_DECIDE && rtol != (PetscReal)PETSC_DEFAULT) ts->rtol = rtol;
4992: if (vrtol) {
4993: PetscCall(PetscObjectReference((PetscObject)vrtol));
4994: PetscCall(VecDestroy(&ts->vrtol));
4995: ts->vrtol = vrtol;
4996: }
4997: PetscFunctionReturn(PETSC_SUCCESS);
4998: }
5000: /*@
5001: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
5003: Logically Collective
5005: Input Parameter:
5006: . ts - time integration context
5008: Output Parameters:
5009: + atol - scalar absolute tolerances, `NULL` to ignore
5010: . vatol - vector of absolute tolerances, `NULL` to ignore
5011: . rtol - scalar relative tolerances, `NULL` to ignore
5012: - vrtol - vector of relative tolerances, `NULL` to ignore
5014: Level: beginner
5016: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
5017: @*/
5018: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
5019: {
5020: PetscFunctionBegin;
5021: if (atol) *atol = ts->atol;
5022: if (vatol) *vatol = ts->vatol;
5023: if (rtol) *rtol = ts->rtol;
5024: if (vrtol) *vrtol = ts->vrtol;
5025: PetscFunctionReturn(PETSC_SUCCESS);
5026: }
5028: /*@
5029: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
5031: Collective
5033: Input Parameters:
5034: + ts - time stepping context
5035: . U - state vector, usually ts->vec_sol
5036: . Y - state vector to be compared to U
5037: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5039: Output Parameters:
5040: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5041: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5042: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5044: Options Database Key:
5045: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5047: Level: developer
5049: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5050: @*/
5051: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5052: {
5053: PetscInt norma_loc, norm_loc, normr_loc;
5055: PetscFunctionBegin;
5057: 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));
5058: if (wnormtype == NORM_2) {
5059: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5060: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5061: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5062: }
5063: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5064: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5065: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5066: PetscFunctionReturn(PETSC_SUCCESS);
5067: }
5069: /*@
5070: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5072: Collective
5074: Input Parameters:
5075: + ts - time stepping context
5076: . E - error vector
5077: . U - state vector, usually ts->vec_sol
5078: . Y - state vector, previous time step
5079: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5081: Output Parameters:
5082: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5083: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5084: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5086: Options Database Key:
5087: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5089: Level: developer
5091: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5092: @*/
5093: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5094: {
5095: PetscInt norma_loc, norm_loc, normr_loc;
5097: PetscFunctionBegin;
5099: 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));
5100: if (wnormtype == NORM_2) {
5101: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5102: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5103: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5104: }
5105: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5106: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5107: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5108: PetscFunctionReturn(PETSC_SUCCESS);
5109: }
5111: /*@
5112: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5114: Logically Collective
5116: Input Parameters:
5117: + ts - time stepping context
5118: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5120: Note:
5121: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5123: Level: intermediate
5125: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5126: @*/
5127: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5128: {
5129: PetscFunctionBegin;
5131: ts->cfltime_local = cfltime;
5132: ts->cfltime = -1.;
5133: PetscFunctionReturn(PETSC_SUCCESS);
5134: }
5136: /*@
5137: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5139: Collective
5141: Input Parameter:
5142: . ts - time stepping context
5144: Output Parameter:
5145: . cfltime - maximum stable time step for forward Euler
5147: Level: advanced
5149: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5150: @*/
5151: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5152: {
5153: PetscFunctionBegin;
5154: if (ts->cfltime < 0) PetscCall(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5155: *cfltime = ts->cfltime;
5156: PetscFunctionReturn(PETSC_SUCCESS);
5157: }
5159: /*@
5160: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5162: Input Parameters:
5163: + ts - the `TS` context.
5164: . xl - lower bound.
5165: - xu - upper bound.
5167: Level: advanced
5169: Note:
5170: If this routine is not called then the lower and upper bounds are set to
5171: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5173: .seealso: [](ch_ts), `TS`
5174: @*/
5175: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5176: {
5177: SNES snes;
5179: PetscFunctionBegin;
5180: PetscCall(TSGetSNES(ts, &snes));
5181: PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5182: PetscFunctionReturn(PETSC_SUCCESS);
5183: }
5185: /*@
5186: TSComputeLinearStability - computes the linear stability function at a point
5188: Collective
5190: Input Parameters:
5191: + ts - the `TS` context
5192: . xr - real part of input argument
5193: - xi - imaginary part of input argument
5195: Output Parameters:
5196: + yr - real part of function value
5197: - yi - imaginary part of function value
5199: Level: developer
5201: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5202: @*/
5203: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5204: {
5205: PetscFunctionBegin;
5207: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5208: PetscFunctionReturn(PETSC_SUCCESS);
5209: }
5211: /*@
5212: TSRestartStep - Flags the solver to restart the next step
5214: Collective
5216: Input Parameter:
5217: . ts - the `TS` context obtained from `TSCreate()`
5219: Level: advanced
5221: Notes:
5222: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5223: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5224: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5225: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5226: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5227: discontinuous source terms).
5229: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5230: @*/
5231: PetscErrorCode TSRestartStep(TS ts)
5232: {
5233: PetscFunctionBegin;
5235: ts->steprestart = PETSC_TRUE;
5236: PetscFunctionReturn(PETSC_SUCCESS);
5237: }
5239: /*@
5240: TSRollBack - Rolls back one time step
5242: Collective
5244: Input Parameter:
5245: . ts - the `TS` context obtained from `TSCreate()`
5247: Level: advanced
5249: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5250: @*/
5251: PetscErrorCode TSRollBack(TS ts)
5252: {
5253: PetscFunctionBegin;
5255: PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5256: PetscTryTypeMethod(ts, rollback);
5257: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
5258: ts->time_step = ts->ptime - ts->ptime_prev;
5259: ts->ptime = ts->ptime_prev;
5260: ts->ptime_prev = ts->ptime_prev_rollback;
5261: ts->steps--;
5262: ts->steprollback = PETSC_TRUE;
5263: PetscFunctionReturn(PETSC_SUCCESS);
5264: }
5266: /*@
5267: TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step
5269: Not collective
5271: Input Parameter:
5272: . ts - the `TS` context obtained from `TSCreate()`
5274: Output Parameter:
5275: . flg - the rollback flag
5277: Level: advanced
5279: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
5280: @*/
5281: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
5282: {
5283: PetscFunctionBegin;
5285: PetscAssertPointer(flg, 2);
5286: *flg = ts->steprollback;
5287: PetscFunctionReturn(PETSC_SUCCESS);
5288: }
5290: /*@
5291: TSGetStages - Get the number of stages and stage values
5293: Input Parameter:
5294: . ts - the `TS` context obtained from `TSCreate()`
5296: Output Parameters:
5297: + ns - the number of stages
5298: - Y - the current stage vectors
5300: Level: advanced
5302: Note:
5303: Both `ns` and `Y` can be `NULL`.
5305: .seealso: [](ch_ts), `TS`, `TSCreate()`
5306: @*/
5307: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5308: {
5309: PetscFunctionBegin;
5311: if (ns) PetscAssertPointer(ns, 2);
5312: if (Y) PetscAssertPointer(Y, 3);
5313: if (!ts->ops->getstages) {
5314: if (ns) *ns = 0;
5315: if (Y) *Y = NULL;
5316: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5317: PetscFunctionReturn(PETSC_SUCCESS);
5318: }
5320: /*@C
5321: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5323: Collective
5325: Input Parameters:
5326: + ts - the `TS` context
5327: . t - current timestep
5328: . U - state vector
5329: . Udot - time derivative of state vector
5330: . shift - shift to apply, see note below
5331: - ctx - an optional user context
5333: Output Parameters:
5334: + J - Jacobian matrix (not altered in this routine)
5335: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
5337: Level: intermediate
5339: Notes:
5340: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5342: dF/dU + shift*dF/dUdot
5344: Most users should not need to explicitly call this routine, as it
5345: is used internally within the nonlinear solvers.
5347: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5348: routine, then it will try to get the coloring from the matrix. This requires that the
5349: matrix have nonzero entries precomputed.
5351: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5352: @*/
5353: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5354: {
5355: SNES snes;
5356: MatFDColoring color;
5357: PetscBool hascolor, matcolor = PETSC_FALSE;
5359: PetscFunctionBegin;
5360: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5361: PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5362: if (!color) {
5363: DM dm;
5364: ISColoring iscoloring;
5366: PetscCall(TSGetDM(ts, &dm));
5367: PetscCall(DMHasColoring(dm, &hascolor));
5368: if (hascolor && !matcolor) {
5369: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5370: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5371: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5372: PetscCall(MatFDColoringSetFromOptions(color));
5373: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5374: PetscCall(ISColoringDestroy(&iscoloring));
5375: } else {
5376: MatColoring mc;
5378: PetscCall(MatColoringCreate(B, &mc));
5379: PetscCall(MatColoringSetDistance(mc, 2));
5380: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5381: PetscCall(MatColoringSetFromOptions(mc));
5382: PetscCall(MatColoringApply(mc, &iscoloring));
5383: PetscCall(MatColoringDestroy(&mc));
5384: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5385: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5386: PetscCall(MatFDColoringSetFromOptions(color));
5387: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5388: PetscCall(ISColoringDestroy(&iscoloring));
5389: }
5390: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5391: PetscCall(PetscObjectDereference((PetscObject)color));
5392: }
5393: PetscCall(TSGetSNES(ts, &snes));
5394: PetscCall(MatFDColoringApply(B, color, U, snes));
5395: if (J != B) {
5396: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5397: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5398: }
5399: PetscFunctionReturn(PETSC_SUCCESS);
5400: }
5402: /*@C
5403: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5405: Input Parameters:
5406: + ts - the `TS` context
5407: - func - function called within `TSFunctionDomainError()`
5409: Calling sequence of `func`:
5410: + ts - the `TS` context
5411: . time - the current time (of the stage)
5412: . state - the state to check if it is valid
5413: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable
5415: Level: intermediate
5417: Notes:
5418: If an implicit ODE solver is being used then, in addition to providing this routine, the
5419: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5420: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5421: Use `TSGetSNES()` to obtain the `SNES` object
5423: Developer Notes:
5424: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5425: since one takes a function pointer and the other does not.
5427: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5428: @*/
5429: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5430: {
5431: PetscFunctionBegin;
5433: ts->functiondomainerror = func;
5434: PetscFunctionReturn(PETSC_SUCCESS);
5435: }
5437: /*@
5438: TSFunctionDomainError - Checks if the current state is valid
5440: Input Parameters:
5441: + ts - the `TS` context
5442: . stagetime - time of the simulation
5443: - Y - state vector to check.
5445: Output Parameter:
5446: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5448: Level: developer
5450: Note:
5451: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5452: to check if the current state is valid.
5454: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5455: @*/
5456: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5457: {
5458: PetscFunctionBegin;
5460: *accept = PETSC_TRUE;
5461: if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5462: PetscFunctionReturn(PETSC_SUCCESS);
5463: }
5465: /*@C
5466: TSClone - This function clones a time step `TS` object.
5468: Collective
5470: Input Parameter:
5471: . tsin - The input `TS`
5473: Output Parameter:
5474: . tsout - The output `TS` (cloned)
5476: Level: developer
5478: Notes:
5479: 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.
5480: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5482: 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
5483: .vb
5484: SNES snes_dup = NULL;
5485: TSGetSNES(ts,&snes_dup);
5486: TSSetSNES(ts,snes_dup);
5487: .ve
5489: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5490: @*/
5491: PetscErrorCode TSClone(TS tsin, TS *tsout)
5492: {
5493: TS t;
5494: SNES snes_start;
5495: DM dm;
5496: TSType type;
5498: PetscFunctionBegin;
5499: PetscAssertPointer(tsin, 1);
5500: *tsout = NULL;
5502: PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));
5504: /* General TS description */
5505: t->numbermonitors = 0;
5506: t->monitorFrequency = 1;
5507: t->setupcalled = 0;
5508: t->ksp_its = 0;
5509: t->snes_its = 0;
5510: t->nwork = 0;
5511: t->rhsjacobian.time = PETSC_MIN_REAL;
5512: t->rhsjacobian.scale = 1.;
5513: t->ijacobian.shift = 1.;
5515: PetscCall(TSGetSNES(tsin, &snes_start));
5516: PetscCall(TSSetSNES(t, snes_start));
5518: PetscCall(TSGetDM(tsin, &dm));
5519: PetscCall(TSSetDM(t, dm));
5521: t->adapt = tsin->adapt;
5522: PetscCall(PetscObjectReference((PetscObject)t->adapt));
5524: t->trajectory = tsin->trajectory;
5525: PetscCall(PetscObjectReference((PetscObject)t->trajectory));
5527: t->event = tsin->event;
5528: if (t->event) t->event->refct++;
5530: t->problem_type = tsin->problem_type;
5531: t->ptime = tsin->ptime;
5532: t->ptime_prev = tsin->ptime_prev;
5533: t->time_step = tsin->time_step;
5534: t->max_time = tsin->max_time;
5535: t->steps = tsin->steps;
5536: t->max_steps = tsin->max_steps;
5537: t->equation_type = tsin->equation_type;
5538: t->atol = tsin->atol;
5539: t->rtol = tsin->rtol;
5540: t->max_snes_failures = tsin->max_snes_failures;
5541: t->max_reject = tsin->max_reject;
5542: t->errorifstepfailed = tsin->errorifstepfailed;
5544: PetscCall(TSGetType(tsin, &type));
5545: PetscCall(TSSetType(t, type));
5547: t->vec_sol = NULL;
5549: t->cfltime = tsin->cfltime;
5550: t->cfltime_local = tsin->cfltime_local;
5551: t->exact_final_time = tsin->exact_final_time;
5553: t->ops[0] = tsin->ops[0];
5555: if (((PetscObject)tsin)->fortran_func_pointers) {
5556: PetscInt i;
5557: PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers));
5558: for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5559: }
5560: *tsout = t;
5561: PetscFunctionReturn(PETSC_SUCCESS);
5562: }
5564: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5565: {
5566: TS ts = (TS)ctx;
5568: PetscFunctionBegin;
5569: PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5570: PetscFunctionReturn(PETSC_SUCCESS);
5571: }
5573: /*@
5574: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5576: Logically Collective
5578: Input Parameter:
5579: . ts - the time stepping routine
5581: Output Parameter:
5582: . flg - `PETSC_TRUE` if the multiply is likely correct
5584: Options Database Key:
5585: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5587: Level: advanced
5589: Note:
5590: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5592: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5593: @*/
5594: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5595: {
5596: Mat J, B;
5597: TSRHSJacobianFn *func;
5598: void *ctx;
5600: PetscFunctionBegin;
5601: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5602: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5603: PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5604: PetscFunctionReturn(PETSC_SUCCESS);
5605: }
5607: /*@C
5608: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5610: Logically Collective
5612: Input Parameter:
5613: . ts - the time stepping routine
5615: Output Parameter:
5616: . flg - `PETSC_TRUE` if the multiply is likely correct
5618: Options Database Key:
5619: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5621: Level: advanced
5623: Notes:
5624: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5626: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5627: @*/
5628: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5629: {
5630: Mat J, B;
5631: void *ctx;
5632: TSRHSJacobianFn *func;
5634: PetscFunctionBegin;
5635: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5636: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5637: PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5638: PetscFunctionReturn(PETSC_SUCCESS);
5639: }
5641: /*@
5642: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
5644: Logically Collective
5646: Input Parameters:
5647: + ts - timestepping context
5648: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5650: Options Database Key:
5651: . -ts_use_splitrhsfunction - <true,false>
5653: Level: intermediate
5655: Note:
5656: This is only for multirate methods
5658: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5659: @*/
5660: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5661: {
5662: PetscFunctionBegin;
5664: ts->use_splitrhsfunction = use_splitrhsfunction;
5665: PetscFunctionReturn(PETSC_SUCCESS);
5666: }
5668: /*@
5669: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
5671: Not Collective
5673: Input Parameter:
5674: . ts - timestepping context
5676: Output Parameter:
5677: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5679: Level: intermediate
5681: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5682: @*/
5683: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5684: {
5685: PetscFunctionBegin;
5687: *use_splitrhsfunction = ts->use_splitrhsfunction;
5688: PetscFunctionReturn(PETSC_SUCCESS);
5689: }
5691: /*@
5692: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
5694: Logically Collective
5696: Input Parameters:
5697: + ts - the time-stepper
5698: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
5700: Level: intermediate
5702: Note:
5703: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
5705: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
5706: @*/
5707: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5708: {
5709: PetscFunctionBegin;
5711: ts->axpy_pattern = str;
5712: PetscFunctionReturn(PETSC_SUCCESS);
5713: }
5715: /*@
5716: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
5718: Collective
5720: Input Parameters:
5721: + ts - the time-stepper
5722: . n - number of the time points (>=2)
5723: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5725: Options Database Key:
5726: . -ts_time_span <t0,...tf> - Sets the time span
5728: Level: intermediate
5730: Notes:
5731: The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
5732: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
5733: The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may
5734: pressure the memory system when using a large number of span points.
5736: .seealso: [](ch_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()`
5737: @*/
5738: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
5739: {
5740: PetscFunctionBegin;
5742: PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
5743: if (ts->tspan && n != ts->tspan->num_span_times) {
5744: PetscCall(PetscFree(ts->tspan->span_times));
5745: PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
5746: PetscCall(PetscMalloc1(n, &ts->tspan->span_times));
5747: }
5748: if (!ts->tspan) {
5749: TSTimeSpan tspan;
5750: PetscCall(PetscNew(&tspan));
5751: PetscCall(PetscMalloc1(n, &tspan->span_times));
5752: tspan->reltol = 1e-6;
5753: tspan->abstol = 10 * PETSC_MACHINE_EPSILON;
5754: tspan->worktol = 0;
5755: ts->tspan = tspan;
5756: }
5757: ts->tspan->num_span_times = n;
5758: PetscCall(PetscArraycpy(ts->tspan->span_times, span_times, n));
5759: PetscCall(TSSetTime(ts, ts->tspan->span_times[0]));
5760: PetscCall(TSSetMaxTime(ts, ts->tspan->span_times[n - 1]));
5761: PetscFunctionReturn(PETSC_SUCCESS);
5762: }
5764: /*@C
5765: TSGetTimeSpan - gets the time span set with `TSSetTimeSpan()`
5767: Not Collective
5769: Input Parameter:
5770: . ts - the time-stepper
5772: Output Parameters:
5773: + n - number of the time points (>=2)
5774: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5776: Level: beginner
5778: Note:
5779: The values obtained are valid until the `TS` object is destroyed.
5781: Both `n` and `span_times` can be `NULL`.
5783: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()`
5784: @*/
5785: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal **span_times)
5786: {
5787: PetscFunctionBegin;
5789: if (n) PetscAssertPointer(n, 2);
5790: if (span_times) PetscAssertPointer(span_times, 3);
5791: if (!ts->tspan) {
5792: if (n) *n = 0;
5793: if (span_times) *span_times = NULL;
5794: } else {
5795: if (n) *n = ts->tspan->num_span_times;
5796: if (span_times) *span_times = ts->tspan->span_times;
5797: }
5798: PetscFunctionReturn(PETSC_SUCCESS);
5799: }
5801: /*@
5802: TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.
5804: Input Parameter:
5805: . ts - the `TS` context obtained from `TSCreate()`
5807: Output Parameters:
5808: + nsol - the number of solutions
5809: - Sols - the solution vectors
5811: Level: intermediate
5813: Notes:
5814: Both `nsol` and `Sols` can be `NULL`.
5816: 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()`.
5817: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span.
5819: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`
5820: @*/
5821: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
5822: {
5823: PetscFunctionBegin;
5825: if (nsol) PetscAssertPointer(nsol, 2);
5826: if (Sols) PetscAssertPointer(Sols, 3);
5827: if (!ts->tspan) {
5828: if (nsol) *nsol = 0;
5829: if (Sols) *Sols = NULL;
5830: } else {
5831: if (nsol) *nsol = ts->tspan->spanctr;
5832: if (Sols) *Sols = ts->tspan->vecs_sol;
5833: }
5834: PetscFunctionReturn(PETSC_SUCCESS);
5835: }
5837: /*@C
5838: TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
5840: Collective
5842: Input Parameters:
5843: + ts - the `TS` context
5844: . J - Jacobian matrix (not altered in this routine)
5845: - B - newly computed Jacobian matrix to use with preconditioner
5847: Level: intermediate
5849: Notes:
5850: This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
5851: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
5852: and multiple fields are involved.
5854: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
5855: structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
5856: usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
5857: `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.
5859: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5860: @*/
5861: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
5862: {
5863: MatColoring mc = NULL;
5864: ISColoring iscoloring = NULL;
5865: MatFDColoring matfdcoloring = NULL;
5867: PetscFunctionBegin;
5868: /* Generate new coloring after eliminating zeros in the matrix */
5869: PetscCall(MatEliminateZeros(B, PETSC_TRUE));
5870: PetscCall(MatColoringCreate(B, &mc));
5871: PetscCall(MatColoringSetDistance(mc, 2));
5872: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5873: PetscCall(MatColoringSetFromOptions(mc));
5874: PetscCall(MatColoringApply(mc, &iscoloring));
5875: PetscCall(MatColoringDestroy(&mc));
5876: /* Replace the old coloring with the new one */
5877: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
5878: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5879: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
5880: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
5881: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
5882: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
5883: PetscCall(ISColoringDestroy(&iscoloring));
5884: PetscFunctionReturn(PETSC_SUCCESS);
5885: }