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 <petscdraw.h>
7: #include <petscconvest.h>
9: /* Logging support */
10: PetscClassId TS_CLASSID, DMTS_CLASSID;
11: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;
13: const char *const TSExactFinalTimeOptions[] = {"UNSPECIFIED", "STEPOVER", "INTERPOLATE", "MATCHSTEP", "TSExactFinalTimeOption", "TS_EXACTFINALTIME_", NULL};
15: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt, TSAdaptType default_type)
16: {
17: PetscFunctionBegin;
19: PetscAssertPointer(default_type, 2);
20: if (!((PetscObject)adapt)->type_name) PetscCall(TSAdaptSetType(adapt, default_type));
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: /*@
25: TSSetFromOptions - Sets various `TS` parameters from the options database
27: Collective
29: Input Parameter:
30: . ts - the `TS` context obtained from `TSCreate()`
32: Options Database Keys:
33: + -ts_type type - see `TSType`
34: . -ts_save_trajectory - checkpoint the solution at each time-step
35: . -ts_max_time time - maximum time to compute to
36: . -ts_time_span t0,...,tf - sets the time span, solutions are computed and stored for each indicated time, init_time and max_time are set
37: . -ts_eval_times t0,...,tn - time points where solutions are computed and stored for each indicated time
38: . -ts_max_steps steps - maximum time-step number to execute until (possibly with nonzero starting value)
39: . -ts_run_steps steps - maximum number of time steps for `TSSolve()` to take on each call
40: . -ts_init_time time - initial time to start computation
41: . -ts_final_time time - final time to compute to (deprecated: use `-ts_max_time`)
42: . -ts_time_step dt - initial time step (only a suggestion, the actual initial time step used differ)
43: . -ts_exact_final_time (stepover,interpolate,matchstep) - whether to stop at the exact given final time and how to compute the solution at that time
44: . -ts_max_snes_failures maxfailures - Maximum number of nonlinear solve failures allowed
45: . -ts_max_step_rejections maxrejects - Maximum number of step rejections before step fails
46: . -ts_error_if_step_fails (true|false) - Error if no step succeeds
47: . -ts_rtol rtol - relative tolerance for local truncation error
48: . -ts_atol atol - Absolute tolerance for local truncation error
49: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function
50: . -ts_rhs_jacobian_test_mult_transpose - test the Jacobian at each iteration against finite difference with RHS function
51: . -ts_adjoint_solve (true|false) - After solving the ODE/DAE solve the adjoint problem (requires `-ts_save_trajectory`)
52: . -ts_fd_color - Use finite differences with coloring to compute IJacobian
53: . -ts_monitor - print information at each timestep
54: . -ts_monitor_cancel - Cancel all monitors
55: . -ts_monitor_wall_clock_time - Monitor wall-clock time, `KSP` iterations, and `SNES` iterations per step
56: . -ts_monitor_lg_solution - Monitor solution graphically
57: . -ts_monitor_lg_error - Monitor error graphically
58: . -ts_monitor_error - Monitors norm of error
59: . -ts_monitor_lg_timestep - Monitor timestep size graphically
60: . -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
61: . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
62: . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
63: . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
64: . -ts_monitor_draw_solution - Monitor solution graphically
65: . -ts_monitor_draw_solution_phase xleft,yleft,xright,yright - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
66: . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
67: . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
68: . -ts_monitor_solution_interval interval - output once every interval (default=1) time steps. Use -1 to only output at the end of the simulation
69: . -ts_monitor_solution_skip_initial - skip writing of initial condition
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_solution_vtk_interval interval - output once every interval (default=1) time steps. Use -1 to only output at the end of the simulation
72: - -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
74: Level: beginner
76: Notes:
77: See `SNESSetFromOptions()` and `KSPSetFromOptions()` for how to control the nonlinear and linear solves used by the time-stepper.
79: Certain `SNES` options get reset for each new nonlinear solver, for example `-snes_lag_jacobian its` and `-snes_lag_preconditioner its`, in order
80: to retain them over the multiple nonlinear solves that `TS` uses you must also provide `-snes_lag_jacobian_persists true` and
81: `-snes_lag_preconditioner_persists true`
83: Developer Notes:
84: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
86: .seealso: [](ch_ts), `TS`, `TSGetType()`
87: @*/
88: PetscErrorCode TSSetFromOptions(TS ts)
89: {
90: PetscBool opt, flg, tflg;
91: char monfilename[PETSC_MAX_PATH_LEN];
92: PetscReal time_step, eval_times[100] = {0};
93: PetscInt num_eval_times = PETSC_STATIC_ARRAY_LENGTH(eval_times);
94: TSExactFinalTimeOption eftopt;
95: char dir[16];
96: TSIFunctionFn *ifun;
97: const char *defaultType;
98: char typeName[256];
100: PetscFunctionBegin;
103: PetscCall(TSRegisterAll());
104: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
106: PetscObjectOptionsBegin((PetscObject)ts);
107: if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
108: else defaultType = ifun ? TSBEULER : TSEULER;
109: PetscCall(PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, sizeof(typeName), &opt));
110: if (opt) PetscCall(TSSetType(ts, typeName));
111: else PetscCall(TSSetType(ts, defaultType));
113: /* Handle generic TS options */
114: PetscCall(PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL));
115: PetscCall(PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL));
116: PetscCall(PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", eval_times, &num_eval_times, &flg));
117: if (flg) PetscCall(TSSetTimeSpan(ts, num_eval_times, eval_times));
118: num_eval_times = PETSC_STATIC_ARRAY_LENGTH(eval_times);
119: PetscCall(PetscOptionsRealArray("-ts_eval_times", "Evaluation time points", "TSSetEvaluationTimes", eval_times, &num_eval_times, &opt));
120: PetscCheck(flg != opt || (!flg && !opt), PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "May not provide -ts_time_span and -ts_eval_times simultaneously");
121: if (opt) PetscCall(TSSetEvaluationTimes(ts, num_eval_times, eval_times));
122: PetscCall(PetscOptionsInt("-ts_max_steps", "Maximum time step number to execute to (possibly with non-zero starting value)", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL));
123: PetscCall(PetscOptionsInt("-ts_run_steps", "Maximum number of time steps to take on each call to TSSolve()", "TSSetRunSteps", ts->run_steps, &ts->run_steps, NULL));
124: PetscCall(PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL));
125: PetscCall(PetscOptionsDeprecated("-ts_dt", "-ts_time_step", "3.25", NULL));
126: PetscCall(PetscOptionsReal("-ts_time_step", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg));
127: if (flg) PetscCall(TSSetTimeStep(ts, time_step));
128: PetscCall(PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg));
129: if (flg) PetscCall(TSSetExactFinalTime(ts, eftopt));
130: PetscCall(PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, &flg));
131: if (flg) PetscCall(TSSetMaxSNESFailures(ts, ts->max_snes_failures));
132: PetscCall(PetscOptionsDeprecated("-ts_max_reject", "-ts_max_step_rejections", "3.25", NULL));
133: PetscCall(PetscOptionsInt("-ts_max_step_rejections", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, &flg));
134: if (flg) PetscCall(TSSetMaxStepRejections(ts, ts->max_reject));
135: PetscCall(PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL));
136: PetscCall(PetscOptionsBoundedReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL, 0));
137: PetscCall(PetscOptionsBoundedReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL, 0));
139: PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL));
140: 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));
141: PetscCall(PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL));
142: #if PetscDefined(HAVE_SAWS)
143: {
144: PetscBool set;
145: flg = PETSC_FALSE;
146: PetscCall(PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set));
147: if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)ts, flg));
148: }
149: #endif
151: /* Monitor options */
152: PetscCall(PetscOptionsDeprecated("-ts_monitor_frequency", "-ts_dmswarm_monitor_moments_interval", "3.24", "Retired in favor of monitor-specific intervals (ts_dmswarm_monitor_moments was the only monitor to use ts_monitor_frequency)"));
153: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL));
154: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_wall_clock_time", "Monitor wall-clock time, KSP iterations, and SNES iterations per step", "TSMonitorWallClockTime", TSMonitorWallClockTime, TSMonitorWallClockTimeSetUp));
155: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL));
156: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, TSMonitorSolutionSetup));
157: PetscCall(TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL));
158: PetscCall(PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg));
159: if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)ts, monfilename));
161: PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt));
162: if (opt) {
163: PetscInt howoften = 1;
164: DM dm;
165: PetscBool net;
167: PetscCall(PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL));
168: PetscCall(TSGetDM(ts, &dm));
169: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net));
170: if (net) {
171: TSMonitorLGCtxNetwork ctx;
172: PetscCall(TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx));
173: PetscCall(TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxNetworkDestroy));
174: PetscCall(PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL));
175: } else {
176: TSMonitorLGCtx ctx;
177: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
178: PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
179: }
180: }
182: PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt));
183: if (opt) {
184: TSMonitorLGCtx ctx;
185: PetscInt howoften = 1;
187: PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL));
188: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
189: PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
190: }
191: PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL));
193: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt));
194: if (opt) {
195: TSMonitorLGCtx ctx;
196: PetscInt howoften = 1;
198: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
199: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
200: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
201: }
202: PetscCall(PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt));
203: if (opt) {
204: TSMonitorLGCtx ctx;
205: PetscInt howoften = 1;
207: PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
208: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
209: PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
210: ctx->semilogy = PETSC_TRUE;
211: }
213: PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt));
214: if (opt) {
215: TSMonitorLGCtx ctx;
216: PetscInt howoften = 1;
218: PetscCall(PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL));
219: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
220: PetscCall(TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
221: }
222: PetscCall(PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt));
223: if (opt) {
224: TSMonitorLGCtx ctx;
225: PetscInt howoften = 1;
227: PetscCall(PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL));
228: PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
229: PetscCall(TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscCtxDestroyFn *)TSMonitorLGCtxDestroy));
230: }
231: PetscCall(PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt));
232: if (opt) {
233: TSMonitorSPEigCtx ctx;
234: PetscInt howoften = 1;
236: PetscCall(PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL));
237: PetscCall(TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
238: PetscCall(TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscCtxDestroyFn *)TSMonitorSPEigCtxDestroy));
239: }
240: PetscCall(PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase space from the DMSwarm", "TSMonitorSPSwarm", &opt));
241: if (opt) {
242: TSMonitorSPCtx ctx;
243: PetscInt howoften = 1, retain = 0;
244: PetscBool phase = PETSC_TRUE, create = PETSC_TRUE, multispecies = PETSC_FALSE;
246: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
247: if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
248: create = PETSC_FALSE;
249: break;
250: }
251: if (create) {
252: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase space from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL));
253: PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL));
254: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL));
255: PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_multi_species", "Color particles by particle species", "TSMonitorSPSwarm", multispecies, &multispecies, NULL));
256: PetscCall(TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, multispecies, &ctx));
257: PetscCall(TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscCtxDestroyFn *)TSMonitorSPCtxDestroy));
258: }
259: }
260: PetscCall(PetscOptionsName("-ts_monitor_hg_swarm", "Display particle histogram from the DMSwarm", "TSMonitorHGSwarm", &opt));
261: if (opt) {
262: TSMonitorHGCtx ctx;
263: PetscInt howoften = 1, Ns = 1;
264: PetscBool velocity = PETSC_FALSE, create = PETSC_TRUE;
266: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
267: if (ts->monitor[i] == TSMonitorHGSwarmSolution) {
268: create = PETSC_FALSE;
269: break;
270: }
271: if (create) {
272: DM sw, dm;
273: PetscInt Nc, Nb;
275: PetscCall(TSGetDM(ts, &sw));
276: PetscCall(DMSwarmGetCellDM(sw, &dm));
277: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc));
278: Nb = PetscMin(20, PetscMax(10, Nc));
279: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL));
280: PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL));
281: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL));
282: PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL));
283: PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx));
284: PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscCtxDestroyFn *)TSMonitorHGCtxDestroy));
285: }
286: }
287: opt = PETSC_FALSE;
288: PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt));
289: if (opt) {
290: TSMonitorDrawCtx ctx;
291: PetscInt howoften = 1;
293: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL));
294: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
295: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
296: }
297: opt = PETSC_FALSE;
298: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt));
299: if (opt) {
300: TSMonitorDrawCtx ctx;
301: PetscReal bounds[4];
302: PetscInt n = 4;
303: PetscDraw draw;
304: PetscDrawAxis axis;
306: PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL));
307: PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field");
308: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx));
309: PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw));
310: PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis));
311: PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]));
312: PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2"));
313: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
314: }
315: opt = PETSC_FALSE;
316: PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt));
317: if (opt) {
318: TSMonitorDrawCtx ctx;
319: PetscInt howoften = 1;
321: PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL));
322: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
323: PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
324: }
325: opt = PETSC_FALSE;
326: PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt));
327: if (opt) {
328: TSMonitorDrawCtx ctx;
329: PetscInt howoften = 1;
331: PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL));
332: PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
333: PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscCtxDestroyFn *)TSMonitorDrawCtxDestroy));
334: }
336: opt = PETSC_FALSE;
337: 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));
338: if (flg) {
339: TSMonitorVTKCtx ctx;
341: PetscCall(TSMonitorSolutionVTKCtxCreate(monfilename, &ctx));
342: PetscCall(PetscOptionsInt("-ts_monitor_solution_vtk_interval", "Save every interval time step (-1 for last step only)", NULL, ctx->interval, &ctx->interval, NULL));
343: PetscCall(TSMonitorSet(ts, (PetscErrorCode (*)(TS, PetscInt, PetscReal, Vec, PetscCtx))TSMonitorSolutionVTK, ctx, (PetscCtxDestroyFn *)TSMonitorSolutionVTKDestroy));
344: }
346: PetscCall(PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg));
347: if (flg) {
348: TSMonitorDMDARayCtx *rayctx;
349: int ray = 0;
350: DMDirection ddir;
351: DM da;
352: PetscMPIInt rank;
354: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
355: if (dir[0] == 'x') ddir = DM_X;
356: else if (dir[0] == 'y') ddir = DM_Y;
357: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
358: sscanf(dir + 2, "%d", &ray);
360: PetscCall(PetscInfo(ts, "Displaying DMDA ray %c = %d\n", dir[0], ray));
361: PetscCall(PetscNew(&rayctx));
362: PetscCall(TSGetDM(ts, &da));
363: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
364: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank));
365: if (rank == 0) PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer));
366: rayctx->lgctx = NULL;
367: PetscCall(TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy));
368: }
369: PetscCall(PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg));
370: if (flg) {
371: TSMonitorDMDARayCtx *rayctx;
372: int ray = 0;
373: DMDirection ddir;
374: DM da;
375: PetscInt howoften = 1;
377: PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
378: if (dir[0] == 'x') ddir = DM_X;
379: else if (dir[0] == 'y') ddir = DM_Y;
380: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
381: sscanf(dir + 2, "%d", &ray);
383: PetscCall(PetscInfo(ts, "Displaying LG DMDA ray %c = %d\n", dir[0], ray));
384: PetscCall(PetscNew(&rayctx));
385: PetscCall(TSGetDM(ts, &da));
386: PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
387: PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx));
388: PetscCall(TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy));
389: }
391: PetscCall(PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt));
392: if (opt) {
393: TSMonitorEnvelopeCtx ctx;
395: PetscCall(TSMonitorEnvelopeCtxCreate(ts, &ctx));
396: PetscCall(TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscCtxDestroyFn *)TSMonitorEnvelopeCtxDestroy));
397: }
398: flg = PETSC_FALSE;
399: PetscCall(PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt));
400: if (opt && flg) PetscCall(TSMonitorCancel(ts));
402: flg = PETSC_FALSE;
403: PetscCall(PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeIJacobianDefaultColor", flg, &flg, NULL));
404: if (flg) {
405: DM dm;
407: PetscCall(TSGetDM(ts, &dm));
408: PetscCall(DMTSUnsetIJacobianContext_Internal(dm));
409: PetscCall(TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL));
410: PetscCall(PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n"));
411: }
413: /* Handle specific TS options */
414: PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);
416: /* Handle TSAdapt options */
417: PetscCall(TSGetAdapt(ts, &ts->adapt));
418: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
419: PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject));
421: /* TS trajectory must be set after TS, since it may use some TS options above */
422: tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
423: PetscCall(PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL));
424: if (tflg) PetscCall(TSSetSaveTrajectory(ts));
426: PetscCall(TSAdjointSetFromOptions(ts, PetscOptionsObject));
428: /* process any options handlers added with PetscObjectAddOptionsHandler() */
429: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject));
430: PetscOptionsEnd();
432: if (ts->trajectory) PetscCall(TSTrajectorySetFromOptions(ts->trajectory, ts));
434: /* why do we have to do this here and not during TSSetUp? */
435: PetscCall(TSGetSNES(ts, &ts->snes));
436: if (ts->problem_type == TS_LINEAR) {
437: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
438: if (!flg) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
439: }
440: PetscCall(SNESSetFromOptions(ts->snes));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: /*@
445: TSGetTrajectory - Gets the trajectory from a `TS` if it exists
447: Collective
449: Input Parameter:
450: . ts - the `TS` context obtained from `TSCreate()`
452: Output Parameter:
453: . tr - the `TSTrajectory` object, if it exists
455: Level: advanced
457: Note:
458: This routine should be called after all `TS` options have been set
460: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectoryCreate()`
461: @*/
462: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
463: {
464: PetscFunctionBegin;
466: *tr = ts->trajectory;
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /*@
471: TSSetSaveTrajectory - Causes the `TS` to save its solutions as it iterates forward in time in a `TSTrajectory` object
473: Collective
475: Input Parameter:
476: . ts - the `TS` context obtained from `TSCreate()`
478: Options Database Keys:
479: + -ts_save_trajectory - saves the trajectory to a file
480: - -ts_trajectory_type (basic|singlefile|memory|visualization) - set trajectory type
482: Level: intermediate
484: Notes:
485: This routine should be called after all `TS` options have been set
487: The `TSTRAJECTORYVISUALIZATION` files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
488: MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m
490: .seealso: [](ch_ts), `TS`, `TSTrajectoryType`, `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`
491: @*/
492: PetscErrorCode TSSetSaveTrajectory(TS ts)
493: {
494: PetscFunctionBegin;
496: if (!ts->trajectory) PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@
501: TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object
503: Collective
505: Input Parameter:
506: . ts - the `TS` context obtained from `TSCreate()`
508: Level: intermediate
510: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
511: @*/
512: PetscErrorCode TSResetTrajectory(TS ts)
513: {
514: PetscFunctionBegin;
516: if (ts->trajectory) {
517: PetscCall(TSTrajectoryDestroy(&ts->trajectory));
518: PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
519: }
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: /*@
524: TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from a `TS`
526: Collective
528: Input Parameter:
529: . ts - the `TS` context obtained from `TSCreate()`
531: Level: intermediate
533: .seealso: [](ch_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
534: @*/
535: PetscErrorCode TSRemoveTrajectory(TS ts)
536: {
537: PetscFunctionBegin;
539: PetscCall(TSTrajectoryDestroy(&ts->trajectory));
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: /*@
544: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
545: set with `TSSetRHSJacobian()`.
547: Collective
549: Input Parameters:
550: + ts - the `TS` context
551: . t - current timestep
552: - U - input vector
554: Output Parameters:
555: + A - Jacobian matrix
556: - B - optional matrix used to compute the preconditioner, often the same as `A`
558: Level: developer
560: Note:
561: Most users should not need to explicitly call this routine, as it
562: is used internally within the ODE integrators.
564: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()`
565: @*/
566: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
567: {
568: VecState Ustate;
569: PetscBool same;
570: DM dm;
571: DMTS tsdm;
572: TSRHSJacobianFn *rhsjacobianfunc;
573: void *ctx;
574: TSRHSFunctionFn *rhsfunction;
576: PetscFunctionBegin;
579: PetscCheckSameComm(ts, 1, U, 3);
580: PetscCall(TSGetDM(ts, &dm));
581: PetscCall(DMGetDMTS(dm, &tsdm));
582: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
583: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx));
584: PetscCall(VecGetState(U, &Ustate));
585: PetscCall(VecStateCompare(Ustate, ts->rhsjacobian.Xstate, &same));
587: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || same) && (rhsfunction != TSComputeRHSFunctionLinear)) PetscFunctionReturn(PETSC_SUCCESS);
589: 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);
590: if (rhsjacobianfunc) {
591: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
592: PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
593: ts->rhsjacs++;
594: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
595: } else {
596: PetscCall(MatZeroEntries(A));
597: if (B && A != B) PetscCall(MatZeroEntries(B));
598: }
599: ts->rhsjacobian.time = t;
600: ts->rhsjacobian.shift = 0;
601: ts->rhsjacobian.scale = 1.;
602: PetscCall(VecGetState(U, &ts->rhsjacobian.Xstate));
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: /*@
607: TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS`
609: Collective
611: Input Parameters:
612: + ts - the `TS` context
613: . t - current time
614: - U - state vector
616: Output Parameter:
617: . y - right-hand side
619: Level: developer
621: Note:
622: Most users should not need to explicitly call this routine, as it
623: is used internally within the nonlinear solvers.
625: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
626: @*/
627: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
628: {
629: TSRHSFunctionFn *rhsfunction;
630: TSIFunctionFn *ifunction;
631: void *ctx;
632: DM dm;
634: PetscFunctionBegin;
638: PetscCall(TSGetDM(ts, &dm));
639: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx));
640: PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
642: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
644: if (rhsfunction) {
645: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, y, 0));
646: PetscCall(VecLockReadPush(U));
647: PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
648: PetscCall(VecLockReadPop(U));
649: ts->rhsfuncs++;
650: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, y, 0));
651: } else PetscCall(VecZeroEntries(y));
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: /*@
656: TSComputeSolutionFunction - Evaluates the solution function.
658: Collective
660: Input Parameters:
661: + ts - the `TS` context
662: - t - current time
664: Output Parameter:
665: . U - the solution
667: Level: developer
669: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
670: @*/
671: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
672: {
673: TSSolutionFn *solutionfunction;
674: void *ctx;
675: DM dm;
677: PetscFunctionBegin;
680: PetscCall(TSGetDM(ts, &dm));
681: PetscCall(DMTSGetSolutionFunction(dm, &solutionfunction, &ctx));
682: if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
685: /*@
686: TSComputeForcingFunction - Evaluates the forcing function.
688: Collective
690: Input Parameters:
691: + ts - the `TS` context
692: - t - current time
694: Output Parameter:
695: . U - the function value
697: Level: developer
699: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
700: @*/
701: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
702: {
703: void *ctx;
704: DM dm;
705: TSForcingFn *forcing;
707: PetscFunctionBegin;
710: PetscCall(TSGetDM(ts, &dm));
711: PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));
713: if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }
717: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
718: {
719: Mat A, B;
720: TSIJacobianFn *ijacobian;
722: PetscFunctionBegin;
723: if (Arhs) *Arhs = NULL;
724: if (Brhs) *Brhs = NULL;
725: PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL));
726: if (Arhs) {
727: if (!ts->Arhs) {
728: if (ijacobian) {
729: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs));
730: PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN));
731: } else {
732: ts->Arhs = A;
733: PetscCall(PetscObjectReference((PetscObject)A));
734: }
735: } else {
736: PetscBool flg;
737: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
738: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
739: if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
740: PetscCall(PetscObjectDereference((PetscObject)ts->Arhs));
741: ts->Arhs = A;
742: PetscCall(PetscObjectReference((PetscObject)A));
743: }
744: }
745: *Arhs = ts->Arhs;
746: }
747: if (Brhs) {
748: if (!ts->Brhs) {
749: if (A != B) {
750: if (ijacobian) {
751: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs));
752: } else {
753: ts->Brhs = B;
754: PetscCall(PetscObjectReference((PetscObject)B));
755: }
756: } else {
757: PetscCall(PetscObjectReference((PetscObject)ts->Arhs));
758: ts->Brhs = ts->Arhs;
759: }
760: }
761: *Brhs = ts->Brhs;
762: }
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
766: /*@
767: TSComputeIFunction - Evaluates the DAE residual written in the implicit form F(t,U,Udot)=0
769: Collective
771: Input Parameters:
772: + ts - the `TS` context
773: . t - current time
774: . U - state vector
775: . Udot - time derivative of state vector
776: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSFunction should be kept separate
778: Output Parameter:
779: . Y - right-hand side
781: Level: developer
783: Note:
784: Most users should not need to explicitly call this routine, as it
785: is used internally within the nonlinear solvers.
787: If the user did not write their equations in implicit form, this
788: function recasts them in implicit form.
790: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
791: @*/
792: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
793: {
794: TSIFunctionFn *ifunction;
795: TSRHSFunctionFn *rhsfunction;
796: void *ctx;
797: DM dm;
799: PetscFunctionBegin;
805: PetscCall(TSGetDM(ts, &dm));
806: PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
807: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
809: PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");
811: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, Udot, Y));
812: if (ifunction) {
813: PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
814: ts->ifuncs++;
815: }
816: if (imex) {
817: if (!ifunction) PetscCall(VecCopy(Udot, Y));
818: } else if (rhsfunction) {
819: if (ifunction) {
820: Vec Frhs;
822: PetscCall(DMGetGlobalVector(dm, &Frhs));
823: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
824: PetscCall(VecAXPY(Y, -1, Frhs));
825: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
826: } else {
827: PetscCall(TSComputeRHSFunction(ts, t, U, Y));
828: PetscCall(VecAYPX(Y, -1, Udot));
829: }
830: }
831: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, Udot, Y));
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: /*
836: TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call `TSComputeRHSJacobian()` on it.
838: Note:
839: This routine is needed when one switches from `TSComputeIJacobian()` to `TSComputeRHSJacobian()` because the Jacobian matrix may be shifted or scaled in `TSComputeIJacobian()`.
841: */
842: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
843: {
844: PetscFunctionBegin;
846: PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
847: PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");
849: if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
850: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
851: if (B && B == ts->Brhs && A != B) {
852: if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
853: if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
854: }
855: ts->rhsjacobian.shift = 0;
856: ts->rhsjacobian.scale = 1.;
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: /*
861: TSComputeIJacobian_Internal - Evaluates the Jacobian of the DAE
863: Collective
865: Input Parameters:
866: + ts - the `TS` context
867: . ijacobian - function to compute LHS Jacobian
868: . rhsjacobian - function to compute RHS Jacobian
869: . ctx - user context for Jacobian functions
870: . t - current timestep
871: . U - state vector
872: . Udot - time derivative of state vector
873: . shift - shift to apply, see note below
874: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSJacobian should be kept separate
876: Output Parameters:
877: + A - Jacobian matrix
878: - B - matrix from which the preconditioner is constructed; often the same as `A`
880: Level: developer
882: Notes:
883: This function exists so that a user can assemble the Jacobian pieces with functions not stores in the `DMTS`. This was necessary in `TSDISCGRAD` since two different representations of the formulation can be stored.
885: If $F(t,U,\dot{U})=0$ is the DAE, the required Jacobian is
886: .vb
887: dF/dU + shift*dF/dUdot
888: .ve
890: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
891: */
892: PetscErrorCode TSComputeIJacobian_Internal(TS ts, TSIJacobianFn *ijacobian, TSRHSJacobianFn *rhsjacobian, void *ctx, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
893: {
894: PetscFunctionBegin;
895: PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
897: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
898: if (ijacobian) {
899: PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
900: ts->ijacs++;
901: }
902: if (imex) {
903: if (!ijacobian) { /* system was written as Udot = G(t,U) */
904: PetscBool assembled;
905: if (rhsjacobian) {
906: Mat Arhs = NULL;
907: PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
908: if (A == Arhs) {
909: 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 */
910: ts->rhsjacobian.time = PETSC_MIN_REAL;
911: }
912: }
913: PetscCall(MatZeroEntries(A));
914: PetscCall(MatAssembled(A, &assembled));
915: if (!assembled) {
916: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
917: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
918: }
919: PetscCall(MatShift(A, shift));
920: if (A != B) {
921: PetscCall(MatZeroEntries(B));
922: PetscCall(MatAssembled(B, &assembled));
923: if (!assembled) {
924: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
925: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
926: }
927: PetscCall(MatShift(B, shift));
928: }
929: }
930: } else {
931: Mat Arhs = NULL, Brhs = NULL;
933: /* RHSJacobian needs to be converted to part of IJacobian if exists */
934: if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
935: if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
936: DM dm;
937: VecState Ustate;
938: PetscBool same;
939: TSRHSFunctionFn *rhsfunction;
941: PetscCall(TSGetDM(ts, &dm));
942: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
943: PetscCall(VecGetState(U, &Ustate));
944: PetscCall(VecStateCompare(Ustate, ts->rhsjacobian.Xstate, &same));
945: if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || same) && rhsfunction != TSComputeRHSFunctionLinear)) && ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */
946: PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
947: if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
948: } else {
949: PetscBool flg;
951: if (ts->rhsjacobian.reuse) { /* Undo the damage */
952: /* MatScale has a short path for this case.
953: However, this code path is taken the first time TSComputeRHSJacobian is called
954: and the matrices have not been assembled yet */
955: PetscCall(TSRecoverRHSJacobian(ts, A, B));
956: }
957: PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
958: PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
959: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
960: if (!flg) {
961: PetscCall(MatScale(A, -1));
962: PetscCall(MatShift(A, shift));
963: }
964: if (A != B) {
965: PetscCall(MatScale(B, -1));
966: PetscCall(MatShift(B, shift));
967: }
968: }
969: ts->rhsjacobian.scale = -1;
970: ts->rhsjacobian.shift = shift;
971: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
972: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
973: PetscCall(MatZeroEntries(A));
974: PetscCall(MatShift(A, shift));
975: if (A != B) {
976: PetscCall(MatZeroEntries(B));
977: PetscCall(MatShift(B, shift));
978: }
979: }
980: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
981: PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
982: if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
983: }
984: }
985: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
986: PetscFunctionReturn(PETSC_SUCCESS);
987: }
989: /*@
990: TSComputeIJacobian - Evaluates the Jacobian of the DAE
992: Collective
994: Input Parameters:
995: + ts - the `TS` context
996: . t - current timestep
997: . U - state vector
998: . Udot - time derivative of state vector
999: . shift - shift to apply, see note below
1000: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSJacobian should be kept separate
1002: Output Parameters:
1003: + A - Jacobian matrix
1004: - B - matrix from which the preconditioner is constructed; often the same as `A`
1006: Level: developer
1008: Notes:
1009: If $F(t,U,\dot{U})=0$ is the DAE, the required Jacobian is
1010: .vb
1011: dF/dU + shift*dF/dUdot
1012: .ve
1013: Most users should not need to explicitly call this routine, as it
1014: is used internally within the nonlinear solvers.
1016: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
1017: @*/
1018: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
1019: {
1020: TSIJacobianFn *ijacobian;
1021: TSRHSJacobianFn *rhsjacobian;
1022: DM dm;
1023: void *ctx;
1025: PetscFunctionBegin;
1032: PetscCall(TSGetDM(ts, &dm));
1033: PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
1034: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1035: PetscCall(TSComputeIJacobian_Internal(ts, ijacobian, rhsjacobian, ctx, t, U, Udot, shift, A, B, imex));
1036: PetscFunctionReturn(PETSC_SUCCESS);
1037: }
1039: /*@
1040: TSSetRHSFunction - Sets the routine for evaluating the function,
1041: where U_t = G(t,u).
1043: Logically Collective
1045: Input Parameters:
1046: + ts - the `TS` context obtained from `TSCreate()`
1047: . r - vector to put the computed right-hand side (or `NULL` to have it created)
1048: . f - routine for evaluating the right-hand-side function
1049: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)
1051: Level: beginner
1053: Note:
1054: You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.
1056: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
1057: @*/
1058: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, PetscCtx ctx)
1059: {
1060: SNES snes;
1061: Vec ralloc = NULL;
1062: DM dm;
1064: PetscFunctionBegin;
1068: PetscCall(TSGetDM(ts, &dm));
1069: PetscCall(DMTSSetRHSFunction(dm, f, ctx));
1070: PetscCall(TSGetSNES(ts, &snes));
1071: if (!r && !ts->dm && ts->vec_sol) {
1072: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1073: r = ralloc;
1074: }
1075: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1076: PetscCall(VecDestroy(&ralloc));
1077: PetscFunctionReturn(PETSC_SUCCESS);
1078: }
1080: /*@
1081: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1083: Logically Collective
1085: Input Parameters:
1086: + ts - the `TS` context obtained from `TSCreate()`
1087: . f - routine for evaluating the solution
1088: - ctx - [optional] user-defined context for private data for the
1089: function evaluation routine (may be `NULL`)
1091: Options Database Keys:
1092: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1093: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
1095: Level: intermediate
1097: Notes:
1098: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1099: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1100: create closed-form solutions with non-physical forcing terms.
1102: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1104: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1105: @*/
1106: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, PetscCtx ctx)
1107: {
1108: DM dm;
1110: PetscFunctionBegin;
1112: PetscCall(TSGetDM(ts, &dm));
1113: PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
1114: PetscFunctionReturn(PETSC_SUCCESS);
1115: }
1117: /*@
1118: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1120: Logically Collective
1122: Input Parameters:
1123: + ts - the `TS` context obtained from `TSCreate()`
1124: . func - routine for evaluating the forcing function
1125: - ctx - [optional] user-defined context for private data for the function evaluation routine
1126: (may be `NULL`)
1128: Level: intermediate
1130: Notes:
1131: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1132: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1133: definition of the problem you are solving and hence possibly introducing bugs.
1135: This replaces the ODE F(u,u_t,t) = 0 the `TS` is solving with F(u,u_t,t) - func(t) = 0
1137: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1138: parameters can be passed in the ctx variable.
1140: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1142: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
1143: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1144: @*/
1145: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, PetscCtx ctx)
1146: {
1147: DM dm;
1149: PetscFunctionBegin;
1151: PetscCall(TSGetDM(ts, &dm));
1152: PetscCall(DMTSSetForcingFunction(dm, func, ctx));
1153: PetscFunctionReturn(PETSC_SUCCESS);
1154: }
1156: /*@
1157: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1158: where U_t = G(U,t), as well as the location to store the matrix.
1160: Logically Collective
1162: Input Parameters:
1163: + ts - the `TS` context obtained from `TSCreate()`
1164: . Amat - (approximate) location to store Jacobian matrix entries computed by `f`
1165: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
1166: . f - the Jacobian evaluation routine
1167: - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1169: Level: beginner
1171: Notes:
1172: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1174: The `TS` solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f()`
1175: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1177: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
1178: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
1179: @*/
1180: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, PetscCtx ctx)
1181: {
1182: SNES snes;
1183: DM dm;
1184: TSIJacobianFn *ijacobian;
1186: PetscFunctionBegin;
1190: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1191: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1193: PetscCall(TSGetDM(ts, &dm));
1194: PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
1195: PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
1196: PetscCall(TSGetSNES(ts, &snes));
1197: if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1198: if (Amat) {
1199: PetscCall(PetscObjectReference((PetscObject)Amat));
1200: PetscCall(MatDestroy(&ts->Arhs));
1201: ts->Arhs = Amat;
1202: }
1203: if (Pmat) {
1204: PetscCall(PetscObjectReference((PetscObject)Pmat));
1205: PetscCall(MatDestroy(&ts->Brhs));
1206: ts->Brhs = Pmat;
1207: }
1208: PetscFunctionReturn(PETSC_SUCCESS);
1209: }
1211: /*@
1212: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1214: Logically Collective
1216: Input Parameters:
1217: + ts - the `TS` context obtained from `TSCreate()`
1218: . r - vector to hold the residual (or `NULL` to have it created internally)
1219: . f - the function evaluation routine
1220: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1222: Level: beginner
1224: Note:
1225: The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function.
1227: .seealso: [](ch_ts), `TS`, `TSIFunctionFn`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`,
1228: `TSSetIJacobian()`
1229: @*/
1230: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunctionFn *f, PetscCtx ctx)
1231: {
1232: SNES snes;
1233: Vec ralloc = NULL;
1234: DM dm;
1236: PetscFunctionBegin;
1240: PetscCall(TSGetDM(ts, &dm));
1241: PetscCall(DMTSSetIFunction(dm, f, ctx));
1243: PetscCall(TSGetSNES(ts, &snes));
1244: if (!r && !ts->dm && ts->vec_sol) {
1245: PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
1246: r = ralloc;
1247: }
1248: PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
1249: PetscCall(VecDestroy(&ralloc));
1250: PetscFunctionReturn(PETSC_SUCCESS);
1251: }
1253: /*@
1254: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.
1256: Not Collective
1258: Input Parameter:
1259: . ts - the `TS` context
1261: Output Parameters:
1262: + r - vector to hold residual (or `NULL`)
1263: . func - the function to compute residual (or `NULL`)
1264: - ctx - the function context (or `NULL`)
1266: Level: advanced
1268: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1269: @*/
1270: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunctionFn **func, PetscCtxRt ctx)
1271: {
1272: SNES snes;
1273: DM dm;
1275: PetscFunctionBegin;
1277: PetscCall(TSGetSNES(ts, &snes));
1278: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1279: PetscCall(TSGetDM(ts, &dm));
1280: PetscCall(DMTSGetIFunction(dm, func, ctx));
1281: PetscFunctionReturn(PETSC_SUCCESS);
1282: }
1284: /*@
1285: TSGetRHSFunction - Returns the vector where the right-hand side is stored and the function/context to compute it.
1287: Not Collective
1289: Input Parameter:
1290: . ts - the `TS` context
1292: Output Parameters:
1293: + r - vector to hold computed right-hand side (or `NULL`)
1294: . func - the function to compute right-hand side (or `NULL`)
1295: - ctx - the function context (or `NULL`)
1297: Level: advanced
1299: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1300: @*/
1301: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunctionFn **func, PetscCtxRt ctx)
1302: {
1303: SNES snes;
1304: DM dm;
1306: PetscFunctionBegin;
1308: PetscCall(TSGetSNES(ts, &snes));
1309: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1310: PetscCall(TSGetDM(ts, &dm));
1311: PetscCall(DMTSGetRHSFunction(dm, func, ctx));
1312: PetscFunctionReturn(PETSC_SUCCESS);
1313: }
1315: /*@
1316: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1317: provided with `TSSetIFunction()`.
1319: Logically Collective
1321: Input Parameters:
1322: + ts - the `TS` context obtained from `TSCreate()`
1323: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
1324: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
1325: . f - the Jacobian evaluation routine
1326: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1328: Level: beginner
1330: Notes:
1331: The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1333: If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
1334: space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.
1336: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1337: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1338: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1339: a and vector W depend on the integration method, step size, and past states. For example with
1340: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1341: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1343: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1345: The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f`
1346: You should not assume the values are the same in the next call to `f` as you set them in the previous call.
1348: In case `TSSetRHSJacobian()` is also used in conjunction with a fully-implicit solver,
1349: multilevel linear solvers, e.g. `PCMG`, will likely not work due to the way `TS` handles rhs matrices.
1351: .seealso: [](ch_ts), `TS`, `TSIJacobianFn`, `TSSetIFunction()`, `TSSetRHSJacobian()`,
1352: `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1353: @*/
1354: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobianFn *f, PetscCtx ctx)
1355: {
1356: SNES snes;
1357: DM dm;
1359: PetscFunctionBegin;
1363: if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
1364: if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);
1366: PetscCall(TSGetDM(ts, &dm));
1367: PetscCall(DMTSSetIJacobian(dm, f, ctx));
1369: PetscCall(TSGetSNES(ts, &snes));
1370: PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
1371: PetscFunctionReturn(PETSC_SUCCESS);
1372: }
1374: /*@
1375: TSRHSJacobianSetReuse - restore the RHS Jacobian before calling the user-provided `TSRHSJacobianFn` function again
1377: Logically Collective
1379: Input Parameters:
1380: + ts - `TS` context obtained from `TSCreate()`
1381: - reuse - `PETSC_TRUE` if the RHS Jacobian
1383: Level: intermediate
1385: Notes:
1386: Without this flag, `TS` will change the sign and shift the RHS Jacobian for a
1387: finite-time-step implicit solve, in which case the user function will need to recompute the
1388: entire Jacobian. The `reuse `flag must be set if the evaluation function assumes that the
1389: matrix entries have not been changed by the `TS`.
1391: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1392: @*/
1393: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1394: {
1395: PetscFunctionBegin;
1396: ts->rhsjacobian.reuse = reuse;
1397: PetscFunctionReturn(PETSC_SUCCESS);
1398: }
1400: /*@
1401: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1403: Logically Collective
1405: Input Parameters:
1406: + ts - the `TS` context obtained from `TSCreate()`
1407: . F - vector to hold the residual (or `NULL` to have it created internally)
1408: . fun - the function evaluation routine
1409: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)
1411: Level: beginner
1413: .seealso: [](ch_ts), `TS`, `TSI2FunctionFn`, `TSSetI2Jacobian()`, `TSSetIFunction()`,
1414: `TSCreate()`, `TSSetRHSFunction()`
1415: @*/
1416: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2FunctionFn *fun, PetscCtx ctx)
1417: {
1418: DM dm;
1420: PetscFunctionBegin;
1423: PetscCall(TSSetIFunction(ts, F, NULL, NULL));
1424: PetscCall(TSGetDM(ts, &dm));
1425: PetscCall(DMTSSetI2Function(dm, fun, ctx));
1426: PetscFunctionReturn(PETSC_SUCCESS);
1427: }
1429: /*@
1430: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.
1432: Not Collective
1434: Input Parameter:
1435: . ts - the `TS` context
1437: Output Parameters:
1438: + r - vector to hold residual (or `NULL`)
1439: . fun - the function to compute residual (or `NULL`)
1440: - ctx - the function context (or `NULL`)
1442: Level: advanced
1444: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1445: @*/
1446: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2FunctionFn **fun, PetscCtxRt ctx)
1447: {
1448: SNES snes;
1449: DM dm;
1451: PetscFunctionBegin;
1453: PetscCall(TSGetSNES(ts, &snes));
1454: PetscCall(SNESGetFunction(snes, r, NULL, NULL));
1455: PetscCall(TSGetDM(ts, &dm));
1456: PetscCall(DMTSGetI2Function(dm, fun, ctx));
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: /*@
1461: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1462: where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.
1464: Logically Collective
1466: Input Parameters:
1467: + ts - the `TS` context obtained from `TSCreate()`
1468: . J - matrix to hold the Jacobian values
1469: . P - matrix for constructing the preconditioner (may be same as `J`)
1470: . jac - the Jacobian evaluation routine, see `TSI2JacobianFn` for the calling sequence
1471: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)
1473: Level: beginner
1475: Notes:
1476: The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve.
1478: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1479: 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.
1480: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1481: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1483: .seealso: [](ch_ts), `TS`, `TSI2JacobianFn`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1484: @*/
1485: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2JacobianFn *jac, PetscCtx ctx)
1486: {
1487: DM dm;
1489: PetscFunctionBegin;
1493: PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
1494: PetscCall(TSGetDM(ts, &dm));
1495: PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
1496: PetscFunctionReturn(PETSC_SUCCESS);
1497: }
1499: /*@
1500: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1502: Not Collective, but parallel objects are returned if `TS` is parallel
1504: Input Parameter:
1505: . ts - The `TS` context obtained from `TSCreate()`
1507: Output Parameters:
1508: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1509: . P - The matrix from which the preconditioner is constructed, often the same as `J`
1510: . jac - The function to compute the Jacobian matrices
1511: - ctx - User-defined context for Jacobian evaluation routine
1513: Level: advanced
1515: Note:
1516: You can pass in `NULL` for any return argument you do not need.
1518: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1519: @*/
1520: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2JacobianFn **jac, PetscCtxRt ctx)
1521: {
1522: SNES snes;
1523: DM dm;
1525: PetscFunctionBegin;
1526: PetscCall(TSGetSNES(ts, &snes));
1527: PetscCall(SNESSetUpMatrices(snes));
1528: PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
1529: PetscCall(TSGetDM(ts, &dm));
1530: PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
1531: PetscFunctionReturn(PETSC_SUCCESS);
1532: }
1534: /*@
1535: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1537: Collective
1539: Input Parameters:
1540: + ts - the `TS` context
1541: . t - current time
1542: . U - state vector
1543: . V - time derivative of state vector (U_t)
1544: - A - second time derivative of state vector (U_tt)
1546: Output Parameter:
1547: . F - the residual vector
1549: Level: developer
1551: Note:
1552: Most users should not need to explicitly call this routine, as it
1553: is used internally within the nonlinear solvers.
1555: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1556: @*/
1557: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1558: {
1559: DM dm;
1560: TSI2FunctionFn *I2Function;
1561: void *ctx;
1562: TSRHSFunctionFn *rhsfunction;
1564: PetscFunctionBegin;
1571: PetscCall(TSGetDM(ts, &dm));
1572: PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
1573: PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
1575: if (!I2Function) {
1576: PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
1577: PetscFunctionReturn(PETSC_SUCCESS);
1578: }
1580: PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, V, F));
1582: PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));
1584: if (rhsfunction) {
1585: Vec Frhs;
1587: PetscCall(DMGetGlobalVector(dm, &Frhs));
1588: PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
1589: PetscCall(VecAXPY(F, -1, Frhs));
1590: PetscCall(DMRestoreGlobalVector(dm, &Frhs));
1591: }
1593: PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, V, F));
1594: PetscFunctionReturn(PETSC_SUCCESS);
1595: }
1597: /*@
1598: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1600: Collective
1602: Input Parameters:
1603: + ts - the `TS` context
1604: . t - current timestep
1605: . U - state vector
1606: . V - time derivative of state vector
1607: . A - second time derivative of state vector
1608: . shiftV - shift to apply, see note below
1609: - shiftA - shift to apply, see note below
1611: Output Parameters:
1612: + J - Jacobian matrix
1613: - P - optional matrix used to construct the preconditioner
1615: Level: developer
1617: Notes:
1618: If $F(t,U,V,A) = 0$ is the DAE, the required Jacobian is
1620: $$
1621: dF/dU + shiftV dF/dV + shiftA dF/dA
1622: $$
1624: Most users should not need to explicitly call this routine, as it
1625: is used internally within the ODE integrators.
1627: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
1628: @*/
1629: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1630: {
1631: DM dm;
1632: TSI2JacobianFn *I2Jacobian;
1633: void *ctx;
1634: TSRHSJacobianFn *rhsjacobian;
1636: PetscFunctionBegin;
1644: PetscCall(TSGetDM(ts, &dm));
1645: PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
1646: PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));
1648: if (!I2Jacobian) {
1649: PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
1650: PetscFunctionReturn(PETSC_SUCCESS);
1651: }
1653: PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, J, P));
1654: PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1655: if (rhsjacobian) {
1656: Mat Jrhs, Prhs;
1657: PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
1658: PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
1659: PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
1660: if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
1661: }
1663: PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, J, P));
1664: PetscFunctionReturn(PETSC_SUCCESS);
1665: }
1667: /*@
1668: TSSetTransientVariable - sets function to transform from state to transient variables
1670: Logically Collective
1672: Input Parameters:
1673: + ts - time stepping context on which to change the transient variable
1674: . tvar - a function that transforms to transient variables, see `TSTransientVariableFn` for the calling sequence
1675: - ctx - a context for tvar
1677: Level: advanced
1679: Notes:
1680: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1681: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1682: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1683: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1684: evaluated via the chain rule, as in
1685: .vb
1686: dF/dP + shift * dF/dCdot dC/dP.
1687: .ve
1689: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSTransientVariableFn`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1690: @*/
1691: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariableFn *tvar, PetscCtx ctx)
1692: {
1693: DM dm;
1695: PetscFunctionBegin;
1697: PetscCall(TSGetDM(ts, &dm));
1698: PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
1699: PetscFunctionReturn(PETSC_SUCCESS);
1700: }
1702: /*@
1703: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1705: Logically Collective
1707: Input Parameters:
1708: + ts - TS on which to compute
1709: - U - state vector to be transformed to transient variables
1711: Output Parameter:
1712: . C - transient (conservative) variable
1714: Level: developer
1716: Developer Notes:
1717: If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed.
1718: This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are
1719: being used.
1721: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1722: @*/
1723: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1724: {
1725: DM dm;
1726: DMTS dmts;
1728: PetscFunctionBegin;
1731: PetscCall(TSGetDM(ts, &dm));
1732: PetscCall(DMGetDMTS(dm, &dmts));
1733: if (dmts->ops->transientvar) {
1735: PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
1736: }
1737: PetscFunctionReturn(PETSC_SUCCESS);
1738: }
1740: /*@
1741: TSHasTransientVariable - determine whether transient variables have been set
1743: Logically Collective
1745: Input Parameter:
1746: . ts - `TS` on which to compute
1748: Output Parameter:
1749: . has - `PETSC_TRUE` if transient variables have been set
1751: Level: developer
1753: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1754: @*/
1755: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1756: {
1757: DM dm;
1758: DMTS dmts;
1760: PetscFunctionBegin;
1762: PetscCall(TSGetDM(ts, &dm));
1763: PetscCall(DMGetDMTS(dm, &dmts));
1764: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1765: PetscFunctionReturn(PETSC_SUCCESS);
1766: }
1768: /*@
1769: TS2SetSolution - Sets the initial solution and time derivative vectors
1770: for use by the `TS` routines handling second order equations.
1772: Logically Collective
1774: Input Parameters:
1775: + ts - the `TS` context obtained from `TSCreate()`
1776: . u - the solution vector
1777: - v - the time derivative vector
1779: Level: beginner
1781: .seealso: [](ch_ts), `TS`
1782: @*/
1783: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1784: {
1785: PetscFunctionBegin;
1789: PetscCall(TSSetSolution(ts, u));
1790: PetscCall(PetscObjectReference((PetscObject)v));
1791: PetscCall(VecDestroy(&ts->vec_dot));
1792: ts->vec_dot = v;
1793: PetscFunctionReturn(PETSC_SUCCESS);
1794: }
1796: /*@
1797: TS2GetSolution - Returns the solution and time derivative at the present timestep
1798: for second order equations.
1800: Not Collective
1802: Input Parameter:
1803: . ts - the `TS` context obtained from `TSCreate()`
1805: Output Parameters:
1806: + u - the vector containing the solution
1807: - v - the vector containing the time derivative
1809: Level: intermediate
1811: Notes:
1812: It is valid to call this routine inside the function
1813: that you are evaluating in order to move to the new timestep. This vector not
1814: changed until the solution at the next timestep has been calculated.
1816: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1817: @*/
1818: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1819: {
1820: PetscFunctionBegin;
1822: if (u) PetscAssertPointer(u, 2);
1823: if (v) PetscAssertPointer(v, 3);
1824: if (u) *u = ts->vec_sol;
1825: if (v) *v = ts->vec_dot;
1826: PetscFunctionReturn(PETSC_SUCCESS);
1827: }
1829: /*@
1830: TSLoad - Loads a `TS` that has been stored in binary with `TSView()`.
1832: Collective
1834: Input Parameters:
1835: + ts - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1836: some related function before a call to `TSLoad()`.
1837: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1839: Level: intermediate
1841: Note:
1842: The type is determined by the data in the file, any type set into the `TS` before this call is ignored.
1844: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1845: @*/
1846: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1847: {
1848: PetscBool isbinary;
1849: PetscInt classid;
1850: char type[256];
1851: DMTS sdm;
1852: DM dm;
1854: PetscFunctionBegin;
1857: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1858: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1860: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
1861: PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
1862: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
1863: PetscCall(TSSetType(ts, type));
1864: PetscTryTypeMethod(ts, load, viewer);
1865: PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
1866: PetscCall(DMLoad(dm, viewer));
1867: PetscCall(TSSetDM(ts, dm));
1868: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
1869: PetscCall(VecLoad(ts->vec_sol, viewer));
1870: PetscCall(DMGetDMTS(ts->dm, &sdm));
1871: PetscCall(DMTSLoad(sdm, viewer));
1872: PetscFunctionReturn(PETSC_SUCCESS);
1873: }
1875: #include <petscdraw.h>
1876: #if PetscDefined(HAVE_SAWS)
1877: #include <petscviewersaws.h>
1878: #endif
1880: /*@
1881: TSViewFromOptions - View a `TS` based on values in the options database
1883: Collective
1885: Input Parameters:
1886: + ts - the `TS` context
1887: . obj - optional object that provides the prefix for the options database keys, pass `NULL` to use the options prefix of `ts`
1888: - name - command line option string to be passed by user
1890: Options Database Key:
1891: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
1893: Level: intermediate
1895: Note:
1896: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
1897: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
1899: .seealso: [](ch_ts), `TS`, `TSView()`, `PetscObjectViewFromOptions()`, `TSCreate()`, `PetscOptionsCreateViewer()`
1900: @*/
1901: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1902: {
1903: PetscFunctionBegin;
1905: PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
1906: PetscFunctionReturn(PETSC_SUCCESS);
1907: }
1909: /*@
1910: TSView - Displays the `TS` data structure.
1912: Collective
1914: Input Parameters:
1915: + ts - the `TS` context obtained from `TSCreate()`
1916: - viewer - visualization context
1918: Options Database Key:
1919: . -ts_view viewer_specification - calls `TSView()` at end of `TSStep()`. See `PetscOptionsCreateViewer()` for the format of `viewer_specification`
1921: Level: beginner
1923: Notes:
1924: The available visualization contexts include
1925: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1926: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1927: output where only the first processor opens
1928: the file. All other processors send their
1929: data to the first processor to print.
1931: The user can open an alternative visualization context with
1932: `PetscViewerASCIIOpen()` - output to a specified file.
1934: In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).
1936: The "initial time step" displayed is the default time step from `TSCreate()` or that set with `TSSetTimeStep()` or `-ts_time_step`
1938: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`, `PetscOptionsCreateViewer()`, `TSViewFromOptions()`
1939: @*/
1940: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1941: {
1942: TSType type;
1943: PetscBool isascii, isstring, issundials, isbinary, isdraw;
1944: DMTS sdm;
1945: #if PetscDefined(HAVE_SAWS)
1946: PetscBool issaws;
1947: #endif
1949: PetscFunctionBegin;
1951: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
1953: PetscCheckSameComm(ts, 1, viewer, 2);
1955: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1956: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
1957: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1958: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1959: #if PetscDefined(HAVE_SAWS)
1960: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
1961: #endif
1962: if (isascii) {
1963: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer));
1964: if (ts->ops->view) {
1965: PetscCall(PetscViewerASCIIPushTab(viewer));
1966: PetscUseTypeMethod(ts, view, viewer);
1967: PetscCall(PetscViewerASCIIPopTab(viewer));
1968: }
1969: PetscCall(PetscViewerASCIIPrintf(viewer, " initial time step=%g\n", (double)ts->initial_time_step));
1970: if (ts->max_steps < PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps));
1971: if (ts->run_steps < PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, " run steps=%" PetscInt_FMT "\n", ts->run_steps));
1972: if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time));
1973: if (ts->max_reject != PETSC_UNLIMITED) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum number of step rejections=%" PetscInt_FMT "\n", ts->max_reject));
1974: if (ts->max_snes_failures != PETSC_UNLIMITED) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum number of SNES failures allowed=%" PetscInt_FMT "\n", ts->max_snes_failures));
1975: if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs));
1976: if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs));
1977: if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs));
1978: if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs));
1979: if (ts->usessnes) {
1980: PetscBool lin;
1981: if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its));
1982: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its));
1983: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
1984: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures));
1985: }
1986: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject));
1987: if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, "));
1988: else PetscCall(PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol));
1989: if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, "using vector of absolute error tolerances\n"));
1990: else PetscCall(PetscViewerASCIIPrintf(viewer, "using absolute error tolerance of %g\n", (double)ts->atol));
1991: PetscCall(PetscViewerASCIIPushTab(viewer));
1992: PetscCall(TSAdaptView(ts->adapt, viewer));
1993: PetscCall(PetscViewerASCIIPopTab(viewer));
1994: } else if (isstring) {
1995: PetscCall(TSGetType(ts, &type));
1996: PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type));
1997: PetscTryTypeMethod(ts, view, viewer);
1998: } else if (isbinary) {
1999: PetscInt classid = TS_FILE_CLASSID;
2000: MPI_Comm comm;
2001: PetscMPIInt rank;
2002: char type[256];
2004: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
2005: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2006: if (rank == 0) {
2007: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
2008: PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
2009: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
2010: }
2011: PetscTryTypeMethod(ts, view, viewer);
2012: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
2013: PetscCall(DMView(ts->dm, viewer));
2014: PetscCall(VecView(ts->vec_sol, viewer));
2015: PetscCall(DMGetDMTS(ts->dm, &sdm));
2016: PetscCall(DMTSView(sdm, viewer));
2017: } else if (isdraw) {
2018: PetscDraw draw;
2019: char str[36];
2020: PetscReal x, y, bottom, h;
2022: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2023: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
2024: PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
2025: PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
2026: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
2027: bottom = y - h;
2028: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
2029: PetscTryTypeMethod(ts, view, viewer);
2030: if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
2031: if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
2032: PetscCall(PetscDrawPopCurrentPoint(draw));
2033: #if PetscDefined(HAVE_SAWS)
2034: } else if (issaws) {
2035: PetscMPIInt rank;
2036: const char *name;
2038: PetscCall(PetscObjectGetName((PetscObject)ts, &name));
2039: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
2040: if (!((PetscObject)ts)->amsmem && rank == 0) {
2041: char dir[1024];
2043: PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
2044: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
2045: PetscCallSAWs(SAWs_Register, dir, &ts->steps, 1, SAWs_READ, SAWs_INT);
2046: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
2047: PetscCallSAWs(SAWs_Register, dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE);
2048: }
2049: PetscTryTypeMethod(ts, view, viewer);
2050: #endif
2051: }
2052: if (ts->snes && ts->usessnes) {
2053: PetscCall(PetscViewerASCIIPushTab(viewer));
2054: PetscCall(SNESView(ts->snes, viewer));
2055: PetscCall(PetscViewerASCIIPopTab(viewer));
2056: }
2057: PetscCall(DMGetDMTS(ts->dm, &sdm));
2058: PetscCall(DMTSView(sdm, viewer));
2060: PetscCall(PetscViewerASCIIPushTab(viewer));
2061: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &issundials));
2062: PetscCall(PetscViewerASCIIPopTab(viewer));
2063: PetscFunctionReturn(PETSC_SUCCESS);
2064: }
2066: /*@
2067: TSSetApplicationContext - Sets an optional user-defined context for the timesteppers that may be accessed, for example inside the user provided
2068: `TS` callbacks with `TSGetApplicationContext()`
2070: Logically Collective
2072: Input Parameters:
2073: + ts - the `TS` context obtained from `TSCreate()`
2074: - ctx - application context
2076: Level: intermediate
2078: Fortran Note:
2079: This only works when `ctx` is a Fortran derived type (it cannot be a `PetscObject`), we recommend writing a Fortran interface definition for this
2080: function that tells the Fortran compiler the derived data type that is passed in as the `ctx` argument. See `TSGetApplicationContext()` for
2081: an example.
2083: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
2084: @*/
2085: PetscErrorCode TSSetApplicationContext(TS ts, PetscCtx ctx)
2086: {
2087: PetscFunctionBegin;
2089: ts->ctx = ctx;
2090: PetscFunctionReturn(PETSC_SUCCESS);
2091: }
2093: /*@
2094: TSGetApplicationContext - Gets the user-defined context for the
2095: timestepper that was set with `TSSetApplicationContext()`
2097: Not Collective
2099: Input Parameter:
2100: . ts - the `TS` context obtained from `TSCreate()`
2102: Output Parameter:
2103: . ctx - a pointer to the application context
2105: Level: intermediate
2107: Fortran Notes:
2108: This only works when the context is a Fortran derived type or a `PetscObject`. Declare `ctx` with
2109: .vb
2110: type(tUsertype), pointer :: ctx
2111: .ve
2113: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
2114: @*/
2115: PetscErrorCode TSGetApplicationContext(TS ts, PetscCtxRt ctx)
2116: {
2117: PetscFunctionBegin;
2119: *(void **)ctx = ts->ctx;
2120: PetscFunctionReturn(PETSC_SUCCESS);
2121: }
2123: /*@
2124: TSGetStepNumber - Gets the number of time steps completed.
2126: Not Collective
2128: Input Parameter:
2129: . ts - the `TS` context obtained from `TSCreate()`
2131: Output Parameter:
2132: . steps - number of steps completed so far
2134: Level: intermediate
2136: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2137: @*/
2138: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2139: {
2140: PetscFunctionBegin;
2142: PetscAssertPointer(steps, 2);
2143: *steps = ts->steps;
2144: PetscFunctionReturn(PETSC_SUCCESS);
2145: }
2147: /*@
2148: TSSetStepNumber - Sets the number of steps completed.
2150: Logically Collective
2152: Input Parameters:
2153: + ts - the `TS` context
2154: - steps - number of steps completed so far
2156: Level: developer
2158: Note:
2159: For most uses of the `TS` solvers the user need not explicitly call
2160: `TSSetStepNumber()`, as the step counter is appropriately updated in
2161: `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2162: reinitialize timestepping by setting the step counter to zero (and time
2163: to the initial time) to solve a similar problem with different initial
2164: conditions or parameters. Other possible use case is to continue
2165: timestepping from a previously interrupted run in such a way that `TS`
2166: monitors will be called with a initial nonzero step counter.
2168: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2169: @*/
2170: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2171: {
2172: PetscFunctionBegin;
2175: PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
2176: ts->steps = steps;
2177: PetscFunctionReturn(PETSC_SUCCESS);
2178: }
2180: /*@
2181: TSSetTimeStep - Allows one to reset the timestep at any time.
2183: Logically Collective
2185: Input Parameters:
2186: + ts - the `TS` context obtained from `TSCreate()`
2187: - time_step - the size of the timestep
2189: Options Database Key:
2190: . -ts_time_step dt - provide the initial time step
2192: Level: intermediate
2194: Notes:
2195: This is only a suggestion, the actual initial time step used may differ
2197: If this is called after `TSSetUp()`, it will not change the initial time step value printed by `TSView()`
2199: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2200: @*/
2201: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2202: {
2203: PetscFunctionBegin;
2206: ts->time_step = time_step;
2207: if (ts->setupcalled == PETSC_FALSE) ts->initial_time_step = time_step;
2208: PetscFunctionReturn(PETSC_SUCCESS);
2209: }
2211: /*@
2212: TSSetExactFinalTime - Determines whether to adapt the final time step to
2213: match the exact final time, to interpolate the solution to the exact final time,
2214: or to just return at the final time `TS` computed (which may be slightly larger
2215: than the requested final time).
2217: Logically Collective
2219: Input Parameters:
2220: + ts - the time-step context
2221: - eftopt - exact final time option
2222: .vb
2223: TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded, just use it
2224: TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time if the final time is exceeded
2225: TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to ensure the computed final time exactly equals the requested final time
2226: .ve
2228: Options Database Key:
2229: . -ts_exact_final_time stepover,interpolate,matchstep - select the final step approach at runtime
2231: Level: beginner
2233: Note:
2234: If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2235: then the final time you selected.
2237: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2238: @*/
2239: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2240: {
2241: PetscFunctionBegin;
2244: ts->exact_final_time = eftopt;
2245: PetscFunctionReturn(PETSC_SUCCESS);
2246: }
2248: /*@
2249: TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`
2251: Not Collective
2253: Input Parameter:
2254: . ts - the `TS` context
2256: Output Parameter:
2257: . eftopt - exact final time option
2259: Level: beginner
2261: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2262: @*/
2263: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2264: {
2265: PetscFunctionBegin;
2267: PetscAssertPointer(eftopt, 2);
2268: *eftopt = ts->exact_final_time;
2269: PetscFunctionReturn(PETSC_SUCCESS);
2270: }
2272: /*@
2273: TSGetTimeStep - Gets the current timestep size.
2275: Not Collective
2277: Input Parameter:
2278: . ts - the `TS` context obtained from `TSCreate()`
2280: Output Parameter:
2281: . dt - the current timestep size
2283: Level: intermediate
2285: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2286: @*/
2287: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2288: {
2289: PetscFunctionBegin;
2291: PetscAssertPointer(dt, 2);
2292: *dt = ts->time_step;
2293: PetscFunctionReturn(PETSC_SUCCESS);
2294: }
2296: /*@
2297: TSGetSolution - Returns the solution at the present timestep. It
2298: is valid to call this routine inside the function that you are evaluating
2299: in order to move to the new timestep. This vector not changed until
2300: the solution at the next timestep has been calculated.
2302: Not Collective, but v returned is parallel if ts is parallel
2304: Input Parameter:
2305: . ts - the `TS` context obtained from `TSCreate()`
2307: Output Parameter:
2308: . v - the vector containing the solution
2310: Level: intermediate
2312: Note:
2313: If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
2314: final time. It returns the solution at the next timestep.
2316: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2317: @*/
2318: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2319: {
2320: PetscFunctionBegin;
2322: PetscAssertPointer(v, 2);
2323: *v = ts->vec_sol;
2324: PetscFunctionReturn(PETSC_SUCCESS);
2325: }
2327: /*@
2328: TSGetSolutionComponents - Returns any solution components at the present
2329: timestep, if available for the time integration method being used.
2330: Solution components are quantities that share the same size and
2331: structure as the solution vector.
2333: Not Collective, but v returned is parallel if ts is parallel
2335: Input Parameters:
2336: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2337: . n - If v is `NULL`, then the number of solution components is
2338: returned through n, else the n-th solution component is
2339: returned in v.
2340: - v - the vector containing the n-th solution component
2341: (may be `NULL` to use this function to find out
2342: the number of solutions components).
2344: Level: advanced
2346: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2347: @*/
2348: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2349: {
2350: PetscFunctionBegin;
2352: if (!ts->ops->getsolutioncomponents) *n = 0;
2353: else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2354: PetscFunctionReturn(PETSC_SUCCESS);
2355: }
2357: /*@
2358: TSGetAuxSolution - Returns an auxiliary solution at the present
2359: timestep, if available for the time integration method being used.
2361: Not Collective, but v returned is parallel if ts is parallel
2363: Input Parameters:
2364: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2365: - v - the vector containing the auxiliary solution
2367: Level: intermediate
2369: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
2370: @*/
2371: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2372: {
2373: PetscFunctionBegin;
2375: if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2376: else PetscCall(VecZeroEntries(*v));
2377: PetscFunctionReturn(PETSC_SUCCESS);
2378: }
2380: /*@
2381: TSGetTimeError - Returns the estimated error vector, if the chosen
2382: `TSType` has an error estimation functionality and `TSSetTimeError()` was called
2384: Not Collective, but v returned is parallel if ts is parallel
2386: Input Parameters:
2387: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2388: . n - current estimate (n=0) or previous one (n=-1)
2389: - v - the vector containing the error (same size as the solution).
2391: Level: intermediate
2393: Note:
2394: MUST call after `TSSetUp()`
2396: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
2397: @*/
2398: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2399: {
2400: PetscFunctionBegin;
2402: if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2403: else PetscCall(VecZeroEntries(*v));
2404: PetscFunctionReturn(PETSC_SUCCESS);
2405: }
2407: /*@
2408: TSSetTimeError - Sets the estimated error vector, if the chosen
2409: `TSType` has an error estimation functionality. This can be used
2410: to restart such a time integrator with a given error vector.
2412: Not Collective, but v returned is parallel if ts is parallel
2414: Input Parameters:
2415: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2416: - v - the vector containing the error (same size as the solution).
2418: Level: intermediate
2420: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
2421: @*/
2422: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2423: {
2424: PetscFunctionBegin;
2426: PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
2427: PetscTryTypeMethod(ts, settimeerror, v);
2428: PetscFunctionReturn(PETSC_SUCCESS);
2429: }
2431: /* ----- Routines to initialize and destroy a timestepper ---- */
2432: /*@
2433: TSSetProblemType - Sets the type of problem to be solved.
2435: Not collective
2437: Input Parameters:
2438: + ts - The `TS`
2439: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2440: .vb
2441: U_t - A U = 0 (linear)
2442: U_t - A(t) U = 0 (linear)
2443: F(t,U,U_t) = 0 (nonlinear)
2444: .ve
2446: Level: beginner
2448: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2449: @*/
2450: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2451: {
2452: PetscFunctionBegin;
2454: ts->problem_type = type;
2455: if (type == TS_LINEAR) {
2456: SNES snes;
2457: PetscCall(TSGetSNES(ts, &snes));
2458: PetscCall(SNESSetType(snes, SNESKSPONLY));
2459: }
2460: PetscFunctionReturn(PETSC_SUCCESS);
2461: }
2463: /*@
2464: TSGetProblemType - Gets the type of problem to be solved.
2466: Not collective
2468: Input Parameter:
2469: . ts - The `TS`
2471: Output Parameter:
2472: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2473: .vb
2474: M U_t = A U
2475: M(t) U_t = A(t) U
2476: F(t,U,U_t)
2477: .ve
2479: Level: beginner
2481: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
2482: @*/
2483: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2484: {
2485: PetscFunctionBegin;
2487: PetscAssertPointer(type, 2);
2488: *type = ts->problem_type;
2489: PetscFunctionReturn(PETSC_SUCCESS);
2490: }
2492: /*
2493: Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2494: */
2495: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2496: {
2497: PetscBool isnone;
2499: PetscFunctionBegin;
2500: PetscCall(TSGetAdapt(ts, &ts->adapt));
2501: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2503: PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
2504: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2505: else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2506: PetscFunctionReturn(PETSC_SUCCESS);
2507: }
2509: /*@
2510: TSSetUp - Sets up the internal data structures for the later use of a timestepper.
2512: Collective
2514: Input Parameter:
2515: . ts - the `TS` context obtained from `TSCreate()`
2517: Level: advanced
2519: Note:
2520: For basic use of the `TS` solvers the user need not explicitly call
2521: `TSSetUp()`, since these actions will automatically occur during
2522: the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this
2523: phase separately, `TSSetUp()` should be called after `TSCreate()`
2524: and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.
2526: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2527: @*/
2528: PetscErrorCode TSSetUp(TS ts)
2529: {
2530: DM dm;
2531: PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2532: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2533: TSIFunctionFn *ifun;
2534: TSIJacobianFn *ijac;
2535: TSI2JacobianFn *i2jac;
2536: TSRHSJacobianFn *rhsjac;
2538: PetscFunctionBegin;
2540: if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
2542: if (!((PetscObject)ts)->type_name) {
2543: PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
2544: PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
2545: }
2547: if (!ts->vec_sol) {
2548: PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
2549: PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
2550: }
2552: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2553: PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
2554: ts->Jacp = ts->Jacprhs;
2555: }
2557: if (ts->quadraturets) {
2558: PetscCall(TSSetUp(ts->quadraturets));
2559: PetscCall(VecDestroy(&ts->vec_costintegrand));
2560: PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
2561: }
2563: PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
2564: if (rhsjac == TSComputeRHSJacobianConstant) {
2565: Mat Amat, Pmat;
2566: SNES snes;
2567: PetscCall(TSGetSNES(ts, &snes));
2568: PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
2569: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2570: * have displaced the RHS matrix */
2571: if (Amat && Amat == ts->Arhs) {
2572: /* 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 */
2573: PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
2574: PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
2575: PetscCall(MatDestroy(&Amat));
2576: }
2577: if (Pmat && Pmat == ts->Brhs) {
2578: PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
2579: PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
2580: PetscCall(MatDestroy(&Pmat));
2581: }
2582: }
2584: PetscCall(TSGetAdapt(ts, &ts->adapt));
2585: PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
2587: PetscTryTypeMethod(ts, setup);
2589: PetscCall(TSSetExactFinalTimeDefault(ts));
2591: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2592: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2593: */
2594: PetscCall(TSGetDM(ts, &dm));
2595: PetscCall(DMSNESGetFunction(dm, &func, NULL));
2596: if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));
2598: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2599: Otherwise, the SNES will use coloring internally to form the Jacobian.
2600: */
2601: PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
2602: PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
2603: PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
2604: PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
2605: if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));
2607: /* if time integration scheme has a starting method, call it */
2608: PetscTryTypeMethod(ts, startingmethod);
2610: ts->setupcalled = PETSC_TRUE;
2611: PetscFunctionReturn(PETSC_SUCCESS);
2612: }
2614: /*@
2615: TSReset - Resets a `TS` context to the state it was in before `TSSetUp()` was called and removes any allocated `Vec` and `Mat` from its data structures
2617: Collective
2619: Input Parameter:
2620: . ts - the `TS` context obtained from `TSCreate()`
2622: Level: developer
2624: Notes:
2625: Any options set on the `TS` object, including those set with `TSSetFromOptions()` remain.
2627: See also `TSSetResize()` to change the size of the system being integrated (for example by adaptive mesh refinement) during the time integration.
2629: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSetResize()`
2630: @*/
2631: PetscErrorCode TSReset(TS ts)
2632: {
2633: TS_RHSSplitLink ilink = ts->tsrhssplit, next;
2635: PetscFunctionBegin;
2638: PetscTryTypeMethod(ts, reset);
2639: if (ts->snes) PetscCall(SNESReset(ts->snes));
2640: if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));
2642: PetscCall(MatDestroy(&ts->Arhs));
2643: PetscCall(MatDestroy(&ts->Brhs));
2644: PetscCall(VecDestroy(&ts->Frhs));
2645: PetscCall(VecDestroy(&ts->vec_sol));
2646: PetscCall(VecDestroy(&ts->vec_sol0));
2647: PetscCall(VecDestroy(&ts->vec_dot));
2648: PetscCall(VecDestroy(&ts->vatol));
2649: PetscCall(VecDestroy(&ts->vrtol));
2650: PetscCall(VecDestroyVecs(ts->nwork, &ts->work));
2652: PetscCall(MatDestroy(&ts->Jacprhs));
2653: PetscCall(MatDestroy(&ts->Jacp));
2654: if (ts->forward_solve) PetscCall(TSForwardReset(ts));
2655: if (ts->quadraturets) {
2656: PetscCall(TSReset(ts->quadraturets));
2657: PetscCall(VecDestroy(&ts->vec_costintegrand));
2658: }
2659: while (ilink) {
2660: next = ilink->next;
2661: PetscCall(TSDestroy(&ilink->ts));
2662: PetscCall(PetscFree(ilink->splitname));
2663: PetscCall(ISDestroy(&ilink->is));
2664: PetscCall(PetscFree(ilink));
2665: ilink = next;
2666: }
2667: ts->tsrhssplit = NULL;
2668: ts->num_rhs_splits = 0;
2669: if (ts->eval_times) {
2670: PetscCall(PetscFree(ts->eval_times->time_points));
2671: PetscCall(PetscFree(ts->eval_times->sol_times));
2672: PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
2673: PetscCall(PetscFree(ts->eval_times));
2674: }
2675: ts->rhsjacobian.time = PETSC_MIN_REAL;
2676: ts->rhsjacobian.scale = 1.0;
2677: ts->ijacobian.shift = 1.0;
2678: ts->setupcalled = PETSC_FALSE;
2679: PetscFunctionReturn(PETSC_SUCCESS);
2680: }
2682: static PetscErrorCode TSResizeReset(TS);
2684: /*@
2685: TSDestroy - Destroys the timestepper context that was created
2686: with `TSCreate()`.
2688: Collective
2690: Input Parameter:
2691: . ts - the `TS` context obtained from `TSCreate()`
2693: Level: beginner
2695: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2696: @*/
2697: PetscErrorCode TSDestroy(TS *ts)
2698: {
2699: PetscFunctionBegin;
2700: if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
2702: if (--((PetscObject)*ts)->refct > 0) {
2703: *ts = NULL;
2704: PetscFunctionReturn(PETSC_SUCCESS);
2705: }
2707: PetscCall(TSReset(*ts));
2708: PetscCall(TSAdjointReset(*ts));
2709: if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
2710: PetscCall(TSResizeReset(*ts));
2712: /* if memory was published with SAWs then destroy it */
2713: PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
2714: PetscTryTypeMethod(*ts, destroy);
2716: PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));
2718: PetscCall(TSAdaptDestroy(&(*ts)->adapt));
2719: PetscCall(TSEventDestroy(&(*ts)->event));
2721: PetscCall(SNESDestroy(&(*ts)->snes));
2722: PetscCall(SNESDestroy(&(*ts)->snesrhssplit));
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: TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the `TS` timestepping context
2735: Collective
2737: Input Parameters:
2738: + ts - the `TS` context obtained from `TSCreate()`
2739: - snes - the nonlinear solver context
2741: Level: developer
2743: Note:
2744: Most users should obtain the `SNES` by calling `TSGetSNES()` rather than setting it with this function.
2746: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetKSP()`, `TSIsImplicit()`, `TSGetSNES()`
2747: @*/
2748: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2749: {
2750: PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);
2752: PetscFunctionBegin;
2755: PetscCall(PetscObjectReference((PetscObject)snes));
2756: PetscCall(SNESDestroy(&ts->snes));
2757: ts->snes = snes;
2758: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2759: PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
2760: if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
2761: PetscFunctionReturn(PETSC_SUCCESS);
2762: }
2764: /*@
2765: TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2766: a `TS` (timestepper) context.
2768: Not Collective, but `snes` is parallel if `ts` is parallel
2770: Input Parameter:
2771: . ts - the `TS` context obtained from `TSCreate()`
2773: Output Parameter:
2774: . snes - the nonlinear solver context
2776: Level: beginner
2778: Notes:
2779: The user can then directly manipulate the `SNES` context to set various
2780: options, etc. Likewise, the user can then extract and manipulate the
2781: `KSP`, and `PC` contexts as well.
2783: For linear problems, use `TSGetKSP()`.
2785: For integrators that do not use `SNES` (that is, explicit methods),
2786: the `snes` exists but is not used. Use `TSIsImplicit()` to determine if the
2787: method is implicit and uses `snes`.
2789: Developer Note:
2790: `TS` manages the life-cycle of the `SNES` object for all `TSType` for the life-time of the `TS` object,
2791: even explicit methods that do not use `SNES`. This is so that `SNES` options are retained between changes to the `TSType` with `TSSetType()`.
2793: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetKSP()`, `TSIsImplicit()`
2794: @*/
2795: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2796: {
2797: PetscFunctionBegin;
2799: PetscAssertPointer(snes, 2);
2800: if (!ts->snes) {
2801: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
2802: PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
2803: PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
2804: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
2805: if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
2806: if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
2807: }
2808: *snes = ts->snes;
2809: PetscFunctionReturn(PETSC_SUCCESS);
2810: }
2812: /*@
2813: TSIsImplicit - Indicates if a `TS` represents an implicit integrator that uses `SNES`
2815: Not Collective
2817: Input Parameter:
2818: . ts - the `TS` context obtained from `TSCreate()`
2820: Output Parameter:
2821: . isimplicit - if the integrator is implicit and uses either `SNES` or `KSP`
2823: Level: beginner
2825: Note:
2826: For integrators that do not use `SNES` (that is, explicit methods), `snes` exists but is not used.
2828: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetKSP()`, `TSGetSNES()`
2829: @*/
2830: PetscErrorCode TSIsImplicit(TS ts, PetscBool *isimplicit)
2831: {
2832: PetscFunctionBegin;
2834: PetscAssertPointer(isimplicit, 2);
2835: *isimplicit = ts->usessnes;
2836: PetscFunctionReturn(PETSC_SUCCESS);
2837: }
2839: /*@
2840: TSGetKSP - Returns the `KSP` (linear solver) associated with
2841: a `TS` (timestepper) context.
2843: Not Collective, but `ksp` is parallel if `ts` is parallel
2845: Input Parameter:
2846: . ts - the `TS` context obtained from `TSCreate()`
2848: Output Parameter:
2849: . ksp - the nonlinear solver context
2851: Level: beginner
2853: Notes:
2854: The user can then directly manipulate the `KSP` context to set various
2855: options, etc. Likewise, the user can then extract and manipulate the
2856: `PC` context as well.
2858: For nonlinear problems (`TS_NONLINEAR`), use `TSGetSNES()` followed by `SNESGetKSP()`.
2860: For integrators that do not use `KSP` (that is, explicit methods),
2861: `TSGetKSP()` returns a `ksp` that is not used. Use `TSIsImplicit()` to determine if
2862: the `ksp` is actually used.
2864: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`, `TSIsImplicit()`
2865: @*/
2866: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2867: {
2868: SNES snes;
2870: PetscFunctionBegin;
2872: PetscAssertPointer(ksp, 2);
2873: PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
2874: PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "For linear problems only; use TSGetSNES() then SNESGetKSP()");
2875: PetscCall(TSGetSNES(ts, &snes));
2876: PetscCall(SNESGetKSP(snes, ksp));
2877: PetscFunctionReturn(PETSC_SUCCESS);
2878: }
2880: /* ----------- Routines to set solver parameters ---------- */
2882: /*@
2883: TSSetMaxSteps - Sets the maximum number of steps to use.
2885: Logically Collective
2887: Input Parameters:
2888: + ts - the `TS` context obtained from `TSCreate()`
2889: - maxsteps - maximum number of steps to use
2891: Options Database Key:
2892: . -ts_max_steps maxsteps - Sets maxsteps
2894: Level: intermediate
2896: Note:
2897: Use `PETSC_DETERMINE` to reset the maximum number of steps to the default from when the object's type was set
2899: The default maximum number of steps is 5,000
2901: Fortran Note:
2902: Use `PETSC_DETERMINE_INTEGER`
2904: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2905: @*/
2906: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2907: {
2908: PetscFunctionBegin;
2911: if (maxsteps == PETSC_DETERMINE) {
2912: ts->max_steps = ts->default_max_steps;
2913: } else {
2914: PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
2915: ts->max_steps = maxsteps;
2916: }
2917: PetscFunctionReturn(PETSC_SUCCESS);
2918: }
2920: /*@
2921: TSGetMaxSteps - Gets the maximum number of steps to use.
2923: Not Collective
2925: Input Parameter:
2926: . ts - the `TS` context obtained from `TSCreate()`
2928: Output Parameter:
2929: . maxsteps - maximum number of steps to use
2931: Level: advanced
2933: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2934: @*/
2935: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2936: {
2937: PetscFunctionBegin;
2939: PetscAssertPointer(maxsteps, 2);
2940: *maxsteps = ts->max_steps;
2941: PetscFunctionReturn(PETSC_SUCCESS);
2942: }
2944: /*@
2945: TSSetRunSteps - Sets the maximum number of steps to take in each call to `TSSolve()`.
2947: If the step count when `TSSolve()` is `start_step`, this will stop the simulation once `current_step - start_step >= run_steps`.
2948: Comparatively, `TSSetMaxSteps()` will stop if `current_step >= max_steps`.
2949: The simulation will stop when either condition is reached.
2951: Logically Collective
2953: Input Parameters:
2954: + ts - the `TS` context obtained from `TSCreate()`
2955: - runsteps - maximum number of steps to take in each call to `TSSolve()`;
2957: Options Database Key:
2958: . -ts_run_steps runsteps - Sets runsteps
2960: Level: intermediate
2962: Note:
2963: The default is `PETSC_UNLIMITED`
2965: .seealso: [](ch_ts), `TS`, `TSGetRunSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`, `TSSetMaxSteps()`
2966: @*/
2967: PetscErrorCode TSSetRunSteps(TS ts, PetscInt runsteps)
2968: {
2969: PetscFunctionBegin;
2972: if (runsteps == PETSC_DETERMINE) {
2973: ts->run_steps = PETSC_UNLIMITED;
2974: } else {
2975: PetscCheck(runsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Max number of steps to take in each call to TSSolve must be non-negative");
2976: ts->run_steps = runsteps;
2977: }
2978: PetscFunctionReturn(PETSC_SUCCESS);
2979: }
2981: /*@
2982: TSGetRunSteps - Gets the maximum number of steps to take in each call to `TSSolve()`.
2984: Not Collective
2986: Input Parameter:
2987: . ts - the `TS` context obtained from `TSCreate()`
2989: Output Parameter:
2990: . runsteps - maximum number of steps to take in each call to `TSSolve`.
2992: Level: advanced
2994: .seealso: [](ch_ts), `TS`, `TSSetRunSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`, `TSGetMaxSteps()`
2995: @*/
2996: PetscErrorCode TSGetRunSteps(TS ts, PetscInt *runsteps)
2997: {
2998: PetscFunctionBegin;
3000: PetscAssertPointer(runsteps, 2);
3001: *runsteps = ts->run_steps;
3002: PetscFunctionReturn(PETSC_SUCCESS);
3003: }
3005: /*@
3006: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
3008: Logically Collective
3010: Input Parameters:
3011: + ts - the `TS` context obtained from `TSCreate()`
3012: - maxtime - final time to step to
3014: Options Database Key:
3015: . -ts_max_time maxtime - Sets maxtime
3017: Level: intermediate
3019: Notes:
3020: Use `PETSC_DETERMINE` to reset the maximum time to the default from when the object's type was set
3022: The default maximum time is 5.0
3024: Fortran Note:
3025: Use `PETSC_DETERMINE_REAL`
3027: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
3028: @*/
3029: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
3030: {
3031: PetscFunctionBegin;
3034: if (maxtime == PETSC_DETERMINE) {
3035: ts->max_time = ts->default_max_time;
3036: } else {
3037: ts->max_time = maxtime;
3038: }
3039: PetscFunctionReturn(PETSC_SUCCESS);
3040: }
3042: /*@
3043: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
3045: Not Collective
3047: Input Parameter:
3048: . ts - the `TS` context obtained from `TSCreate()`
3050: Output Parameter:
3051: . maxtime - final time to step to
3053: Level: advanced
3055: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
3056: @*/
3057: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
3058: {
3059: PetscFunctionBegin;
3061: PetscAssertPointer(maxtime, 2);
3062: *maxtime = ts->max_time;
3063: PetscFunctionReturn(PETSC_SUCCESS);
3064: }
3066: // PetscClangLinter pragma disable: -fdoc-*
3067: /*@
3068: TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.
3070: Level: deprecated
3072: @*/
3073: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
3074: {
3075: PetscFunctionBegin;
3077: PetscCall(TSSetTime(ts, initial_time));
3078: PetscCall(TSSetTimeStep(ts, time_step));
3079: PetscFunctionReturn(PETSC_SUCCESS);
3080: }
3082: // PetscClangLinter pragma disable: -fdoc-*
3083: /*@
3084: TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.
3086: Level: deprecated
3088: @*/
3089: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
3090: {
3091: PetscFunctionBegin;
3093: if (maxsteps) {
3094: PetscAssertPointer(maxsteps, 2);
3095: *maxsteps = ts->max_steps;
3096: }
3097: if (maxtime) {
3098: PetscAssertPointer(maxtime, 3);
3099: *maxtime = ts->max_time;
3100: }
3101: PetscFunctionReturn(PETSC_SUCCESS);
3102: }
3104: // PetscClangLinter pragma disable: -fdoc-*
3105: /*@
3106: TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.
3108: Level: deprecated
3110: @*/
3111: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
3112: {
3113: PetscFunctionBegin;
3114: if (maxsteps != PETSC_CURRENT) PetscCall(TSSetMaxSteps(ts, maxsteps));
3115: if (maxtime != (PetscReal)PETSC_CURRENT) PetscCall(TSSetMaxTime(ts, maxtime));
3116: PetscFunctionReturn(PETSC_SUCCESS);
3117: }
3119: // PetscClangLinter pragma disable: -fdoc-*
3120: /*@
3121: TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.
3123: Level: deprecated
3125: @*/
3126: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
3127: {
3128: return TSGetStepNumber(ts, steps);
3129: }
3131: // PetscClangLinter pragma disable: -fdoc-*
3132: /*@
3133: TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.
3135: Level: deprecated
3137: @*/
3138: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
3139: {
3140: return TSGetStepNumber(ts, steps);
3141: }
3143: /*@
3144: TSSetSolution - Sets the initial solution vector
3145: for use by the `TS` routines.
3147: Logically Collective
3149: Input Parameters:
3150: + ts - the `TS` context obtained from `TSCreate()`
3151: - u - the solution vector
3153: Level: beginner
3155: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
3156: @*/
3157: PetscErrorCode TSSetSolution(TS ts, Vec u)
3158: {
3159: DM dm;
3161: PetscFunctionBegin;
3164: PetscCall(PetscObjectReference((PetscObject)u));
3165: PetscCall(VecDestroy(&ts->vec_sol));
3166: ts->vec_sol = u;
3168: PetscCall(TSGetDM(ts, &dm));
3169: PetscCall(DMShellSetGlobalVector(dm, u));
3170: PetscFunctionReturn(PETSC_SUCCESS);
3171: }
3173: /*@
3174: TSSetPreStep - Sets the general-purpose function
3175: called once at the beginning of each time step.
3177: Logically Collective
3179: Input Parameters:
3180: + ts - The `TS` context obtained from `TSCreate()`
3181: - func - The function
3183: Calling sequence of `func`:
3184: . ts - the `TS` context
3186: Level: intermediate
3188: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
3189: @*/
3190: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
3191: {
3192: PetscFunctionBegin;
3194: ts->prestep = func;
3195: PetscFunctionReturn(PETSC_SUCCESS);
3196: }
3198: /*@
3199: TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`
3201: Collective
3203: Input Parameter:
3204: . ts - The `TS` context obtained from `TSCreate()`
3206: Level: developer
3208: Note:
3209: `TSPreStep()` is typically used within time stepping implementations,
3210: so most users would not generally call this routine themselves.
3212: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3213: @*/
3214: PetscErrorCode TSPreStep(TS ts)
3215: {
3216: PetscFunctionBegin;
3218: if (ts->prestep) {
3219: Vec U;
3220: VecState Ustate;
3221: PetscBool same;
3223: PetscCall(TSGetSolution(ts, &U));
3224: PetscCall(VecGetState(U, &Ustate));
3225: PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3226: PetscCall(TSGetSolution(ts, &U));
3227: PetscCall(VecStateCompareUpdate(U, &Ustate, &same));
3228: if (!same) PetscCall(TSRestartStep(ts));
3229: }
3230: PetscFunctionReturn(PETSC_SUCCESS);
3231: }
3233: /*@
3234: TSSetPreStage - Sets the general-purpose function
3235: called once at the beginning of each stage.
3237: Logically Collective
3239: Input Parameters:
3240: + ts - The `TS` context obtained from `TSCreate()`
3241: - func - The function
3243: Calling sequence of `func`:
3244: + ts - the `TS` context
3245: - stagetime - the stage time
3247: Level: intermediate
3249: Note:
3250: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3251: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3252: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3254: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3255: @*/
3256: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime))
3257: {
3258: PetscFunctionBegin;
3260: ts->prestage = func;
3261: PetscFunctionReturn(PETSC_SUCCESS);
3262: }
3264: /*@
3265: TSSetPostStage - Sets the general-purpose function
3266: called once at the end of each stage.
3268: Logically Collective
3270: Input Parameters:
3271: + ts - The `TS` context obtained from `TSCreate()`
3272: - func - The function
3274: Calling sequence of `func`:
3275: + ts - the `TS` context
3276: . stagetime - the stage time
3277: . stageindex - the stage index
3278: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3280: Level: intermediate
3282: Note:
3283: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3284: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3285: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3287: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3288: @*/
3289: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
3290: {
3291: PetscFunctionBegin;
3293: ts->poststage = func;
3294: PetscFunctionReturn(PETSC_SUCCESS);
3295: }
3297: /*@
3298: TSSetPostEvaluate - Sets the general-purpose function
3299: called at the end of each step evaluation.
3301: Logically Collective
3303: Input Parameters:
3304: + ts - The `TS` context obtained from `TSCreate()`
3305: - func - The function
3307: Calling sequence of `func`:
3308: . ts - the `TS` context
3310: Level: intermediate
3312: Note:
3313: The function set by `TSSetPostEvaluate()` is called after the solution is evaluated, or after the step rollback.
3314: Inside the `func` callback, the solution vector can be obtained with `TSGetSolution()`, and modified, if need be.
3315: The time step can be obtained with `TSGetTimeStep()`, and the time at the start of the step - via `TSGetTime()`.
3316: The potential changes to the solution vector introduced by event handling (`postevent()`) are not relevant for `TSSetPostEvaluate()`,
3317: but are relevant for `TSSetPostStep()`, according to the function call scheme in `TSSolve()`, as shown below
3318: .vb
3319: ...
3320: Step()
3321: PostEvaluate()
3322: EventHandling()
3323: step_rollback ? PostEvaluate() : PostStep()
3324: ...
3325: .ve
3326: where EventHandling() may result in one of the following three outcomes
3327: .vb
3328: (1) | successful step | solution intact
3329: (2) | successful step | solution modified by `postevent()`
3330: (3) | step_rollback | solution rolled back
3331: .ve
3333: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3334: @*/
3335: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
3336: {
3337: PetscFunctionBegin;
3339: ts->postevaluate = func;
3340: PetscFunctionReturn(PETSC_SUCCESS);
3341: }
3343: /*@
3344: TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`
3346: Collective
3348: Input Parameters:
3349: + ts - The `TS` context obtained from `TSCreate()`
3350: - stagetime - The absolute time of the current stage
3352: Level: developer
3354: Note:
3355: `TSPreStage()` is typically used within time stepping implementations,
3356: most users would not generally call this routine themselves.
3358: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3359: @*/
3360: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3361: {
3362: PetscFunctionBegin;
3364: if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3365: PetscFunctionReturn(PETSC_SUCCESS);
3366: }
3368: /*@
3369: TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`
3371: Collective
3373: Input Parameters:
3374: + ts - The `TS` context obtained from `TSCreate()`
3375: . stagetime - The absolute time of the current stage
3376: . stageindex - Stage number
3377: - Y - Array of vectors (of size = total number of stages) with the stage solutions
3379: Level: developer
3381: Note:
3382: `TSPostStage()` is typically used within time stepping implementations,
3383: most users would not generally call this routine themselves.
3385: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3386: @*/
3387: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec Y[])
3388: {
3389: PetscFunctionBegin;
3391: if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3392: PetscFunctionReturn(PETSC_SUCCESS);
3393: }
3395: /*@
3396: TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`
3398: Collective
3400: Input Parameter:
3401: . ts - The `TS` context obtained from `TSCreate()`
3403: Level: developer
3405: Note:
3406: `TSPostEvaluate()` is typically used within time stepping implementations,
3407: most users would not generally call this routine themselves.
3409: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3410: @*/
3411: PetscErrorCode TSPostEvaluate(TS ts)
3412: {
3413: PetscFunctionBegin;
3415: if (ts->postevaluate) {
3416: Vec U;
3417: VecState Ustate;
3418: PetscBool same;
3420: PetscCall(TSGetSolution(ts, &U));
3421: PetscCall(VecGetState(U, &Ustate));
3422: PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3423: PetscCall(VecStateCompareUpdate(U, &Ustate, &same));
3424: if (!same) PetscCall(TSRestartStep(ts));
3425: }
3426: PetscFunctionReturn(PETSC_SUCCESS);
3427: }
3429: /*@
3430: TSSetPostStep - Sets the general-purpose function
3431: called once at the end of each successful time step.
3433: Logically Collective
3435: Input Parameters:
3436: + ts - The `TS` context obtained from `TSCreate()`
3437: - func - The function
3439: Calling sequence of `func`:
3440: . ts - the `TS` context
3442: Level: intermediate
3444: Note:
3445: The function set by `TSSetPostStep()` is called after each successful step. If the event handler locates an event at the
3446: given step, and `postevent()` modifies the solution vector, the solution vector obtained by `TSGetSolution()` inside `func` will
3447: contain the changes. To get the solution without these changes, use `TSSetPostEvaluate()` to set the appropriate callback.
3448: The scheme of the relevant function calls in `TSSolve()` is shown below
3449: .vb
3450: ...
3451: Step()
3452: PostEvaluate()
3453: EventHandling()
3454: step_rollback ? PostEvaluate() : PostStep()
3455: ...
3456: .ve
3457: where EventHandling() may result in one of the following three outcomes
3458: .vb
3459: (1) | successful step | solution intact
3460: (2) | successful step | solution modified by `postevent()`
3461: (3) | step_rollback | solution rolled back
3462: .ve
3464: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3465: @*/
3466: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
3467: {
3468: PetscFunctionBegin;
3470: ts->poststep = func;
3471: PetscFunctionReturn(PETSC_SUCCESS);
3472: }
3474: /*@
3475: TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`
3477: Collective
3479: Input Parameter:
3480: . ts - The `TS` context obtained from `TSCreate()`
3482: Note:
3483: `TSPostStep()` is typically used within time stepping implementations,
3484: so most users would not generally call this routine themselves.
3486: Level: developer
3488: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPostStep()`
3489: @*/
3490: PetscErrorCode TSPostStep(TS ts)
3491: {
3492: PetscFunctionBegin;
3494: if (ts->poststep) {
3495: Vec U;
3496: VecState Ustate;
3497: PetscBool same;
3499: PetscCall(TSGetSolution(ts, &U));
3500: PetscCall(VecGetState(U, &Ustate));
3501: PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3502: PetscCall(TSGetSolution(ts, &U));
3503: PetscCall(VecStateCompareUpdate(U, &Ustate, &same));
3504: if (!same) PetscCall(TSRestartStep(ts));
3505: }
3506: PetscFunctionReturn(PETSC_SUCCESS);
3507: }
3509: /*@
3510: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3512: Collective
3514: Input Parameters:
3515: + ts - time stepping context
3516: - t - time to interpolate to
3518: Output Parameter:
3519: . U - state at given time
3521: Level: intermediate
3523: Developer Notes:
3524: `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3526: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3527: @*/
3528: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3529: {
3530: PetscFunctionBegin;
3533: 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);
3534: PetscUseTypeMethod(ts, interpolate, t, U);
3535: PetscFunctionReturn(PETSC_SUCCESS);
3536: }
3538: /*@
3539: TSStep - Steps one time step
3541: Collective
3543: Input Parameter:
3544: . ts - the `TS` context obtained from `TSCreate()`
3546: Level: developer
3548: Notes:
3549: The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.
3551: The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
3552: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3554: This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
3555: time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.
3557: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3558: @*/
3559: PetscErrorCode TSStep(TS ts)
3560: {
3561: static PetscBool cite = PETSC_FALSE;
3562: PetscReal ptime;
3564: PetscFunctionBegin;
3566: PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3567: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3568: " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3569: " journal = {arXiv e-preprints},\n"
3570: " eprint = {1806.01437},\n"
3571: " archivePrefix = {arXiv},\n"
3572: " year = {2018}\n}\n",
3573: &cite));
3574: PetscCall(TSSetUp(ts));
3575: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
3576: if (ts->eval_times)
3577: ts->eval_times->worktol = 0; /* In each step of TSSolve() 'eval_times->worktol' will be meaningfully defined (later) only once:
3578: in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */
3580: PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->run_steps != PETSC_INT_MAX || ts->max_steps != PETSC_INT_MAX, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime(), TSSetMaxSteps(), or TSSetRunSteps() or use -ts_max_time <time>, -ts_max_steps <steps>, -ts_run_steps <steps>");
3581: 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()");
3582: 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");
3584: if (!ts->vec_sol0) PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sol0));
3585: PetscCall(VecCopy(ts->vec_sol, ts->vec_sol0));
3586: ts->time_step0 = ts->time_step;
3588: if (!ts->steps) ts->ptime_prev = ts->ptime;
3589: ptime = ts->ptime;
3591: ts->ptime_prev_rollback = ts->ptime_prev;
3592: ts->reason = TS_CONVERGED_ITERATING;
3594: PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
3595: PetscUseTypeMethod(ts, step);
3596: PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));
3598: if (ts->reason >= 0) {
3599: ts->ptime_prev = ptime;
3600: ts->steps++;
3601: ts->steprollback = PETSC_FALSE;
3602: ts->steprestart = PETSC_FALSE;
3603: ts->stepresize = PETSC_FALSE;
3604: }
3606: if (ts->reason < 0 && ts->errorifstepfailed) {
3607: PetscCall(TSMonitorCancel(ts));
3608: if (ts->usessnes && ts->snes) PetscCall(SNESMonitorCancel(ts->snes));
3609: PetscCheck(ts->reason != TS_DIVERGED_NONLINEAR_SOLVE, PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s, increase -ts_max_snes_failures or use unlimited to attempt recovery", TSConvergedReasons[ts->reason]);
3610: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
3611: }
3612: PetscFunctionReturn(PETSC_SUCCESS);
3613: }
3615: /*@
3616: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3617: at the end of a time step with a given order of accuracy.
3619: Collective
3621: Input Parameters:
3622: + ts - time stepping context
3623: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3625: Input/Output Parameter:
3626: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3627: on output, the actual order of the error evaluation
3629: Output Parameter:
3630: . wlte - the weighted local truncation error norm
3632: Level: advanced
3634: Note:
3635: If the timestepper cannot evaluate the error in a particular step
3636: (eg. in the first step or restart steps after event handling),
3637: this routine returns wlte=-1.0 .
3639: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3640: @*/
3641: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3642: {
3643: PetscFunctionBegin;
3647: if (order) PetscAssertPointer(order, 3);
3649: PetscAssertPointer(wlte, 4);
3650: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
3651: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3652: PetscFunctionReturn(PETSC_SUCCESS);
3653: }
3655: /*@
3656: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3658: Collective
3660: Input Parameters:
3661: + ts - time stepping context
3662: . order - desired order of accuracy
3663: - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)
3665: Output Parameter:
3666: . U - state at the end of the current step
3668: Level: advanced
3670: Notes:
3671: This function cannot be called until all stages have been evaluated.
3673: 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.
3675: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
3676: @*/
3677: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3678: {
3679: PetscFunctionBegin;
3683: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3684: PetscFunctionReturn(PETSC_SUCCESS);
3685: }
3687: /*@
3688: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3690: Not collective
3692: Input Parameter:
3693: . ts - time stepping context
3695: Output Parameter:
3696: . initCondition - The function which computes an initial condition
3698: Calling sequence of `initCondition`:
3699: + ts - The timestepping context
3700: - u - The input vector in which the initial condition is stored
3702: Level: advanced
3704: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3705: @*/
3706: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
3707: {
3708: PetscFunctionBegin;
3710: PetscAssertPointer(initCondition, 2);
3711: *initCondition = ts->ops->initcondition;
3712: PetscFunctionReturn(PETSC_SUCCESS);
3713: }
3715: /*@
3716: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3718: Logically collective
3720: Input Parameters:
3721: + ts - time stepping context
3722: - initCondition - The function which computes an initial condition
3724: Calling sequence of `initCondition`:
3725: + ts - The timestepping context
3726: - e - The input vector in which the initial condition is to be stored
3728: Level: advanced
3730: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3731: @*/
3732: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
3733: {
3734: PetscFunctionBegin;
3737: ts->ops->initcondition = initCondition;
3738: PetscFunctionReturn(PETSC_SUCCESS);
3739: }
3741: /*@
3742: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3744: Collective
3746: Input Parameters:
3747: + ts - time stepping context
3748: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3750: Level: advanced
3752: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3753: @*/
3754: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3755: {
3756: PetscFunctionBegin;
3759: PetscTryTypeMethod(ts, initcondition, u);
3760: PetscFunctionReturn(PETSC_SUCCESS);
3761: }
3763: /*@
3764: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3766: Not collective
3768: Input Parameter:
3769: . ts - time stepping context
3771: Output Parameter:
3772: . exactError - The function which computes the solution error
3774: Calling sequence of `exactError`:
3775: + ts - The timestepping context
3776: . u - The approximate solution vector
3777: - e - The vector in which the error is stored
3779: Level: advanced
3781: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
3782: @*/
3783: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
3784: {
3785: PetscFunctionBegin;
3787: PetscAssertPointer(exactError, 2);
3788: *exactError = ts->ops->exacterror;
3789: PetscFunctionReturn(PETSC_SUCCESS);
3790: }
3792: /*@
3793: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3795: Logically collective
3797: Input Parameters:
3798: + ts - time stepping context
3799: - exactError - The function which computes the solution error
3801: Calling sequence of `exactError`:
3802: + ts - The timestepping context
3803: . u - The approximate solution vector
3804: - e - The vector in which the error is stored
3806: Level: advanced
3808: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3809: @*/
3810: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
3811: {
3812: PetscFunctionBegin;
3815: ts->ops->exacterror = exactError;
3816: PetscFunctionReturn(PETSC_SUCCESS);
3817: }
3819: /*@
3820: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3822: Collective
3824: Input Parameters:
3825: + ts - time stepping context
3826: . u - The approximate solution
3827: - e - The `Vec` used to store the error
3829: Level: advanced
3831: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3832: @*/
3833: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3834: {
3835: PetscFunctionBegin;
3839: PetscTryTypeMethod(ts, exacterror, u, e);
3840: PetscFunctionReturn(PETSC_SUCCESS);
3841: }
3843: /*@
3844: TSSetResize - Sets the resize callbacks.
3846: Logically Collective
3848: Input Parameters:
3849: + ts - The `TS` context obtained from `TSCreate()`
3850: . rollback - Whether a resize will restart the step
3851: . setup - The setup function
3852: . transfer - The transfer function
3853: - ctx - [optional] The user-defined context
3855: Calling sequence of `setup`:
3856: + ts - the `TS` context
3857: . step - the current step
3858: . time - the current time
3859: . state - the current vector of state
3860: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
3861: - ctx - user defined context
3863: Calling sequence of `transfer`:
3864: + ts - the `TS` context
3865: . nv - the number of vectors to be transferred
3866: . vecsin - array of vectors to be transferred
3867: . vecsout - array of transferred vectors
3868: - ctx - user defined context
3870: Notes:
3871: The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
3872: depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
3873: and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
3874: that the size of the discrete problem has changed.
3875: In both cases, the solver will collect the needed vectors that will be
3876: transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
3877: current solution vector, and other vectors needed by the specific solver used.
3878: For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
3879: Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
3880: will be automatically reset if the sizes are changed and they must be specified again by the user
3881: inside the `transfer` function.
3882: The input and output arrays passed to `transfer` are allocated by PETSc.
3883: Vectors in `vecsout` must be created by the user.
3884: Ownership of vectors in `vecsout` is transferred to PETSc.
3886: Level: advanced
3888: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
3889: @*/
3890: PetscErrorCode TSSetResize(TS ts, PetscBool rollback, PetscErrorCode (*setup)(TS ts, PetscInt step, PetscReal time, Vec state, PetscBool *resize, PetscCtx ctx), PetscErrorCode (*transfer)(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], PetscCtx ctx), PetscCtx ctx)
3891: {
3892: PetscFunctionBegin;
3894: ts->resizerollback = rollback;
3895: ts->resizesetup = setup;
3896: ts->resizetransfer = transfer;
3897: ts->resizectx = ctx;
3898: PetscFunctionReturn(PETSC_SUCCESS);
3899: }
3901: /*
3902: TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.
3904: Collective
3906: Input Parameters:
3907: + ts - The `TS` context obtained from `TSCreate()`
3908: - flg - If `PETSC_TRUE` each TS implementation (e.g. `TSBDF`) will register vectors to be transferred, if `PETSC_FALSE` vectors will be imported from transferred vectors.
3910: Level: developer
3912: Note:
3913: `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
3914: used within time stepping implementations,
3915: so most users would not generally call this routine themselves.
3917: .seealso: [](ch_ts), `TS`, `TSSetResize()`
3918: @*/
3919: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
3920: {
3921: PetscFunctionBegin;
3923: PetscTryTypeMethod(ts, resizeregister, flg);
3924: /* PetscTryTypeMethod(adapt, resizeregister, flg); */
3925: PetscFunctionReturn(PETSC_SUCCESS);
3926: }
3928: static PetscErrorCode TSResizeReset(TS ts)
3929: {
3930: PetscFunctionBegin;
3932: PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
3933: PetscFunctionReturn(PETSC_SUCCESS);
3934: }
3936: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
3937: {
3938: PetscFunctionBegin;
3941: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
3942: if (ts->resizetransfer) {
3943: PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
3944: PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
3945: }
3946: for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
3947: PetscFunctionReturn(PETSC_SUCCESS);
3948: }
3950: /*@
3951: TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.
3953: Collective
3955: Input Parameters:
3956: + ts - The `TS` context obtained from `TSCreate()`
3957: . name - A string identifying the vector
3958: - vec - The vector
3960: Level: developer
3962: Note:
3963: `TSResizeRegisterVec()` is typically used within time stepping implementations,
3964: so most users would not generally call this routine themselves.
3966: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
3967: @*/
3968: PetscErrorCode TSResizeRegisterVec(TS ts, const char name[], Vec vec)
3969: {
3970: PetscFunctionBegin;
3972: PetscAssertPointer(name, 2);
3974: PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
3975: PetscFunctionReturn(PETSC_SUCCESS);
3976: }
3978: /*@
3979: TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.
3981: Collective
3983: Input Parameters:
3984: + ts - The `TS` context obtained from `TSCreate()`
3985: . name - A string identifying the vector
3986: - vec - The vector
3988: Level: developer
3990: Note:
3991: `TSResizeRetrieveVec()` is typically used within time stepping implementations,
3992: so most users would not generally call this routine themselves.
3994: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
3995: @*/
3996: PetscErrorCode TSResizeRetrieveVec(TS ts, const char name[], Vec *vec)
3997: {
3998: PetscFunctionBegin;
4000: PetscAssertPointer(name, 2);
4001: PetscAssertPointer(vec, 3);
4002: PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
4003: PetscFunctionReturn(PETSC_SUCCESS);
4004: }
4006: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
4007: {
4008: PetscInt cnt;
4009: PetscObjectList tmp;
4010: Vec *vecsin = NULL;
4011: const char **namesin = NULL;
4013: PetscFunctionBegin;
4014: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
4015: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
4016: if (names) PetscCall(PetscMalloc1(cnt, &namesin));
4017: if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
4018: for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
4019: if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
4020: if (vecs) vecsin[cnt] = (Vec)tmp->obj;
4021: if (names) namesin[cnt] = tmp->name;
4022: cnt++;
4023: }
4024: }
4025: if (nv) *nv = cnt;
4026: if (names) *names = namesin;
4027: if (vecs) *vecs = vecsin;
4028: PetscFunctionReturn(PETSC_SUCCESS);
4029: }
4031: /*@
4032: TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`
4034: Collective
4036: Input Parameter:
4037: . ts - The `TS` context obtained from `TSCreate()`
4039: Level: developer
4041: Note:
4042: `TSResize()` is typically used within time stepping implementations,
4043: so most users would not generally call this routine themselves.
4045: .seealso: [](ch_ts), `TS`, `TSSetResize()`
4046: @*/
4047: PetscErrorCode TSResize(TS ts)
4048: {
4049: PetscInt nv = 0;
4050: const char **names = NULL;
4051: Vec *vecsin = NULL;
4052: const char *solname = "ts:vec_sol";
4054: PetscFunctionBegin;
4056: if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
4057: if (ts->resizesetup) {
4058: PetscCall(VecLockReadPush(ts->vec_sol));
4059: PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &ts->stepresize, ts->resizectx));
4060: PetscCall(VecLockReadPop(ts->vec_sol));
4061: if (ts->stepresize) {
4062: if (ts->resizerollback) {
4063: PetscCall(TSRollBack(ts));
4064: ts->time_step = ts->time_step0;
4065: }
4066: PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
4067: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
4068: }
4069: }
4071: PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
4072: if (nv) {
4073: Vec *vecsout, vecsol;
4075: /* Reset internal objects */
4076: PetscCall(TSReset(ts));
4078: /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
4079: PetscCall(PetscCalloc1(nv, &vecsout));
4080: PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
4081: for (PetscInt i = 0; i < nv; i++) {
4082: const char *name;
4083: char *oname;
4085: PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
4086: PetscCall(PetscStrallocpy(name, &oname));
4087: PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
4088: if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
4089: PetscCall(PetscFree(oname));
4090: PetscCall(VecDestroy(&vecsout[i]));
4091: }
4092: PetscCall(PetscFree(vecsout));
4093: PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */
4095: PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
4096: if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
4097: PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
4098: }
4100: PetscCall(PetscFree(names));
4101: PetscCall(PetscFree(vecsin));
4102: PetscCall(TSResizeReset(ts));
4103: PetscFunctionReturn(PETSC_SUCCESS);
4104: }
4106: /*@
4107: TSSolve - Steps the requested number of timesteps.
4109: Collective
4111: Input Parameters:
4112: + ts - the `TS` context obtained from `TSCreate()`
4113: - u - the solution vector (can be `NULL` if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
4114: otherwise it must contain the initial conditions and will contain the solution at the final requested time
4116: Level: beginner
4118: Note:
4119: The final time returned by this function may be different from the time of the internally
4120: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
4121: stepped over the final time.
4123: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
4124: @*/
4125: PetscErrorCode TSSolve(TS ts, Vec u)
4126: {
4127: Vec solution;
4129: PetscFunctionBegin;
4133: PetscCall(TSSetExactFinalTimeDefault(ts));
4134: 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 */
4135: if (!ts->vec_sol || u == ts->vec_sol) {
4136: PetscCall(VecDuplicate(u, &solution));
4137: PetscCall(TSSetSolution(ts, solution));
4138: PetscCall(VecDestroy(&solution)); /* grant ownership */
4139: }
4140: PetscCall(VecCopy(u, ts->vec_sol));
4141: PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
4142: } else if (u) PetscCall(TSSetSolution(ts, u));
4143: PetscCall(TSSetUp(ts));
4144: PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
4146: PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->run_steps != PETSC_INT_MAX || ts->max_steps != PETSC_INT_MAX, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime(), TSSetMaxSteps(), or TSSetRunSteps() or use -ts_max_time <time>, -ts_max_steps <steps>, -ts_run_steps <steps>");
4147: 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()");
4148: 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");
4149: PetscCheck(!(ts->eval_times && ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP), PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "You must use TS_EXACTFINALTIME_MATCHSTEP when using time span or evaluation times");
4151: if (ts->eval_times) {
4152: if (!ts->eval_times->sol_vecs) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
4153: for (PetscInt i = 0; i < ts->eval_times->num_time_points; i++) {
4154: PetscBool is_close = PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[i], ts->eval_times->reltol * ts->time_step + ts->eval_times->abstol, 0);
4155: if (ts->ptime <= ts->eval_times->time_points[i] || is_close) {
4156: ts->eval_times->time_point_idx = i;
4158: PetscBool is_ptime_in_sol_times = PETSC_FALSE; // If current solution has already been saved, we should not save it again
4159: if (ts->eval_times->sol_idx > 0) is_ptime_in_sol_times = PetscIsCloseAtTol(ts->ptime, ts->eval_times->sol_times[ts->eval_times->sol_idx - 1], ts->eval_times->reltol * ts->time_step + ts->eval_times->abstol, 0);
4160: if (is_close && !is_ptime_in_sol_times) {
4161: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4162: ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4163: ts->eval_times->sol_idx++;
4164: ts->eval_times->time_point_idx++;
4165: }
4166: break;
4167: }
4168: }
4169: }
4171: if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));
4173: /* reset number of steps only when the step is not restarted. ARKIMEX
4174: restarts the step after an event. Resetting these counters in such case causes
4175: TSTrajectory to incorrectly save the output files
4176: */
4177: /* reset time step and iteration counters */
4178: if (!ts->steps) {
4179: ts->ksp_its = 0;
4180: ts->snes_its = 0;
4181: ts->num_snes_failures = 0;
4182: ts->reject = 0;
4183: ts->steprestart = PETSC_TRUE;
4184: ts->steprollback = PETSC_FALSE;
4185: ts->stepresize = PETSC_FALSE;
4186: ts->rhsjacobian.time = PETSC_MIN_REAL;
4187: }
4189: /* make sure initial time step does not overshoot final time or the next point in evaluation times */
4190: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
4191: PetscReal maxdt;
4192: PetscReal dt = ts->time_step;
4194: if (ts->eval_times) maxdt = ts->eval_times->time_points[ts->eval_times->time_point_idx] - ts->ptime;
4195: else maxdt = ts->max_time - ts->ptime;
4196: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
4197: }
4198: ts->reason = TS_CONVERGED_ITERATING;
4200: {
4201: PetscViewer viewer;
4202: PetscViewerFormat format;
4203: PetscBool flg;
4204: static PetscBool incall = PETSC_FALSE;
4206: if (!incall) {
4207: /* Estimate the convergence rate of the time discretization */
4208: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
4209: if (flg) {
4210: PetscConvEst conv;
4211: DM dm;
4212: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
4213: PetscInt Nf;
4214: PetscBool checkTemporal = PETSC_TRUE;
4216: incall = PETSC_TRUE;
4217: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
4218: PetscCall(TSGetDM(ts, &dm));
4219: PetscCall(DMGetNumFields(dm, &Nf));
4220: PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
4221: PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
4222: PetscCall(PetscConvEstUseTS(conv, checkTemporal));
4223: PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
4224: PetscCall(PetscConvEstSetFromOptions(conv));
4225: PetscCall(PetscConvEstSetUp(conv));
4226: PetscCall(PetscConvEstGetConvRate(conv, alpha));
4227: PetscCall(PetscViewerPushFormat(viewer, format));
4228: PetscCall(PetscConvEstRateView(conv, alpha, viewer));
4229: PetscCall(PetscViewerPopFormat(viewer));
4230: PetscCall(PetscViewerDestroy(&viewer));
4231: PetscCall(PetscConvEstDestroy(&conv));
4232: PetscCall(PetscFree(alpha));
4233: incall = PETSC_FALSE;
4234: }
4235: }
4236: }
4238: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));
4240: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
4241: PetscUseTypeMethod(ts, solve);
4242: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4243: ts->solvetime = ts->ptime;
4244: solution = ts->vec_sol;
4245: } else { /* Step the requested number of timesteps. */
4246: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4247: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4249: if (!ts->steps) {
4250: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4251: PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
4252: }
4254: ts->start_step = ts->steps; // records starting step
4255: while (!ts->reason) {
4256: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4257: if (!ts->steprollback || (ts->stepresize && ts->resizerollback)) PetscCall(TSPreStep(ts));
4258: PetscCall(TSStep(ts));
4259: if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
4260: if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
4261: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
4262: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4263: PetscCall(TSForwardCostIntegral(ts));
4264: if (ts->reason >= 0) ts->steps++;
4265: }
4266: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
4267: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
4268: PetscCall(TSForwardStep(ts));
4269: if (ts->reason >= 0) ts->steps++;
4270: }
4271: PetscCall(TSPostEvaluate(ts));
4272: 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. */
4273: if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
4274: if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
4275: /* check convergence */
4276: if (!ts->reason) {
4277: if ((ts->steps - ts->start_step) >= ts->run_steps) ts->reason = TS_CONVERGED_ITS;
4278: else if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
4279: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
4280: }
4281: if (!ts->steprollback) {
4282: PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
4283: PetscCall(TSPostStep(ts));
4284: if (!ts->resizerollback) PetscCall(TSResize(ts));
4286: if (ts->eval_times && ts->eval_times->time_point_idx < ts->eval_times->num_time_points && ts->reason >= 0) {
4287: PetscCheck(ts->eval_times->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(eval_times->worktol > 0) in TSSolve()");
4288: if (PetscIsCloseAtTol(ts->ptime, ts->eval_times->time_points[ts->eval_times->time_point_idx], ts->eval_times->worktol, 0)) {
4289: ts->eval_times->sol_times[ts->eval_times->sol_idx] = ts->ptime;
4290: PetscCall(VecCopy(ts->vec_sol, ts->eval_times->sol_vecs[ts->eval_times->sol_idx]));
4291: ts->eval_times->sol_idx++;
4292: ts->eval_times->time_point_idx++;
4293: }
4294: }
4295: }
4296: }
4297: PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
4299: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
4300: if (!u) u = ts->vec_sol;
4301: PetscCall(TSInterpolate(ts, ts->max_time, u));
4302: ts->solvetime = ts->max_time;
4303: solution = u;
4304: PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
4305: } else {
4306: if (u) PetscCall(VecCopy(ts->vec_sol, u));
4307: ts->solvetime = ts->ptime;
4308: solution = ts->vec_sol;
4309: }
4310: }
4312: PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
4313: PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
4314: PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
4315: if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
4316: PetscFunctionReturn(PETSC_SUCCESS);
4317: }
4319: /*@
4320: TSGetTime - Gets the time of the most recently completed step.
4322: Not Collective
4324: Input Parameter:
4325: . ts - the `TS` context obtained from `TSCreate()`
4327: Output Parameter:
4328: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
4330: Level: beginner
4332: Note:
4333: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
4334: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
4336: .seealso: [](ch_ts), `TS`, `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
4337: @*/
4338: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
4339: {
4340: PetscFunctionBegin;
4342: PetscAssertPointer(t, 2);
4343: *t = ts->ptime;
4344: PetscFunctionReturn(PETSC_SUCCESS);
4345: }
4347: /*@
4348: TSGetPrevTime - Gets the starting time of the previously completed step.
4350: Not Collective
4352: Input Parameter:
4353: . ts - the `TS` context obtained from `TSCreate()`
4355: Output Parameter:
4356: . t - the previous time
4358: Level: beginner
4360: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
4361: @*/
4362: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
4363: {
4364: PetscFunctionBegin;
4366: PetscAssertPointer(t, 2);
4367: *t = ts->ptime_prev;
4368: PetscFunctionReturn(PETSC_SUCCESS);
4369: }
4371: /*@
4372: TSSetTime - Allows one to reset the time.
4374: Logically Collective
4376: Input Parameters:
4377: + ts - the `TS` context obtained from `TSCreate()`
4378: - t - the time
4380: Level: intermediate
4382: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
4383: @*/
4384: PetscErrorCode TSSetTime(TS ts, PetscReal t)
4385: {
4386: PetscFunctionBegin;
4389: ts->ptime = t;
4390: PetscFunctionReturn(PETSC_SUCCESS);
4391: }
4393: /*@
4394: TSSetOptionsPrefix - Sets the prefix used for searching for all
4395: TS options in the database.
4397: Logically Collective
4399: Input Parameters:
4400: + ts - The `TS` context
4401: - prefix - The prefix to prepend to all option names
4403: Level: advanced
4405: Note:
4406: A hyphen (-) must NOT be given at the beginning of the prefix name.
4407: The first character of all runtime options is AUTOMATICALLY the
4408: hyphen.
4410: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
4411: @*/
4412: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
4413: {
4414: SNES snes;
4416: PetscFunctionBegin;
4418: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
4419: PetscCall(TSGetSNES(ts, &snes));
4420: PetscCall(SNESSetOptionsPrefix(snes, prefix));
4421: PetscFunctionReturn(PETSC_SUCCESS);
4422: }
4424: /*@
4425: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
4426: TS options in the database.
4428: Logically Collective
4430: Input Parameters:
4431: + ts - The `TS` context
4432: - prefix - The prefix to prepend to all option names
4434: Level: advanced
4436: Note:
4437: A hyphen (-) must NOT be given at the beginning of the prefix name.
4438: The first character of all runtime options is AUTOMATICALLY the
4439: hyphen.
4441: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
4442: @*/
4443: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
4444: {
4445: SNES snes;
4447: PetscFunctionBegin;
4449: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
4450: PetscCall(TSGetSNES(ts, &snes));
4451: PetscCall(SNESAppendOptionsPrefix(snes, prefix));
4452: PetscFunctionReturn(PETSC_SUCCESS);
4453: }
4455: /*@
4456: TSGetOptionsPrefix - Sets the prefix used for searching for all
4457: `TS` options in the database.
4459: Not Collective
4461: Input Parameter:
4462: . ts - The `TS` context
4464: Output Parameter:
4465: . prefix - A pointer to the prefix string used
4467: Level: intermediate
4469: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
4470: @*/
4471: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
4472: {
4473: PetscFunctionBegin;
4475: PetscAssertPointer(prefix, 2);
4476: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
4477: PetscFunctionReturn(PETSC_SUCCESS);
4478: }
4480: /*@
4481: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
4483: Not Collective, but parallel objects are returned if ts is parallel
4485: Input Parameter:
4486: . ts - The `TS` context obtained from `TSCreate()`
4488: Output Parameters:
4489: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`)
4490: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`)
4491: . func - Function to compute the Jacobian of the RHS (or `NULL`)
4492: - ctx - User-defined context for Jacobian evaluation routine (or `NULL`)
4494: Level: intermediate
4496: Note:
4497: You can pass in `NULL` for any return argument you do not need.
4499: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4500: @*/
4501: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, PetscCtxRt ctx)
4502: {
4503: DM dm;
4505: PetscFunctionBegin;
4506: if (Amat || Pmat) {
4507: SNES snes;
4508: PetscCall(TSGetSNES(ts, &snes));
4509: PetscCall(SNESSetUpMatrices(snes));
4510: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4511: }
4512: PetscCall(TSGetDM(ts, &dm));
4513: PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
4514: PetscFunctionReturn(PETSC_SUCCESS);
4515: }
4517: /*@
4518: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
4520: Not Collective, but parallel objects are returned if ts is parallel
4522: Input Parameter:
4523: . ts - The `TS` context obtained from `TSCreate()`
4525: Output Parameters:
4526: + Amat - The (approximate) Jacobian of F(t,U,U_t)
4527: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
4528: . f - The function to compute the matrices
4529: - ctx - User-defined context for Jacobian evaluation routine
4531: Level: advanced
4533: Note:
4534: You can pass in `NULL` for any return argument you do not need.
4536: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
4537: @*/
4538: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, PetscCtxRt ctx)
4539: {
4540: DM dm;
4542: PetscFunctionBegin;
4543: if (Amat || Pmat) {
4544: SNES snes;
4545: PetscCall(TSGetSNES(ts, &snes));
4546: PetscCall(SNESSetUpMatrices(snes));
4547: PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
4548: }
4549: PetscCall(TSGetDM(ts, &dm));
4550: PetscCall(DMTSGetIJacobian(dm, f, ctx));
4551: PetscFunctionReturn(PETSC_SUCCESS);
4552: }
4554: #include <petsc/private/dmimpl.h>
4555: /*@
4556: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4558: Logically Collective
4560: Input Parameters:
4561: + ts - the `TS` integrator object
4562: - dm - the dm, cannot be `NULL`
4564: Level: intermediate
4566: Notes:
4567: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4568: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4569: different problems using the same function space.
4571: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4572: @*/
4573: PetscErrorCode TSSetDM(TS ts, DM dm)
4574: {
4575: SNES snes;
4576: DMTS tsdm;
4578: PetscFunctionBegin;
4581: PetscCall(PetscObjectReference((PetscObject)dm));
4582: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4583: if (ts->dm->dmts && !dm->dmts) {
4584: PetscCall(DMCopyDMTS(ts->dm, dm));
4585: PetscCall(DMGetDMTS(ts->dm, &tsdm));
4586: /* Grant write privileges to the replacement DM */
4587: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4588: }
4589: PetscCall(DMDestroy(&ts->dm));
4590: }
4591: ts->dm = dm;
4593: PetscCall(TSGetSNES(ts, &snes));
4594: PetscCall(SNESSetDM(snes, dm));
4595: PetscFunctionReturn(PETSC_SUCCESS);
4596: }
4598: /*@
4599: TSGetDM - Gets the `DM` that may be used by some preconditioners
4601: Not Collective
4603: Input Parameter:
4604: . ts - the `TS`
4606: Output Parameter:
4607: . dm - the `DM`
4609: Level: intermediate
4611: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4612: @*/
4613: PetscErrorCode TSGetDM(TS ts, DM *dm)
4614: {
4615: PetscFunctionBegin;
4617: if (!ts->dm) {
4618: PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
4619: if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
4620: }
4621: *dm = ts->dm;
4622: PetscFunctionReturn(PETSC_SUCCESS);
4623: }
4625: /*@
4626: SNESTSFormFunction - Function to evaluate nonlinear residual defined by an ODE solver algorithm implemented within `TS`
4628: Logically Collective
4630: Input Parameters:
4631: + snes - nonlinear solver
4632: . U - the current state at which to evaluate the residual
4633: - ctx - application context, must be a `TS`
4635: Output Parameter:
4636: . F - the nonlinear residual
4638: Level: developer
4640: Note:
4641: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4642: It is most frequently passed to `MatFDColoringSetFunction()`.
4644: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4645: @*/
4646: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, PetscCtx ctx)
4647: {
4648: TS ts = (TS)ctx;
4650: PetscFunctionBegin;
4655: PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
4656: PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
4657: PetscFunctionReturn(PETSC_SUCCESS);
4658: }
4660: /*@
4661: SNESTSFormJacobian - Function to evaluate the Jacobian defined by an ODE solver algorithm implemented within `TS`
4663: Collective
4665: Input Parameters:
4666: + snes - nonlinear solver
4667: . U - the current state at which to evaluate the residual
4668: - ctx - application context, must be a `TS`
4670: Output Parameters:
4671: + A - the Jacobian
4672: - B - the matrix used to construct the preconditioner (often the same as `A`)
4674: Level: developer
4676: Note:
4677: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4679: .seealso: [](ch_ts), `SNESSetJacobian()`
4680: @*/
4681: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, PetscCtx ctx)
4682: {
4683: TS ts = (TS)ctx;
4685: PetscFunctionBegin;
4691: PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
4692: PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
4693: PetscFunctionReturn(PETSC_SUCCESS);
4694: }
4696: /*@
4697: TSComputeRHSFunctionLinear - Evaluate the right-hand side via the user-provided Jacobian, for linear problems Udot = A U only
4699: Collective
4701: Input Parameters:
4702: + ts - time stepping context
4703: . t - time at which to evaluate
4704: . U - state at which to evaluate
4705: - ctx - context
4707: Output Parameter:
4708: . F - right-hand side
4710: Level: intermediate
4712: Note:
4713: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right-hand side for linear problems.
4714: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4716: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4717: @*/
4718: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, PetscCtx ctx)
4719: {
4720: Mat Arhs, Brhs;
4722: PetscFunctionBegin;
4723: PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
4724: /* undo the damage caused by shifting */
4725: PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
4726: PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
4727: PetscCall(MatMult(Arhs, U, F));
4728: PetscFunctionReturn(PETSC_SUCCESS);
4729: }
4731: /*@
4732: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
4734: Collective
4736: Input Parameters:
4737: + ts - time stepping context
4738: . t - time at which to evaluate
4739: . U - state at which to evaluate
4740: - ctx - context
4742: Output Parameters:
4743: + A - Jacobian
4744: - B - matrix used to construct the preconditioner, often the same as `A`
4746: Level: intermediate
4748: Note:
4749: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4751: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4752: @*/
4753: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, PetscCtx ctx)
4754: {
4755: PetscFunctionBegin;
4756: PetscFunctionReturn(PETSC_SUCCESS);
4757: }
4759: /*@
4760: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4762: Collective
4764: Input Parameters:
4765: + ts - time stepping context
4766: . t - time at which to evaluate
4767: . U - state at which to evaluate
4768: . Udot - time derivative of state vector
4769: - ctx - context
4771: Output Parameter:
4772: . F - left hand side
4774: Level: intermediate
4776: Notes:
4777: 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
4778: user is required to write their own `TSComputeIFunction()`.
4779: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4780: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4782: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4784: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4785: @*/
4786: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
4787: {
4788: Mat A, B;
4790: PetscFunctionBegin;
4791: PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
4792: PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
4793: PetscCall(MatMult(A, Udot, F));
4794: PetscFunctionReturn(PETSC_SUCCESS);
4795: }
4797: /*@
4798: TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE
4800: Collective
4802: Input Parameters:
4803: + ts - time stepping context
4804: . t - time at which to evaluate
4805: . U - state at which to evaluate
4806: . Udot - time derivative of state vector
4807: . shift - shift to apply
4808: - ctx - context
4810: Output Parameters:
4811: + A - pointer to operator
4812: - B - pointer to matrix from which the preconditioner is built (often `A`)
4814: Level: advanced
4816: Notes:
4817: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4819: It is only appropriate for problems of the form
4821: $$
4822: M \dot{U} = F(U,t)
4823: $$
4825: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4826: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4827: an implicit operator of the form
4829: $$
4830: shift M + J
4831: $$
4833: 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
4834: a copy of M or reassemble it when requested.
4836: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4837: @*/
4838: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscCtx ctx)
4839: {
4840: PetscFunctionBegin;
4841: PetscCall(MatScale(A, shift / ts->ijacobian.shift));
4842: ts->ijacobian.shift = shift;
4843: PetscFunctionReturn(PETSC_SUCCESS);
4844: }
4846: /*@
4847: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4849: Not Collective
4851: Input Parameter:
4852: . ts - the `TS` context
4854: Output Parameter:
4855: . equation_type - see `TSEquationType`
4857: Level: beginner
4859: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4860: @*/
4861: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4862: {
4863: PetscFunctionBegin;
4865: PetscAssertPointer(equation_type, 2);
4866: *equation_type = ts->equation_type;
4867: PetscFunctionReturn(PETSC_SUCCESS);
4868: }
4870: /*@
4871: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4873: Not Collective
4875: Input Parameters:
4876: + ts - the `TS` context
4877: - equation_type - see `TSEquationType`
4879: Level: advanced
4881: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4882: @*/
4883: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4884: {
4885: PetscFunctionBegin;
4887: ts->equation_type = equation_type;
4888: PetscFunctionReturn(PETSC_SUCCESS);
4889: }
4891: /*@
4892: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4894: Not Collective
4896: Input Parameter:
4897: . ts - the `TS` context
4899: Output Parameter:
4900: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4901: manual pages for the individual convergence tests for complete lists
4903: Level: beginner
4905: Note:
4906: Can only be called after the call to `TSSolve()` is complete.
4908: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4909: @*/
4910: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4911: {
4912: PetscFunctionBegin;
4914: PetscAssertPointer(reason, 2);
4915: *reason = ts->reason;
4916: PetscFunctionReturn(PETSC_SUCCESS);
4917: }
4919: /*@
4920: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4922: Logically Collective; reason must contain common value
4924: Input Parameters:
4925: + ts - the `TS` context
4926: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4927: manual pages for the individual convergence tests for complete lists
4929: Level: advanced
4931: Note:
4932: Can only be called while `TSSolve()` is active.
4934: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4935: @*/
4936: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4937: {
4938: PetscFunctionBegin;
4940: ts->reason = reason;
4941: PetscFunctionReturn(PETSC_SUCCESS);
4942: }
4944: /*@
4945: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4947: Not Collective
4949: Input Parameter:
4950: . ts - the `TS` context
4952: Output Parameter:
4953: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4955: Level: beginner
4957: Note:
4958: Can only be called after the call to `TSSolve()` is complete.
4960: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4961: @*/
4962: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4963: {
4964: PetscFunctionBegin;
4966: PetscAssertPointer(ftime, 2);
4967: *ftime = ts->solvetime;
4968: PetscFunctionReturn(PETSC_SUCCESS);
4969: }
4971: /*@
4972: TSGetSNESIterations - Gets the total number of nonlinear iterations
4973: used by the time integrator.
4975: Not Collective
4977: Input Parameter:
4978: . ts - `TS` context
4980: Output Parameter:
4981: . nits - number of nonlinear iterations
4983: Level: intermediate
4985: Note:
4986: This counter is reset to zero for each successive call to `TSSolve()`.
4988: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4989: @*/
4990: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4991: {
4992: PetscFunctionBegin;
4994: PetscAssertPointer(nits, 2);
4995: *nits = ts->snes_its;
4996: PetscFunctionReturn(PETSC_SUCCESS);
4997: }
4999: /*@
5000: TSGetKSPIterations - Gets the total number of linear iterations
5001: used by the time integrator.
5003: Not Collective
5005: Input Parameter:
5006: . ts - `TS` context
5008: Output Parameter:
5009: . lits - number of linear iterations
5011: Level: intermediate
5013: Note:
5014: This counter is reset to zero for each successive call to `TSSolve()`.
5016: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`
5017: @*/
5018: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
5019: {
5020: PetscFunctionBegin;
5022: PetscAssertPointer(lits, 2);
5023: *lits = ts->ksp_its;
5024: PetscFunctionReturn(PETSC_SUCCESS);
5025: }
5027: /*@
5028: TSGetStepRejections - Gets the total number of rejected steps.
5030: Not Collective
5032: Input Parameter:
5033: . ts - `TS` context
5035: Output Parameter:
5036: . rejects - number of steps rejected
5038: Level: intermediate
5040: Note:
5041: This counter is reset to zero for each successive call to `TSSolve()`.
5043: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
5044: @*/
5045: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
5046: {
5047: PetscFunctionBegin;
5049: PetscAssertPointer(rejects, 2);
5050: *rejects = ts->reject;
5051: PetscFunctionReturn(PETSC_SUCCESS);
5052: }
5054: /*@
5055: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
5057: Not Collective
5059: Input Parameter:
5060: . ts - `TS` context
5062: Output Parameter:
5063: . fails - number of failed nonlinear solves
5065: Level: intermediate
5067: Note:
5068: This counter is reset to zero for each successive call to `TSSolve()`.
5070: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
5071: @*/
5072: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
5073: {
5074: PetscFunctionBegin;
5076: PetscAssertPointer(fails, 2);
5077: *fails = ts->num_snes_failures;
5078: PetscFunctionReturn(PETSC_SUCCESS);
5079: }
5081: /*@
5082: TSSetMaxStepRejections - Sets the maximum number of step rejections allowed in a single time-step attempt before a time step fails in `TSSolve()` with `TS_DIVERGED_STEP_REJECTED`
5084: Not Collective
5086: Input Parameters:
5087: + ts - `TS` context
5088: - rejects - maximum number of rejected steps, pass `PETSC_UNLIMITED` for unlimited
5090: Options Database Key:
5091: . -ts_max_step_rejections - Maximum number of step rejections before a step fails
5093: Level: intermediate
5095: Developer Note:
5096: The options database name is incorrect.
5098: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`,
5099: `TSGetConvergedReason()`, `TSSolve()`, `TS_DIVERGED_STEP_REJECTED`
5100: @*/
5101: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
5102: {
5103: PetscFunctionBegin;
5105: if (rejects == PETSC_UNLIMITED || rejects == -1) {
5106: ts->max_reject = PETSC_UNLIMITED;
5107: } else {
5108: PetscCheck(rejects >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of rejections");
5109: ts->max_reject = rejects;
5110: }
5111: PetscFunctionReturn(PETSC_SUCCESS);
5112: }
5114: /*@
5115: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves allowed before `TSSolve()` is ended with a `TSConvergedReason` of `TS_DIVERGED_NONLINEAR_SOLVE`
5117: Not Collective
5119: Input Parameters:
5120: + ts - `TS` context
5121: - fails - maximum number of failed nonlinear solves, pass `PETSC_UNLIMITED` to allow any number of failures.
5123: Options Database Key:
5124: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
5126: Level: intermediate
5128: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`,
5129: `TSGetConvergedReason()`, `TS_DIVERGED_NONLINEAR_SOLVE`, `TSConvergedReason`
5130: @*/
5131: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
5132: {
5133: PetscFunctionBegin;
5135: if (fails == PETSC_UNLIMITED || fails == -1) {
5136: ts->max_snes_failures = PETSC_UNLIMITED;
5137: } else {
5138: PetscCheck(fails >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures");
5139: ts->max_snes_failures = fails;
5140: }
5141: PetscFunctionReturn(PETSC_SUCCESS);
5142: }
5144: /*@
5145: TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`
5147: Not Collective
5149: Input Parameters:
5150: + ts - `TS` context
5151: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
5153: Options Database Key:
5154: . -ts_error_if_step_fails - Error if no step succeeds
5156: Level: intermediate
5158: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
5159: @*/
5160: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
5161: {
5162: PetscFunctionBegin;
5164: ts->errorifstepfailed = err;
5165: PetscFunctionReturn(PETSC_SUCCESS);
5166: }
5168: /*@
5169: TSGetAdapt - Get the adaptive controller context for the current method
5171: Collective if controller has not yet been created
5173: Input Parameter:
5174: . ts - time stepping context
5176: Output Parameter:
5177: . adapt - adaptive controller
5179: Level: intermediate
5181: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
5182: @*/
5183: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
5184: {
5185: PetscFunctionBegin;
5187: PetscAssertPointer(adapt, 2);
5188: if (!ts->adapt) {
5189: PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
5190: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
5191: }
5192: *adapt = ts->adapt;
5193: PetscFunctionReturn(PETSC_SUCCESS);
5194: }
5196: /*@
5197: TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller
5199: Logically Collective
5201: Input Parameters:
5202: + ts - time integration context
5203: . atol - scalar absolute tolerances
5204: . vatol - vector of absolute tolerances or `NULL`, used in preference to `atol` if present
5205: . rtol - scalar relative tolerances
5206: - vrtol - vector of relative tolerances or `NULL`, used in preference to `rtol` if present
5208: Options Database Keys:
5209: + -ts_rtol rtol - relative tolerance for local truncation error
5210: - -ts_atol atol - Absolute tolerance for local truncation error
5212: Level: beginner
5214: Notes:
5215: `PETSC_CURRENT` or `PETSC_DETERMINE` may be used for `atol` or `rtol` to indicate the current value
5216: or the default value from when the object's type was set.
5218: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
5219: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
5220: computed only for the differential or the algebraic part then this can be done using the vector of
5221: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
5222: differential part and infinity for the algebraic part, the LTE calculation will include only the
5223: differential variables.
5225: Fortran Note:
5226: Use `PETSC_CURRENT_INTEGER` or `PETSC_DETERMINE_INTEGER`.
5228: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
5229: @*/
5230: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
5231: {
5232: PetscFunctionBegin;
5233: if (atol == (PetscReal)PETSC_DETERMINE) {
5234: ts->atol = ts->default_atol;
5235: } else if (atol != (PetscReal)PETSC_CURRENT) {
5236: PetscCheck(atol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)atol);
5237: ts->atol = atol;
5238: }
5240: if (vatol) {
5241: PetscCall(PetscObjectReference((PetscObject)vatol));
5242: PetscCall(VecDestroy(&ts->vatol));
5243: ts->vatol = vatol;
5244: }
5246: if (rtol == (PetscReal)PETSC_DETERMINE) {
5247: ts->rtol = ts->default_rtol;
5248: } else if (rtol != (PetscReal)PETSC_CURRENT) {
5249: PetscCheck(rtol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative", (double)rtol);
5250: ts->rtol = rtol;
5251: }
5253: if (vrtol) {
5254: PetscCall(PetscObjectReference((PetscObject)vrtol));
5255: PetscCall(VecDestroy(&ts->vrtol));
5256: ts->vrtol = vrtol;
5257: }
5258: PetscFunctionReturn(PETSC_SUCCESS);
5259: }
5261: /*@
5262: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
5264: Logically Collective
5266: Input Parameter:
5267: . ts - time integration context
5269: Output Parameters:
5270: + atol - scalar absolute tolerances, `NULL` to ignore
5271: . vatol - vector of absolute tolerances, `NULL` to ignore
5272: . rtol - scalar relative tolerances, `NULL` to ignore
5273: - vrtol - vector of relative tolerances, `NULL` to ignore
5275: Level: beginner
5277: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
5278: @*/
5279: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
5280: {
5281: PetscFunctionBegin;
5282: if (atol) *atol = ts->atol;
5283: if (vatol) *vatol = ts->vatol;
5284: if (rtol) *rtol = ts->rtol;
5285: if (vrtol) *vrtol = ts->vrtol;
5286: PetscFunctionReturn(PETSC_SUCCESS);
5287: }
5289: /*@
5290: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
5292: Collective
5294: Input Parameters:
5295: + ts - time stepping context
5296: . U - state vector, usually ts->vec_sol
5297: . Y - state vector to be compared to U
5298: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5300: Output Parameters:
5301: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5302: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5303: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5305: Options Database Key:
5306: . -ts_adapt_wnormtype wnormtype - 2, INFINITY
5308: Level: developer
5310: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
5311: @*/
5312: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5313: {
5314: PetscInt norma_loc, norm_loc, normr_loc;
5316: PetscFunctionBegin;
5321: PetscAssertPointer(norm, 5);
5322: PetscAssertPointer(norma, 6);
5323: PetscAssertPointer(normr, 7);
5324: PetscCall(VecErrorWeightedNorms(U, Y, NULL, wnormtype, ts->atol, ts->vatol, ts->rtol, ts->vrtol, ts->adapt->ignore_max, norm, &norm_loc, norma, &norma_loc, normr, &normr_loc));
5325: if (wnormtype == NORM_2) {
5326: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5327: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5328: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5329: }
5330: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5331: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5332: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5333: PetscFunctionReturn(PETSC_SUCCESS);
5334: }
5336: /*@
5337: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5339: Collective
5341: Input Parameters:
5342: + ts - time stepping context
5343: . E - error vector
5344: . U - state vector, usually ts->vec_sol
5345: . Y - state vector, previous time step
5346: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5348: Output Parameters:
5349: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5350: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5351: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5353: Options Database Key:
5354: . -ts_adapt_wnormtype wnormtype - 2, INFINITY
5356: Level: developer
5358: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
5359: @*/
5360: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5361: {
5362: PetscInt norma_loc, norm_loc, normr_loc;
5364: PetscFunctionBegin;
5366: PetscCall(VecErrorWeightedNorms(U, Y, E, wnormtype, ts->atol, ts->vatol, ts->rtol, ts->vrtol, ts->adapt->ignore_max, norm, &norm_loc, norma, &norma_loc, normr, &normr_loc));
5367: if (wnormtype == NORM_2) {
5368: if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
5369: if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
5370: if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
5371: }
5372: PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
5373: PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
5374: PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
5375: PetscFunctionReturn(PETSC_SUCCESS);
5376: }
5378: /*@
5379: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5381: Logically Collective
5383: Input Parameters:
5384: + ts - time stepping context
5385: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5387: Note:
5388: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5390: Level: intermediate
5392: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5393: @*/
5394: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5395: {
5396: PetscFunctionBegin;
5398: ts->cfltime_local = cfltime;
5399: ts->cfltime = -1.;
5400: PetscFunctionReturn(PETSC_SUCCESS);
5401: }
5403: /*@
5404: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5406: Collective
5408: Input Parameter:
5409: . ts - time stepping context
5411: Output Parameter:
5412: . cfltime - maximum stable time step for forward Euler
5414: Level: advanced
5416: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
5417: @*/
5418: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5419: {
5420: PetscFunctionBegin;
5421: if (ts->cfltime < 0) PetscCallMPI(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
5422: *cfltime = ts->cfltime;
5423: PetscFunctionReturn(PETSC_SUCCESS);
5424: }
5426: /*@
5427: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5429: Input Parameters:
5430: + ts - the `TS` context.
5431: . xl - lower bound.
5432: - xu - upper bound.
5434: Level: advanced
5436: Note:
5437: If this routine is not called then the lower and upper bounds are set to
5438: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5440: .seealso: [](ch_ts), `TS`
5441: @*/
5442: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5443: {
5444: SNES snes;
5446: PetscFunctionBegin;
5447: PetscCall(TSGetSNES(ts, &snes));
5448: PetscCall(SNESVISetVariableBounds(snes, xl, xu));
5449: PetscFunctionReturn(PETSC_SUCCESS);
5450: }
5452: /*@
5453: TSComputeLinearStability - computes the linear stability function at a point
5455: Collective
5457: Input Parameters:
5458: + ts - the `TS` context
5459: . xr - real part of input argument
5460: - xi - imaginary part of input argument
5462: Output Parameters:
5463: + yr - real part of function value
5464: - yi - imaginary part of function value
5466: Level: developer
5468: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5469: @*/
5470: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5471: {
5472: PetscFunctionBegin;
5474: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5475: PetscFunctionReturn(PETSC_SUCCESS);
5476: }
5478: /*@
5479: TSRestartStep - Flags the solver to restart the next step
5481: Collective
5483: Input Parameter:
5484: . ts - the `TS` context obtained from `TSCreate()`
5486: Level: advanced
5488: Notes:
5489: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5490: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5491: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5492: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5493: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5494: discontinuous source terms).
5496: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5497: @*/
5498: PetscErrorCode TSRestartStep(TS ts)
5499: {
5500: PetscFunctionBegin;
5502: ts->steprestart = PETSC_TRUE;
5503: PetscFunctionReturn(PETSC_SUCCESS);
5504: }
5506: /*@
5507: TSRollBack - Rolls back one time step
5509: Collective
5511: Input Parameter:
5512: . ts - the `TS` context obtained from `TSCreate()`
5514: Level: advanced
5516: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5517: @*/
5518: PetscErrorCode TSRollBack(TS ts)
5519: {
5520: PetscFunctionBegin;
5522: PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
5523: PetscTryTypeMethod(ts, rollback);
5524: PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
5525: ts->time_step = ts->ptime - ts->ptime_prev;
5526: ts->ptime = ts->ptime_prev;
5527: ts->ptime_prev = ts->ptime_prev_rollback;
5528: ts->steps--;
5529: ts->steprollback = PETSC_TRUE;
5530: PetscFunctionReturn(PETSC_SUCCESS);
5531: }
5533: /*@
5534: TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step
5536: Not collective
5538: Input Parameter:
5539: . ts - the `TS` context obtained from `TSCreate()`
5541: Output Parameter:
5542: . flg - the rollback flag
5544: Level: advanced
5546: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
5547: @*/
5548: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
5549: {
5550: PetscFunctionBegin;
5552: PetscAssertPointer(flg, 2);
5553: *flg = ts->steprollback;
5554: PetscFunctionReturn(PETSC_SUCCESS);
5555: }
5557: /*@
5558: TSGetStepResize - Get the internal flag indicating if the current step is after a resize.
5560: Not collective
5562: Input Parameter:
5563: . ts - the `TS` context obtained from `TSCreate()`
5565: Output Parameter:
5566: . flg - the resize flag
5568: Level: advanced
5570: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetResize()`
5571: @*/
5572: PetscErrorCode TSGetStepResize(TS ts, PetscBool *flg)
5573: {
5574: PetscFunctionBegin;
5576: PetscAssertPointer(flg, 2);
5577: *flg = ts->stepresize;
5578: PetscFunctionReturn(PETSC_SUCCESS);
5579: }
5581: /*@
5582: TSGetStages - Get the number of stages and stage values
5584: Input Parameter:
5585: . ts - the `TS` context obtained from `TSCreate()`
5587: Output Parameters:
5588: + ns - the number of stages
5589: - Y - the current stage vectors
5591: Level: advanced
5593: Note:
5594: Both `ns` and `Y` can be `NULL`.
5596: .seealso: [](ch_ts), `TS`, `TSCreate()`
5597: @*/
5598: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5599: {
5600: PetscFunctionBegin;
5602: if (ns) PetscAssertPointer(ns, 2);
5603: if (Y) PetscAssertPointer(Y, 3);
5604: if (!ts->ops->getstages) {
5605: if (ns) *ns = 0;
5606: if (Y) *Y = NULL;
5607: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5608: PetscFunctionReturn(PETSC_SUCCESS);
5609: }
5611: /*@
5612: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5614: Collective
5616: Input Parameters:
5617: + ts - the `TS` context
5618: . t - current timestep
5619: . U - state vector
5620: . Udot - time derivative of state vector
5621: . shift - shift to apply, see note below
5622: - ctx - an optional application context
5624: Output Parameters:
5625: + J - Jacobian matrix (not altered in this routine)
5626: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
5628: Level: intermediate
5630: Notes:
5631: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5633: dF/dU + shift*dF/dUdot
5635: Most users should not need to explicitly call this routine, as it
5636: is used internally within the nonlinear solvers.
5638: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5639: routine, then it will try to get the coloring from the matrix. This requires that the
5640: matrix have nonzero entries precomputed.
5642: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5643: @*/
5644: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, PetscCtx ctx)
5645: {
5646: SNES snes;
5647: MatFDColoring color;
5648: PetscBool hascolor, matcolor = PETSC_FALSE;
5650: PetscFunctionBegin;
5651: PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
5652: PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
5653: if (!color) {
5654: DM dm;
5655: ISColoring iscoloring;
5657: PetscCall(TSGetDM(ts, &dm));
5658: PetscCall(DMHasColoring(dm, &hascolor));
5659: if (hascolor && !matcolor) {
5660: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
5661: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5662: PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5663: PetscCall(MatFDColoringSetFromOptions(color));
5664: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5665: PetscCall(ISColoringDestroy(&iscoloring));
5666: } else {
5667: MatColoring mc;
5669: PetscCall(MatColoringCreate(B, &mc));
5670: PetscCall(MatColoringSetDistance(mc, 2));
5671: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
5672: PetscCall(MatColoringSetFromOptions(mc));
5673: PetscCall(MatColoringApply(mc, &iscoloring));
5674: PetscCall(MatColoringDestroy(&mc));
5675: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
5676: PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
5677: PetscCall(MatFDColoringSetFromOptions(color));
5678: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
5679: PetscCall(ISColoringDestroy(&iscoloring));
5680: }
5681: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
5682: PetscCall(PetscObjectDereference((PetscObject)color));
5683: }
5684: PetscCall(TSGetSNES(ts, &snes));
5685: PetscCall(MatFDColoringApply(B, color, U, snes));
5686: if (J != B) {
5687: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
5688: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
5689: }
5690: PetscFunctionReturn(PETSC_SUCCESS);
5691: }
5693: /*@
5694: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5696: Logically collective
5698: Input Parameters:
5699: + ts - the `TS` context
5700: - func - function called within `TSFunctionDomainError()`
5702: Calling sequence of `func`:
5703: + ts - the `TS` context
5704: . time - the current time (of the stage)
5705: . state - the state to check if it is valid
5706: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable
5708: Level: intermediate
5710: Notes:
5711: `accept` must be collectively specified.
5712: If an implicit ODE solver is being used then, in addition to providing this routine, the
5713: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5714: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5715: Use `TSGetSNES()` to obtain the `SNES` object
5717: Developer Notes:
5718: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5719: since one takes a function pointer and the other does not.
5721: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5722: @*/
5723: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
5724: {
5725: PetscFunctionBegin;
5727: ts->functiondomainerror = func;
5728: PetscFunctionReturn(PETSC_SUCCESS);
5729: }
5731: /*@
5732: TSFunctionDomainError - Checks if the current state is valid
5734: Collective
5736: Input Parameters:
5737: + ts - the `TS` context
5738: . stagetime - time of the simulation
5739: - Y - state vector to check.
5741: Output Parameter:
5742: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5744: Level: developer
5746: Note:
5747: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5748: to check if the current state is valid.
5750: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
5751: @*/
5752: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5753: {
5754: PetscFunctionBegin;
5758: PetscAssertPointer(accept, 4);
5759: *accept = PETSC_TRUE;
5760: if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
5761: PetscFunctionReturn(PETSC_SUCCESS);
5762: }
5764: /*@
5765: TSClone - This function clones a time step `TS` object.
5767: Collective
5769: Input Parameter:
5770: . tsin - The input `TS`
5772: Output Parameter:
5773: . tsout - The output `TS` (cloned)
5775: Level: developer
5777: Notes:
5778: 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.
5779: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5781: 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
5782: .vb
5783: SNES snes_dup = NULL;
5784: TSGetSNES(ts,&snes_dup);
5785: TSSetSNES(ts,snes_dup);
5786: .ve
5788: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5789: @*/
5790: PetscErrorCode TSClone(TS tsin, TS *tsout)
5791: {
5792: TS t;
5793: SNES snes_start;
5794: DM dm;
5795: TSType type;
5797: PetscFunctionBegin;
5798: PetscAssertPointer(tsin, 1);
5799: *tsout = NULL;
5801: PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));
5803: /* General TS description */
5804: t->numbermonitors = 0;
5805: t->setupcalled = PETSC_FALSE;
5806: t->ksp_its = 0;
5807: t->snes_its = 0;
5808: t->nwork = 0;
5809: t->rhsjacobian.time = PETSC_MIN_REAL;
5810: t->rhsjacobian.scale = 1.;
5811: t->ijacobian.shift = 1.;
5813: PetscCall(TSGetSNES(tsin, &snes_start));
5814: PetscCall(TSSetSNES(t, snes_start));
5816: PetscCall(TSGetDM(tsin, &dm));
5817: PetscCall(TSSetDM(t, dm));
5819: t->adapt = tsin->adapt;
5820: PetscCall(PetscObjectReference((PetscObject)t->adapt));
5822: t->trajectory = tsin->trajectory;
5823: PetscCall(PetscObjectReference((PetscObject)t->trajectory));
5825: t->event = tsin->event;
5826: if (t->event) t->event->refct++;
5828: t->problem_type = tsin->problem_type;
5829: t->ptime = tsin->ptime;
5830: t->ptime_prev = tsin->ptime_prev;
5831: t->time_step = tsin->time_step;
5832: t->max_time = tsin->max_time;
5833: t->steps = tsin->steps;
5834: t->max_steps = tsin->max_steps;
5835: t->equation_type = tsin->equation_type;
5836: t->atol = tsin->atol;
5837: t->rtol = tsin->rtol;
5838: t->max_snes_failures = tsin->max_snes_failures;
5839: t->max_reject = tsin->max_reject;
5840: t->errorifstepfailed = tsin->errorifstepfailed;
5842: PetscCall(TSGetType(tsin, &type));
5843: PetscCall(TSSetType(t, type));
5845: t->vec_sol = NULL;
5847: t->cfltime = tsin->cfltime;
5848: t->cfltime_local = tsin->cfltime_local;
5849: t->exact_final_time = tsin->exact_final_time;
5851: t->ops[0] = tsin->ops[0];
5853: if (((PetscObject)tsin)->fortran_func_pointers) {
5854: PetscCall(PetscMalloc((10) * sizeof(PetscFortranCallbackFn *), &((PetscObject)t)->fortran_func_pointers));
5855: for (PetscInt i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5856: }
5857: *tsout = t;
5858: PetscFunctionReturn(PETSC_SUCCESS);
5859: }
5861: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(PetscCtx ctx, Vec x, Vec y)
5862: {
5863: TS ts = (TS)ctx;
5865: PetscFunctionBegin;
5866: PetscCall(TSComputeRHSFunction(ts, 0, x, y));
5867: PetscFunctionReturn(PETSC_SUCCESS);
5868: }
5870: /*@
5871: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5873: Logically Collective
5875: Input Parameter:
5876: . ts - the time stepping routine
5878: Output Parameter:
5879: . flg - `PETSC_TRUE` if the multiply is likely correct
5881: Options Database Key:
5882: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5884: Level: advanced
5886: Note:
5887: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5889: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5890: @*/
5891: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5892: {
5893: Mat J, B;
5894: TSRHSJacobianFn *func;
5895: void *ctx;
5897: PetscFunctionBegin;
5898: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5899: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5900: PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5901: PetscFunctionReturn(PETSC_SUCCESS);
5902: }
5904: /*@
5905: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5907: Logically Collective
5909: Input Parameter:
5910: . ts - the time stepping routine
5912: Output Parameter:
5913: . flg - `PETSC_TRUE` if the multiply is likely correct
5915: Options Database Key:
5916: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5918: Level: advanced
5920: Notes:
5921: This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian
5923: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5924: @*/
5925: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5926: {
5927: Mat J, B;
5928: void *ctx;
5929: TSRHSJacobianFn *func;
5931: PetscFunctionBegin;
5932: PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
5933: PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
5934: PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
5935: PetscFunctionReturn(PETSC_SUCCESS);
5936: }
5938: /*@
5939: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
5941: Logically Collective
5943: Input Parameters:
5944: + ts - timestepping context
5945: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5947: Options Database Key:
5948: . -ts_use_splitrhsfunction (true|false) - use the split RHS function for multirate solvers
5950: Level: intermediate
5952: Note:
5953: This is only for multirate methods
5955: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
5956: @*/
5957: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5958: {
5959: PetscFunctionBegin;
5961: ts->use_splitrhsfunction = use_splitrhsfunction;
5962: PetscFunctionReturn(PETSC_SUCCESS);
5963: }
5965: /*@
5966: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
5968: Not Collective
5970: Input Parameter:
5971: . ts - timestepping context
5973: Output Parameter:
5974: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5976: Level: intermediate
5978: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
5979: @*/
5980: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5981: {
5982: PetscFunctionBegin;
5984: *use_splitrhsfunction = ts->use_splitrhsfunction;
5985: PetscFunctionReturn(PETSC_SUCCESS);
5986: }
5988: /*@
5989: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
5991: Logically Collective
5993: Input Parameters:
5994: + ts - the time-stepper
5995: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
5997: Level: intermediate
5999: Note:
6000: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
6002: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
6003: @*/
6004: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
6005: {
6006: PetscFunctionBegin;
6008: ts->axpy_pattern = str;
6009: PetscFunctionReturn(PETSC_SUCCESS);
6010: }
6012: /*@
6013: TSSetEvaluationTimes - sets the evaluation points. The solution will be computed and stored for each time requested
6015: Collective
6017: Input Parameters:
6018: + ts - the time-stepper
6019: . n - number of the time points
6020: - time_points - array of the time points, must be increasing
6022: Options Database Key:
6023: . -ts_eval_times t0,...,tn - Sets the evaluation times
6025: Level: intermediate
6027: Notes:
6028: The elements in `time_points` must be all increasing. They correspond to the intermediate points to be saved.
6030: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
6032: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using evaluation times may
6033: pressure the memory system when using a large number of time points.
6035: .seealso: [](ch_ts), `TS`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`, `TSSetTimeSpan()`
6036: @*/
6037: PetscErrorCode TSSetEvaluationTimes(TS ts, PetscInt n, PetscReal time_points[])
6038: {
6039: PetscBool is_sorted;
6041: PetscFunctionBegin;
6043: if (ts->eval_times) { // Reset eval_times
6044: ts->eval_times->sol_idx = 0;
6045: ts->eval_times->time_point_idx = 0;
6046: if (n != ts->eval_times->num_time_points) {
6047: PetscCall(PetscFree(ts->eval_times->time_points));
6048: PetscCall(PetscFree(ts->eval_times->sol_times));
6049: PetscCall(VecDestroyVecs(ts->eval_times->num_time_points, &ts->eval_times->sol_vecs));
6050: } else {
6051: PetscCall(PetscArrayzero(ts->eval_times->sol_times, n));
6052: for (PetscInt i = 0; i < n; i++) PetscCall(VecZeroEntries(ts->eval_times->sol_vecs[i]));
6053: }
6054: } else { // Create/initialize eval_times
6055: TSEvaluationTimes eval_times;
6056: PetscCall(PetscNew(&eval_times));
6057: PetscCall(PetscMalloc1(n, &eval_times->time_points));
6058: PetscCall(PetscMalloc1(n, &eval_times->sol_times));
6059: eval_times->reltol = 1e-6;
6060: eval_times->abstol = 10 * PETSC_MACHINE_EPSILON;
6061: eval_times->worktol = 0;
6062: ts->eval_times = eval_times;
6063: }
6064: ts->eval_times->num_time_points = n;
6065: PetscCall(PetscSortedReal(n, time_points, &is_sorted));
6066: PetscCheck(is_sorted, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "time_points array must be sorted");
6067: PetscCall(PetscArraycpy(ts->eval_times->time_points, time_points, n));
6068: // Note: ts->vec_sol not guaranteed to exist, so ts->eval_times->sol_vecs allocated at TSSolve time
6069: PetscFunctionReturn(PETSC_SUCCESS);
6070: }
6072: /*@
6073: TSGetEvaluationTimes - gets the evaluation times set with `TSSetEvaluationTimes()`
6075: Not Collective
6077: Input Parameter:
6078: . ts - the time-stepper
6080: Output Parameters:
6081: + n - number of the time points
6082: - time_points - array of the time points
6084: Level: beginner
6086: Note:
6087: The values obtained are valid until the `TS` object is destroyed.
6089: Both `n` and `time_points` can be `NULL`.
6091: Also used to see time points set by `TSSetTimeSpan()`.
6093: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6094: @*/
6095: PetscErrorCode TSGetEvaluationTimes(TS ts, PetscInt *n, const PetscReal *time_points[])
6096: {
6097: PetscFunctionBegin;
6099: if (n) PetscAssertPointer(n, 2);
6100: if (time_points) PetscAssertPointer(time_points, 3);
6101: if (!ts->eval_times) {
6102: if (n) *n = 0;
6103: if (time_points) *time_points = NULL;
6104: } else {
6105: if (n) *n = ts->eval_times->num_time_points;
6106: if (time_points) *time_points = ts->eval_times->time_points;
6107: }
6108: PetscFunctionReturn(PETSC_SUCCESS);
6109: }
6111: /*@
6112: TSGetEvaluationSolutions - Get the number of solutions and the solutions at the evaluation time points specified
6114: Input Parameter:
6115: . ts - the `TS` context obtained from `TSCreate()`
6117: Output Parameters:
6118: + nsol - the number of solutions
6119: . sol_times - array of solution times corresponding to the solution vectors. See note below
6120: - Sols - the solution vectors
6122: Level: intermediate
6124: Notes:
6125: Both `nsol` and `Sols` can be `NULL`.
6127: Some time points in the evaluation points may be skipped by `TS` so that `nsol` is less than the number of points specified by `TSSetEvaluationTimes()`.
6128: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain evaluation times.
6130: Also used to see view solutions requested by `TSSetTimeSpan()`.
6132: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`
6133: @*/
6134: PetscErrorCode TSGetEvaluationSolutions(TS ts, PetscInt *nsol, const PetscReal *sol_times[], Vec *Sols[])
6135: {
6136: PetscFunctionBegin;
6138: if (nsol) PetscAssertPointer(nsol, 2);
6139: if (sol_times) PetscAssertPointer(sol_times, 3);
6140: if (Sols) PetscAssertPointer(Sols, 4);
6141: if (!ts->eval_times) {
6142: if (nsol) *nsol = 0;
6143: if (sol_times) *sol_times = NULL;
6144: if (Sols) *Sols = NULL;
6145: } else {
6146: if (nsol) *nsol = ts->eval_times->sol_idx;
6147: if (sol_times) *sol_times = ts->eval_times->sol_times;
6148: if (Sols) *Sols = ts->eval_times->sol_vecs;
6149: }
6150: PetscFunctionReturn(PETSC_SUCCESS);
6151: }
6153: /*@
6154: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
6156: Collective
6158: Input Parameters:
6159: + ts - the time-stepper
6160: . n - number of the time points (>=2)
6161: - span_times - array of the time points, must be increasing. The first element and the last element are the initial time and the final time respectively.
6163: Options Database Key:
6164: . -ts_time_span t0,...,tf - Sets the time span
6166: Level: intermediate
6168: Notes:
6169: This function is identical to `TSSetEvaluationTimes()`, except that it also sets the initial time and final time for the `ts` to the first and last `span_times` entries.
6171: The elements in `span_times` must be all increasing. They correspond to the intermediate points to be saved.
6173: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
6175: The intermediate solutions are saved in a vector array that can be accessed with `TSGetEvaluationSolutions()`. Thus using time span may
6176: pressure the memory system when using a large number of span points.
6178: .seealso: [](ch_ts), `TS`, `TSSetEvaluationTimes()`, `TSGetEvaluationTimes()`, `TSGetEvaluationSolutions()`
6179: @*/
6180: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal span_times[])
6181: {
6182: PetscFunctionBegin;
6184: PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
6185: PetscCall(TSSetEvaluationTimes(ts, n, span_times));
6186: PetscCall(TSSetTime(ts, span_times[0]));
6187: PetscCall(TSSetMaxTime(ts, span_times[n - 1]));
6188: PetscFunctionReturn(PETSC_SUCCESS);
6189: }
6191: /*@
6192: TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
6194: Collective
6196: Input Parameters:
6197: + ts - the `TS` context
6198: . J - Jacobian matrix (not altered in this routine)
6199: - B - newly computed Jacobian matrix to use with preconditioner
6201: Level: intermediate
6203: Notes:
6204: This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
6205: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
6206: and multiple fields are involved.
6208: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
6209: structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
6210: usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
6211: `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.
6213: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
6214: @*/
6215: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
6216: {
6217: MatColoring mc = NULL;
6218: ISColoring iscoloring = NULL;
6219: MatFDColoring matfdcoloring = NULL;
6221: PetscFunctionBegin;
6222: /* Generate new coloring after eliminating zeros in the matrix */
6223: PetscCall(MatEliminateZeros(B, PETSC_TRUE));
6224: PetscCall(MatColoringCreate(B, &mc));
6225: PetscCall(MatColoringSetDistance(mc, 2));
6226: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
6227: PetscCall(MatColoringSetFromOptions(mc));
6228: PetscCall(MatColoringApply(mc, &iscoloring));
6229: PetscCall(MatColoringDestroy(&mc));
6230: /* Replace the old coloring with the new one */
6231: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
6232: PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESTSFormFunction, (void *)ts));
6233: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
6234: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
6235: PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
6236: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
6237: PetscCall(ISColoringDestroy(&iscoloring));
6238: PetscFunctionReturn(PETSC_SUCCESS);
6239: }