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: #define SkipSmallValue(a, b, tol) \
11: if (PetscAbsScalar(a) < tol || PetscAbsScalar(b) < tol) continue;
13: /* Logging support */
14: PetscClassId TS_CLASSID, DMTS_CLASSID;
15: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;
17: const char *const TSExactFinalTimeOptions[] = {"UNSPECIFIED", "STEPOVER", "INTERPOLATE", "MATCHSTEP", "TSExactFinalTimeOption", "TS_EXACTFINALTIME_", NULL};
19: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt, TSAdaptType default_type)
20: {
21: PetscFunctionBegin;
24: if (!((PetscObject)adapt)->type_name) PetscCall(TSAdaptSetType(adapt, default_type));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: /*@
29: TSSetFromOptions - Sets various `TS` parameters from the options database
31: Collective
33: Input Parameter:
34: . ts - the `TS` context obtained from `TSCreate()`
36: Options Database Keys:
37: + -ts_type <type> - EULER, BEULER, SUNDIALS, PSEUDO, CN, RK, THETA, ALPHA, GLLE, SSP, GLEE, BSYMP, IRK, see `TSType`
38: . -ts_save_trajectory - checkpoint the solution at each time-step
39: . -ts_max_time <time> - maximum time to compute to
40: . -ts_time_span <t0,...tf> - sets the time span, solutions are computed and stored for each indicated time
41: . -ts_max_steps <steps> - maximum number of time-steps to take
42: . -ts_init_time <time> - initial time to start computation
43: . -ts_final_time <time> - final time to compute to (deprecated: use `-ts_max_time`)
44: . -ts_dt <dt> - initial time step
45: . -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
46: . -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
47: . -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
48: . -ts_error_if_step_fails <true,false> - Error if no step succeeds
49: . -ts_rtol <rtol> - relative tolerance for local truncation error
50: . -ts_atol <atol> - Absolute tolerance for local truncation error
51: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function
52: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - test the Jacobian at each iteration against finite difference with RHS function
53: . -ts_adjoint_solve <yes,no> - After solving the ODE/DAE solve the adjoint problem (requires `-ts_save_trajectory`)
54: . -ts_fd_color - Use finite differences with coloring to compute IJacobian
55: . -ts_monitor - print information at each timestep
56: . -ts_monitor_cancel - Cancel all monitors
57: . -ts_monitor_lg_solution - Monitor solution graphically
58: . -ts_monitor_lg_error - Monitor error graphically
59: . -ts_monitor_error - Monitors norm of error
60: . -ts_monitor_lg_timestep - Monitor timestep size graphically
61: . -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
62: . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
63: . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
64: . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
65: . -ts_monitor_draw_solution - Monitor solution graphically
66: . -ts_monitor_draw_solution_phase <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
67: . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
68: . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
69: -ts_monitor_solution_interval <interval> - output once every interval (default=1) time steps; used with -ts_monitor_solution
70: . -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)
71: - -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
73: Level: beginner
75: Notes:
76: See `SNESSetFromOptions()` and `KSPSetFromOptions()` for how to control the nonlinear and linear solves used by the time-stepper.
78: Certain `SNES` options get reset for each new nonlinear solver, for example `-snes_lag_jacobian its` and `-snes_lag_preconditioner its`, in order
79: to retain them over the multiple nonlinear solves that `TS` uses you mush also provide `-snes_lag_jacobian_persists true` and
80: `-snes_lag_preconditioner_persists true`
82: Developer Note:
83: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
85: .seealso: [](ch_ts), `TS`, `TSGetType()`
86: @*/
87: PetscErrorCode TSSetFromOptions(TS ts)
88: {
89: PetscBool opt, flg, tflg;
90: char monfilename[PETSC_MAX_PATH_LEN];
91: PetscReal time_step, tspan[100];
92: PetscInt nt = PETSC_STATIC_ARRAY_LENGTH(tspan);
93: TSExactFinalTimeOption eftopt;
94: char dir[16];
95: TSIFunction ifun;
96: const char *defaultType;
97: char typeName[256];
99: PetscFunctionBegin;
102: PetscCall(TSRegisterAll());
103: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
105: PetscObjectOptionsBegin((PetscObject)ts);
106: if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
107: else defaultType = ifun ? TSBEULER : TSEULER;
108: PetscCall(PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, 256, &opt));
109: if (opt) PetscCall(TSSetType(ts, typeName));
110: else PetscCall(TSSetType(ts, defaultType));
112: /* Handle generic TS options */
113: PetscCall(PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL));
114: PetscCall(PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL));
115: PetscCall(PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", tspan, &nt, &flg));
116: if (flg) PetscCall(TSSetTimeSpan(ts, nt, tspan));
117: PetscCall(PetscOptionsInt("-ts_max_steps", "Maximum number of time steps", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL));
118: PetscCall(PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL));
119: PetscCall(PetscOptionsReal("-ts_dt", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg));
120: if (flg) PetscCall(TSSetTimeStep(ts, time_step));
121: PetscCall(PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg));
122: if (flg) PetscCall(TSSetExactFinalTime(ts, eftopt));
123: PetscCall(PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, NULL));
124: PetscCall(PetscOptionsInt("-ts_max_reject", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, NULL));
125: PetscCall(PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL));
126: PetscCall(PetscOptionsReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL));
127: PetscCall(PetscOptionsReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL));
129: PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL));
130: PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult_transpose", "Test the RHS Jacobian transpose for consistency with RHS at each solve ", "None", ts->testjacobiantranspose, &ts->testjacobiantranspose, NULL));
131: PetscCall(PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL));
132: #if defined(PETSC_HAVE_SAWS)
133: {
134: PetscBool set;
135: flg = PETSC_FALSE;
136: PetscCall(PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set));
137: if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)ts, flg));
138: }
139: #endif
141: /* Monitor options */
142: PetscCall(PetscOptionsInt("-ts_monitor_frequency", "Number of time steps between monitor output", "TSMonitorSetFrequency", ts->monitorFrequency, &ts->monitorFrequency, NULL));
143: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL));
144: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL));
145: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, NULL));
146: PetscCall(TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL));
148: PetscCall(PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg));
149: if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)ts, monfilename));
151: PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt));
152: if (opt) {
153: PetscInt howoften = 1;
154: DM dm;
155: PetscBool net;
157: PetscCall(PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL));
158: PetscCall(TSGetDM(ts, &dm));
159: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net));
160: if (net) {
161: TSMonitorLGCtxNetwork ctx;
162: PetscCall(TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx));
163: PetscCall(TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxNetworkDestroy));
164: PetscCall(PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL));
165: } else {
166: TSMonitorLGCtx ctx;
167: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
168: PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
169: }
170: }
172: PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt));
173: if (opt) {
174: TSMonitorLGCtx ctx;
175: PetscInt howoften = 1;
177: PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL));
178: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
179: PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
180: }
181: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL));
183: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt));
184: if (opt) {
185: TSMonitorLGCtx ctx;
186: PetscInt howoften = 1;
188: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
189: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
190: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
191: }
192: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt));
193: if (opt) {
194: TSMonitorLGCtx ctx;
195: PetscInt howoften = 1;
197: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
198: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
199: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
200: ctx->semilogy = PETSC_TRUE;
201: }
203: PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt));
204: if (opt) {
205: TSMonitorLGCtx ctx;
206: PetscInt howoften = 1;
208: PetscCall(PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL));
209: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
210: PetscCall(TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
211: }
212: PetscCall(PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt));
213: if (opt) {
214: TSMonitorLGCtx ctx;
215: PetscInt howoften = 1;
217: PetscCall(PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL));
218: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
219: PetscCall(TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy));
220: }
221: PetscCall(PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt));
222: if (opt) {
223: TSMonitorSPEigCtx ctx;
224: PetscInt howoften = 1;
226: PetscCall(PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL));
227: PetscCall(TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
228: PetscCall(TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscErrorCode(*)(void **))TSMonitorSPEigCtxDestroy));
229: }
230: PetscCall(PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase space from the DMSwarm", "TSMonitorSPSwarm", &opt));
231: if (opt) {
232: TSMonitorSPCtx ctx;
233: PetscInt howoften = 1, retain = 0;
234: PetscBool phase = PETSC_TRUE, create = PETSC_TRUE, multispecies = PETSC_FALSE;
236: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
237: if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
238: create = PETSC_FALSE;
239: break;
240: }
241: if (create) {
242: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase space from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL));
243: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL));
244: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL));
245: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_multi_species", "Color particles by particle species", "TSMonitorSPSwarm", multispecies, &multispecies, NULL));
246: PetscCall(TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, multispecies, &ctx));
247: PetscCall(TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorSPCtxDestroy));
248: }
249: }
250: PetscCall(PetscOptionsName("-ts_monitor_hg_swarm", "Display particle histogram from the DMSwarm", "TSMonitorHGSwarm", &opt));
251: if (opt) {
252: TSMonitorHGCtx ctx;
253: PetscInt howoften = 1, Ns = 1;
254: PetscBool velocity = PETSC_FALSE, create = PETSC_TRUE;
256: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
257: if (ts->monitor[i] == TSMonitorHGSwarmSolution) {
258: create = PETSC_FALSE;
259: break;
260: }
261: if (create) {
262: DM sw, dm;
263: PetscInt Nc, Nb;
265: PetscCall(TSGetDM(ts, &sw));
266: PetscCall(DMSwarmGetCellDM(sw, &dm));
267: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc));
268: Nb = PetscMin(20, PetscMax(10, Nc));
269: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL));
270: PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL));
271: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL));
272: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL));
273: PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx));
274: PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorHGCtxDestroy));
275: }
276: }
277: opt = PETSC_FALSE;
278: PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt));
279: if (opt) {
280: TSMonitorDrawCtx ctx;
281: PetscInt howoften = 1;
283: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL));
284: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
285: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
286: }
287: opt = PETSC_FALSE;
288: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt));
289: if (opt) {
290: TSMonitorDrawCtx ctx;
291: PetscReal bounds[4];
292: PetscInt n = 4;
293: PetscDraw draw;
294: PetscDrawAxis axis;
296: PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL));
297: PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field");
298: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx));
299: PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw));
300: PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis));
301: PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]));
302: PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2"));
303: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
304: }
305: opt = PETSC_FALSE;
306: PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt));
307: if (opt) {
308: TSMonitorDrawCtx ctx;
309: PetscInt howoften = 1;
311: PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL));
312: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
313: PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
314: }
315: opt = PETSC_FALSE;
316: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt));
317: if (opt) {
318: TSMonitorDrawCtx ctx;
319: PetscInt howoften = 1;
321: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL));
322: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
323: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy));
324: }
326: opt = PETSC_FALSE;
327: PetscCall(PetscOptionsString("-ts_monitor_solution_vtk", "Save each time step to a binary file, use filename-%%03" PetscInt_FMT ".vts", "TSMonitorSolutionVTK", NULL, monfilename, sizeof(monfilename), &flg));
328: if (flg) {
329: const char *ptr = NULL, *ptr2 = NULL;
330: char *filetemplate;
331: PetscCheck(monfilename[0], PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03" PetscInt_FMT ".vts");
332: /* Do some cursory validation of the input. */
333: PetscCall(PetscStrstr(monfilename, "%", (char **)&ptr));
334: PetscCheck(ptr, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03" PetscInt_FMT ".vts");
335: for (ptr++; ptr && *ptr; ptr++) {
336: PetscCall(PetscStrchr("DdiouxX", *ptr, (char **)&ptr2));
337: 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");
338: if (ptr2) break;
339: }
340: PetscCall(PetscStrallocpy(monfilename, &filetemplate));
341: PetscCall(TSMonitorSet(ts, TSMonitorSolutionVTK, filetemplate, (PetscErrorCode(*)(void **))TSMonitorSolutionVTKDestroy));
342: }
344: PetscCall(PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg));
345: if (flg) {
346: TSMonitorDMDARayCtx *rayctx;
347: int ray = 0;
348: DMDirection ddir;
349: DM da;
350: PetscMPIInt rank;
352: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
353: if (dir[0] == 'x') ddir = DM_X;
354: else if (dir[0] == 'y') ddir = DM_Y;
355: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
356: sscanf(dir + 2, "%d", &ray);
358: PetscCall(PetscInfo(((PetscObject)ts), "Displaying DMDA ray %c = %d\n", dir[0], ray));
359: PetscCall(PetscNew(&rayctx));
360: PetscCall(TSGetDM(ts, &da));
361: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
362: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank));
363: if (rank == 0) PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer));
364: rayctx->lgctx = NULL;
365: PetscCall(TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy));
366: }
367: PetscCall(PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg));
368: if (flg) {
369: TSMonitorDMDARayCtx *rayctx;
370: int ray = 0;
371: DMDirection ddir;
372: DM da;
373: PetscInt howoften = 1;
375: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
376: if (dir[0] == 'x') ddir = DM_X;
377: else if (dir[0] == 'y') ddir = DM_Y;
378: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
379: sscanf(dir + 2, "%d", &ray);
381: PetscCall(PetscInfo(((PetscObject)ts), "Displaying LG DMDA ray %c = %d\n", dir[0], ray));
382: PetscCall(PetscNew(&rayctx));
383: PetscCall(TSGetDM(ts, &da));
384: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
385: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx));
386: PetscCall(TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy));
387: }
389: PetscCall(PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt));
390: if (opt) {
391: TSMonitorEnvelopeCtx ctx;
393: PetscCall(TSMonitorEnvelopeCtxCreate(ts, &ctx));
394: PetscCall(TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscErrorCode(*)(void **))TSMonitorEnvelopeCtxDestroy));
395: }
396: flg = PETSC_FALSE;
397: PetscCall(PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt));
398: if (opt && flg) PetscCall(TSMonitorCancel(ts));
400: flg = PETSC_FALSE;
401: PetscCall(PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeIJacobianDefaultColor", flg, &flg, NULL));
402: if (flg) {
403: DM dm;
405: PetscCall(TSGetDM(ts, &dm));
406: PetscCall(DMTSUnsetIJacobianContext_Internal(dm));
407: PetscCall(TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL));
408: PetscCall(PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n"));
409: }
411: /* Handle specific TS options */
412: PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);
414: /* Handle TSAdapt options */
415: PetscCall(TSGetAdapt(ts, &ts->adapt));
416: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
417: PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject));
419: /* TS trajectory must be set after TS, since it may use some TS options above */
420: tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
421: PetscCall(PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL));
422: if (tflg) PetscCall(TSSetSaveTrajectory(ts));
424: PetscCall(TSAdjointSetFromOptions(ts, PetscOptionsObject));
426: /* process any options handlers added with PetscObjectAddOptionsHandler() */
427: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject));
428: PetscOptionsEnd();
430: if (ts->trajectory) PetscCall(TSTrajectorySetFromOptions(ts->trajectory, ts));
432: /* why do we have to do this here and not during TSSetUp? */
433: PetscCall(TSGetSNES(ts, &ts->snes));
434: if (ts->problem_type == TS_LINEAR) {
435: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
436: if (!flg) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
437: }
438: PetscCall(SNESSetFromOptions(ts->snes));
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: /*@
443: TSGetTrajectory - Gets the trajectory from a `TS` if it exists
445: Collective
447: Input Parameter:
448: . ts - the `TS` context obtained from `TSCreate()`
450: Output Parameter:
451: . tr - the `TSTrajectory` object, if it exists
453: Level: advanced
455: Note:
456: This routine should be called after all `TS` options have been set
458: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectory`, `TSTrajectoryCreate()`
459: @*/
460: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
461: {
462: PetscFunctionBegin;
464: *tr = ts->trajectory;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@
469: TSSetSaveTrajectory - Causes the `TS` to save its solutions as it iterates forward in time in a `TSTrajectory` object
471: Collective
473: Input Parameter:
474: . ts - the `TS` context obtained from `TSCreate()`
476: Options Database Keys:
477: + -ts_save_trajectory - saves the trajectory to a file
478: - -ts_trajectory_type type - set trajectory type
480: Level: intermediate
482: Notes:
483: This routine should be called after all `TS` options have been set
485: The `TSTRAJECTORYVISUALIZATION` files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
486: MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m
488: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`
489: @*/
490: PetscErrorCode TSSetSaveTrajectory(TS ts)
491: {
492: PetscFunctionBegin;
494: if (!ts->trajectory) PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /*@
499: TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object
501: Collective
503: Input Parameter:
504: . ts - the `TS` context obtained from `TSCreate()`
506: Level: intermediate
508: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
509: @*/
510: PetscErrorCode TSResetTrajectory(TS ts)
511: {
512: PetscFunctionBegin;
514: if (ts->trajectory) {
515: PetscCall(TSTrajectoryDestroy(&ts->trajectory));
516: PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
517: }
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: /*@
522: TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from a `TS`
524: Collective
526: Input Parameter:
527: . ts - the `TS` context obtained from `TSCreate()`
529: Level: intermediate
531: .seealso: [](ch_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
532: @*/
533: PetscErrorCode TSRemoveTrajectory(TS ts)
534: {
535: PetscFunctionBegin;
537: if (ts->trajectory) PetscCall(TSTrajectoryDestroy(&ts->trajectory));
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: /*@
542: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
543: set with `TSSetRHSJacobian()`.
545: Collective
547: Input Parameters:
548: + ts - the `TS` context
549: . t - current timestep
550: - U - input vector
552: Output Parameters:
553: + A - Jacobian matrix
554: - B - optional preconditioning matrix
556: Level: developer
558: Note:
559: Most users should not need to explicitly call this routine, as it
560: is used internally within the nonlinear solvers.
562: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()`
563: @*/
564: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
565: {
566: PetscObjectState Ustate;
567: PetscObjectId Uid;
568: DM dm;
569: DMTS tsdm;
570: TSRHSJacobian rhsjacobianfunc;
571: void *ctx;
572: TSRHSFunction rhsfunction;
574: PetscFunctionBegin;
577: PetscCheckSameComm(ts, 1, U, 3);
578: PetscCall(TSGetDM(ts, &dm));
579: PetscCall(DMGetDMTS(dm, &tsdm));
580: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
581: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx));
582: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
583: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
585: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) PetscFunctionReturn(PETSC_SUCCESS);
587: 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);
588: if (rhsjacobianfunc) {
589: PetscCall(PetscLogEventBegin(TS_JacobianEval, ts, U, A, B));
590: PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
591: ts->rhsjacs++;
592: PetscCall(PetscLogEventEnd(TS_JacobianEval, ts, U, A, B));
593: } else {
594: PetscCall(MatZeroEntries(A));
595: if (B && A != B) PetscCall(MatZeroEntries(B));
596: }
597: ts->rhsjacobian.time = t;
598: ts->rhsjacobian.shift = 0;
599: ts->rhsjacobian.scale = 1.;
600: PetscCall(PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid));
601: PetscCall(PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: /*@
606: TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS`
608: Collective
610: Input Parameters:
611: + ts - the `TS` context
612: . t - current time
613: - U - state vector
615: Output Parameter:
616: . y - right hand side
618: Level: developer
620: Note:
621: Most users should not need to explicitly call this routine, as it
622: is used internally within the nonlinear solvers.
624: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
625: @*/
626: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
627: {
628: TSRHSFunction rhsfunction;
629: TSIFunction ifunction;
630: void *ctx;
631: DM dm;
633: PetscFunctionBegin;
637: PetscCall(TSGetDM(ts, &dm));
638: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx));
639: PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
641: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
643: if (rhsfunction) {
644: PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, y, 0));
645: PetscCall(VecLockReadPush(U));
646: PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
647: PetscCall(VecLockReadPop(U));
648: ts->rhsfuncs++;
649: PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, y, 0));
650: } else PetscCall(VecZeroEntries(y));
651: PetscFunctionReturn(PETSC_SUCCESS);
652: }
654: /*@
655: TSComputeSolutionFunction - Evaluates the solution function.
657: Collective
659: Input Parameters:
660: + ts - the `TS` context
661: - t - current time
663: Output Parameter:
664: . U - the solution
666: Level: developer
668: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
669: @*/
670: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
671: {
672: TSSolutionFunction solutionfunction;
673: void *ctx;
674: DM dm;
676: PetscFunctionBegin;
679: PetscCall(TSGetDM(ts, &dm));
680: PetscCall(DMTSGetSolutionFunction(dm, &solutionfunction, &ctx));
681: if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
684: /*@
685: TSComputeForcingFunction - Evaluates the forcing function.
687: Collective
689: Input Parameters:
690: + ts - the `TS` context
691: - t - current time
693: Output Parameter:
694: . U - the function value
696: Level: developer
698: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
699: @*/
700: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
701: {
702: void *ctx;
703: DM dm;
704: TSForcingFunction forcing;
706: PetscFunctionBegin;
709: PetscCall(TSGetDM(ts, &dm));
710: PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));
712: if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: static PetscErrorCode TSGetRHSVec_Private(TS ts, Vec *Frhs)
717: {
718: Vec F;
720: PetscFunctionBegin;
721: *Frhs = NULL;
722: PetscCall(TSGetIFunction(ts, &F, NULL, NULL));
723: if (!ts->Frhs) PetscCall(VecDuplicate(F, &ts->Frhs));
724: *Frhs = ts->Frhs;
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
729: {
730: Mat A, B;
731: TSIJacobian ijacobian;
733: PetscFunctionBegin;
734: if (Arhs) *Arhs = NULL;
735: if (Brhs) *Brhs = NULL;
736: PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL));
737: if (Arhs) {
738: if (!ts->Arhs) {
739: if (ijacobian) {
740: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs));
741: PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN));
742: } else {
743: ts->Arhs = A;
744: PetscCall(PetscObjectReference((PetscObject)A));
745: }
746: } else {
747: PetscBool flg;
748: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
749: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
750: if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
751: PetscCall(PetscObjectDereference((PetscObject)ts->Arhs));
752: ts->Arhs = A;
753: PetscCall(PetscObjectReference((PetscObject)A));
754: }
755: }
756: *Arhs = ts->Arhs;
757: }
758: if (Brhs) {
759: if (!ts->Brhs) {
760: if (A != B) {
761: if (ijacobian) {
762: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs));
763: } else {
764: ts->Brhs = B;
765: PetscCall(PetscObjectReference((PetscObject)B));
766: }
767: } else {
768: PetscCall(PetscObjectReference((PetscObject)ts->Arhs));
769: ts->Brhs = ts->Arhs;
770: }
771: }
772: *Brhs = ts->Brhs;
773: }
774: PetscFunctionReturn(PETSC_SUCCESS);
775: }
777: /*@
778: TSComputeIFunction - Evaluates the DAE residual written in the implicit form F(t,U,Udot)=0
780: Collective
782: Input Parameters:
783: + ts - the `TS` context
784: . t - current time
785: . U - state vector
786: . Udot - time derivative of state vector
787: - imex - flag indicates if the method is `TSIMEX` so that the RHSFunction should be kept separate
789: Output Parameter:
790: . Y - right hand side
792: Level: developer
794: Note:
795: Most users should not need to explicitly call this routine, as it
796: is used internally within the nonlinear solvers.
798: If the user did did not write their equations in implicit form, this
799: function recasts them in implicit form.
801: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
802: @*/
803: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
804: {
805: TSIFunction ifunction;
806: TSRHSFunction rhsfunction;
807: void *ctx;
808: DM dm;
810: PetscFunctionBegin;
816: PetscCall(TSGetDM(ts, &dm));
817: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
818: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
820: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
822: PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, Udot, Y));
823: if (ifunction) {
824: PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
825: ts->ifuncs++;
826: }
827: if (imex) {
828: if (!ifunction) PetscCall(VecCopy(Udot, Y));
829: } else if (rhsfunction) {
830: if (ifunction) {
831: Vec Frhs;
832: PetscCall(TSGetRHSVec_Private(ts, &Frhs));
833: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
834: PetscCall(VecAXPY(Y, -1, Frhs));
835: } else {
836: PetscCall(TSComputeRHSFunction(ts, t, U, Y));
837: PetscCall(VecAYPX(Y, -1, Udot));
838: }
839: }
840: PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, Udot, Y));
841: PetscFunctionReturn(PETSC_SUCCESS);
842: }
844: /*
845: TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call `TSComputeRHSJacobian()` on it.
847: Note:
848: This routine is needed when one switches from `TSComputeIJacobian()` to `TSComputeRHSJacobian()` because the Jacobian matrix may be shifted or scaled in `TSComputeIJacobian()`.
850: */
851: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
852: {
853: PetscFunctionBegin;
855: PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
856: PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");
858: if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
859: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
860: if (B && B == ts->Brhs && A != B) {
861: if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
862: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
863: }
864: ts->rhsjacobian.shift = 0;
865: ts->rhsjacobian.scale = 1.;
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: /*@
870: TSComputeIJacobian - Evaluates the Jacobian of the DAE
872: Collective
874: Input
875: Input Parameters:
876: + ts - the `TS` context
877: . t - current timestep
878: . U - state vector
879: . Udot - time derivative of state vector
880: . shift - shift to apply, see note below
881: - imex - flag indicates if the method is `TSIMEX` so that the RHSJacobian should be kept separate
883: Output Parameters:
884: + A - Jacobian matrix
885: - B - matrix from which the preconditioner is constructed; often the same as `A`
887: Level: developer
889: Notes:
890: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
891: .vb
892: dF/dU + shift*dF/dUdot
893: .ve
894: Most users should not need to explicitly call this routine, as it
895: is used internally within the nonlinear solvers.
897: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
898: @*/
899: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
900: {
901: TSIJacobian ijacobian;
902: TSRHSJacobian rhsjacobian;
903: DM dm;
904: void *ctx;
906: PetscFunctionBegin;
915: PetscCall(TSGetDM(ts, &dm));
916: PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
917: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
919: PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
921: PetscCall(PetscLogEventBegin(TS_JacobianEval, ts, U, A, B));
922: if (ijacobian) {
923: PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
924: ts->ijacs++;
925: }
926: if (imex) {
927: if (!ijacobian) { /* system was written as Udot = G(t,U) */
928: PetscBool assembled;
929: if (rhsjacobian) {
930: Mat Arhs = NULL;
931: PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
932: if (A == Arhs) {
933: 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 */
934: ts->rhsjacobian.time = PETSC_MIN_REAL;
935: }
936: }
937: PetscCall(MatZeroEntries(A));
938: PetscCall(MatAssembled(A, &assembled));
939: if (!assembled) {
940: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
941: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
942: }
943: PetscCall(MatShift(A, shift));
944: if (A != B) {
945: PetscCall(MatZeroEntries(B));
946: PetscCall(MatAssembled(B, &assembled));
947: if (!assembled) {
948: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
949: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
950: }
951: PetscCall(MatShift(B, shift));
952: }
953: }
954: } else {
955: Mat Arhs = NULL, Brhs = NULL;
957: /* RHSJacobian needs to be converted to part of IJacobian if exists */
958: if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
959: if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
960: PetscObjectState Ustate;
961: PetscObjectId Uid;
962: TSRHSFunction rhsfunction;
964: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
965: PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
966: PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
967: if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
968: ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */
969: PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
970: if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
971: } else {
972: PetscBool flg;
974: if (ts->rhsjacobian.reuse) { /* Undo the damage */
975: /* MatScale has a short path for this case.
976: However, this code path is taken the first time TSComputeRHSJacobian is called
977: and the matrices have not been assembled yet */
978: PetscCall(TSRecoverRHSJacobian(ts, A, B));
979: }
980: PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
981: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
982: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
983: if (!flg) {
984: PetscCall(MatScale(A, -1));
985: PetscCall(MatShift(A, shift));
986: }
987: if (A != B) {
988: PetscCall(MatScale(B, -1));
989: PetscCall(MatShift(B, shift));
990: }
991: }
992: ts->rhsjacobian.scale = -1;
993: ts->rhsjacobian.shift = shift;
994: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
995: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
996: PetscCall(MatZeroEntries(A));
997: PetscCall(MatShift(A, shift));
998: if (A != B) {
999: PetscCall(MatZeroEntries(B));
1000: PetscCall(MatShift(B, shift));
1001: }
1002: }
1003: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
1004: PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
1005: if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
1006: }
1007: }
1008: PetscCall(PetscLogEventEnd(TS_JacobianEval, ts, U, A, B));
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }
1012: /*@C
1013: TSSetRHSFunction - Sets the routine for evaluating the function,
1014: where U_t = G(t,u).
1016: Logically Collective
1018: Input Parameters:
1019: + ts - the `TS` context obtained from `TSCreate()`
1020: . r - vector to put the computed right hand side (or `NULL` to have it created)
1021: . f - routine for evaluating the right-hand-side function
1022: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
1024: Calling sequence of `f`:
1025: $ PetscErrorCode f(TS ts, PetscReal t, Vec u, Vec F, void *ctx)
1026: + ts - timestep context
1027: . t - current timestep
1028: . u - input vector
1029: . F - function vector
1030: - ctx - [optional] user-defined function context
1032: Level: beginner
1034: Note:
1035: You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.
1037: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1038: @*/
1039: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, PetscErrorCode (*f)(TS, PetscReal, Vec, Vec, void *), void *ctx)
1040: {
1041: SNES snes;
1042: Vec ralloc = NULL;
1043: DM dm;
1045: PetscFunctionBegin;
1049: PetscCall(TSGetDM(ts, &dm));
1050: PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1051: PetscCall(TSGetSNES(ts, &snes));
1052: if (!r && !ts->dm && ts->vec_sol) {
1053: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1054: r = ralloc;
1055: }
1056: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1057: PetscCall(VecDestroy(&ralloc));
1058: PetscFunctionReturn(PETSC_SUCCESS);
1059: }
1061: /*@C
1062: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1064: Logically Collective
1066: Input Parameters:
1067: + ts - the `TS` context obtained from `TSCreate()`
1068: . f - routine for evaluating the solution
1069: - ctx - [optional] user-defined context for private data for the
1070: function evaluation routine (may be `NULL`)
1072: Calling sequence of `f`:
1073: $ PetscErrorCode f(TS ts, PetscReal t, Vec u, void *ctx)
1074: + ts - the integrator context
1075: . t - current timestep
1076: . u - output vector
1077: - ctx - [optional] user-defined function context
1079: Options Database Keys:
1080: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1081: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
1083: Level: intermediate
1085: Notes:
1086: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1087: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1088: create closed-form solutions with non-physical forcing terms.
1090: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1092: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1093: @*/
1094: PetscErrorCode TSSetSolutionFunction(TS ts, PetscErrorCode (*f)(TS, PetscReal, Vec, void *), void *ctx)
1095: {
1096: DM dm;
1098: PetscFunctionBegin;
1100: PetscCall(TSGetDM(ts, &dm));
1101: PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1102: PetscFunctionReturn(PETSC_SUCCESS);
1103: }
1105: /*@C
1106: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1108: Logically Collective
1110: Input Parameters:
1111: + ts - the `TS` context obtained from `TSCreate()`
1112: . func - routine for evaluating the forcing function
1113: - ctx - [optional] user-defined context for private data for the
1114: function evaluation routine (may be `NULL`)
1116: Calling sequence of func:
1117: $ PetscErrorCode func(TS ts, PetscReal t, Vec f, void *ctx)
1118: + ts - the integrator context
1119: . t - current timestep
1120: . f - output vector
1121: - ctx - [optional] user-defined function context
1123: Level: intermediate
1125: Notes:
1126: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1127: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1128: definition of the problem you are solving and hence possibly introducing bugs.
1130: This replaces the ODE F(u,u_t,t) = 0 the TS is solving with F(u,u_t,t) - func(t) = 0
1132: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1133: parameters can be passed in the ctx variable.
1135: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1137: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1138: @*/
1139: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFunction func, void *ctx)
1140: {
1141: DM dm;
1143: PetscFunctionBegin;
1145: PetscCall(TSGetDM(ts, &dm));
1146: PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1147: PetscFunctionReturn(PETSC_SUCCESS);
1148: }
1150: /*@C
1151: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1152: where U_t = G(U,t), as well as the location to store the matrix.
1154: Logically Collective
1156: Input Parameters:
1157: + ts - the `TS` context obtained from `TSCreate()`
1158: . Amat - (approximate) location to store Jacobian matrix entries computed by `f`
1159: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1160: . f - the Jacobian evaluation routine
1161: - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1163: Calling sequence of `f`:
1164: $ PetscErrorCode f(TS ts, PetscReal t, Vec u, Mat A, Mat B, void *ctx)
1165: + ts - the `TS` context obtained from `TSCreate()`
1166: . t - current timestep
1167: . u - input vector
1168: . Amat - (approximate) Jacobian matrix
1169: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1170: - ctx - [optional] user-defined context for matrix evaluation routine
1172: Level: beginner
1174: Notes:
1175: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1177: The `TS` solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f()
1178: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1180: .seealso: [](ch_ts), `TS`, `SNESComputeJacobianDefaultColor()`, `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`
1181: @*/
1182: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobian f, void *ctx)
1183: {
1184: SNES snes;
1185: DM dm;
1186: TSIJacobian ijacobian;
1188: PetscFunctionBegin;
1192: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1193: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1195: PetscCall(TSGetDM(ts, &dm));
1196: PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1197: PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1198: PetscCall(TSGetSNES(ts, &snes));
1199: if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1200: if (Amat) {
1201: PetscCall(PetscObjectReference((PetscObject)Amat));
1202: PetscCall(MatDestroy(&ts->Arhs));
1203: ts->Arhs = Amat;
1204: }
1205: if (Pmat) {
1206: PetscCall(PetscObjectReference((PetscObject)Pmat));
1207: PetscCall(MatDestroy(&ts->Brhs));
1208: ts->Brhs = Pmat;
1209: }
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: /*@C
1214: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1216: Logically Collective
1218: Input Parameters:
1219: + ts - the `TS` context obtained from `TSCreate()`
1220: . r - vector to hold the residual (or `NULL` to have it created internally)
1221: . f - the function evaluation routine
1222: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1224: Calling sequence of `f`:
1225: $ PetscErrorCode f(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
1226: + ts - the `TS` context obtained from `TSCreate()`
1227: . t - time at step/stage being solved
1228: . u - state vector
1229: . u_t - time derivative of state vector
1230: . F - function vector
1231: - ctx - [optional] user-defined context for matrix evaluation routine
1233: Level: beginner
1235: Note:
1236: The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function.
1238: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`, `TSSetIJacobian()`
1239: @*/
1240: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunction f, void *ctx)
1241: {
1242: SNES snes;
1243: Vec ralloc = NULL;
1244: DM dm;
1246: PetscFunctionBegin;
1250: PetscCall(TSGetDM(ts, &dm));
1251: PetscCall(DMTSSetIFunction(dm, f, ctx));
1253: PetscCall(TSGetSNES(ts, &snes));
1254: if (!r && !ts->dm && ts->vec_sol) {
1255: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1256: r = ralloc;
1257: }
1258: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1259: PetscCall(VecDestroy(&ralloc));
1260: PetscFunctionReturn(PETSC_SUCCESS);
1261: }
1263: /*@C
1264: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.
1266: Not Collective
1268: Input Parameter:
1269: . ts - the `TS` context
1271: Output Parameters:
1272: + r - vector to hold residual (or `NULL`)
1273: . func - the function to compute residual (or `NULL`)
1274: - ctx - the function context (or `NULL`)
1276: Level: advanced
1278: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1279: @*/
1280: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunction *func, void **ctx)
1281: {
1282: SNES snes;
1283: DM dm;
1285: PetscFunctionBegin;
1287: PetscCall(TSGetSNES(ts, &snes));
1288: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1289: PetscCall(TSGetDM(ts, &dm));
1290: PetscCall(DMTSGetIFunction(dm, func, ctx));
1291: PetscFunctionReturn(PETSC_SUCCESS);
1292: }
1294: /*@C
1295: TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it.
1297: Not Collective
1299: Input Parameter:
1300: . ts - the `TS` context
1302: Output Parameters:
1303: + r - vector to hold computed right hand side (or `NULL`)
1304: . func - the function to compute right hand side (or `NULL`)
1305: - ctx - the function context (or `NULL`)
1307: Level: advanced
1309: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1310: @*/
1311: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunction *func, void **ctx)
1312: {
1313: SNES snes;
1314: DM dm;
1316: PetscFunctionBegin;
1318: PetscCall(TSGetSNES(ts, &snes));
1319: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1320: PetscCall(TSGetDM(ts, &dm));
1321: PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1322: PetscFunctionReturn(PETSC_SUCCESS);
1323: }
1325: /*@C
1326: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1327: provided with `TSSetIFunction()`.
1329: Logically Collective
1331: Input Parameters:
1332: + ts - the `TS` context obtained from `TSCreate()`
1333: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
1334: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1335: . f - the Jacobian evaluation routine
1336: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1338: Calling sequence of `f`:
1339: $ PetscErrorCode f(TS ts, PetscReal t, Vec U, Vec U_t, PetscReal a, Mat Amat, Mat Pmat, void *ctx)
1340: + ts - the `TS` context obtained from `TSCreate()`
1341: . t - time at step/stage being solved
1342: . U - state vector
1343: . U_t - time derivative of state vector
1344: . a - shift
1345: . Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1346: . Pmat - matrix used for constructing preconditioner, usually the same as `Amat`
1347: - ctx - [optional] user-defined context for matrix evaluation routine
1349: Level: beginner
1351: Notes:
1352: The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1354: If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
1355: space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.
1357: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1358: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1359: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1360: a and vector W depend on the integration method, step size, and past states. For example with
1361: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1362: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1364: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1366: The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f`
1367: You should not assume the values are the same in the next call to `f` as you set them in the previous call.
1369: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetRHSJacobian()`, `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1370: @*/
1371: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobian f, void *ctx)
1372: {
1373: SNES snes;
1374: DM dm;
1376: PetscFunctionBegin;
1380: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1381: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1383: PetscCall(TSGetDM(ts, &dm));
1384: PetscCall(DMTSSetIJacobian(dm, f, ctx));
1386: PetscCall(TSGetSNES(ts, &snes));
1387: PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1388: PetscFunctionReturn(PETSC_SUCCESS);
1389: }
1391: /*@
1392: TSRHSJacobianSetReuse - restore RHS Jacobian before re-evaluating. Without this flag, `TS` will change the sign and
1393: shift the RHS Jacobian for a finite-time-step implicit solve, in which case the user function will need to recompute
1394: the entire Jacobian. The reuse flag must be set if the evaluation function will assume that the matrix entries have
1395: not been changed by the `TS`.
1397: Logically Collective
1399: Input Parameters:
1400: + ts - `TS` context obtained from `TSCreate()`
1401: - reuse - `PETSC_TRUE` if the RHS Jacobian
1403: Level: intermediate
1405: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1406: @*/
1407: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1408: {
1409: PetscFunctionBegin;
1410: ts->rhsjacobian.reuse = reuse;
1411: PetscFunctionReturn(PETSC_SUCCESS);
1412: }
1414: /*@C
1415: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1417: Logically Collective
1419: Input Parameters:
1420: + ts - the `TS` context obtained from `TSCreate()`
1421: . F - vector to hold the residual (or `NULL` to have it created internally)
1422: . fun - the function evaluation routine
1423: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1425: Calling sequence of `fun`:
1426: $ PetscErrorCode fun(TS ts, PetscReal t, Vec U, Vec U_t, Vec U_tt, Vec F,void *ctx);
1427: + ts - the `TS` context obtained from `TSCreate()`
1428: . t - time at step/stage being solved
1429: . U - state vector
1430: . U_t - time derivative of state vector
1431: . U_tt - second time derivative of state vector
1432: . F - function vector
1433: - ctx - [optional] user-defined context for matrix evaluation routine (may be `NULL`)
1435: Level: beginner
1437: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`, `TSSetIFunction()`, `TSCreate()`, `TSSetRHSFunction()`
1438: @*/
1439: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2Function fun, void *ctx)
1440: {
1441: DM dm;
1443: PetscFunctionBegin;
1446: PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1447: PetscCall(TSGetDM(ts, &dm));
1448: PetscCall(DMTSSetI2Function(dm, fun, ctx));
1449: PetscFunctionReturn(PETSC_SUCCESS);
1450: }
1452: /*@C
1453: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.
1455: Not Collective
1457: Input Parameter:
1458: . ts - the `TS` context
1460: Output Parameters:
1461: + r - vector to hold residual (or `NULL`)
1462: . fun - the function to compute residual (or `NULL`)
1463: - ctx - the function context (or `NULL`)
1465: Level: advanced
1467: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1468: @*/
1469: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2Function *fun, void **ctx)
1470: {
1471: SNES snes;
1472: DM dm;
1474: PetscFunctionBegin;
1476: PetscCall(TSGetSNES(ts, &snes));
1477: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1478: PetscCall(TSGetDM(ts, &dm));
1479: PetscCall(DMTSGetI2Function(dm, fun, ctx));
1480: PetscFunctionReturn(PETSC_SUCCESS);
1481: }
1483: /*@C
1484: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1485: where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.
1487: Logically Collective
1489: Input Parameters:
1490: + ts - the `TS` context obtained from `TSCreate()`
1491: . J - matrix to hold the Jacobian values
1492: . P - matrix for constructing the preconditioner (may be same as `J`)
1493: . jac - the Jacobian evaluation routine
1494: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1496: Calling sequence of `jac`:
1497: $ PetscErrorCode jac(TS ts, PetscReal t, Vec U, Vec U_t, Vec U_tt, PetscReal v, PetscReal a, Mat J, Mat P, void *ctx)
1498: + ts - the `TS` context obtained from `TSCreate()`
1499: . t - time at step/stage being solved
1500: . U - state vector
1501: . U_t - time derivative of state vector
1502: . U_tt - second time derivative of state vector
1503: . v - shift for U_t
1504: . a - shift for U_tt
1505: . J - Jacobian of G(U) = F(t,U,W+v*U,W'+a*U), equivalent to dF/dU + v*dF/dU_t + a*dF/dU_tt
1506: . P - preconditioning matrix for J, may be same as J
1507: - ctx - [optional] user-defined context for matrix evaluation routine
1509: Level: beginner
1511: Notes:
1512: The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1514: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1515: 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.
1516: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1517: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1519: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1520: @*/
1521: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2Jacobian jac, void *ctx)
1522: {
1523: DM dm;
1525: PetscFunctionBegin;
1529: PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1530: PetscCall(TSGetDM(ts, &dm));
1531: PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1532: PetscFunctionReturn(PETSC_SUCCESS);
1533: }
1535: /*@C
1536: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1538: Not Collective, but parallel objects are returned if `TS` is parallel
1540: Input Parameter:
1541: . ts - The `TS` context obtained from `TSCreate()`
1543: Output Parameters:
1544: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1545: . P - The matrix from which the preconditioner is constructed, often the same as `J`
1546: . jac - The function to compute the Jacobian matrices
1547: - ctx - User-defined context for Jacobian evaluation routine
1549: Level: advanced
1551: Note:
1552: You can pass in `NULL` for any return argument you do not need.
1554: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1555: @*/
1556: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2Jacobian *jac, void **ctx)
1557: {
1558: SNES snes;
1559: DM dm;
1561: PetscFunctionBegin;
1562: PetscCall(TSGetSNES(ts, &snes));
1563: PetscCall(SNESSetUpMatrices(snes));
1564: PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1565: PetscCall(TSGetDM(ts, &dm));
1566: PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1567: PetscFunctionReturn(PETSC_SUCCESS);
1568: }
1570: /*@
1571: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1573: Collective
1575: Input Parameters:
1576: + ts - the `TS` context
1577: . t - current time
1578: . U - state vector
1579: . V - time derivative of state vector (U_t)
1580: - A - second time derivative of state vector (U_tt)
1582: Output Parameter:
1583: . F - the residual vector
1585: Level: developer
1587: Note:
1588: Most users should not need to explicitly call this routine, as it
1589: is used internally within the nonlinear solvers.
1591: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1592: @*/
1593: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1594: {
1595: DM dm;
1596: TSI2Function I2Function;
1597: void *ctx;
1598: TSRHSFunction rhsfunction;
1600: PetscFunctionBegin;
1607: PetscCall(TSGetDM(ts, &dm));
1608: PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1609: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
1611: if (!I2Function) {
1612: PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, V, F));
1618: PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));
1620: if (rhsfunction) {
1621: Vec Frhs;
1622: PetscCall(TSGetRHSVec_Private(ts, &Frhs));
1623: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1624: PetscCall(VecAXPY(F, -1, Frhs));
1625: }
1627: PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, V, F));
1628: PetscFunctionReturn(PETSC_SUCCESS);
1629: }
1631: /*@
1632: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1634: Collective
1636: Input Parameters:
1637: + ts - the `TS` context
1638: . t - current timestep
1639: . U - state vector
1640: . V - time derivative of state vector
1641: . A - second time derivative of state vector
1642: . shiftV - shift to apply, see note below
1643: - shiftA - shift to apply, see note below
1645: Output Parameters:
1646: + J - Jacobian matrix
1647: - P - optional preconditioning matrix
1649: Level: developer
1651: Notes:
1652: If F(t,U,V,A)=0 is the DAE, the required Jacobian is
1654: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1656: Most users should not need to explicitly call this routine, as it
1657: is used internally within the nonlinear solvers.
1659: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1660: @*/
1661: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1662: {
1663: DM dm;
1664: TSI2Jacobian I2Jacobian;
1665: void *ctx;
1666: TSRHSJacobian rhsjacobian;
1668: PetscFunctionBegin;
1676: PetscCall(TSGetDM(ts, &dm));
1677: PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1678: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1680: if (!I2Jacobian) {
1681: PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1682: PetscFunctionReturn(PETSC_SUCCESS);
1683: }
1685: PetscCall(PetscLogEventBegin(TS_JacobianEval, ts, U, J, P));
1686: PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1687: if (rhsjacobian) {
1688: Mat Jrhs, Prhs;
1689: PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1690: PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1691: PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1692: if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1693: }
1695: PetscCall(PetscLogEventEnd(TS_JacobianEval, ts, U, J, P));
1696: PetscFunctionReturn(PETSC_SUCCESS);
1697: }
1699: /*@C
1700: TSSetTransientVariable - sets function to transform from state to transient variables
1702: Logically Collective
1704: Input Parameters:
1705: + ts - time stepping context on which to change the transient variable
1706: . tvar - a function that transforms to transient variables
1707: - ctx - a context for tvar
1709: Calling sequence of `tvar`:
1710: $ PetscErrorCode tvar(TS ts, Vec p, Vec c, void *ctx)
1711: + ts - timestep context
1712: . p - input vector (primitive form)
1713: . c - output vector, transient variables (conservative form)
1714: - ctx - [optional] user-defined function context
1716: Level: advanced
1718: Notes:
1719: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1720: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1721: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1722: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1723: evaluated via the chain rule, as in
1724: .vb
1725: dF/dP + shift * dF/dCdot dC/dP.
1726: .ve
1728: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1729: @*/
1730: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariable tvar, void *ctx)
1731: {
1732: DM dm;
1734: PetscFunctionBegin;
1736: PetscCall(TSGetDM(ts, &dm));
1737: PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1738: PetscFunctionReturn(PETSC_SUCCESS);
1739: }
1741: /*@
1742: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1744: Logically Collective
1746: Input Parameters:
1747: + ts - TS on which to compute
1748: - U - state vector to be transformed to transient variables
1750: Output Parameter:
1751: . C - transient (conservative) variable
1753: Level: developer
1755: Developer Note:
1756: If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed.
1757: This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are
1758: being used.
1760: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1761: @*/
1762: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1763: {
1764: DM dm;
1765: DMTS dmts;
1767: PetscFunctionBegin;
1770: PetscCall(TSGetDM(ts, &dm));
1771: PetscCall(DMGetDMTS(dm, &dmts));
1772: if (dmts->ops->transientvar) {
1774: PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1775: }
1776: PetscFunctionReturn(PETSC_SUCCESS);
1777: }
1779: /*@
1780: TSHasTransientVariable - determine whether transient variables have been set
1782: Logically Collective
1784: Input Parameter:
1785: . ts - `TS` on which to compute
1787: Output Parameter:
1788: . has - `PETSC_TRUE` if transient variables have been set
1790: Level: developer
1792: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1793: @*/
1794: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1795: {
1796: DM dm;
1797: DMTS dmts;
1799: PetscFunctionBegin;
1801: PetscCall(TSGetDM(ts, &dm));
1802: PetscCall(DMGetDMTS(dm, &dmts));
1803: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1804: PetscFunctionReturn(PETSC_SUCCESS);
1805: }
1807: /*@
1808: TS2SetSolution - Sets the initial solution and time derivative vectors
1809: for use by the `TS` routines handling second order equations.
1811: Logically Collective
1813: Input Parameters:
1814: + ts - the `TS` context obtained from `TSCreate()`
1815: . u - the solution vector
1816: - v - the time derivative vector
1818: Level: beginner
1820: .seealso: [](ch_ts), `TS`
1821: @*/
1822: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1823: {
1824: PetscFunctionBegin;
1828: PetscCall(TSSetSolution(ts, u));
1829: PetscCall(PetscObjectReference((PetscObject)v));
1830: PetscCall(VecDestroy(&ts->vec_dot));
1831: ts->vec_dot = v;
1832: PetscFunctionReturn(PETSC_SUCCESS);
1833: }
1835: /*@
1836: TS2GetSolution - Returns the solution and time derivative at the present timestep
1837: for second order equations. It is valid to call this routine inside the function
1838: that you are evaluating in order to move to the new timestep. This vector not
1839: changed until the solution at the next timestep has been calculated.
1841: Not Collective, but `u` returned is parallel if `TS` is parallel
1843: Input Parameter:
1844: . ts - the `TS` context obtained from `TSCreate()`
1846: Output Parameters:
1847: + u - the vector containing the solution
1848: - v - the vector containing the time derivative
1850: Level: intermediate
1852: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1853: @*/
1854: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1855: {
1856: PetscFunctionBegin;
1860: if (u) *u = ts->vec_sol;
1861: if (v) *v = ts->vec_dot;
1862: PetscFunctionReturn(PETSC_SUCCESS);
1863: }
1865: /*@C
1866: TSLoad - Loads a `TS` that has been stored in binary with `TSView()`.
1868: Collective
1870: Input Parameters:
1871: + newdm - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1872: some related function before a call to `TSLoad()`.
1873: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1875: Level: intermediate
1877: Note:
1878: The type is determined by the data in the file, any type set into the `TS` before this call is ignored.
1880: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1881: @*/
1882: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1883: {
1884: PetscBool isbinary;
1885: PetscInt classid;
1886: char type[256];
1887: DMTS sdm;
1888: DM dm;
1890: PetscFunctionBegin;
1893: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1894: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1896: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1897: PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1898: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1899: PetscCall(TSSetType(ts, type));
1900: PetscTryTypeMethod(ts, load, viewer);
1901: PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1902: PetscCall(DMLoad(dm, viewer));
1903: PetscCall(TSSetDM(ts, dm));
1904: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1905: PetscCall(VecLoad(ts->vec_sol, viewer));
1906: PetscCall(DMGetDMTS(ts->dm, &sdm));
1907: PetscCall(DMTSLoad(sdm, viewer));
1908: PetscFunctionReturn(PETSC_SUCCESS);
1909: }
1911: #include <petscdraw.h>
1912: #if defined(PETSC_HAVE_SAWS)
1913: #include <petscviewersaws.h>
1914: #endif
1916: /*@C
1917: TSViewFromOptions - View a `TS` based on values in the options database
1919: Collective
1921: Input Parameters:
1922: + ts - the `TS` context
1923: . obj - Optional object that provides the prefix for the options database keys
1924: - name - command line option string to be passed by user
1926: Level: intermediate
1928: .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1929: @*/
1930: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1931: {
1932: PetscFunctionBegin;
1934: PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1935: PetscFunctionReturn(PETSC_SUCCESS);
1936: }
1938: /*@C
1939: TSView - Prints the `TS` data structure.
1941: Collective
1943: Input Parameters:
1944: + ts - the `TS` context obtained from `TSCreate()`
1945: - viewer - visualization context
1947: Options Database Key:
1948: . -ts_view - calls `TSView()` at end of `TSStep()`
1950: Level: beginner
1952: Notes:
1953: The available visualization contexts include
1954: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1955: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1956: output where only the first processor opens
1957: the file. All other processors send their
1958: data to the first processor to print.
1960: The user can open an alternative visualization context with
1961: `PetscViewerASCIIOpen()` - output to a specified file.
1963: In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).
1965: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1966: @*/
1967: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1968: {
1969: TSType type;
1970: PetscBool iascii, isstring, isundials, isbinary, isdraw;
1971: DMTS sdm;
1972: #if defined(PETSC_HAVE_SAWS)
1973: PetscBool issaws;
1974: #endif
1976: PetscFunctionBegin;
1978: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1980: PetscCheckSameComm(ts, 1, viewer, 2);
1982: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1983: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1984: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1985: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1986: #if defined(PETSC_HAVE_SAWS)
1987: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1988: #endif
1989: if (iascii) {
1990: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer));
1991: if (ts->ops->view) {
1992: PetscCall(PetscViewerASCIIPushTab(viewer));
1993: PetscUseTypeMethod(ts, view, viewer);
1994: PetscCall(PetscViewerASCIIPopTab(viewer));
1995: }
1996: if (ts->max_steps < PETSC_MAX_INT) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps));
1997: if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time));
1998: if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs));
1999: if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs));
2000: if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs));
2001: if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs));
2002: if (ts->usessnes) {
2003: PetscBool lin;
2004: if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its));
2005: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its));
2006: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
2007: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures));
2008: }
2009: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject));
2010: if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, "));
2011: else PetscCall(PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol));
2012: if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of absolute error tolerances\n"));
2013: else PetscCall(PetscViewerASCIIPrintf(viewer, " using absolute error tolerance of %g\n", (double)ts->atol));
2014: PetscCall(PetscViewerASCIIPushTab(viewer));
2015: PetscCall(TSAdaptView(ts->adapt, viewer));
2016: PetscCall(PetscViewerASCIIPopTab(viewer));
2017: } else if (isstring) {
2018: PetscCall(TSGetType(ts, &type));
2019: PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type));
2020: PetscTryTypeMethod(ts, view, viewer);
2021: } else if (isbinary) {
2022: PetscInt classid = TS_FILE_CLASSID;
2023: MPI_Comm comm;
2024: PetscMPIInt rank;
2025: char type[256];
2027: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
2028: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2029: if (rank == 0) {
2030: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
2031: PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
2032: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
2033: }
2034: PetscTryTypeMethod(ts, view, viewer);
2035: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
2036: PetscCall(DMView(ts->dm, viewer));
2037: PetscCall(VecView(ts->vec_sol, viewer));
2038: PetscCall(DMGetDMTS(ts->dm, &sdm));
2039: PetscCall(DMTSView(sdm, viewer));
2040: } else if (isdraw) {
2041: PetscDraw draw;
2042: char str[36];
2043: PetscReal x, y, bottom, h;
2045: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2046: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
2047: PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
2048: PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
2049: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
2050: bottom = y - h;
2051: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
2052: PetscTryTypeMethod(ts, view, viewer);
2053: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
2054: if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
2055: PetscCall(PetscDrawPopCurrentPoint(draw));
2056: #if defined(PETSC_HAVE_SAWS)
2057: } else if (issaws) {
2058: PetscMPIInt rank;
2059: const char *name;
2061: PetscCall(PetscObjectGetName((PetscObject)ts, &name));
2062: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
2063: if (!((PetscObject)ts)->amsmem && rank == 0) {
2064: char dir[1024];
2066: PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
2067: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
2068: PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
2069: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
2070: PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
2071: }
2072: PetscTryTypeMethod(ts, view, viewer);
2073: #endif
2074: }
2075: if (ts->snes && ts->usessnes) {
2076: PetscCall(PetscViewerASCIIPushTab(viewer));
2077: PetscCall(SNESView(ts->snes, viewer));
2078: PetscCall(PetscViewerASCIIPopTab(viewer));
2079: }
2080: PetscCall(DMGetDMTS(ts->dm, &sdm));
2081: PetscCall(DMTSView(sdm, viewer));
2083: PetscCall(PetscViewerASCIIPushTab(viewer));
2084: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials));
2085: PetscCall(PetscViewerASCIIPopTab(viewer));
2086: PetscFunctionReturn(PETSC_SUCCESS);
2087: }
2089: /*@
2090: TSSetApplicationContext - Sets an optional user-defined context for
2091: the timesteppers.
2093: Logically Collective
2095: Input Parameters:
2096: + ts - the `TS` context obtained from `TSCreate()`
2097: - usrP - user context
2099: Level: intermediate
2101: Fortran Note:
2102: You must write a Fortran interface definition for this
2103: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
2105: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2106: @*/
2107: PetscErrorCode TSSetApplicationContext(TS ts, void *usrP)
2108: {
2109: PetscFunctionBegin;
2111: ts->user = usrP;
2112: PetscFunctionReturn(PETSC_SUCCESS);
2113: }
2115: /*@
2116: TSGetApplicationContext - Gets the user-defined context for the
2117: timestepper that was set with `TSSetApplicationContext()`
2119: Not Collective
2121: Input Parameter:
2122: . ts - the `TS` context obtained from `TSCreate()`
2124: Output Parameter:
2125: . usrP - user context
2127: Level: intermediate
2129: Fortran Note:
2130: You must write a Fortran interface definition for this
2131: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
2133: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2134: @*/
2135: PetscErrorCode TSGetApplicationContext(TS ts, void *usrP)
2136: {
2137: PetscFunctionBegin;
2139: *(void **)usrP = ts->user;
2140: PetscFunctionReturn(PETSC_SUCCESS);
2141: }
2143: /*@
2144: TSGetStepNumber - Gets the number of time steps completed.
2146: Not Collective
2148: Input Parameter:
2149: . ts - the `TS` context obtained from `TSCreate()`
2151: Output Parameter:
2152: . steps - number of steps completed so far
2154: Level: intermediate
2156: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2157: @*/
2158: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2159: {
2160: PetscFunctionBegin;
2163: *steps = ts->steps;
2164: PetscFunctionReturn(PETSC_SUCCESS);
2165: }
2167: /*@
2168: TSSetStepNumber - Sets the number of steps completed.
2170: Logically Collective
2172: Input Parameters:
2173: + ts - the `TS` context
2174: - steps - number of steps completed so far
2176: Level: developer
2178: Note:
2179: For most uses of the `TS` solvers the user need not explicitly call
2180: `TSSetStepNumber()`, as the step counter is appropriately updated in
2181: `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2182: reinitialize timestepping by setting the step counter to zero (and time
2183: to the initial time) to solve a similar problem with different initial
2184: conditions or parameters. Other possible use case is to continue
2185: timestepping from a previously interrupted run in such a way that `TS`
2186: monitors will be called with a initial nonzero step counter.
2188: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2189: @*/
2190: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2191: {
2192: PetscFunctionBegin;
2195: PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2196: ts->steps = steps;
2197: PetscFunctionReturn(PETSC_SUCCESS);
2198: }
2200: /*@
2201: TSSetTimeStep - Allows one to reset the timestep at any time,
2202: useful for simple pseudo-timestepping codes.
2204: Logically Collective
2206: Input Parameters:
2207: + ts - the `TS` context obtained from `TSCreate()`
2208: - time_step - the size of the timestep
2210: Level: intermediate
2212: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2213: @*/
2214: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2215: {
2216: PetscFunctionBegin;
2219: ts->time_step = time_step;
2220: PetscFunctionReturn(PETSC_SUCCESS);
2221: }
2223: /*@
2224: TSSetExactFinalTime - Determines whether to adapt the final time step to
2225: match the exact final time, interpolate solution to the exact final time,
2226: or just return at the final time `TS` computed.
2228: Logically Collective
2230: Input Parameters:
2231: + ts - the time-step context
2232: - eftopt - exact final time option
2233: .vb
2234: TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded
2235: TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2236: TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2237: .ve
2239: Options Database Key:
2240: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step at runtime
2242: Level: beginner
2244: Note:
2245: If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2246: then the final time you selected.
2248: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2249: @*/
2250: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2251: {
2252: PetscFunctionBegin;
2255: ts->exact_final_time = eftopt;
2256: PetscFunctionReturn(PETSC_SUCCESS);
2257: }
2259: /*@
2260: TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`
2262: Not Collective
2264: Input Parameter:
2265: . ts - the `TS` context
2267: Output Parameter:
2268: . eftopt - exact final time option
2270: Level: beginner
2272: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2273: @*/
2274: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2275: {
2276: PetscFunctionBegin;
2279: *eftopt = ts->exact_final_time;
2280: PetscFunctionReturn(PETSC_SUCCESS);
2281: }
2283: /*@
2284: TSGetTimeStep - Gets the current timestep size.
2286: Not Collective
2288: Input Parameter:
2289: . ts - the `TS` context obtained from `TSCreate()`
2291: Output Parameter:
2292: . dt - the current timestep size
2294: Level: intermediate
2296: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2297: @*/
2298: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2299: {
2300: PetscFunctionBegin;
2303: *dt = ts->time_step;
2304: PetscFunctionReturn(PETSC_SUCCESS);
2305: }
2307: /*@
2308: TSGetSolution - Returns the solution at the present timestep. It
2309: is valid to call this routine inside the function that you are evaluating
2310: in order to move to the new timestep. This vector not changed until
2311: the solution at the next timestep has been calculated.
2313: Not Collective, but v returned is parallel if ts is parallel
2315: Input Parameter:
2316: . ts - the `TS` context obtained from `TSCreate()`
2318: Output Parameter:
2319: . v - the vector containing the solution
2321: Level: intermediate
2323: Note:
2324: If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
2325: final time. It returns the solution at the next timestep.
2327: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2328: @*/
2329: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2330: {
2331: PetscFunctionBegin;
2334: *v = ts->vec_sol;
2335: PetscFunctionReturn(PETSC_SUCCESS);
2336: }
2338: /*@
2339: TSGetSolutionComponents - Returns any solution components at the present
2340: timestep, if available for the time integration method being used.
2341: Solution components are quantities that share the same size and
2342: structure as the solution vector.
2344: Not Collective, but v returned is parallel if ts is parallel
2346: Parameters :
2347: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2348: . n - If v is `NULL`, then the number of solution components is
2349: returned through n, else the n-th solution component is
2350: returned in v.
2351: - v - the vector containing the n-th solution component
2352: (may be `NULL` to use this function to find out
2353: the number of solutions components).
2355: Level: advanced
2357: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2358: @*/
2359: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2360: {
2361: PetscFunctionBegin;
2363: if (!ts->ops->getsolutioncomponents) *n = 0;
2364: else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2365: PetscFunctionReturn(PETSC_SUCCESS);
2366: }
2368: /*@
2369: TSGetAuxSolution - Returns an auxiliary solution at the present
2370: timestep, if available for the time integration method being used.
2372: Not Collective, but v returned is parallel if ts is parallel
2374: Parameters :
2375: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2376: - v - the vector containing the auxiliary solution
2378: Level: intermediate
2380: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2381: @*/
2382: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2383: {
2384: PetscFunctionBegin;
2386: if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2387: else PetscCall(VecZeroEntries(*v));
2388: PetscFunctionReturn(PETSC_SUCCESS);
2389: }
2391: /*@
2392: TSGetTimeError - Returns the estimated error vector, if the chosen
2393: `TSType` has an error estimation functionality and `TSSetTimeError()` was called
2395: Not Collective, but v returned is parallel if ts is parallel
2397: Parameters :
2398: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2399: . n - current estimate (n=0) or previous one (n=-1)
2400: - v - the vector containing the error (same size as the solution).
2402: Level: intermediate
2404: Note:
2405: MUST call after `TSSetUp()`
2407: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2408: @*/
2409: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2410: {
2411: PetscFunctionBegin;
2413: if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2414: else PetscCall(VecZeroEntries(*v));
2415: PetscFunctionReturn(PETSC_SUCCESS);
2416: }
2418: /*@
2419: TSSetTimeError - Sets the estimated error vector, if the chosen
2420: `TSType` has an error estimation functionality. This can be used
2421: to restart such a time integrator with a given error vector.
2423: Not Collective, but v returned is parallel if ts is parallel
2425: Parameters :
2426: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2427: - v - the vector containing the error (same size as the solution).
2429: Level: intermediate
2431: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError)`
2432: @*/
2433: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2434: {
2435: PetscFunctionBegin;
2437: PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2438: PetscTryTypeMethod(ts, settimeerror, v);
2439: PetscFunctionReturn(PETSC_SUCCESS);
2440: }
2442: /* ----- Routines to initialize and destroy a timestepper ---- */
2443: /*@
2444: TSSetProblemType - Sets the type of problem to be solved.
2446: Not collective
2448: Input Parameters:
2449: + ts - The `TS`
2450: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2451: .vb
2452: U_t - A U = 0 (linear)
2453: U_t - A(t) U = 0 (linear)
2454: F(t,U,U_t) = 0 (nonlinear)
2455: .ve
2457: Level: beginner
2459: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2460: @*/
2461: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2462: {
2463: PetscFunctionBegin;
2465: ts->problem_type = type;
2466: if (type == TS_LINEAR) {
2467: SNES snes;
2468: PetscCall(TSGetSNES(ts, &snes));
2469: PetscCall(SNESSetType(snes, SNESKSPONLY));
2470: }
2471: PetscFunctionReturn(PETSC_SUCCESS);
2472: }
2474: /*@C
2475: TSGetProblemType - Gets the type of problem to be solved.
2477: Not collective
2479: Input Parameter:
2480: . ts - The `TS`
2482: Output Parameter:
2483: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2484: .vb
2485: M U_t = A U
2486: M(t) U_t = A(t) U
2487: F(t,U,U_t)
2488: .ve
2490: Level: beginner
2492: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2493: @*/
2494: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2495: {
2496: PetscFunctionBegin;
2499: *type = ts->problem_type;
2500: PetscFunctionReturn(PETSC_SUCCESS);
2501: }
2503: /*
2504: Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2505: */
2506: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2507: {
2508: PetscBool isnone;
2510: PetscFunctionBegin;
2511: PetscCall(TSGetAdapt(ts, &ts->adapt));
2512: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2514: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2515: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2516: else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2517: PetscFunctionReturn(PETSC_SUCCESS);
2518: }
2520: /*@
2521: TSSetUp - Sets up the internal data structures for the later use of a timestepper.
2523: Collective
2525: Input Parameter:
2526: . ts - the `TS` context obtained from `TSCreate()`
2528: Level: advanced
2530: Note:
2531: For basic use of the `TS` solvers the user need not explicitly call
2532: `TSSetUp()`, since these actions will automatically occur during
2533: the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this
2534: phase separately, `TSSetUp()` should be called after `TSCreate()`
2535: and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.
2537: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2538: @*/
2539: PetscErrorCode TSSetUp(TS ts)
2540: {
2541: DM dm;
2542: PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2543: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2544: TSIFunction ifun;
2545: TSIJacobian ijac;
2546: TSI2Jacobian i2jac;
2547: TSRHSJacobian rhsjac;
2549: PetscFunctionBegin;
2551: if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
2553: if (!((PetscObject)ts)->type_name) {
2554: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2555: PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2556: }
2558: if (!ts->vec_sol) {
2559: PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2560: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2561: }
2563: if (ts->tspan) {
2564: if (!ts->tspan->vecs_sol) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2565: }
2566: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2567: PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2568: ts->Jacp = ts->Jacprhs;
2569: }
2571: if (ts->quadraturets) {
2572: PetscCall(TSSetUp(ts->quadraturets));
2573: PetscCall(VecDestroy(&ts->vec_costintegrand));
2574: PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2575: }
2577: PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
2578: if (rhsjac == TSComputeRHSJacobianConstant) {
2579: Mat Amat, Pmat;
2580: SNES snes;
2581: PetscCall(TSGetSNES(ts, &snes));
2582: PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
2583: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2584: * have displaced the RHS matrix */
2585: if (Amat && Amat == ts->Arhs) {
2586: /* 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 */
2587: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
2588: PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
2589: PetscCall(MatDestroy(&Amat));
2590: }
2591: if (Pmat && Pmat == ts->Brhs) {
2592: PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
2593: PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
2594: PetscCall(MatDestroy(&Pmat));
2595: }
2596: }
2598: PetscCall(TSGetAdapt(ts, &ts->adapt));
2599: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2601: PetscTryTypeMethod(ts, setup);
2603: PetscCall(TSSetExactFinalTimeDefault(ts));
2605: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2606: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2607: */
2608: PetscCall(TSGetDM(ts, &dm));
2609: PetscCall(DMSNESGetFunction(dm, &func, NULL));
2610: if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));
2612: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2613: Otherwise, the SNES will use coloring internally to form the Jacobian.
2614: */
2615: PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2616: PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2617: PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2618: PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2619: if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));
2621: /* if time integration scheme has a starting method, call it */
2622: PetscTryTypeMethod(ts, startingmethod);
2624: ts->setupcalled = PETSC_TRUE;
2625: PetscFunctionReturn(PETSC_SUCCESS);
2626: }
2628: /*@
2629: TSReset - Resets a `TS` context and removes any allocated `Vec`s and `Mat`s.
2631: Collective
2633: Input Parameter:
2634: . ts - the `TS` context obtained from `TSCreate()`
2636: Level: beginner
2638: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()`
2639: @*/
2640: PetscErrorCode TSReset(TS ts)
2641: {
2642: TS_RHSSplitLink ilink = ts->tsrhssplit, next;
2644: PetscFunctionBegin;
2647: PetscTryTypeMethod(ts, reset);
2648: if (ts->snes) PetscCall(SNESReset(ts->snes));
2649: if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));
2651: PetscCall(MatDestroy(&ts->Arhs));
2652: PetscCall(MatDestroy(&ts->Brhs));
2653: PetscCall(VecDestroy(&ts->Frhs));
2654: PetscCall(VecDestroy(&ts->vec_sol));
2655: PetscCall(VecDestroy(&ts->vec_dot));
2656: PetscCall(VecDestroy(&ts->vatol));
2657: PetscCall(VecDestroy(&ts->vrtol));
2658: PetscCall(VecDestroyVecs(ts->nwork, &ts->work));
2660: PetscCall(MatDestroy(&ts->Jacprhs));
2661: PetscCall(MatDestroy(&ts->Jacp));
2662: if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2663: if (ts->quadraturets) {
2664: PetscCall(TSReset(ts->quadraturets));
2665: PetscCall(VecDestroy(&ts->vec_costintegrand));
2666: }
2667: while (ilink) {
2668: next = ilink->next;
2669: PetscCall(TSDestroy(&ilink->ts));
2670: PetscCall(PetscFree(ilink->splitname));
2671: PetscCall(ISDestroy(&ilink->is));
2672: PetscCall(PetscFree(ilink));
2673: ilink = next;
2674: }
2675: ts->tsrhssplit = NULL;
2676: ts->num_rhs_splits = 0;
2677: if (ts->tspan) {
2678: PetscCall(PetscFree(ts->tspan->span_times));
2679: PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
2680: PetscCall(PetscFree(ts->tspan));
2681: }
2682: ts->setupcalled = PETSC_FALSE;
2683: PetscFunctionReturn(PETSC_SUCCESS);
2684: }
2686: /*@C
2687: TSDestroy - Destroys the timestepper context that was created
2688: with `TSCreate()`.
2690: Collective
2692: Input Parameter:
2693: . ts - the `TS` context obtained from `TSCreate()`
2695: Level: beginner
2697: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2698: @*/
2699: PetscErrorCode TSDestroy(TS *ts)
2700: {
2701: PetscFunctionBegin;
2702: if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2704: if (--((PetscObject)(*ts))->refct > 0) {
2705: *ts = NULL;
2706: PetscFunctionReturn(PETSC_SUCCESS);
2707: }
2709: PetscCall(TSReset(*ts));
2710: PetscCall(TSAdjointReset(*ts));
2711: if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2713: /* if memory was published with SAWs then destroy it */
2714: PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2715: PetscTryTypeMethod((*ts), destroy);
2717: PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));
2719: PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2720: PetscCall(TSEventDestroy(&(*ts)->event));
2722: PetscCall(SNESDestroy(&(*ts)->snes));
2723: PetscCall(DMDestroy(&(*ts)->dm));
2724: PetscCall(TSMonitorCancel((*ts)));
2725: PetscCall(TSAdjointMonitorCancel((*ts)));
2727: PetscCall(TSDestroy(&(*ts)->quadraturets));
2728: PetscCall(PetscHeaderDestroy(ts));
2729: PetscFunctionReturn(PETSC_SUCCESS);
2730: }
2732: /*@
2733: TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2734: a `TS` (timestepper) context. Valid only for nonlinear problems.
2736: Not Collective, but snes is parallel if ts is parallel
2738: Input Parameter:
2739: . ts - the `TS` context obtained from `TSCreate()`
2741: Output Parameter:
2742: . snes - the nonlinear solver context
2744: Level: beginner
2746: Notes:
2747: The user can then directly manipulate the `SNES` context to set various
2748: options, etc. Likewise, the user can then extract and manipulate the
2749: `KSP`, and `PC` contexts as well.
2751: `TSGetSNES()` does not work for integrators that do not use `SNES`; in
2752: this case `TSGetSNES()` returns `NULL` in `snes`.
2754: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2755: @*/
2756: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2757: {
2758: PetscFunctionBegin;
2761: if (!ts->snes) {
2762: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2763: PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2764: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2765: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2766: if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2767: if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2768: }
2769: *snes = ts->snes;
2770: PetscFunctionReturn(PETSC_SUCCESS);
2771: }
2773: /*@
2774: TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the timestepping context
2776: Collective
2778: Input Parameters:
2779: + ts - the `TS` context obtained from `TSCreate()`
2780: - snes - the nonlinear solver context
2782: Level: developer
2784: Note:
2785: Most users should have the `TS` created by calling `TSGetSNES()`
2787: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2788: @*/
2789: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2790: {
2791: PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);
2793: PetscFunctionBegin;
2796: PetscCall(PetscObjectReference((PetscObject)snes));
2797: PetscCall(SNESDestroy(&ts->snes));
2799: ts->snes = snes;
2801: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2802: PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2803: if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2804: PetscFunctionReturn(PETSC_SUCCESS);
2805: }
2807: /*@
2808: TSGetKSP - Returns the `KSP` (linear solver) associated with
2809: a `TS` (timestepper) context.
2811: Not Collective, but `ksp` is parallel if `ts` is parallel
2813: Input Parameter:
2814: . ts - the `TS` context obtained from `TSCreate()`
2816: Output Parameter:
2817: . ksp - the nonlinear solver context
2819: Level: beginner
2821: Notes:
2822: The user can then directly manipulate the `KSP` context to set various
2823: options, etc. Likewise, the user can then extract and manipulate the
2824: `PC` context as well.
2826: `TSGetKSP()` does not work for integrators that do not use `KSP`;
2827: in this case `TSGetKSP()` returns `NULL` in `ksp`.
2829: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2830: @*/
2831: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2832: {
2833: SNES snes;
2835: PetscFunctionBegin;
2838: PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2839: PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
2840: PetscCall(TSGetSNES(ts, &snes));
2841: PetscCall(SNESGetKSP(snes, ksp));
2842: PetscFunctionReturn(PETSC_SUCCESS);
2843: }
2845: /* ----------- Routines to set solver parameters ---------- */
2847: /*@
2848: TSSetMaxSteps - Sets the maximum number of steps to use.
2850: Logically Collective
2852: Input Parameters:
2853: + ts - the `TS` context obtained from `TSCreate()`
2854: - maxsteps - maximum number of steps to use
2856: Options Database Key:
2857: . -ts_max_steps <maxsteps> - Sets maxsteps
2859: Level: intermediate
2861: Note:
2862: The default maximum number of steps is 5000
2864: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2865: @*/
2866: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2867: {
2868: PetscFunctionBegin;
2871: PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2872: ts->max_steps = maxsteps;
2873: PetscFunctionReturn(PETSC_SUCCESS);
2874: }
2876: /*@
2877: TSGetMaxSteps - Gets the maximum number of steps to use.
2879: Not Collective
2881: Input Parameter:
2882: . ts - the `TS` context obtained from `TSCreate()`
2884: Output Parameter:
2885: . maxsteps - maximum number of steps to use
2887: Level: advanced
2889: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2890: @*/
2891: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2892: {
2893: PetscFunctionBegin;
2896: *maxsteps = ts->max_steps;
2897: PetscFunctionReturn(PETSC_SUCCESS);
2898: }
2900: /*@
2901: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
2903: Logically Collective
2905: Input Parameters:
2906: + ts - the `TS` context obtained from `TSCreate()`
2907: - maxtime - final time to step to
2909: Options Database Key:
2910: . -ts_max_time <maxtime> - Sets maxtime
2912: Level: intermediate
2914: Notes:
2915: The default maximum time is 5.0
2917: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2918: @*/
2919: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2920: {
2921: PetscFunctionBegin;
2924: ts->max_time = maxtime;
2925: PetscFunctionReturn(PETSC_SUCCESS);
2926: }
2928: /*@
2929: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
2931: Not Collective
2933: Input Parameter:
2934: . ts - the `TS` context obtained from `TSCreate()`
2936: Output Parameter:
2937: . maxtime - final time to step to
2939: Level: advanced
2941: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2942: @*/
2943: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2944: {
2945: PetscFunctionBegin;
2948: *maxtime = ts->max_time;
2949: PetscFunctionReturn(PETSC_SUCCESS);
2950: }
2952: /*@
2953: TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.
2955: Level: deprecated
2957: @*/
2958: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2959: {
2960: PetscFunctionBegin;
2962: PetscCall(TSSetTime(ts, initial_time));
2963: PetscCall(TSSetTimeStep(ts, time_step));
2964: PetscFunctionReturn(PETSC_SUCCESS);
2965: }
2967: /*@
2968: TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.
2970: Level: deprecated
2972: @*/
2973: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2974: {
2975: PetscFunctionBegin;
2977: if (maxsteps) {
2979: *maxsteps = ts->max_steps;
2980: }
2981: if (maxtime) {
2983: *maxtime = ts->max_time;
2984: }
2985: PetscFunctionReturn(PETSC_SUCCESS);
2986: }
2988: /*@
2989: TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.
2991: Level: deprecated
2993: @*/
2994: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2995: {
2996: PetscFunctionBegin;
3000: if (maxsteps >= 0) ts->max_steps = maxsteps;
3001: if (maxtime != (PetscReal)PETSC_DEFAULT) ts->max_time = maxtime;
3002: PetscFunctionReturn(PETSC_SUCCESS);
3003: }
3005: /*@
3006: TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.
3008: Level: deprecated
3010: @*/
3011: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
3012: {
3013: return TSGetStepNumber(ts, steps);
3014: }
3016: /*@
3017: TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.
3019: Level: deprecated
3021: @*/
3022: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
3023: {
3024: return TSGetStepNumber(ts, steps);
3025: }
3027: /*@
3028: TSSetSolution - Sets the initial solution vector
3029: for use by the `TS` routines.
3031: Logically Collective
3033: Input Parameters:
3034: + ts - the `TS` context obtained from `TSCreate()`
3035: - u - the solution vector
3037: Level: beginner
3039: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
3040: @*/
3041: PetscErrorCode TSSetSolution(TS ts, Vec u)
3042: {
3043: DM dm;
3045: PetscFunctionBegin;
3048: PetscCall(PetscObjectReference((PetscObject)u));
3049: PetscCall(VecDestroy(&ts->vec_sol));
3050: ts->vec_sol = u;
3052: PetscCall(TSGetDM(ts, &dm));
3053: PetscCall(DMShellSetGlobalVector(dm, u));
3054: PetscFunctionReturn(PETSC_SUCCESS);
3055: }
3057: /*@C
3058: TSSetPreStep - Sets the general-purpose function
3059: called once at the beginning of each time step.
3061: Logically Collective
3063: Input Parameters:
3064: + ts - The `TS` context obtained from `TSCreate()`
3065: - func - The function
3067: Calling sequence of `func`:
3068: .vb
3069: PetscErrorCode func (TS ts)
3070: .ve
3072: Level: intermediate
3074: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3075: @*/
3076: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
3077: {
3078: PetscFunctionBegin;
3080: ts->prestep = func;
3081: PetscFunctionReturn(PETSC_SUCCESS);
3082: }
3084: /*@
3085: TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`
3087: Collective
3089: Input Parameter:
3090: . ts - The `TS` context obtained from `TSCreate()`
3092: Level: developer
3094: Note:
3095: `TSPreStep()` is typically used within time stepping implementations,
3096: so most users would not generally call this routine themselves.
3098: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3099: @*/
3100: PetscErrorCode TSPreStep(TS ts)
3101: {
3102: PetscFunctionBegin;
3104: if (ts->prestep) {
3105: Vec U;
3106: PetscObjectId idprev;
3107: PetscBool sameObject;
3108: PetscObjectState sprev, spost;
3110: PetscCall(TSGetSolution(ts, &U));
3111: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3112: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3113: PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3114: PetscCall(TSGetSolution(ts, &U));
3115: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3116: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3117: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3118: }
3119: PetscFunctionReturn(PETSC_SUCCESS);
3120: }
3122: /*@C
3123: TSSetPreStage - Sets the general-purpose function
3124: called once at the beginning of each stage.
3126: Logically Collective
3128: Input Parameters:
3129: + ts - The `TS` context obtained from `TSCreate()`
3130: - func - The function
3132: Calling sequence of `func`:
3133: .vb
3134: PetscErrorCode func(TS ts, PetscReal stagetime)
3135: .ve
3137: Level: intermediate
3139: Note:
3140: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3141: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3142: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3144: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3145: @*/
3146: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS, PetscReal))
3147: {
3148: PetscFunctionBegin;
3150: ts->prestage = func;
3151: PetscFunctionReturn(PETSC_SUCCESS);
3152: }
3154: /*@C
3155: TSSetPostStage - Sets the general-purpose function, provided with `TSSetPostStep()`,
3156: called once at the end of each stage.
3158: Logically Collective
3160: Input Parameters:
3161: + ts - The `TS` context obtained from `TSCreate()`
3162: - func - The function
3164: Calling sequence of `func`:
3165: .vb
3166: PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y)
3167: .ve
3169: Level: intermediate
3171: Note:
3172: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3173: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3174: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3176: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3177: @*/
3178: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS, PetscReal, PetscInt, Vec *))
3179: {
3180: PetscFunctionBegin;
3182: ts->poststage = func;
3183: PetscFunctionReturn(PETSC_SUCCESS);
3184: }
3186: /*@C
3187: TSSetPostEvaluate - Sets the general-purpose function
3188: called once at the end of each step evaluation.
3190: Logically Collective
3192: Input Parameters:
3193: + ts - The `TS` context obtained from `TSCreate()`
3194: - func - The function
3196: Calling sequence of `func`:
3197: .vb
3198: PetscErrorCode func(TS ts)
3199: .ve
3201: Level: intermediate
3203: Note:
3204: Semantically, `TSSetPostEvaluate()` differs from `TSSetPostStep()` since the function it sets is called before event-handling
3205: thus guaranteeing the same solution (computed by the time-stepper) will be passed to it. On the other hand, `TSPostStep()`
3206: may be passed a different solution, possibly changed by the event handler. `TSPostEvaluate()` is called after the next step
3207: solution is evaluated allowing to modify it, if need be. The solution can be obtained with `TSGetSolution()`, the time step
3208: with `TSGetTimeStep()`, and the time at the start of the step is available via `TSGetTime()`
3210: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3211: @*/
3212: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS))
3213: {
3214: PetscFunctionBegin;
3216: ts->postevaluate = func;
3217: PetscFunctionReturn(PETSC_SUCCESS);
3218: }
3220: /*@
3221: TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`
3223: Collective
3225: Input Parameters:
3226: . ts - The `TS` context obtained from `TSCreate()`
3227: stagetime - The absolute time of the current stage
3229: Level: developer
3231: Note:
3232: `TSPreStage()` is typically used within time stepping implementations,
3233: most users would not generally call this routine themselves.
3235: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3236: @*/
3237: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3238: {
3239: PetscFunctionBegin;
3241: if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3242: PetscFunctionReturn(PETSC_SUCCESS);
3243: }
3245: /*@
3246: TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`
3248: Collective
3250: Input Parameters:
3251: . ts - The `TS` context obtained from `TSCreate()`
3252: stagetime - The absolute time of the current stage
3253: stageindex - Stage number
3254: Y - Array of vectors (of size = total number
3255: of stages) with the stage solutions
3257: Level: developer
3259: Note:
3260: `TSPostStage()` is typically used within time stepping implementations,
3261: most users would not generally call this routine themselves.
3263: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3264: @*/
3265: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3266: {
3267: PetscFunctionBegin;
3269: if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3270: PetscFunctionReturn(PETSC_SUCCESS);
3271: }
3273: /*@
3274: TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`
3276: Collective
3278: Input Parameter:
3279: . ts - The `TS` context obtained from `TSCreate()`
3281: Level: developer
3283: Note:
3284: `TSPostEvaluate()` is typically used within time stepping implementations,
3285: most users would not generally call this routine themselves.
3287: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3288: @*/
3289: PetscErrorCode TSPostEvaluate(TS ts)
3290: {
3291: PetscFunctionBegin;
3293: if (ts->postevaluate) {
3294: Vec U;
3295: PetscObjectState sprev, spost;
3297: PetscCall(TSGetSolution(ts, &U));
3298: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3299: PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3300: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3301: if (sprev != spost) PetscCall(TSRestartStep(ts));
3302: }
3303: PetscFunctionReturn(PETSC_SUCCESS);
3304: }
3306: /*@C
3307: TSSetPostStep - Sets the general-purpose function
3308: called once at the end of each time step.
3310: Logically Collective
3312: Input Parameters:
3313: + ts - The `TS` context obtained from `TSCreate()`
3314: - func - The function
3316: Calling sequence of `func`:
3317: $ PetscErrorCode func(TS ts)
3319: Level: intermediate
3321: Note:
3322: The function set by `TSSetPostStep()` is called after each successful step. The solution vector
3323: obtained by `TSGetSolution()` may be different than that computed at the step end if the event handler
3324: locates an event and `TSPostEvent()` modifies it. Use `TSSetPostEvaluate()` if an unmodified solution is needed instead.
3326: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3327: @*/
3328: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
3329: {
3330: PetscFunctionBegin;
3332: ts->poststep = func;
3333: PetscFunctionReturn(PETSC_SUCCESS);
3334: }
3336: /*@
3337: TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`
3339: Collective
3341: Input Parameter:
3342: . ts - The `TS` context obtained from `TSCreate()`
3344: Note:
3345: `TSPostStep()` is typically used within time stepping implementations,
3346: so most users would not generally call this routine themselves.
3348: Level: developer
3350: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()`
3351: @*/
3352: PetscErrorCode TSPostStep(TS ts)
3353: {
3354: PetscFunctionBegin;
3356: if (ts->poststep) {
3357: Vec U;
3358: PetscObjectId idprev;
3359: PetscBool sameObject;
3360: PetscObjectState sprev, spost;
3362: PetscCall(TSGetSolution(ts, &U));
3363: PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
3364: PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
3365: PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3366: PetscCall(TSGetSolution(ts, &U));
3367: PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
3368: PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
3369: if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
3370: }
3371: PetscFunctionReturn(PETSC_SUCCESS);
3372: }
3374: /*@
3375: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3377: Collective
3379: Input Parameters:
3380: + ts - time stepping context
3381: - t - time to interpolate to
3383: Output Parameter:
3384: . U - state at given time
3386: Level: intermediate
3388: Developer Note:
3389: `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3391: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3392: @*/
3393: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3394: {
3395: PetscFunctionBegin;
3398: 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);
3399: PetscUseTypeMethod(ts, interpolate, t, U);
3400: PetscFunctionReturn(PETSC_SUCCESS);
3401: }
3403: /*@
3404: TSStep - Steps one time step
3406: Collective
3408: Input Parameter:
3409: . ts - the `TS` context obtained from `TSCreate()`
3411: Level: developer
3413: Notes:
3414: The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.
3416: The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
3417: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3419: This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
3420: time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.
3422: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3423: @*/
3424: PetscErrorCode TSStep(TS ts)
3425: {
3426: static PetscBool cite = PETSC_FALSE;
3427: PetscReal ptime;
3429: PetscFunctionBegin;
3431: PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3432: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3433: " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3434: " journal = {arXiv e-preprints},\n"
3435: " eprint = {1806.01437},\n"
3436: " archivePrefix = {arXiv},\n"
3437: " year = {2018}\n}\n",
3438: &cite));
3439: PetscCall(TSSetUp(ts));
3440: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3442: 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>");
3443: 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()");
3444: 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");
3446: if (!ts->steps) ts->ptime_prev = ts->ptime;
3447: ptime = ts->ptime;
3448: ts->ptime_prev_rollback = ts->ptime_prev;
3449: ts->reason = TS_CONVERGED_ITERATING;
3451: PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3452: PetscUseTypeMethod(ts, step);
3453: PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));
3455: if (ts->tspan && PetscIsCloseAtTol(ts->ptime, ts->tspan->span_times[ts->tspan->spanctr], ts->tspan->reltol * ts->time_step + ts->tspan->abstol, 0) && ts->tspan->spanctr < ts->tspan->num_span_times)
3456: PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[ts->tspan->spanctr++]));
3457: if (ts->reason >= 0) {
3458: ts->ptime_prev = ptime;
3459: ts->steps++;
3460: ts->steprollback = PETSC_FALSE;
3461: ts->steprestart = PETSC_FALSE;
3462: }
3463: if (!ts->reason) {
3464: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3465: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3466: }
3468: if (ts->reason < 0 && ts->errorifstepfailed) {
3469: PetscCall(TSMonitorCancel(ts));
3470: 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]);
3471: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3472: }
3473: PetscFunctionReturn(PETSC_SUCCESS);
3474: }
3476: /*@
3477: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3478: at the end of a time step with a given order of accuracy.
3480: Collective
3482: Input Parameters:
3483: + ts - time stepping context
3484: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3486: Input/Output Parameter:
3487: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3488: on output, the actual order of the error evaluation
3490: Output Parameter:
3491: . wlte - the weighted local truncation error norm
3493: Level: advanced
3495: Note:
3496: If the timestepper cannot evaluate the error in a particular step
3497: (eg. in the first step or restart steps after event handling),
3498: this routine returns wlte=-1.0 .
3500: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3501: @*/
3502: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3503: {
3504: PetscFunctionBegin;
3511: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3512: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3513: PetscFunctionReturn(PETSC_SUCCESS);
3514: }
3516: /*@
3517: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3519: Collective
3521: Input Parameters:
3522: + ts - time stepping context
3523: . order - desired order of accuracy
3524: - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)
3526: Output Parameter:
3527: . U - state at the end of the current step
3529: Level: advanced
3531: Notes:
3532: This function cannot be called until all stages have been evaluated.
3534: 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.
3536: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3537: @*/
3538: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3539: {
3540: PetscFunctionBegin;
3544: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3545: PetscFunctionReturn(PETSC_SUCCESS);
3546: }
3548: /*@C
3549: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3551: Not collective
3553: Input Parameter:
3554: . ts - time stepping context
3556: Output Parameter:
3557: . initConditions - The function which computes an initial condition
3559: Calling sequence of `initCondition`:
3560: .vb
3561: PetscErrorCode initCondition(TS ts, Vec u)
3562: .ve
3563: + ts - The timestepping context
3564: - u - The input vector in which the initial condition is stored
3566: Level: advanced
3568: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3569: @*/
3570: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS, Vec))
3571: {
3572: PetscFunctionBegin;
3575: *initCondition = ts->ops->initcondition;
3576: PetscFunctionReturn(PETSC_SUCCESS);
3577: }
3579: /*@C
3580: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3582: Logically collective
3584: Input Parameters:
3585: + ts - time stepping context
3586: - initCondition - The function which computes an initial condition
3588: Calling sequence of `initCondition`:
3589: $ PetscErrorCode initCondition(TS ts, Vec u)
3590: + ts - The timestepping context
3591: - u - The input vector in which the initial condition is to be stored
3593: Level: advanced
3595: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3596: @*/
3597: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS, Vec))
3598: {
3599: PetscFunctionBegin;
3602: ts->ops->initcondition = initCondition;
3603: PetscFunctionReturn(PETSC_SUCCESS);
3604: }
3606: /*@
3607: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3609: Collective
3611: Input Parameters:
3612: + ts - time stepping context
3613: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3615: Level: advanced
3617: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3618: @*/
3619: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3620: {
3621: PetscFunctionBegin;
3624: PetscTryTypeMethod(ts, initcondition, u);
3625: PetscFunctionReturn(PETSC_SUCCESS);
3626: }
3628: /*@C
3629: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3631: Not collective
3633: Input Parameter:
3634: . ts - time stepping context
3636: Output Parameter:
3637: . exactError - The function which computes the solution error
3639: Calling sequence of `exactError`:
3640: $ PetscErrorCode exactError(TS ts, Vec u, Vec e)
3641: + ts - The timestepping context
3642: . u - The approximate solution vector
3643: - e - The vector in which the error is stored
3645: Level: advanced
3647: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3648: @*/
3649: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS, Vec, Vec))
3650: {
3651: PetscFunctionBegin;
3654: *exactError = ts->ops->exacterror;
3655: PetscFunctionReturn(PETSC_SUCCESS);
3656: }
3658: /*@C
3659: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3661: Logically collective
3663: Input Parameters:
3664: + ts - time stepping context
3665: - exactError - The function which computes the solution error
3667: Calling sequence of `exactError`:
3668: $ PetscErrorCode exactError(TS ts, Vec u)
3669: + ts - The timestepping context
3670: . u - The approximate solution vector
3671: - e - The vector in which the error is stored
3673: Level: advanced
3675: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3676: @*/
3677: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS, Vec, Vec))
3678: {
3679: PetscFunctionBegin;
3682: ts->ops->exacterror = exactError;
3683: PetscFunctionReturn(PETSC_SUCCESS);
3684: }
3686: /*@
3687: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3689: Collective
3691: Input Parameters:
3692: + ts - time stepping context
3693: . u - The approximate solution
3694: - e - The `Vec` used to store the error
3696: Level: advanced
3698: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3699: @*/
3700: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3701: {
3702: PetscFunctionBegin;
3706: PetscTryTypeMethod(ts, exacterror, u, e);
3707: PetscFunctionReturn(PETSC_SUCCESS);
3708: }
3710: /*@
3711: TSSolve - Steps the requested number of timesteps.
3713: Collective
3715: Input Parameters:
3716: + ts - the `TS` context obtained from `TSCreate()`
3717: - u - the solution vector (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
3718: otherwise must contain the initial conditions and will contain the solution at the final requested time
3720: Level: beginner
3722: Notes:
3723: The final time returned by this function may be different from the time of the internally
3724: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
3725: stepped over the final time.
3727: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3728: @*/
3729: PetscErrorCode TSSolve(TS ts, Vec u)
3730: {
3731: Vec solution;
3733: PetscFunctionBegin;
3737: PetscCall(TSSetExactFinalTimeDefault(ts));
3738: 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 */
3739: if (!ts->vec_sol || u == ts->vec_sol) {
3740: PetscCall(VecDuplicate(u, &solution));
3741: PetscCall(TSSetSolution(ts, solution));
3742: PetscCall(VecDestroy(&solution)); /* grant ownership */
3743: }
3744: PetscCall(VecCopy(u, ts->vec_sol));
3745: PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
3746: } else if (u) PetscCall(TSSetSolution(ts, u));
3747: PetscCall(TSSetUp(ts));
3748: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3750: 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>");
3751: 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()");
3752: 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");
3753: 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");
3755: 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 */
3756: PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]));
3757: ts->tspan->spanctr = 1;
3758: }
3760: if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));
3762: /* reset number of steps only when the step is not restarted. ARKIMEX
3763: restarts the step after an event. Resetting these counters in such case causes
3764: TSTrajectory to incorrectly save the output files
3765: */
3766: /* reset time step and iteration counters */
3767: if (!ts->steps) {
3768: ts->ksp_its = 0;
3769: ts->snes_its = 0;
3770: ts->num_snes_failures = 0;
3771: ts->reject = 0;
3772: ts->steprestart = PETSC_TRUE;
3773: ts->steprollback = PETSC_FALSE;
3774: ts->rhsjacobian.time = PETSC_MIN_REAL;
3775: }
3777: /* make sure initial time step does not overshoot final time or the next point in tspan */
3778: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
3779: PetscReal maxdt;
3780: PetscReal dt = ts->time_step;
3782: if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
3783: else maxdt = ts->max_time - ts->ptime;
3784: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
3785: }
3786: ts->reason = TS_CONVERGED_ITERATING;
3788: {
3789: PetscViewer viewer;
3790: PetscViewerFormat format;
3791: PetscBool flg;
3792: static PetscBool incall = PETSC_FALSE;
3794: if (!incall) {
3795: /* Estimate the convergence rate of the time discretization */
3796: PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
3797: if (flg) {
3798: PetscConvEst conv;
3799: DM dm;
3800: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
3801: PetscInt Nf;
3802: PetscBool checkTemporal = PETSC_TRUE;
3804: incall = PETSC_TRUE;
3805: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
3806: PetscCall(TSGetDM(ts, &dm));
3807: PetscCall(DMGetNumFields(dm, &Nf));
3808: PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
3809: PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
3810: PetscCall(PetscConvEstUseTS(conv, checkTemporal));
3811: PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
3812: PetscCall(PetscConvEstSetFromOptions(conv));
3813: PetscCall(PetscConvEstSetUp(conv));
3814: PetscCall(PetscConvEstGetConvRate(conv, alpha));
3815: PetscCall(PetscViewerPushFormat(viewer, format));
3816: PetscCall(PetscConvEstRateView(conv, alpha, viewer));
3817: PetscCall(PetscViewerPopFormat(viewer));
3818: PetscCall(PetscViewerDestroy(&viewer));
3819: PetscCall(PetscConvEstDestroy(&conv));
3820: PetscCall(PetscFree(alpha));
3821: incall = PETSC_FALSE;
3822: }
3823: }
3824: }
3826: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));
3828: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
3829: PetscUseTypeMethod(ts, solve);
3830: if (u) PetscCall(VecCopy(ts->vec_sol, u));
3831: ts->solvetime = ts->ptime;
3832: solution = ts->vec_sol;
3833: } else { /* Step the requested number of timesteps. */
3834: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3835: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3837: if (!ts->steps) {
3838: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
3839: PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
3840: }
3842: while (!ts->reason) {
3843: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
3844: if (!ts->steprollback) PetscCall(TSPreStep(ts));
3845: PetscCall(TSStep(ts));
3846: if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
3847: if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
3848: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
3849: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
3850: PetscCall(TSForwardCostIntegral(ts));
3851: if (ts->reason >= 0) ts->steps++;
3852: }
3853: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
3854: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
3855: PetscCall(TSForwardStep(ts));
3856: if (ts->reason >= 0) ts->steps++;
3857: }
3858: PetscCall(TSPostEvaluate(ts));
3859: 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. */
3860: if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
3861: if (!ts->steprollback) {
3862: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
3863: PetscCall(TSPostStep(ts));
3864: }
3865: }
3866: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
3868: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
3869: PetscCall(TSInterpolate(ts, ts->max_time, u));
3870: ts->solvetime = ts->max_time;
3871: solution = u;
3872: PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
3873: } else {
3874: if (u) PetscCall(VecCopy(ts->vec_sol, u));
3875: ts->solvetime = ts->ptime;
3876: solution = ts->vec_sol;
3877: }
3878: }
3880: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
3881: PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
3882: PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
3883: if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
3884: PetscFunctionReturn(PETSC_SUCCESS);
3885: }
3887: /*@
3888: TSGetTime - Gets the time of the most recently completed step.
3890: Not Collective
3892: Input Parameter:
3893: . ts - the `TS` context obtained from `TSCreate()`
3895: Output Parameter:
3896: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
3898: Level: beginner
3900: Note:
3901: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
3902: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
3904: .seealso: [](ch_ts), TS`, ``TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
3905: @*/
3906: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
3907: {
3908: PetscFunctionBegin;
3911: *t = ts->ptime;
3912: PetscFunctionReturn(PETSC_SUCCESS);
3913: }
3915: /*@
3916: TSGetPrevTime - Gets the starting time of the previously completed step.
3918: Not Collective
3920: Input Parameter:
3921: . ts - the `TS` context obtained from `TSCreate()`
3923: Output Parameter:
3924: . t - the previous time
3926: Level: beginner
3928: .seealso: [](ch_ts), TS`, ``TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
3929: @*/
3930: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
3931: {
3932: PetscFunctionBegin;
3935: *t = ts->ptime_prev;
3936: PetscFunctionReturn(PETSC_SUCCESS);
3937: }
3939: /*@
3940: TSSetTime - Allows one to reset the time.
3942: Logically Collective
3944: Input Parameters:
3945: + ts - the `TS` context obtained from `TSCreate()`
3946: - time - the time
3948: Level: intermediate
3950: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
3951: @*/
3952: PetscErrorCode TSSetTime(TS ts, PetscReal t)
3953: {
3954: PetscFunctionBegin;
3957: ts->ptime = t;
3958: PetscFunctionReturn(PETSC_SUCCESS);
3959: }
3961: /*@C
3962: TSSetOptionsPrefix - Sets the prefix used for searching for all
3963: TS options in the database.
3965: Logically Collective
3967: Input Parameters:
3968: + ts - The `TS` context
3969: - prefix - The prefix to prepend to all option names
3971: Level: advanced
3973: Note:
3974: A hyphen (-) must NOT be given at the beginning of the prefix name.
3975: The first character of all runtime options is AUTOMATICALLY the
3976: hyphen.
3978: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
3979: @*/
3980: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
3981: {
3982: SNES snes;
3984: PetscFunctionBegin;
3986: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
3987: PetscCall(TSGetSNES(ts, &snes));
3988: PetscCall(SNESSetOptionsPrefix(snes, prefix));
3989: PetscFunctionReturn(PETSC_SUCCESS);
3990: }
3992: /*@C
3993: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
3994: TS options in the database.
3996: Logically Collective
3998: Input Parameters:
3999: + ts - The `TS` context
4000: - prefix - The prefix to prepend to all option names
4002: Level: advanced
4004: Note:
4005: A hyphen (-) must NOT be given at the beginning of the prefix name.
4006: The first character of all runtime options is AUTOMATICALLY the
4007: hyphen.
4009: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4010: @*/
4011: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4012: {
4013: SNES snes;
4015: PetscFunctionBegin;
4017: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4018: PetscCall(TSGetSNES(ts, &snes));
4019: PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4020: PetscFunctionReturn(PETSC_SUCCESS);
4021: }
4023: /*@C
4024: TSGetOptionsPrefix - Sets the prefix used for searching for all
4025: `TS` options in the database.
4027: Not Collective
4029: Input Parameter:
4030: . ts - The `TS` context
4032: Output Parameter:
4033: . prefix - A pointer to the prefix string used
4035: Level: intermediate
4037: Fortran Note:
4038: The user should pass in a string 'prefix' of
4039: sufficient length to hold the prefix.
4041: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4042: @*/
4043: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4044: {
4045: PetscFunctionBegin;
4048: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4049: PetscFunctionReturn(PETSC_SUCCESS);
4050: }
4052: /*@C
4053: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4055: Not Collective, but parallel objects are returned if ts is parallel
4057: Input Parameter:
4058: . ts - The `TS` context obtained from `TSCreate()`
4060: Output Parameters:
4061: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`)
4062: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`)
4063: . func - Function to compute the Jacobian of the RHS (or `NULL`)
4064: - ctx - User-defined context for Jacobian evaluation routine (or `NULL`)
4066: Level: intermediate
4068: Note:
4069: You can pass in `NULL` for any return argument you do not need.
4071: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4073: @*/
4074: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobian *func, void **ctx)
4075: {
4076: DM dm;
4078: PetscFunctionBegin;
4079: if (Amat || Pmat) {
4080: SNES snes;
4081: PetscCall(TSGetSNES(ts, &snes));
4082: PetscCall(SNESSetUpMatrices(snes));
4083: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4084: }
4085: PetscCall(TSGetDM(ts, &dm));
4086: PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4087: PetscFunctionReturn(PETSC_SUCCESS);
4088: }
4090: /*@C
4091: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4093: Not Collective, but parallel objects are returned if ts is parallel
4095: Input Parameter:
4096: . ts - The `TS` context obtained from `TSCreate()`
4098: Output Parameters:
4099: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4100: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4101: . f - The function to compute the matrices
4102: - ctx - User-defined context for Jacobian evaluation routine
4104: Level: advanced
4106: Note:
4107: You can pass in `NULL` for any return argument you do not need.
4109: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4110: @*/
4111: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobian *f, void **ctx)
4112: {
4113: DM dm;
4115: PetscFunctionBegin;
4116: if (Amat || Pmat) {
4117: SNES snes;
4118: PetscCall(TSGetSNES(ts, &snes));
4119: PetscCall(SNESSetUpMatrices(snes));
4120: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4121: }
4122: PetscCall(TSGetDM(ts, &dm));
4123: PetscCall(DMTSGetIJacobian(dm, f, ctx));
4124: PetscFunctionReturn(PETSC_SUCCESS);
4125: }
4127: #include <petsc/private/dmimpl.h>
4128: /*@
4129: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4131: Logically Collective
4133: Input Parameters:
4134: + ts - the `TS` integrator object
4135: - dm - the dm, cannot be `NULL`
4137: Level: intermediate
4139: Notes:
4140: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4141: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4142: different problems using the same function space.
4144: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4145: @*/
4146: PetscErrorCode TSSetDM(TS ts, DM dm)
4147: {
4148: SNES snes;
4149: DMTS tsdm;
4151: PetscFunctionBegin;
4154: PetscCall(PetscObjectReference((PetscObject)dm));
4155: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4156: if (ts->dm->dmts && !dm->dmts) {
4157: PetscCall(DMCopyDMTS(ts->dm, dm));
4158: PetscCall(DMGetDMTS(ts->dm, &tsdm));
4159: /* Grant write privileges to the replacement DM */
4160: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4161: }
4162: PetscCall(DMDestroy(&ts->dm));
4163: }
4164: ts->dm = dm;
4166: PetscCall(TSGetSNES(ts, &snes));
4167: PetscCall(SNESSetDM(snes, dm));
4168: PetscFunctionReturn(PETSC_SUCCESS);
4169: }
4171: /*@
4172: TSGetDM - Gets the `DM` that may be used by some preconditioners
4174: Not Collective
4176: Input Parameter:
4177: . ts - the `TS`
4179: Output Parameter:
4180: . dm - the `DM`
4182: Level: intermediate
4184: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4185: @*/
4186: PetscErrorCode TSGetDM(TS ts, DM *dm)
4187: {
4188: PetscFunctionBegin;
4190: if (!ts->dm) {
4191: PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4192: if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4193: }
4194: *dm = ts->dm;
4195: PetscFunctionReturn(PETSC_SUCCESS);
4196: }
4198: /*@
4199: SNESTSFormFunction - Function to evaluate nonlinear residual
4201: Logically Collective
4203: Input Parameters:
4204: + snes - nonlinear solver
4205: . U - the current state at which to evaluate the residual
4206: - ctx - user context, must be a TS
4208: Output Parameter:
4209: . F - the nonlinear residual
4211: Level: advanced
4213: Note:
4214: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4215: It is most frequently passed to `MatFDColoringSetFunction()`.
4217: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4218: @*/
4219: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4220: {
4221: TS ts = (TS)ctx;
4223: PetscFunctionBegin;
4228: PetscCall((ts->ops->snesfunction)(snes, U, F, ts));
4229: PetscFunctionReturn(PETSC_SUCCESS);
4230: }
4232: /*@
4233: SNESTSFormJacobian - Function to evaluate the Jacobian
4235: Collective
4237: Input Parameters:
4238: + snes - nonlinear solver
4239: . U - the current state at which to evaluate the residual
4240: - ctx - user context, must be a `TS`
4242: Output Parameters:
4243: + A - the Jacobian
4244: - B - the preconditioning matrix (may be the same as A)
4246: Level: developer
4248: Note:
4249: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4251: .seealso: [](ch_ts), `SNESSetJacobian()`
4252: @*/
4253: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4254: {
4255: TS ts = (TS)ctx;
4257: PetscFunctionBegin;
4265: PetscCall((ts->ops->snesjacobian)(snes, U, A, B, ts));
4266: PetscFunctionReturn(PETSC_SUCCESS);
4267: }
4269: /*@C
4270: TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems Udot = A U only
4272: Collective
4274: Input Parameters:
4275: + ts - time stepping context
4276: . t - time at which to evaluate
4277: . U - state at which to evaluate
4278: - ctx - context
4280: Output Parameter:
4281: . F - right hand side
4283: Level: intermediate
4285: Note:
4286: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right hand side for linear problems.
4287: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4289: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4290: @*/
4291: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4292: {
4293: Mat Arhs, Brhs;
4295: PetscFunctionBegin;
4296: PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4297: /* undo the damage caused by shifting */
4298: PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4299: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4300: PetscCall(MatMult(Arhs, U, F));
4301: PetscFunctionReturn(PETSC_SUCCESS);
4302: }
4304: /*@C
4305: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
4307: Collective
4309: Input Parameters:
4310: + ts - time stepping context
4311: . t - time at which to evaluate
4312: . U - state at which to evaluate
4313: - ctx - context
4315: Output Parameters:
4316: + A - pointer to operator
4317: - B - pointer to preconditioning matrix
4319: Level: intermediate
4321: Note:
4322: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4324: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4325: @*/
4326: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4327: {
4328: PetscFunctionBegin;
4329: PetscFunctionReturn(PETSC_SUCCESS);
4330: }
4332: /*@C
4333: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4335: Collective
4337: Input Parameters:
4338: + ts - time stepping context
4339: . t - time at which to evaluate
4340: . U - state at which to evaluate
4341: . Udot - time derivative of state vector
4342: - ctx - context
4344: Output Parameter:
4345: . F - left hand side
4347: Level: intermediate
4349: Notes:
4350: 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
4351: user is required to write their own `TSComputeIFunction()`.
4352: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4353: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4355: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4357: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4358: @*/
4359: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4360: {
4361: Mat A, B;
4363: PetscFunctionBegin;
4364: PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4365: PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4366: PetscCall(MatMult(A, Udot, F));
4367: PetscFunctionReturn(PETSC_SUCCESS);
4368: }
4370: /*@C
4371: TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE
4373: Collective
4375: Input Parameters:
4376: + ts - time stepping context
4377: . t - time at which to evaluate
4378: . U - state at which to evaluate
4379: . Udot - time derivative of state vector
4380: . shift - shift to apply
4381: - ctx - context
4383: Output Parameters:
4384: + A - pointer to operator
4385: - B - pointer to preconditioning matrix
4387: Level: advanced
4389: Notes:
4390: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4392: It is only appropriate for problems of the form
4394: $ M Udot = F(U,t)
4396: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4397: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4398: an implicit operator of the form
4400: $ shift*M + J
4402: 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
4403: a copy of M or reassemble it when requested.
4405: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4406: @*/
4407: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4408: {
4409: PetscFunctionBegin;
4410: PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4411: ts->ijacobian.shift = shift;
4412: PetscFunctionReturn(PETSC_SUCCESS);
4413: }
4415: /*@
4416: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4418: Not Collective
4420: Input Parameter:
4421: . ts - the `TS` context
4423: Output Parameter:
4424: . equation_type - see `TSEquationType`
4426: Level: beginner
4428: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4429: @*/
4430: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4431: {
4432: PetscFunctionBegin;
4435: *equation_type = ts->equation_type;
4436: PetscFunctionReturn(PETSC_SUCCESS);
4437: }
4439: /*@
4440: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4442: Not Collective
4444: Input Parameters:
4445: + ts - the `TS` context
4446: - equation_type - see `TSEquationType`
4448: Level: advanced
4450: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4451: @*/
4452: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4453: {
4454: PetscFunctionBegin;
4456: ts->equation_type = equation_type;
4457: PetscFunctionReturn(PETSC_SUCCESS);
4458: }
4460: /*@
4461: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4463: Not Collective
4465: Input Parameter:
4466: . ts - the `TS` context
4468: Output Parameter:
4469: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4470: manual pages for the individual convergence tests for complete lists
4472: Level: beginner
4474: Note:
4475: Can only be called after the call to `TSSolve()` is complete.
4477: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4478: @*/
4479: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4480: {
4481: PetscFunctionBegin;
4484: *reason = ts->reason;
4485: PetscFunctionReturn(PETSC_SUCCESS);
4486: }
4488: /*@
4489: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4491: Logically Collective; reason must contain common value
4493: Input Parameters:
4494: + ts - the `TS` context
4495: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4496: manual pages for the individual convergence tests for complete lists
4498: Level: advanced
4500: Note:
4501: Can only be called while `TSSolve()` is active.
4503: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4504: @*/
4505: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4506: {
4507: PetscFunctionBegin;
4509: ts->reason = reason;
4510: PetscFunctionReturn(PETSC_SUCCESS);
4511: }
4513: /*@
4514: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4516: Not Collective
4518: Input Parameter:
4519: . ts - the `TS` context
4521: Output Parameter:
4522: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4524: Level: beginner
4526: Note:
4527: Can only be called after the call to `TSSolve()` is complete.
4529: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4530: @*/
4531: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4532: {
4533: PetscFunctionBegin;
4536: *ftime = ts->solvetime;
4537: PetscFunctionReturn(PETSC_SUCCESS);
4538: }
4540: /*@
4541: TSGetSNESIterations - Gets the total number of nonlinear iterations
4542: used by the time integrator.
4544: Not Collective
4546: Input Parameter:
4547: . ts - `TS` context
4549: Output Parameter:
4550: . nits - number of nonlinear iterations
4552: Level: intermediate
4554: Note:
4555: This counter is reset to zero for each successive call to `TSSolve()`.
4557: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4558: @*/
4559: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4560: {
4561: PetscFunctionBegin;
4564: *nits = ts->snes_its;
4565: PetscFunctionReturn(PETSC_SUCCESS);
4566: }
4568: /*@
4569: TSGetKSPIterations - Gets the total number of linear iterations
4570: used by the time integrator.
4572: Not Collective
4574: Input Parameter:
4575: . ts - `TS` context
4577: Output Parameter:
4578: . lits - number of linear iterations
4580: Level: intermediate
4582: Note:
4583: This counter is reset to zero for each successive call to `TSSolve()`.
4585: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()`
4586: @*/
4587: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4588: {
4589: PetscFunctionBegin;
4592: *lits = ts->ksp_its;
4593: PetscFunctionReturn(PETSC_SUCCESS);
4594: }
4596: /*@
4597: TSGetStepRejections - Gets the total number of rejected steps.
4599: Not Collective
4601: Input Parameter:
4602: . ts - `TS` context
4604: Output Parameter:
4605: . rejects - number of steps rejected
4607: Level: intermediate
4609: Note:
4610: This counter is reset to zero for each successive call to `TSSolve()`.
4612: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4613: @*/
4614: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4615: {
4616: PetscFunctionBegin;
4619: *rejects = ts->reject;
4620: PetscFunctionReturn(PETSC_SUCCESS);
4621: }
4623: /*@
4624: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
4626: Not Collective
4628: Input Parameter:
4629: . ts - `TS` context
4631: Output Parameter:
4632: . fails - number of failed nonlinear solves
4634: Level: intermediate
4636: Note:
4637: This counter is reset to zero for each successive call to `TSSolve()`.
4639: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4640: @*/
4641: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4642: {
4643: PetscFunctionBegin;
4646: *fails = ts->num_snes_failures;
4647: PetscFunctionReturn(PETSC_SUCCESS);
4648: }
4650: /*@
4651: TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails
4653: Not Collective
4655: Input Parameters:
4656: + ts - `TS` context
4657: - rejects - maximum number of rejected steps, pass -1 for unlimited
4659: Options Database Key:
4660: . -ts_max_reject - Maximum number of step rejections before a step fails
4662: Level: intermediate
4664: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4665: @*/
4666: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4667: {
4668: PetscFunctionBegin;
4670: ts->max_reject = rejects;
4671: PetscFunctionReturn(PETSC_SUCCESS);
4672: }
4674: /*@
4675: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves
4677: Not Collective
4679: Input Parameters:
4680: + ts - `TS` context
4681: - fails - maximum number of failed nonlinear solves, pass -1 for unlimited
4683: Options Database Key:
4684: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
4686: Level: intermediate
4688: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4689: @*/
4690: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4691: {
4692: PetscFunctionBegin;
4694: ts->max_snes_failures = fails;
4695: PetscFunctionReturn(PETSC_SUCCESS);
4696: }
4698: /*@
4699: TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`
4701: Not Collective
4703: Input Parameters:
4704: + ts - `TS` context
4705: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
4707: Options Database Key:
4708: . -ts_error_if_step_fails - Error if no step succeeds
4710: Level: intermediate
4712: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4713: @*/
4714: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
4715: {
4716: PetscFunctionBegin;
4718: ts->errorifstepfailed = err;
4719: PetscFunctionReturn(PETSC_SUCCESS);
4720: }
4722: /*@
4723: TSGetAdapt - Get the adaptive controller context for the current method
4725: Collective on `ts` if controller has not been created yet
4727: Input Parameter:
4728: . ts - time stepping context
4730: Output Parameter:
4731: . adapt - adaptive controller
4733: Level: intermediate
4735: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
4736: @*/
4737: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
4738: {
4739: PetscFunctionBegin;
4742: if (!ts->adapt) {
4743: PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
4744: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
4745: }
4746: *adapt = ts->adapt;
4747: PetscFunctionReturn(PETSC_SUCCESS);
4748: }
4750: /*@
4751: TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller
4753: Logically Collective
4755: Input Parameters:
4756: + ts - time integration context
4757: . atol - scalar absolute tolerances, `PETSC_DECIDE` to leave current value
4758: . vatol - vector of absolute tolerances or `NULL`, used in preference to atol if present
4759: . rtol - scalar relative tolerances, `PETSC_DECIDE` to leave current value
4760: - vrtol - vector of relative tolerances or `NULL`, used in preference to atol if present
4762: Options Database Keys:
4763: + -ts_rtol <rtol> - relative tolerance for local truncation error
4764: - -ts_atol <atol> - Absolute tolerance for local truncation error
4766: Level: beginner
4768: Notes:
4769: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
4770: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
4771: computed only for the differential or the algebraic part then this can be done using the vector of
4772: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
4773: differential part and infinity for the algebraic part, the LTE calculation will include only the
4774: differential variables.
4776: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
4777: @*/
4778: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
4779: {
4780: PetscFunctionBegin;
4781: if (atol != (PetscReal)PETSC_DECIDE && atol != (PetscReal)PETSC_DEFAULT) ts->atol = atol;
4782: if (vatol) {
4783: PetscCall(PetscObjectReference((PetscObject)vatol));
4784: PetscCall(VecDestroy(&ts->vatol));
4785: ts->vatol = vatol;
4786: }
4787: if (rtol != (PetscReal)PETSC_DECIDE && rtol != (PetscReal)PETSC_DEFAULT) ts->rtol = rtol;
4788: if (vrtol) {
4789: PetscCall(PetscObjectReference((PetscObject)vrtol));
4790: PetscCall(VecDestroy(&ts->vrtol));
4791: ts->vrtol = vrtol;
4792: }
4793: PetscFunctionReturn(PETSC_SUCCESS);
4794: }
4796: /*@
4797: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
4799: Logically Collective
4801: Input Parameter:
4802: . ts - time integration context
4804: Output Parameters:
4805: + atol - scalar absolute tolerances, `NULL` to ignore
4806: . vatol - vector of absolute tolerances, `NULL` to ignore
4807: . rtol - scalar relative tolerances, `NULL` to ignore
4808: - vrtol - vector of relative tolerances, `NULL` to ignore
4810: Level: beginner
4812: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
4813: @*/
4814: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
4815: {
4816: PetscFunctionBegin;
4817: if (atol) *atol = ts->atol;
4818: if (vatol) *vatol = ts->vatol;
4819: if (rtol) *rtol = ts->rtol;
4820: if (vrtol) *vrtol = ts->vrtol;
4821: PetscFunctionReturn(PETSC_SUCCESS);
4822: }
4824: /*@
4825: TSErrorWeightedNorm2 - compute a weighted 2-norm of the difference between two state vectors
4827: Collective
4829: Input Parameters:
4830: + ts - time stepping context
4831: . U - state vector, usually ts->vec_sol
4832: - Y - state vector to be compared to U
4834: Output Parameters:
4835: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
4836: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
4837: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
4839: Level: developer
4841: .seealso: [](ch_ts), `TS`, `TSErrorWeightedNorm()`, `TSErrorWeightedNormInfinity()`
4842: @*/
4843: PetscErrorCode TSErrorWeightedNorm2(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4844: {
4845: PetscInt i, n, N, rstart;
4846: PetscInt n_loc, na_loc, nr_loc;
4847: PetscReal n_glb, na_glb, nr_glb;
4848: const PetscScalar *u, *y;
4849: PetscReal sum, suma, sumr, gsum, gsuma, gsumr, diff;
4850: PetscReal tol, tola, tolr;
4851: PetscReal err_loc[6], err_glb[6];
4853: PetscFunctionBegin;
4859: PetscCheckSameComm(U, 2, Y, 3);
4863: PetscCheck(U != Y, PetscObjectComm((PetscObject)U), PETSC_ERR_ARG_IDN, "U and Y cannot be the same vector");
4865: PetscCall(VecGetSize(U, &N));
4866: PetscCall(VecGetLocalSize(U, &n));
4867: PetscCall(VecGetOwnershipRange(U, &rstart, NULL));
4868: PetscCall(VecGetArrayRead(U, &u));
4869: PetscCall(VecGetArrayRead(Y, &y));
4870: sum = 0.;
4871: n_loc = 0;
4872: suma = 0.;
4873: na_loc = 0;
4874: sumr = 0.;
4875: nr_loc = 0;
4876: if (ts->vatol && ts->vrtol) {
4877: const PetscScalar *atol, *rtol;
4878: PetscCall(VecGetArrayRead(ts->vatol, &atol));
4879: PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
4880: for (i = 0; i < n; i++) {
4881: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4882: diff = PetscAbsScalar(y[i] - u[i]);
4883: tola = PetscRealPart(atol[i]);
4884: if (tola > 0.) {
4885: suma += PetscSqr(diff / tola);
4886: na_loc++;
4887: }
4888: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4889: if (tolr > 0.) {
4890: sumr += PetscSqr(diff / tolr);
4891: nr_loc++;
4892: }
4893: tol = tola + tolr;
4894: if (tol > 0.) {
4895: sum += PetscSqr(diff / tol);
4896: n_loc++;
4897: }
4898: }
4899: PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
4900: PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
4901: } else if (ts->vatol) { /* vector atol, scalar rtol */
4902: const PetscScalar *atol;
4903: PetscCall(VecGetArrayRead(ts->vatol, &atol));
4904: for (i = 0; i < n; i++) {
4905: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4906: diff = PetscAbsScalar(y[i] - u[i]);
4907: tola = PetscRealPart(atol[i]);
4908: if (tola > 0.) {
4909: suma += PetscSqr(diff / tola);
4910: na_loc++;
4911: }
4912: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4913: if (tolr > 0.) {
4914: sumr += PetscSqr(diff / tolr);
4915: nr_loc++;
4916: }
4917: tol = tola + tolr;
4918: if (tol > 0.) {
4919: sum += PetscSqr(diff / tol);
4920: n_loc++;
4921: }
4922: }
4923: PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
4924: } else if (ts->vrtol) { /* scalar atol, vector rtol */
4925: const PetscScalar *rtol;
4926: PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
4927: for (i = 0; i < n; i++) {
4928: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4929: diff = PetscAbsScalar(y[i] - u[i]);
4930: tola = ts->atol;
4931: if (tola > 0.) {
4932: suma += PetscSqr(diff / tola);
4933: na_loc++;
4934: }
4935: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4936: if (tolr > 0.) {
4937: sumr += PetscSqr(diff / tolr);
4938: nr_loc++;
4939: }
4940: tol = tola + tolr;
4941: if (tol > 0.) {
4942: sum += PetscSqr(diff / tol);
4943: n_loc++;
4944: }
4945: }
4946: PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
4947: } else { /* scalar atol, scalar rtol */
4948: for (i = 0; i < n; i++) {
4949: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4950: diff = PetscAbsScalar(y[i] - u[i]);
4951: tola = ts->atol;
4952: if (tola > 0.) {
4953: suma += PetscSqr(diff / tola);
4954: na_loc++;
4955: }
4956: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4957: if (tolr > 0.) {
4958: sumr += PetscSqr(diff / tolr);
4959: nr_loc++;
4960: }
4961: tol = tola + tolr;
4962: if (tol > 0.) {
4963: sum += PetscSqr(diff / tol);
4964: n_loc++;
4965: }
4966: }
4967: }
4968: PetscCall(VecRestoreArrayRead(U, &u));
4969: PetscCall(VecRestoreArrayRead(Y, &y));
4971: err_loc[0] = sum;
4972: err_loc[1] = suma;
4973: err_loc[2] = sumr;
4974: err_loc[3] = (PetscReal)n_loc;
4975: err_loc[4] = (PetscReal)na_loc;
4976: err_loc[5] = (PetscReal)nr_loc;
4978: PetscCall(MPIU_Allreduce(err_loc, err_glb, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
4980: gsum = err_glb[0];
4981: gsuma = err_glb[1];
4982: gsumr = err_glb[2];
4983: n_glb = err_glb[3];
4984: na_glb = err_glb[4];
4985: nr_glb = err_glb[5];
4987: *norm = 0.;
4988: if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
4989: *norma = 0.;
4990: if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
4991: *normr = 0.;
4992: if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);
4994: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
4995: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
4996: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
4997: PetscFunctionReturn(PETSC_SUCCESS);
4998: }
5000: /*@
5001: TSErrorWeightedNormInfinity - compute a weighted infinity-norm of the difference between two state vectors
5003: Collective
5005: Input Parameters:
5006: + ts - time stepping context
5007: . U - state vector, usually ts->vec_sol
5008: - Y - state vector to be compared to U
5010: Output Parameters:
5011: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5012: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5013: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
5015: Level: developer
5017: .seealso: [](ch_ts), `TS`, `TSErrorWeightedNorm()`, `TSErrorWeightedNorm2()`
5018: @*/
5019: PetscErrorCode TSErrorWeightedNormInfinity(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5020: {
5021: PetscInt i, n, N, rstart;
5022: const PetscScalar *u, *y;
5023: PetscReal max, gmax, maxa, gmaxa, maxr, gmaxr;
5024: PetscReal tol, tola, tolr, diff;
5025: PetscReal err_loc[3], err_glb[3];
5027: PetscFunctionBegin;
5033: PetscCheckSameComm(U, 2, Y, 3);
5037: PetscCheck(U != Y, PetscObjectComm((PetscObject)U), PETSC_ERR_ARG_IDN, "U and Y cannot be the same vector");
5039: PetscCall(VecGetSize(U, &N));
5040: PetscCall(VecGetLocalSize(U, &n));
5041: PetscCall(VecGetOwnershipRange(U, &rstart, NULL));
5042: PetscCall(VecGetArrayRead(U, &u));
5043: PetscCall(VecGetArrayRead(Y, &y));
5045: max = 0.;
5046: maxa = 0.;
5047: maxr = 0.;
5049: if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
5050: const PetscScalar *atol, *rtol;
5051: PetscCall(VecGetArrayRead(ts->vatol, &atol));
5052: PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
5054: for (i = 0; i < n; i++) {
5055: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5056: diff = PetscAbsScalar(y[i] - u[i]);
5057: tola = PetscRealPart(atol[i]);
5058: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5059: tol = tola + tolr;
5060: if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
5061: if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
5062: if (tol > 0.) max = PetscMax(max, diff / tol);
5063: }
5064: PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5065: PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5066: } else if (ts->vatol) { /* vector atol, scalar rtol */
5067: const PetscScalar *atol;
5068: PetscCall(VecGetArrayRead(ts->vatol, &atol));
5069: for (i = 0; i < n; i++) {
5070: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5071: diff = PetscAbsScalar(y[i] - u[i]);
5072: tola = PetscRealPart(atol[i]);
5073: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5074: tol = tola + tolr;
5075: if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
5076: if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
5077: if (tol > 0.) max = PetscMax(max, diff / tol);
5078: }
5079: PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5080: } else if (ts->vrtol) { /* scalar atol, vector rtol */
5081: const PetscScalar *rtol;
5082: PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
5084: for (i = 0; i < n; i++) {
5085: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5086: diff = PetscAbsScalar(y[i] - u[i]);
5087: tola = ts->atol;
5088: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5089: tol = tola + tolr;
5090: if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
5091: if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
5092: if (tol > 0.) max = PetscMax(max, diff / tol);
5093: }
5094: PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5095: } else { /* scalar atol, scalar rtol */
5097: for (i = 0; i < n; i++) {
5098: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5099: diff = PetscAbsScalar(y[i] - u[i]);
5100: tola = ts->atol;
5101: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5102: tol = tola + tolr;
5103: if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
5104: if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
5105: if (tol > 0.) max = PetscMax(max, diff / tol);
5106: }
5107: }
5108: PetscCall(VecRestoreArrayRead(U, &u));
5109: PetscCall(VecRestoreArrayRead(Y, &y));
5110: err_loc[0] = max;
5111: err_loc[1] = maxa;
5112: err_loc[2] = maxr;
5113: PetscCall(MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
5114: gmax = err_glb[0];
5115: gmaxa = err_glb[1];
5116: gmaxr = err_glb[2];
5118: *norm = gmax;
5119: *norma = gmaxa;
5120: *normr = gmaxr;
5121: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5122: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5123: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5124: PetscFunctionReturn(PETSC_SUCCESS);
5125: }
5127: /*@
5128: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
5130: Collective
5132: Input Parameters:
5133: + ts - time stepping context
5134: . U - state vector, usually ts->vec_sol
5135: . Y - state vector to be compared to U
5136: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5138: Output Parameters:
5139: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5140: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5141: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5143: Options Database Key:
5144: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5146: Level: developer
5148: .seealso: [](ch_ts), `TS`, `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`, `TSErrorWeightedENorm`
5149: @*/
5150: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5151: {
5152: PetscFunctionBegin;
5153: if (wnormtype == NORM_2) PetscCall(TSErrorWeightedNorm2(ts, U, Y, norm, norma, normr));
5154: else if (wnormtype == NORM_INFINITY) PetscCall(TSErrorWeightedNormInfinity(ts, U, Y, norm, norma, normr));
5155: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5156: PetscFunctionReturn(PETSC_SUCCESS);
5157: }
5159: /*@
5160: TSErrorWeightedENorm2 - compute a weighted 2 error norm based on supplied absolute and relative tolerances
5162: Collective
5164: Input Parameters:
5165: + ts - time stepping context
5166: . E - error vector
5167: . U - state vector, usually ts->vec_sol
5168: - Y - state vector, previous time step
5170: Output Parameters:
5171: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5172: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5173: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
5175: Level: developer
5177: .seealso: [](ch_ts), `TS`, `TSErrorWeightedENorm()`, `TSErrorWeightedENormInfinity()`
5178: @*/
5179: PetscErrorCode TSErrorWeightedENorm2(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5180: {
5181: PetscInt i, n, N, rstart;
5182: PetscInt n_loc, na_loc, nr_loc;
5183: PetscReal n_glb, na_glb, nr_glb;
5184: const PetscScalar *e, *u, *y;
5185: PetscReal err, sum, suma, sumr, gsum, gsuma, gsumr;
5186: PetscReal tol, tola, tolr;
5187: PetscReal err_loc[6], err_glb[6];
5189: PetscFunctionBegin;
5197: PetscCheckSameComm(E, 2, U, 3);
5198: PetscCheckSameComm(U, 3, Y, 4);
5203: PetscCall(VecGetSize(E, &N));
5204: PetscCall(VecGetLocalSize(E, &n));
5205: PetscCall(VecGetOwnershipRange(E, &rstart, NULL));
5206: PetscCall(VecGetArrayRead(E, &e));
5207: PetscCall(VecGetArrayRead(U, &u));
5208: PetscCall(VecGetArrayRead(Y, &y));
5209: sum = 0.;
5210: n_loc = 0;
5211: suma = 0.;
5212: na_loc = 0;
5213: sumr = 0.;
5214: nr_loc = 0;
5215: if (ts->vatol && ts->vrtol) {
5216: const PetscScalar *atol, *rtol;
5217: PetscCall(VecGetArrayRead(ts->vatol, &atol));
5218: PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
5219: for (i = 0; i < n; i++) {
5220: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5221: err = PetscAbsScalar(e[i]);
5222: tola = PetscRealPart(atol[i]);
5223: if (tola > 0.) {
5224: suma += PetscSqr(err / tola);
5225: na_loc++;
5226: }
5227: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5228: if (tolr > 0.) {
5229: sumr += PetscSqr(err / tolr);
5230: nr_loc++;
5231: }
5232: tol = tola + tolr;
5233: if (tol > 0.) {
5234: sum += PetscSqr(err / tol);
5235: n_loc++;
5236: }
5237: }
5238: PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5239: PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5240: } else if (ts->vatol) { /* vector atol, scalar rtol */
5241: const PetscScalar *atol;
5242: PetscCall(VecGetArrayRead(ts->vatol, &atol));
5243: for (i = 0; i < n; i++) {
5244: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5245: err = PetscAbsScalar(e[i]);
5246: tola = PetscRealPart(atol[i]);
5247: if (tola > 0.) {
5248: suma += PetscSqr(err / tola);
5249: na_loc++;
5250: }
5251: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5252: if (tolr > 0.) {
5253: sumr += PetscSqr(err / tolr);
5254: nr_loc++;
5255: }
5256: tol = tola + tolr;
5257: if (tol > 0.) {
5258: sum += PetscSqr(err / tol);
5259: n_loc++;
5260: }
5261: }
5262: PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5263: } else if (ts->vrtol) { /* scalar atol, vector rtol */
5264: const PetscScalar *rtol;
5265: PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
5266: for (i = 0; i < n; i++) {
5267: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5268: err = PetscAbsScalar(e[i]);
5269: tola = ts->atol;
5270: if (tola > 0.) {
5271: suma += PetscSqr(err / tola);
5272: na_loc++;
5273: }
5274: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5275: if (tolr > 0.) {
5276: sumr += PetscSqr(err / tolr);
5277: nr_loc++;
5278: }
5279: tol = tola + tolr;
5280: if (tol > 0.) {
5281: sum += PetscSqr(err / tol);
5282: n_loc++;
5283: }
5284: }
5285: PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5286: } else { /* scalar atol, scalar rtol */
5287: for (i = 0; i < n; i++) {
5288: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5289: err = PetscAbsScalar(e[i]);
5290: tola = ts->atol;
5291: if (tola > 0.) {
5292: suma += PetscSqr(err / tola);
5293: na_loc++;
5294: }
5295: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5296: if (tolr > 0.) {
5297: sumr += PetscSqr(err / tolr);
5298: nr_loc++;
5299: }
5300: tol = tola + tolr;
5301: if (tol > 0.) {
5302: sum += PetscSqr(err / tol);
5303: n_loc++;
5304: }
5305: }
5306: }
5307: PetscCall(VecRestoreArrayRead(E, &e));
5308: PetscCall(VecRestoreArrayRead(U, &u));
5309: PetscCall(VecRestoreArrayRead(Y, &y));
5311: err_loc[0] = sum;
5312: err_loc[1] = suma;
5313: err_loc[2] = sumr;
5314: err_loc[3] = (PetscReal)n_loc;
5315: err_loc[4] = (PetscReal)na_loc;
5316: err_loc[5] = (PetscReal)nr_loc;
5318: PetscCall(MPIU_Allreduce(err_loc, err_glb, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
5320: gsum = err_glb[0];
5321: gsuma = err_glb[1];
5322: gsumr = err_glb[2];
5323: n_glb = err_glb[3];
5324: na_glb = err_glb[4];
5325: nr_glb = err_glb[5];
5327: *norm = 0.;
5328: if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
5329: *norma = 0.;
5330: if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
5331: *normr = 0.;
5332: if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);
5334: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5335: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5336: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5337: PetscFunctionReturn(PETSC_SUCCESS);
5338: }
5340: /*@
5341: TSErrorWeightedENormInfinity - compute a weighted infinity error norm based on supplied absolute and relative tolerances
5343: Collective
5345: Input Parameters:
5346: + ts - time stepping context
5347: . E - error vector
5348: . U - state vector, usually ts->vec_sol
5349: - Y - state vector, previous time step
5351: Output Parameters:
5352: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5353: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5354: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
5356: Level: developer
5358: .seealso: [](ch_ts), `TS`, `TSErrorWeightedENorm()`, `TSErrorWeightedENorm2()`
5359: @*/
5360: PetscErrorCode TSErrorWeightedENormInfinity(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5361: {
5362: PetscInt i, n, N, rstart;
5363: const PetscScalar *e, *u, *y;
5364: PetscReal err, max, gmax, maxa, gmaxa, maxr, gmaxr;
5365: PetscReal tol, tola, tolr;
5366: PetscReal err_loc[3], err_glb[3];
5368: PetscFunctionBegin;
5376: PetscCheckSameComm(E, 2, U, 3);
5377: PetscCheckSameComm(U, 3, Y, 4);
5382: PetscCall(VecGetSize(E, &N));
5383: PetscCall(VecGetLocalSize(E, &n));
5384: PetscCall(VecGetOwnershipRange(E, &rstart, NULL));
5385: PetscCall(VecGetArrayRead(E, &e));
5386: PetscCall(VecGetArrayRead(U, &u));
5387: PetscCall(VecGetArrayRead(Y, &y));
5389: max = 0.;
5390: maxa = 0.;
5391: maxr = 0.;
5393: if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
5394: const PetscScalar *atol, *rtol;
5395: PetscCall(VecGetArrayRead(ts->vatol, &atol));
5396: PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
5398: for (i = 0; i < n; i++) {
5399: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5400: err = PetscAbsScalar(e[i]);
5401: tola = PetscRealPart(atol[i]);
5402: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5403: tol = tola + tolr;
5404: if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5405: if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5406: if (tol > 0.) max = PetscMax(max, err / tol);
5407: }
5408: PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5409: PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5410: } else if (ts->vatol) { /* vector atol, scalar rtol */
5411: const PetscScalar *atol;
5412: PetscCall(VecGetArrayRead(ts->vatol, &atol));
5413: for (i = 0; i < n; i++) {
5414: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5415: err = PetscAbsScalar(e[i]);
5416: tola = PetscRealPart(atol[i]);
5417: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5418: tol = tola + tolr;
5419: if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5420: if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5421: if (tol > 0.) max = PetscMax(max, err / tol);
5422: }
5423: PetscCall(VecRestoreArrayRead(ts->vatol, &atol));
5424: } else if (ts->vrtol) { /* scalar atol, vector rtol */
5425: const PetscScalar *rtol;
5426: PetscCall(VecGetArrayRead(ts->vrtol, &rtol));
5428: for (i = 0; i < n; i++) {
5429: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5430: err = PetscAbsScalar(e[i]);
5431: tola = ts->atol;
5432: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5433: tol = tola + tolr;
5434: if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5435: if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5436: if (tol > 0.) max = PetscMax(max, err / tol);
5437: }
5438: PetscCall(VecRestoreArrayRead(ts->vrtol, &rtol));
5439: } else { /* scalar atol, scalar rtol */
5441: for (i = 0; i < n; i++) {
5442: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5443: err = PetscAbsScalar(e[i]);
5444: tola = ts->atol;
5445: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5446: tol = tola + tolr;
5447: if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5448: if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5449: if (tol > 0.) max = PetscMax(max, err / tol);
5450: }
5451: }
5452: PetscCall(VecRestoreArrayRead(E, &e));
5453: PetscCall(VecRestoreArrayRead(U, &u));
5454: PetscCall(VecRestoreArrayRead(Y, &y));
5455: err_loc[0] = max;
5456: err_loc[1] = maxa;
5457: err_loc[2] = maxr;
5458: PetscCall(MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
5459: gmax = err_glb[0];
5460: gmaxa = err_glb[1];
5461: gmaxr = err_glb[2];
5463: *norm = gmax;
5464: *norma = gmaxa;
5465: *normr = gmaxr;
5466: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5467: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5468: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5469: PetscFunctionReturn(PETSC_SUCCESS);
5470: }
5472: /*@
5473: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5475: Collective
5477: Input Parameters:
5478: + ts - time stepping context
5479: . E - error vector
5480: . U - state vector, usually ts->vec_sol
5481: . Y - state vector, previous time step
5482: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5484: Output Parameters:
5485: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5486: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5487: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5489: Options Database Key:
5490: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5492: Level: developer
5494: .seealso: [](ch_ts), `TS`, `TSErrorWeightedENormInfinity()`, `TSErrorWeightedENorm2()`, `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`
5495: @*/
5496: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5497: {
5498: PetscFunctionBegin;
5499: if (wnormtype == NORM_2) PetscCall(TSErrorWeightedENorm2(ts, E, U, Y, norm, norma, normr));
5500: else if (wnormtype == NORM_INFINITY) PetscCall(TSErrorWeightedENormInfinity(ts, E, U, Y, norm, norma, normr));
5501: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5502: PetscFunctionReturn(PETSC_SUCCESS);
5503: }
5505: /*@
5506: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5508: Logically Collective
5510: Input Parameters:
5511: + ts - time stepping context
5512: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5514: Note:
5515: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5517: Level: intermediate
5519: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5520: @*/
5521: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5522: {
5523: PetscFunctionBegin;
5525: ts->cfltime_local = cfltime;
5526: ts->cfltime = -1.;
5527: PetscFunctionReturn(PETSC_SUCCESS);
5528: }
5530: /*@
5531: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5533: Collective
5535: Input Parameter:
5536: . ts - time stepping context
5538: Output Parameter:
5539: . cfltime - maximum stable time step for forward Euler
5541: Level: advanced
5543: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5544: @*/
5545: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5546: {
5547: PetscFunctionBegin;
5548: if (ts->cfltime < 0) PetscCall(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5549: *cfltime = ts->cfltime;
5550: PetscFunctionReturn(PETSC_SUCCESS);
5551: }
5553: /*@
5554: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5556: Input Parameters:
5557: + ts - the `TS` context.
5558: . xl - lower bound.
5559: - xu - upper bound.
5561: Level: advanced
5563: Note:
5564: If this routine is not called then the lower and upper bounds are set to
5565: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5567: .seealso: [](ch_ts), `TS`
5568: @*/
5569: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5570: {
5571: SNES snes;
5573: PetscFunctionBegin;
5574: PetscCall(TSGetSNES(ts, &snes));
5575: PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5576: PetscFunctionReturn(PETSC_SUCCESS);
5577: }
5579: /*@
5580: TSComputeLinearStability - computes the linear stability function at a point
5582: Collective
5584: Input Parameters:
5585: + ts - the `TS` context
5586: . xr - real part of input argument
5587: - xi - imaginary part of input argument
5589: Output Parameters:
5590: + yr - real part of function value
5591: - yi - imaginary part of function value
5593: Level: developer
5595: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5596: @*/
5597: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5598: {
5599: PetscFunctionBegin;
5601: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5602: PetscFunctionReturn(PETSC_SUCCESS);
5603: }
5605: /*@
5606: TSRestartStep - Flags the solver to restart the next step
5608: Collective
5610: Input Parameter:
5611: . ts - the `TS` context obtained from `TSCreate()`
5613: Level: advanced
5615: Notes:
5616: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5617: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5618: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5619: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5620: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5621: discontinuous source terms).
5623: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5624: @*/
5625: PetscErrorCode TSRestartStep(TS ts)
5626: {
5627: PetscFunctionBegin;
5629: ts->steprestart = PETSC_TRUE;
5630: PetscFunctionReturn(PETSC_SUCCESS);
5631: }
5633: /*@
5634: TSRollBack - Rolls back one time step
5636: Collective
5638: Input Parameter:
5639: . ts - the `TS` context obtained from `TSCreate()`
5641: Level: advanced
5643: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5644: @*/
5645: PetscErrorCode TSRollBack(TS ts)
5646: {
5647: PetscFunctionBegin;
5649: PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5650: PetscUseTypeMethod(ts, rollback);
5651: ts->time_step = ts->ptime - ts->ptime_prev;
5652: ts->ptime = ts->ptime_prev;
5653: ts->ptime_prev = ts->ptime_prev_rollback;
5654: ts->steps--;
5655: ts->steprollback = PETSC_TRUE;
5656: PetscFunctionReturn(PETSC_SUCCESS);
5657: }
5659: /*@
5660: TSGetStages - Get the number of stages and stage values
5662: Input Parameter:
5663: . ts - the `TS` context obtained from `TSCreate()`
5665: Output Parameters:
5666: + ns - the number of stages
5667: - Y - the current stage vectors
5669: Level: advanced
5671: Note:
5672: Both `ns` and `Y` can be `NULL`.
5674: .seealso: [](ch_ts), `TS`, `TSCreate()`
5675: @*/
5676: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5677: {
5678: PetscFunctionBegin;
5682: if (!ts->ops->getstages) {
5683: if (ns) *ns = 0;
5684: if (Y) *Y = NULL;
5685: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5686: PetscFunctionReturn(PETSC_SUCCESS);
5687: }
5689: /*@C
5690: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5692: Collective
5694: Input Parameters:
5695: + ts - the `TS` context
5696: . t - current timestep
5697: . U - state vector
5698: . Udot - time derivative of state vector
5699: . shift - shift to apply, see note below
5700: - ctx - an optional user context
5702: Output Parameters:
5703: + J - Jacobian matrix (not altered in this routine)
5704: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
5706: Level: intermediate
5708: Notes:
5709: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5711: dF/dU + shift*dF/dUdot
5713: Most users should not need to explicitly call this routine, as it
5714: is used internally within the nonlinear solvers.
5716: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5717: routine, then it will try to get the coloring from the matrix. This requires that the
5718: matrix have nonzero entries precomputed.
5720: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5721: @*/
5722: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5723: {
5724: SNES snes;
5725: MatFDColoring color;
5726: PetscBool hascolor, matcolor = PETSC_FALSE;
5728: PetscFunctionBegin;
5729: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5730: PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5731: if (!color) {
5732: DM dm;
5733: ISColoring iscoloring;
5735: PetscCall(TSGetDM(ts, &dm));
5736: PetscCall(DMHasColoring(dm, &hascolor));
5737: if (hascolor && !matcolor) {
5738: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5739: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5740: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5741: PetscCall(MatFDColoringSetFromOptions(color));
5742: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5743: PetscCall(ISColoringDestroy(&iscoloring));
5744: } else {
5745: MatColoring mc;
5747: PetscCall(MatColoringCreate(B, &mc));
5748: PetscCall(MatColoringSetDistance(mc, 2));
5749: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5750: PetscCall(MatColoringSetFromOptions(mc));
5751: PetscCall(MatColoringApply(mc, &iscoloring));
5752: PetscCall(MatColoringDestroy(&mc));
5753: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5754: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
5755: PetscCall(MatFDColoringSetFromOptions(color));
5756: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5757: PetscCall(ISColoringDestroy(&iscoloring));
5758: }
5759: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5760: PetscCall(PetscObjectDereference((PetscObject)color));
5761: }
5762: PetscCall(TSGetSNES(ts, &snes));
5763: PetscCall(MatFDColoringApply(B, color, U, snes));
5764: if (J != B) {
5765: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5766: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5767: }
5768: PetscFunctionReturn(PETSC_SUCCESS);
5769: }
5771: /*@
5772: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5774: Input Parameters:
5775: + ts - the `TS` context
5776: - func - function called within `TSFunctionDomainError()`
5778: Calling sequence of `func`:
5779: $ PetscErrorCode func(TS ts, PetscReal time, Vec state, PetscBool reject)
5780: + ts - the TS context
5781: . time - the current time (of the stage)
5782: . state - the state to check if it is valid
5783: - reject - (output parameter) `PETSC_FALSE` if the state is acceptable, `PETSC_TRUE` if not acceptable
5785: Level: intermediate
5787: Notes:
5788: If an implicit ODE solver is being used then, in addition to providing this routine, the
5789: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5790: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5791: Use `TSGetSNES()` to obtain the `SNES` object
5793: Developer Note:
5794: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5795: since one takes a function pointer and the other does not.
5797: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5798: @*/
5800: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS, PetscReal, Vec, PetscBool *))
5801: {
5802: PetscFunctionBegin;
5804: ts->functiondomainerror = func;
5805: PetscFunctionReturn(PETSC_SUCCESS);
5806: }
5808: /*@
5809: TSFunctionDomainError - Checks if the current state is valid
5811: Input Parameters:
5812: + ts - the `TS` context
5813: . stagetime - time of the simulation
5814: - Y - state vector to check.
5816: Output Parameter:
5817: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5819: Level: developer
5821: Note:
5822: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5823: to check if the current state is valid.
5825: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5826: @*/
5827: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5828: {
5829: PetscFunctionBegin;
5831: *accept = PETSC_TRUE;
5832: if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5833: PetscFunctionReturn(PETSC_SUCCESS);
5834: }
5836: /*@C
5837: TSClone - This function clones a time step `TS` object.
5839: Collective
5841: Input Parameter:
5842: . tsin - The input `TS`
5844: Output Parameter:
5845: . tsout - The output `TS` (cloned)
5847: Level: developer
5849: Notes:
5850: 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.
5851: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5853: 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
5854: .vb
5855: SNES snes_dup = NULL;
5856: TSGetSNES(ts,&snes_dup);
5857: TSSetSNES(ts,snes_dup);
5858: .ve
5860: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5861: @*/
5862: PetscErrorCode TSClone(TS tsin, TS *tsout)
5863: {
5864: TS t;
5865: SNES snes_start;
5866: DM dm;
5867: TSType type;
5869: PetscFunctionBegin;
5871: *tsout = NULL;
5873: PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));
5875: /* General TS description */
5876: t->numbermonitors = 0;
5877: t->monitorFrequency = 1;
5878: t->setupcalled = 0;
5879: t->ksp_its = 0;
5880: t->snes_its = 0;
5881: t->nwork = 0;
5882: t->rhsjacobian.time = PETSC_MIN_REAL;
5883: t->rhsjacobian.scale = 1.;
5884: t->ijacobian.shift = 1.;
5886: PetscCall(TSGetSNES(tsin, &snes_start));
5887: PetscCall(TSSetSNES(t, snes_start));
5889: PetscCall(TSGetDM(tsin, &dm));
5890: PetscCall(TSSetDM(t, dm));
5892: t->adapt = tsin->adapt;
5893: PetscCall(PetscObjectReference((PetscObject)t->adapt));
5895: t->trajectory = tsin->trajectory;
5896: PetscCall(PetscObjectReference((PetscObject)t->trajectory));
5898: t->event = tsin->event;
5899: if (t->event) t->event->refct++;
5901: t->problem_type = tsin->problem_type;
5902: t->ptime = tsin->ptime;
5903: t->ptime_prev = tsin->ptime_prev;
5904: t->time_step = tsin->time_step;
5905: t->max_time = tsin->max_time;
5906: t->steps = tsin->steps;
5907: t->max_steps = tsin->max_steps;
5908: t->equation_type = tsin->equation_type;
5909: t->atol = tsin->atol;
5910: t->rtol = tsin->rtol;
5911: t->max_snes_failures = tsin->max_snes_failures;
5912: t->max_reject = tsin->max_reject;
5913: t->errorifstepfailed = tsin->errorifstepfailed;
5915: PetscCall(TSGetType(tsin, &type));
5916: PetscCall(TSSetType(t, type));
5918: t->vec_sol = NULL;
5920: t->cfltime = tsin->cfltime;
5921: t->cfltime_local = tsin->cfltime_local;
5922: t->exact_final_time = tsin->exact_final_time;
5924: PetscCall(PetscMemcpy(t->ops, tsin->ops, sizeof(struct _TSOps)));
5926: if (((PetscObject)tsin)->fortran_func_pointers) {
5927: PetscInt i;
5928: PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers));
5929: for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5930: }
5931: *tsout = t;
5932: PetscFunctionReturn(PETSC_SUCCESS);
5933: }
5935: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5936: {
5937: TS ts = (TS)ctx;
5939: PetscFunctionBegin;
5940: PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5941: PetscFunctionReturn(PETSC_SUCCESS);
5942: }
5944: /*@
5945: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5947: Logically Collective
5949: Input Parameter:
5950: . TS - the time stepping routine
5952: Output Parameter:
5953: . flg - `PETSC_TRUE` if the multiply is likely correct
5955: Options Database Key:
5956: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5958: Level: advanced
5960: Note:
5961: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5963: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5964: @*/
5965: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5966: {
5967: Mat J, B;
5968: TSRHSJacobian func;
5969: void *ctx;
5971: PetscFunctionBegin;
5972: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5973: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5974: PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5975: PetscFunctionReturn(PETSC_SUCCESS);
5976: }
5978: /*@C
5979: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5981: Logically Collective
5983: Input Parameter:
5984: . TS - the time stepping routine
5986: Output Parameter:
5987: . flg - `PETSC_TRUE` if the multiply is likely correct
5989: Options Database Key:
5990: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5992: Level: advanced
5994: Notes:
5995: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5997: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5998: @*/
5999: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
6000: {
6001: Mat J, B;
6002: void *ctx;
6003: TSRHSJacobian func;
6005: PetscFunctionBegin;
6006: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
6007: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
6008: PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
6009: PetscFunctionReturn(PETSC_SUCCESS);
6010: }
6012: /*@
6013: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
6015: Logically Collective
6017: Input Parameters:
6018: + ts - timestepping context
6019: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
6021: Options Database Key:
6022: . -ts_use_splitrhsfunction - <true,false>
6024: Level: intermediate
6026: Note:
6027: This is only for multirate methods
6029: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
6030: @*/
6031: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
6032: {
6033: PetscFunctionBegin;
6035: ts->use_splitrhsfunction = use_splitrhsfunction;
6036: PetscFunctionReturn(PETSC_SUCCESS);
6037: }
6039: /*@
6040: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
6042: Not Collective
6044: Input Parameter:
6045: . ts - timestepping context
6047: Output Parameter:
6048: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
6050: Level: intermediate
6052: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
6053: @*/
6054: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
6055: {
6056: PetscFunctionBegin;
6058: *use_splitrhsfunction = ts->use_splitrhsfunction;
6059: PetscFunctionReturn(PETSC_SUCCESS);
6060: }
6062: /*@
6063: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
6065: Logically Collective
6067: Input Parameters:
6068: + ts - the time-stepper
6069: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
6071: Level: intermediate
6073: Note:
6074: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
6076: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
6077: @*/
6078: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
6079: {
6080: PetscFunctionBegin;
6082: ts->axpy_pattern = str;
6083: PetscFunctionReturn(PETSC_SUCCESS);
6084: }
6086: /*@
6087: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
6089: Collective
6091: Input Parameters:
6092: + ts - the time-stepper
6093: . n - number of the time points (>=2)
6094: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
6096: Options Database Key:
6097: . -ts_time_span <t0,...tf> - Sets the time span
6099: Level: intermediate
6101: Notes:
6102: The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
6103: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
6104: The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may
6105: pressure the memory system when using a large number of span points.
6107: .seealso: [](ch_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()`
6108: @*/
6109: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
6110: {
6111: PetscFunctionBegin;
6113: PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
6114: if (ts->tspan && n != ts->tspan->num_span_times) {
6115: PetscCall(PetscFree(ts->tspan->span_times));
6116: PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
6117: PetscCall(PetscMalloc1(n, &ts->tspan->span_times));
6118: }
6119: if (!ts->tspan) {
6120: TSTimeSpan tspan;
6121: PetscCall(PetscNew(&tspan));
6122: PetscCall(PetscMalloc1(n, &tspan->span_times));
6123: tspan->reltol = 1e-6;
6124: tspan->abstol = 10 * PETSC_MACHINE_EPSILON;
6125: ts->tspan = tspan;
6126: }
6127: ts->tspan->num_span_times = n;
6128: PetscCall(PetscArraycpy(ts->tspan->span_times, span_times, n));
6129: PetscCall(TSSetTime(ts, ts->tspan->span_times[0]));
6130: PetscCall(TSSetMaxTime(ts, ts->tspan->span_times[n - 1]));
6131: PetscFunctionReturn(PETSC_SUCCESS);
6132: }
6134: /*@C
6135: TSGetTimeSpan - gets the time span set with `TSSetTimeSpan()`
6137: Not Collective
6139: Input Parameter:
6140: . ts - the time-stepper
6142: Output Parameters:
6143: + n - number of the time points (>=2)
6144: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
6146: Level: beginner
6148: Note:
6149: The values obtained are valid until the `TS` object is destroyed.
6151: Both `n` and `span_times` can be `NULL`.
6153: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()`
6154: @*/
6155: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal **span_times)
6156: {
6157: PetscFunctionBegin;
6161: if (!ts->tspan) {
6162: if (n) *n = 0;
6163: if (span_times) *span_times = NULL;
6164: } else {
6165: if (n) *n = ts->tspan->num_span_times;
6166: if (span_times) *span_times = ts->tspan->span_times;
6167: }
6168: PetscFunctionReturn(PETSC_SUCCESS);
6169: }
6171: /*@
6172: TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.
6174: Input Parameter:
6175: . ts - the `TS` context obtained from `TSCreate()`
6177: Output Parameters:
6178: + nsol - the number of solutions
6179: - Sols - the solution vectors
6181: Level: intermediate
6183: Notes:
6184: Both `nsol` and `Sols` can be `NULL`.
6186: 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()`.
6187: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span.
6189: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`
6190: @*/
6191: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
6192: {
6193: PetscFunctionBegin;
6197: if (!ts->tspan) {
6198: if (nsol) *nsol = 0;
6199: if (Sols) *Sols = NULL;
6200: } else {
6201: if (nsol) *nsol = ts->tspan->spanctr;
6202: if (Sols) *Sols = ts->tspan->vecs_sol;
6203: }
6204: PetscFunctionReturn(PETSC_SUCCESS);
6205: }
6207: /*@C
6208: TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
6210: Collective
6212: Input Parameters:
6213: + ts - the `TS` context
6214: . J - Jacobian matrix (not altered in this routine)
6215: - B - newly computed Jacobian matrix to use with preconditioner
6217: Level: intermediate
6219: Notes:
6220: This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
6221: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
6222: and multiple fields are involved.
6224: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
6225: structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
6226: usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
6227: `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.
6229: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
6230: @*/
6231: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
6232: {
6233: MatColoring mc = NULL;
6234: ISColoring iscoloring = NULL;
6235: MatFDColoring matfdcoloring = NULL;
6237: PetscFunctionBegin;
6238: /* Generate new coloring after eliminating zeros in the matrix */
6239: PetscCall(MatEliminateZeros(B));
6240: PetscCall(MatColoringCreate(B, &mc));
6241: PetscCall(MatColoringSetDistance(mc, 2));
6242: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
6243: PetscCall(MatColoringSetFromOptions(mc));
6244: PetscCall(MatColoringApply(mc, &iscoloring));
6245: PetscCall(MatColoringDestroy(&mc));
6246: /* Replace the old coloring with the new one */
6247: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
6248: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts));
6249: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
6250: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
6251: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
6252: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
6253: PetscCall(ISColoringDestroy(&iscoloring));
6254: PetscFunctionReturn(PETSC_SUCCESS);
6255: }