Actual source code: ts.c
1: #include <petsc/private/tsimpl.h>
2: #include <petscdmshell.h>
3: #include <petscdmda.h>
4: #include <petscviewer.h>
5: #include <petscdraw.h>
6: #include <petscconvest.h>
8: #define SkipSmallValue(a, b, tol) \
9: if (PetscAbsScalar(a) < tol || PetscAbsScalar(b) < tol) continue;
11: /* Logging support */
12: PetscClassId TS_CLASSID, DMTS_CLASSID;
13: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;
15: const char *const TSExactFinalTimeOptions[] = {"UNSPECIFIED", "STEPOVER", "INTERPOLATE", "MATCHSTEP", "TSExactFinalTimeOption", "TS_EXACTFINALTIME_", NULL};
17: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt, TSAdaptType default_type)
18: {
21: if (!((PetscObject)adapt)->type_name) TSAdaptSetType(adapt, default_type);
22: return 0;
23: }
25: /*@
26: TSSetFromOptions - Sets various `TS` parameters from user options.
28: Collective
30: Input Parameter:
31: . ts - the `TS` context obtained from `TSCreate()`
33: Options Database Keys:
34: + -ts_type <type> - EULER, BEULER, SUNDIALS, PSEUDO, CN, RK, THETA, ALPHA, GLLE, SSP, GLEE, BSYMP, IRK
35: . -ts_save_trajectory - checkpoint the solution at each time-step
36: . -ts_max_time <time> - maximum time to compute to
37: . -ts_time_span <t0,...tf> - sets the time span, solutions are computed and stored for each indicated time
38: . -ts_max_steps <steps> - maximum number of time-steps to take
39: . -ts_init_time <time> - initial time to start computation
40: . -ts_final_time <time> - final time to compute to (deprecated: use -ts_max_time)
41: . -ts_dt <dt> - initial time step
42: . -ts_exact_final_time <stepover,interpolate,matchstep> - whether to stop at the exact given final time and how to compute the solution at that time
43: . -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed
44: . -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails
45: . -ts_error_if_step_fails <true,false> - Error if no step succeeds
46: . -ts_rtol <rtol> - relative tolerance for local truncation error
47: . -ts_atol <atol> - Absolute tolerance for local truncation error
48: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function
49: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - test the Jacobian at each iteration against finite difference with RHS function
50: . -ts_adjoint_solve <yes,no> - After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory)
51: . -ts_fd_color - Use finite differences with coloring to compute IJacobian
52: . -ts_monitor - print information at each timestep
53: . -ts_monitor_cancel - Cancel all monitors
54: . -ts_monitor_lg_solution - Monitor solution graphically
55: . -ts_monitor_lg_error - Monitor error graphically
56: . -ts_monitor_error - Monitors norm of error
57: . -ts_monitor_lg_timestep - Monitor timestep size graphically
58: . -ts_monitor_lg_timestep_log - Monitor log timestep size graphically
59: . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically
60: . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically
61: . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically
62: . -ts_monitor_draw_solution - Monitor solution graphically
63: . -ts_monitor_draw_solution_phase <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
64: . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
65: . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
66: . -ts_monitor_solution_vtk <filename.vts,filename.vtu> - Save each time step to a binary file, use filename-%%03" PetscInt_FMT ".vts (filename-%%03" PetscInt_FMT ".vtu)
67: - -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time
69: Level: beginner
71: Notes:
72: See `SNESSetFromOptions()` and `KSPSetFromOptions()` for how to control the nonlinear and linear solves used by the time-stepper.
74: Certain `SNES` options get reset for each new nonlinear solver, for example -snes_lag_jacobian <its> and -snes_lag_preconditioner <its>, in order
75: to retain them over the multiple nonlinear solves that TS uses you mush also provide -snes_lag_jacobian_persists true and
76: -snes_lag_preconditioner_persists true
78: Developer Note:
79: We should unify all the -ts_monitor options in the way that -xxx_view has been unified
81: .seealso: [](chapter_ts), `TS`, `TSGetType()`
82: @*/
83: PetscErrorCode TSSetFromOptions(TS ts)
84: {
85: PetscBool opt, flg, tflg;
86: char monfilename[PETSC_MAX_PATH_LEN];
87: PetscReal time_step, tspan[100];
88: PetscInt nt = PETSC_STATIC_ARRAY_LENGTH(tspan);
89: TSExactFinalTimeOption eftopt;
90: char dir[16];
91: TSIFunction ifun;
92: const char *defaultType;
93: char typeName[256];
97: TSRegisterAll();
98: TSGetIFunction(ts, NULL, &ifun, NULL);
100: PetscObjectOptionsBegin((PetscObject)ts);
101: if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
102: else defaultType = ifun ? TSBEULER : TSEULER;
103: PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, 256, &opt);
104: if (opt) TSSetType(ts, typeName);
105: else TSSetType(ts, defaultType);
107: /* Handle generic TS options */
108: PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL);
109: PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL);
110: PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", tspan, &nt, &flg);
111: if (flg) TSSetTimeSpan(ts, nt, tspan);
112: PetscOptionsInt("-ts_max_steps", "Maximum number of time steps", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL);
113: PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL);
114: PetscOptionsReal("-ts_dt", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg);
115: if (flg) TSSetTimeStep(ts, time_step);
116: PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg);
117: if (flg) TSSetExactFinalTime(ts, eftopt);
118: PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, NULL);
119: PetscOptionsInt("-ts_max_reject", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, NULL);
120: PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL);
121: PetscOptionsReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL);
122: PetscOptionsReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL);
124: PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL);
125: 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);
126: PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL);
127: #if defined(PETSC_HAVE_SAWS)
128: {
129: PetscBool set;
130: flg = PETSC_FALSE;
131: PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set);
132: if (set) PetscObjectSAWsSetBlock((PetscObject)ts, flg);
133: }
134: #endif
136: /* Monitor options */
137: PetscOptionsInt("-ts_monitor_frequency", "Number of time steps between monitor output", "TSMonitorSetFrequency", ts->monitorFrequency, &ts->monitorFrequency, NULL);
138: TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL);
139: TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL);
140: TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, NULL);
141: TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL);
143: PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg);
144: if (flg) PetscPythonMonitorSet((PetscObject)ts, monfilename);
146: PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt);
147: if (opt) {
148: PetscInt howoften = 1;
149: DM dm;
150: PetscBool net;
152: PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL);
153: TSGetDM(ts, &dm);
154: PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net);
155: if (net) {
156: TSMonitorLGCtxNetwork ctx;
157: TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx);
158: TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxNetworkDestroy);
159: PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL);
160: } else {
161: TSMonitorLGCtx ctx;
162: TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
163: TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
164: }
165: }
167: PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt);
168: if (opt) {
169: TSMonitorLGCtx ctx;
170: PetscInt howoften = 1;
172: PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL);
173: TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
174: TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
175: }
176: TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL);
178: PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt);
179: if (opt) {
180: TSMonitorLGCtx ctx;
181: PetscInt howoften = 1;
183: PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL);
184: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
185: TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
186: }
187: PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt);
188: if (opt) {
189: TSMonitorLGCtx ctx;
190: PetscInt howoften = 1;
192: PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL);
193: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
194: TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
195: ctx->semilogy = PETSC_TRUE;
196: }
198: PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt);
199: if (opt) {
200: TSMonitorLGCtx ctx;
201: PetscInt howoften = 1;
203: PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL);
204: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
205: TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
206: }
207: PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt);
208: if (opt) {
209: TSMonitorLGCtx ctx;
210: PetscInt howoften = 1;
212: PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL);
213: TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx);
214: TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy);
215: }
216: PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt);
217: if (opt) {
218: TSMonitorSPEigCtx ctx;
219: PetscInt howoften = 1;
221: PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL);
222: TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
223: TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscErrorCode(*)(void **))TSMonitorSPEigCtxDestroy);
224: }
225: PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase from the DMSwarm", "TSMonitorSPSwarm", &opt);
226: if (opt) {
227: TSMonitorSPCtx ctx;
228: PetscInt howoften = 1, retain = 0;
229: PetscBool phase = PETSC_TRUE, create = PETSC_TRUE;
231: for (PetscInt i = 0; i < ts->numbermonitors; ++i)
232: if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
233: create = PETSC_FALSE;
234: break;
235: }
236: if (create) {
237: PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL);
238: PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL);
239: PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL);
240: TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, &ctx);
241: TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorSPCtxDestroy);
242: }
243: }
244: opt = PETSC_FALSE;
245: PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt);
246: if (opt) {
247: TSMonitorDrawCtx ctx;
248: PetscInt howoften = 1;
250: PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL);
251: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
252: TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
253: }
254: opt = PETSC_FALSE;
255: PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt);
256: if (opt) {
257: TSMonitorDrawCtx ctx;
258: PetscReal bounds[4];
259: PetscInt n = 4;
260: PetscDraw draw;
261: PetscDrawAxis axis;
263: PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL);
265: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx);
266: PetscViewerDrawGetDraw(ctx->viewer, 0, &draw);
267: PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis);
268: PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]);
269: PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2");
270: TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
271: }
272: opt = PETSC_FALSE;
273: PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt);
274: if (opt) {
275: TSMonitorDrawCtx ctx;
276: PetscInt howoften = 1;
278: PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL);
279: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
280: TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
281: }
282: opt = PETSC_FALSE;
283: PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt);
284: if (opt) {
285: TSMonitorDrawCtx ctx;
286: PetscInt howoften = 1;
288: PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL);
289: TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx);
290: TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy);
291: }
293: opt = PETSC_FALSE;
294: 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);
295: if (flg) {
296: const char *ptr, *ptr2;
297: char *filetemplate;
299: /* Do some cursory validation of the input. */
300: PetscStrstr(monfilename, "%", (char **)&ptr);
302: for (ptr++; ptr && *ptr; ptr++) {
303: PetscStrchr("DdiouxX", *ptr, (char **)&ptr2);
305: if (ptr2) break;
306: }
307: PetscStrallocpy(monfilename, &filetemplate);
308: TSMonitorSet(ts, TSMonitorSolutionVTK, filetemplate, (PetscErrorCode(*)(void **))TSMonitorSolutionVTKDestroy);
309: }
311: PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg);
312: if (flg) {
313: TSMonitorDMDARayCtx *rayctx;
314: int ray = 0;
315: DMDirection ddir;
316: DM da;
317: PetscMPIInt rank;
320: if (dir[0] == 'x') ddir = DM_X;
321: else if (dir[0] == 'y') ddir = DM_Y;
322: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
323: sscanf(dir + 2, "%d", &ray);
325: PetscInfo(((PetscObject)ts), "Displaying DMDA ray %c = %d\n", dir[0], ray);
326: PetscNew(&rayctx);
327: TSGetDM(ts, &da);
328: DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);
329: MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank);
330: if (rank == 0) PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer);
331: rayctx->lgctx = NULL;
332: TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy);
333: }
334: PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg);
335: if (flg) {
336: TSMonitorDMDARayCtx *rayctx;
337: int ray = 0;
338: DMDirection ddir;
339: DM da;
340: PetscInt howoften = 1;
343: if (dir[0] == 'x') ddir = DM_X;
344: else if (dir[0] == 'y') ddir = DM_Y;
345: else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
346: sscanf(dir + 2, "%d", &ray);
348: PetscInfo(((PetscObject)ts), "Displaying LG DMDA ray %c = %d\n", dir[0], ray);
349: PetscNew(&rayctx);
350: TSGetDM(ts, &da);
351: DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);
352: TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx);
353: TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy);
354: }
356: PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt);
357: if (opt) {
358: TSMonitorEnvelopeCtx ctx;
360: TSMonitorEnvelopeCtxCreate(ts, &ctx);
361: TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscErrorCode(*)(void **))TSMonitorEnvelopeCtxDestroy);
362: }
363: flg = PETSC_FALSE;
364: PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt);
365: if (opt && flg) TSMonitorCancel(ts);
367: flg = PETSC_FALSE;
368: PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeJacobianDefaultColor", flg, &flg, NULL);
369: if (flg) {
370: DM dm;
372: TSGetDM(ts, &dm);
373: DMTSUnsetIJacobianContext_Internal(dm);
374: TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL);
375: PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n");
376: }
378: /* Handle specific TS options */
379: PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);
381: /* Handle TSAdapt options */
382: TSGetAdapt(ts, &ts->adapt);
383: TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type);
384: TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject);
386: /* TS trajectory must be set after TS, since it may use some TS options above */
387: tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
388: PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL);
389: if (tflg) TSSetSaveTrajectory(ts);
391: TSAdjointSetFromOptions(ts, PetscOptionsObject);
393: /* process any options handlers added with PetscObjectAddOptionsHandler() */
394: PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject);
395: PetscOptionsEnd();
397: if (ts->trajectory) TSTrajectorySetFromOptions(ts->trajectory, ts);
399: /* why do we have to do this here and not during TSSetUp? */
400: TSGetSNES(ts, &ts->snes);
401: if (ts->problem_type == TS_LINEAR) {
402: PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, "");
403: if (!flg) SNESSetType(ts->snes, SNESKSPONLY);
404: }
405: SNESSetFromOptions(ts->snes);
406: return 0;
407: }
409: /*@
410: TSGetTrajectory - Gets the trajectory from a `TS` if it exists
412: Collective
414: Input Parameters:
415: . ts - the `TS` context obtained from `TSCreate()`
417: Output Parameters:
418: . tr - the `TSTrajectory` object, if it exists
420: Level: advanced
422: Note:
423: This routine should be called after all `TS` options have been set
425: .seealso: [](chapter_ts), `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectory`, `TSTrajectoryCreate()`
426: @*/
427: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
428: {
430: *tr = ts->trajectory;
431: return 0;
432: }
434: /*@
435: TSSetSaveTrajectory - Causes the `TS` to save its solutions as it iterates forward in time in a `TSTrajectory` object
437: Collective
439: Input Parameter:
440: . ts - the `TS` context obtained from `TSCreate()`
442: Options Database Keys:
443: + -ts_save_trajectory - saves the trajectory to a file
444: - -ts_trajectory_type type - set trajectory type
446: Level: intermediate
448: Notes:
449: This routine should be called after all `TS` options have been set
451: The `TSTRAJECTORYVISUALIZATION` files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
452: MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m
454: .seealso: [](chapter_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`
455: @*/
456: PetscErrorCode TSSetSaveTrajectory(TS ts)
457: {
459: if (!ts->trajectory) TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory);
460: return 0;
461: }
463: /*@
464: TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object
466: Collective
468: Input Parameters:
469: . ts - the `TS` context obtained from `TSCreate()`
471: Level: intermediate
473: .seealso: [](chapter_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
474: @*/
475: PetscErrorCode TSResetTrajectory(TS ts)
476: {
478: if (ts->trajectory) {
479: TSTrajectoryDestroy(&ts->trajectory);
480: TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory);
481: }
482: return 0;
483: }
485: /*@
486: TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from TS
488: Collective
490: Input Parameters:
491: . ts - the `TS` context obtained from `TSCreate()`
493: Level: intermediate
495: .seealso: [](chapter_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
496: @*/
497: PetscErrorCode TSRemoveTrajectory(TS ts)
498: {
500: if (ts->trajectory) TSTrajectoryDestroy(&ts->trajectory);
501: return 0;
502: }
504: /*@
505: TSComputeRHSJacobian - Computes the Jacobian matrix that has been
506: set with `TSSetRHSJacobian()`.
508: Collective
510: Input Parameters:
511: + ts - the `TS` context
512: . t - current timestep
513: - U - input vector
515: Output Parameters:
516: + A - Jacobian matrix
517: - B - optional preconditioning matrix
519: Level: developer
521: Note:
522: Most users should not need to explicitly call this routine, as it
523: is used internally within the nonlinear solvers.
525: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()`
526: @*/
527: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
528: {
529: PetscObjectState Ustate;
530: PetscObjectId Uid;
531: DM dm;
532: DMTS tsdm;
533: TSRHSJacobian rhsjacobianfunc;
534: void *ctx;
535: TSRHSFunction rhsfunction;
540: TSGetDM(ts, &dm);
541: DMGetDMTS(dm, &tsdm);
542: DMTSGetRHSFunction(dm, &rhsfunction, NULL);
543: DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx);
544: PetscObjectStateGet((PetscObject)U, &Ustate);
545: PetscObjectGetId((PetscObject)U, &Uid);
547: if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) return 0;
550: if (rhsjacobianfunc) {
551: PetscLogEventBegin(TS_JacobianEval, ts, U, A, B);
552: PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
553: ts->rhsjacs++;
554: PetscLogEventEnd(TS_JacobianEval, ts, U, A, B);
555: } else {
556: MatZeroEntries(A);
557: if (B && A != B) MatZeroEntries(B);
558: }
559: ts->rhsjacobian.time = t;
560: ts->rhsjacobian.shift = 0;
561: ts->rhsjacobian.scale = 1.;
562: PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid);
563: PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate);
564: return 0;
565: }
567: /*@
568: TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS`
570: Collective
572: Input Parameters:
573: + ts - the `TS` context
574: . t - current time
575: - U - state vector
577: Output Parameter:
578: . y - right hand side
580: Level: developer
582: Note:
583: Most users should not need to explicitly call this routine, as it
584: is used internally within the nonlinear solvers.
586: .seealso: [](chapter_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
587: @*/
588: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
589: {
590: TSRHSFunction rhsfunction;
591: TSIFunction ifunction;
592: void *ctx;
593: DM dm;
598: TSGetDM(ts, &dm);
599: DMTSGetRHSFunction(dm, &rhsfunction, &ctx);
600: DMTSGetIFunction(dm, &ifunction, NULL);
604: if (rhsfunction) {
605: PetscLogEventBegin(TS_FunctionEval, ts, U, y, 0);
606: VecLockReadPush(U);
607: PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
608: VecLockReadPop(U);
609: ts->rhsfuncs++;
610: PetscLogEventEnd(TS_FunctionEval, ts, U, y, 0);
611: } else VecZeroEntries(y);
612: return 0;
613: }
615: /*@
616: TSComputeSolutionFunction - Evaluates the solution function.
618: Collective
620: Input Parameters:
621: + ts - the `TS` context
622: - t - current time
624: Output Parameter:
625: . U - the solution
627: Level: developer
629: .seealso: [](chapter_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
630: @*/
631: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
632: {
633: TSSolutionFunction solutionfunction;
634: void *ctx;
635: DM dm;
639: TSGetDM(ts, &dm);
640: DMTSGetSolutionFunction(dm, &solutionfunction, &ctx);
641: if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
642: return 0;
643: }
644: /*@
645: TSComputeForcingFunction - Evaluates the forcing function.
647: Collective
649: Input Parameters:
650: + ts - the `TS` context
651: - t - current time
653: Output Parameter:
654: . U - the function value
656: Level: developer
658: .seealso: [](chapter_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
659: @*/
660: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
661: {
662: void *ctx;
663: DM dm;
664: TSForcingFunction forcing;
668: TSGetDM(ts, &dm);
669: DMTSGetForcingFunction(dm, &forcing, &ctx);
671: if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
672: return 0;
673: }
675: static PetscErrorCode TSGetRHSVec_Private(TS ts, Vec *Frhs)
676: {
677: Vec F;
679: *Frhs = NULL;
680: TSGetIFunction(ts, &F, NULL, NULL);
681: if (!ts->Frhs) VecDuplicate(F, &ts->Frhs);
682: *Frhs = ts->Frhs;
683: return 0;
684: }
686: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
687: {
688: Mat A, B;
689: TSIJacobian ijacobian;
691: if (Arhs) *Arhs = NULL;
692: if (Brhs) *Brhs = NULL;
693: TSGetIJacobian(ts, &A, &B, &ijacobian, NULL);
694: if (Arhs) {
695: if (!ts->Arhs) {
696: if (ijacobian) {
697: MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs);
698: TSSetMatStructure(ts, SAME_NONZERO_PATTERN);
699: } else {
700: ts->Arhs = A;
701: PetscObjectReference((PetscObject)A);
702: }
703: } else {
704: PetscBool flg;
705: SNESGetUseMatrixFree(ts->snes, NULL, &flg);
706: /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
707: if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
708: PetscObjectDereference((PetscObject)ts->Arhs);
709: ts->Arhs = A;
710: PetscObjectReference((PetscObject)A);
711: }
712: }
713: *Arhs = ts->Arhs;
714: }
715: if (Brhs) {
716: if (!ts->Brhs) {
717: if (A != B) {
718: if (ijacobian) {
719: MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs);
720: } else {
721: ts->Brhs = B;
722: PetscObjectReference((PetscObject)B);
723: }
724: } else {
725: PetscObjectReference((PetscObject)ts->Arhs);
726: ts->Brhs = ts->Arhs;
727: }
728: }
729: *Brhs = ts->Brhs;
730: }
731: return 0;
732: }
734: /*@
735: TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,U,Udot)=0
737: Collective
739: Input Parameters:
740: + ts - the `TS` context
741: . t - current time
742: . U - state vector
743: . Udot - time derivative of state vector
744: - imex - flag indicates if the method is `TSIMEX` so that the RHSFunction should be kept separate
746: Output Parameter:
747: . Y - right hand side
749: Level: developer
751: Note:
752: Most users should not need to explicitly call this routine, as it
753: is used internally within the nonlinear solvers.
755: If the user did did not write their equations in implicit form, this
756: function recasts them in implicit form.
758: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
759: @*/
760: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
761: {
762: TSIFunction ifunction;
763: TSRHSFunction rhsfunction;
764: void *ctx;
765: DM dm;
772: TSGetDM(ts, &dm);
773: DMTSGetIFunction(dm, &ifunction, &ctx);
774: DMTSGetRHSFunction(dm, &rhsfunction, NULL);
778: PetscLogEventBegin(TS_FunctionEval, ts, U, Udot, Y);
779: if (ifunction) {
780: PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
781: ts->ifuncs++;
782: }
783: if (imex) {
784: if (!ifunction) VecCopy(Udot, Y);
785: } else if (rhsfunction) {
786: if (ifunction) {
787: Vec Frhs;
788: TSGetRHSVec_Private(ts, &Frhs);
789: TSComputeRHSFunction(ts, t, U, Frhs);
790: VecAXPY(Y, -1, Frhs);
791: } else {
792: TSComputeRHSFunction(ts, t, U, Y);
793: VecAYPX(Y, -1, Udot);
794: }
795: }
796: PetscLogEventEnd(TS_FunctionEval, ts, U, Udot, Y);
797: return 0;
798: }
800: /*
801: TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call TSComputeRHSJacobian() on it.
803: Note:
804: This routine is needed when one switches from TSComputeIJacobian() to TSComputeRHSJacobian() because the Jacobian matrix may be shifted or scaled in TSComputeIJacobian().
806: */
807: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
808: {
813: if (ts->rhsjacobian.shift) MatShift(A, -ts->rhsjacobian.shift);
814: if (ts->rhsjacobian.scale == -1.) MatScale(A, -1);
815: if (B && B == ts->Brhs && A != B) {
816: if (ts->rhsjacobian.shift) MatShift(B, -ts->rhsjacobian.shift);
817: if (ts->rhsjacobian.scale == -1.) MatScale(B, -1);
818: }
819: ts->rhsjacobian.shift = 0;
820: ts->rhsjacobian.scale = 1.;
821: return 0;
822: }
824: /*@
825: TSComputeIJacobian - Evaluates the Jacobian of the DAE
827: Collective
829: Input
830: Input Parameters:
831: + ts - the `TS` context
832: . t - current timestep
833: . U - state vector
834: . Udot - time derivative of state vector
835: . shift - shift to apply, see note below
836: - imex - flag indicates if the method is `TSIMEX` so that the RHSJacobian should be kept separate
838: Output Parameters:
839: + A - Jacobian matrix
840: - B - matrix from which the preconditioner is constructed; often the same as A
842: Level: developer
844: Notes:
845: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
847: dF/dU + shift*dF/dUdot
849: Most users should not need to explicitly call this routine, as it
850: is used internally within the nonlinear solvers.
852: .seealso: [](chapter_ts), `TS`, `TSSetIJacobian()`
853: @*/
854: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
855: {
856: TSIJacobian ijacobian;
857: TSRHSJacobian rhsjacobian;
858: DM dm;
859: void *ctx;
869: TSGetDM(ts, &dm);
870: DMTSGetIJacobian(dm, &ijacobian, &ctx);
871: DMTSGetRHSJacobian(dm, &rhsjacobian, NULL);
875: PetscLogEventBegin(TS_JacobianEval, ts, U, A, B);
876: if (ijacobian) {
877: PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
878: ts->ijacs++;
879: }
880: if (imex) {
881: if (!ijacobian) { /* system was written as Udot = G(t,U) */
882: PetscBool assembled;
883: if (rhsjacobian) {
884: Mat Arhs = NULL;
885: TSGetRHSMats_Private(ts, &Arhs, NULL);
886: if (A == Arhs) {
888: ts->rhsjacobian.time = PETSC_MIN_REAL;
889: }
890: }
891: MatZeroEntries(A);
892: MatAssembled(A, &assembled);
893: if (!assembled) {
894: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
895: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
896: }
897: MatShift(A, shift);
898: if (A != B) {
899: MatZeroEntries(B);
900: MatAssembled(B, &assembled);
901: if (!assembled) {
902: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
903: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
904: }
905: MatShift(B, shift);
906: }
907: }
908: } else {
909: Mat Arhs = NULL, Brhs = NULL;
911: /* RHSJacobian needs to be converted to part of IJacobian if exists */
912: if (rhsjacobian) TSGetRHSMats_Private(ts, &Arhs, &Brhs);
913: if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
914: PetscObjectState Ustate;
915: PetscObjectId Uid;
916: TSRHSFunction rhsfunction;
918: DMTSGetRHSFunction(dm, &rhsfunction, NULL);
919: PetscObjectStateGet((PetscObject)U, &Ustate);
920: PetscObjectGetId((PetscObject)U, &Uid);
921: if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
922: ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */
923: MatShift(A, shift - ts->rhsjacobian.shift); /* revert the old shift and add the new shift with a single call to MatShift */
924: if (A != B) MatShift(B, shift - ts->rhsjacobian.shift);
925: } else {
926: PetscBool flg;
928: if (ts->rhsjacobian.reuse) { /* Undo the damage */
929: /* MatScale has a short path for this case.
930: However, this code path is taken the first time TSComputeRHSJacobian is called
931: and the matrices have not been assembled yet */
932: TSRecoverRHSJacobian(ts, A, B);
933: }
934: TSComputeRHSJacobian(ts, t, U, A, B);
935: SNESGetUseMatrixFree(ts->snes, NULL, &flg);
936: /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
937: if (!flg) {
938: MatScale(A, -1);
939: MatShift(A, shift);
940: }
941: if (A != B) {
942: MatScale(B, -1);
943: MatShift(B, shift);
944: }
945: }
946: ts->rhsjacobian.scale = -1;
947: ts->rhsjacobian.shift = shift;
948: } else if (Arhs) { /* Both IJacobian and RHSJacobian */
949: if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
950: MatZeroEntries(A);
951: MatShift(A, shift);
952: if (A != B) {
953: MatZeroEntries(B);
954: MatShift(B, shift);
955: }
956: }
957: TSComputeRHSJacobian(ts, t, U, Arhs, Brhs);
958: MatAXPY(A, -1, Arhs, ts->axpy_pattern);
959: if (A != B) MatAXPY(B, -1, Brhs, ts->axpy_pattern);
960: }
961: }
962: PetscLogEventEnd(TS_JacobianEval, ts, U, A, B);
963: return 0;
964: }
966: /*@C
967: TSSetRHSFunction - Sets the routine for evaluating the function,
968: where U_t = G(t,u).
970: Logically Collective
972: Input Parameters:
973: + ts - the `TS` context obtained from `TSCreate()`
974: . r - vector to put the computed right hand side (or NULL to have it created)
975: . f - routine for evaluating the right-hand-side function
976: - ctx - [optional] user-defined context for private data for the
977: function evaluation routine (may be NULL)
979: Calling sequence of f:
980: $ PetscErrorCode f(TS ts,PetscReal t,Vec u,Vec F,void *ctx);
982: + ts - timestep context
983: . t - current timestep
984: . u - input vector
985: . F - function vector
986: - ctx - [optional] user-defined function context
988: Level: beginner
990: Note:
991: You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.
993: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
994: @*/
995: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, PetscErrorCode (*f)(TS, PetscReal, Vec, Vec, void *), void *ctx)
996: {
997: SNES snes;
998: Vec ralloc = NULL;
999: DM dm;
1004: TSGetDM(ts, &dm);
1005: DMTSSetRHSFunction(dm, f, ctx);
1006: TSGetSNES(ts, &snes);
1007: if (!r && !ts->dm && ts->vec_sol) {
1008: VecDuplicate(ts->vec_sol, &ralloc);
1009: r = ralloc;
1010: }
1011: SNESSetFunction(snes, r, SNESTSFormFunction, ts);
1012: VecDestroy(&ralloc);
1013: return 0;
1014: }
1016: /*@C
1017: TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE
1019: Logically Collective
1021: Input Parameters:
1022: + ts - the `TS` context obtained from `TSCreate()`
1023: . f - routine for evaluating the solution
1024: - ctx - [optional] user-defined context for private data for the
1025: function evaluation routine (may be NULL)
1027: Calling sequence of f:
1028: $ PetscErrorCode f(TS ts,PetscReal t,Vec u,void *ctx);
1030: + t - current timestep
1031: . u - output vector
1032: - ctx - [optional] user-defined function context
1034: Options Database Keys:
1035: + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
1036: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`
1038: Level: intermediate
1040: Notes:
1041: This routine is used for testing accuracy of time integration schemes when you already know the solution.
1042: If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
1043: create closed-form solutions with non-physical forcing terms.
1045: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1047: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
1048: @*/
1049: PetscErrorCode TSSetSolutionFunction(TS ts, PetscErrorCode (*f)(TS, PetscReal, Vec, void *), void *ctx)
1050: {
1051: DM dm;
1054: TSGetDM(ts, &dm);
1055: DMTSSetSolutionFunction(dm, f, ctx);
1056: return 0;
1057: }
1059: /*@C
1060: TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE
1062: Logically Collective
1064: Input Parameters:
1065: + ts - the `TS` context obtained from `TSCreate()`
1066: . func - routine for evaluating the forcing function
1067: - ctx - [optional] user-defined context for private data for the
1068: function evaluation routine (may be NULL)
1070: Calling sequence of func:
1071: $ PetscErrorCode func (TS ts,PetscReal t,Vec f,void *ctx);
1073: + t - current timestep
1074: . f - output vector
1075: - ctx - [optional] user-defined function context
1077: Level: intermediate
1079: Notes:
1080: This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
1081: create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
1082: definition of the problem you are solving and hence possibly introducing bugs.
1084: This replaces the ODE F(u,u_t,t) = 0 the TS is solving with F(u,u_t,t) - func(t) = 0
1086: This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
1087: parameters can be passed in the ctx variable.
1089: For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.
1091: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
1092: @*/
1093: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFunction func, void *ctx)
1094: {
1095: DM dm;
1098: TSGetDM(ts, &dm);
1099: DMTSSetForcingFunction(dm, func, ctx);
1100: return 0;
1101: }
1103: /*@C
1104: TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
1105: where U_t = G(U,t), as well as the location to store the matrix.
1107: Logically Collective
1109: Input Parameters:
1110: + ts - the `TS` context obtained from `TSCreate()`
1111: . Amat - (approximate) Jacobian matrix
1112: . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1113: . f - the Jacobian evaluation routine
1114: - ctx - [optional] user-defined context for private data for the
1115: Jacobian evaluation routine (may be NULL)
1117: Calling sequence of f:
1118: $ PetscErrorCode f(TS ts,PetscReal t,Vec u,Mat A,Mat B,void *ctx);
1120: + t - current timestep
1121: . u - input vector
1122: . Amat - (approximate) Jacobian matrix
1123: . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat)
1124: - ctx - [optional] user-defined context for matrix evaluation routine
1126: Level: beginner
1128: Notes:
1129: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1131: The `TS` solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f()
1132: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1134: .seealso: [](chapter_ts), `TS`, `SNESComputeJacobianDefaultColor()`, `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`
1135: @*/
1136: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobian f, void *ctx)
1137: {
1138: SNES snes;
1139: DM dm;
1140: TSIJacobian ijacobian;
1148: TSGetDM(ts, &dm);
1149: DMTSSetRHSJacobian(dm, f, ctx);
1150: DMTSGetIJacobian(dm, &ijacobian, NULL);
1151: TSGetSNES(ts, &snes);
1152: if (!ijacobian) SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts);
1153: if (Amat) {
1154: PetscObjectReference((PetscObject)Amat);
1155: MatDestroy(&ts->Arhs);
1156: ts->Arhs = Amat;
1157: }
1158: if (Pmat) {
1159: PetscObjectReference((PetscObject)Pmat);
1160: MatDestroy(&ts->Brhs);
1161: ts->Brhs = Pmat;
1162: }
1163: return 0;
1164: }
1166: /*@C
1167: TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.
1169: Logically Collective
1171: Input Parameters:
1172: + ts - the `TS` context obtained from `TSCreate()`
1173: . r - vector to hold the residual (or NULL to have it created internally)
1174: . f - the function evaluation routine
1175: - ctx - user-defined context for private data for the function evaluation routine (may be NULL)
1177: Calling sequence of f:
1178: $ PetscErrorCode f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx);
1180: + t - time at step/stage being solved
1181: . u - state vector
1182: . u_t - time derivative of state vector
1183: . F - function vector
1184: - ctx - [optional] user-defined context for matrix evaluation routine
1186: Level: beginner
1188: Note:
1189: The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function.
1191: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`, `TSSetIJacobian()`
1192: @*/
1193: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunction f, void *ctx)
1194: {
1195: SNES snes;
1196: Vec ralloc = NULL;
1197: DM dm;
1202: TSGetDM(ts, &dm);
1203: DMTSSetIFunction(dm, f, ctx);
1205: TSGetSNES(ts, &snes);
1206: if (!r && !ts->dm && ts->vec_sol) {
1207: VecDuplicate(ts->vec_sol, &ralloc);
1208: r = ralloc;
1209: }
1210: SNESSetFunction(snes, r, SNESTSFormFunction, ts);
1211: VecDestroy(&ralloc);
1212: return 0;
1213: }
1215: /*@C
1216: TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.
1218: Not Collective
1220: Input Parameter:
1221: . ts - the `TS` context
1223: Output Parameters:
1224: + r - vector to hold residual (or NULL)
1225: . func - the function to compute residual (or NULL)
1226: - ctx - the function context (or NULL)
1228: Level: advanced
1230: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
1231: @*/
1232: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunction *func, void **ctx)
1233: {
1234: SNES snes;
1235: DM dm;
1238: TSGetSNES(ts, &snes);
1239: SNESGetFunction(snes, r, NULL, NULL);
1240: TSGetDM(ts, &dm);
1241: DMTSGetIFunction(dm, func, ctx);
1242: return 0;
1243: }
1245: /*@C
1246: TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it.
1248: Not Collective
1250: Input Parameter:
1251: . ts - the `TS` context
1253: Output Parameters:
1254: + r - vector to hold computed right hand side (or NULL)
1255: . func - the function to compute right hand side (or NULL)
1256: - ctx - the function context (or NULL)
1258: Level: advanced
1260: .seealso: [](chapter_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
1261: @*/
1262: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunction *func, void **ctx)
1263: {
1264: SNES snes;
1265: DM dm;
1268: TSGetSNES(ts, &snes);
1269: SNESGetFunction(snes, r, NULL, NULL);
1270: TSGetDM(ts, &dm);
1271: DMTSGetRHSFunction(dm, func, ctx);
1272: return 0;
1273: }
1275: /*@C
1276: TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
1277: provided with `TSSetIFunction()`.
1279: Logically Collective
1281: Input Parameters:
1282: + ts - the `TS` context obtained from `TSCreate()`
1283: . Amat - (approximate) Jacobian matrix
1284: . Pmat - matrix used to compute preconditioner (usually the same as Amat)
1285: . f - the Jacobian evaluation routine
1286: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
1288: Calling sequence of f:
1289: $ PetscErrorCode f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat Amat,Mat Pmat,void *ctx);
1291: + t - time at step/stage being solved
1292: . U - state vector
1293: . U_t - time derivative of state vector
1294: . a - shift
1295: . Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t
1296: . Pmat - matrix used for constructing preconditioner, usually the same as Amat
1297: - ctx - [optional] user-defined context for matrix evaluation routine
1299: Level: beginner
1301: Notes:
1302: The matrices Amat and Pmat are exactly the matrices that are used by `SNES` for the nonlinear solve.
1304: If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
1305: space to Amat and the `KSP` solvers will automatically use that null space as needed during the solution process.
1307: The matrix dF/dU + a*dF/dU_t you provide turns out to be
1308: the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
1309: The time integrator internally approximates U_t by W+a*U where the positive "shift"
1310: a and vector W depend on the integration method, step size, and past states. For example with
1311: the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
1312: W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt
1314: You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value
1316: The TS solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f()
1317: You should not assume the values are the same in the next call to f() as you set them in the previous call.
1319: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `TSSetRHSJacobian()`, `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
1320: @*/
1321: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobian f, void *ctx)
1322: {
1323: SNES snes;
1324: DM dm;
1332: TSGetDM(ts, &dm);
1333: DMTSSetIJacobian(dm, f, ctx);
1335: TSGetSNES(ts, &snes);
1336: SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts);
1337: return 0;
1338: }
1340: /*@
1341: TSRHSJacobianSetReuse - restore RHS Jacobian before re-evaluating. Without this flag, `TS` will change the sign and
1342: shift the RHS Jacobian for a finite-time-step implicit solve, in which case the user function will need to recompute
1343: the entire Jacobian. The reuse flag must be set if the evaluation function will assume that the matrix entries have
1344: not been changed by the TS.
1346: Logically Collective
1348: Input Parameters:
1349: + ts - `TS` context obtained from `TSCreate()`
1350: - reuse - `PETSC_TRUE` if the RHS Jacobian
1352: Level: intermediate
1354: .seealso: [](chapter_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
1355: @*/
1356: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
1357: {
1358: ts->rhsjacobian.reuse = reuse;
1359: return 0;
1360: }
1362: /*@C
1363: TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.
1365: Logically Collective
1367: Input Parameters:
1368: + ts - the `TS` context obtained from `TSCreate()`
1369: . F - vector to hold the residual (or NULL to have it created internally)
1370: . fun - the function evaluation routine
1371: - ctx - user-defined context for private data for the function evaluation routine (may be NULL)
1373: Calling sequence of fun:
1374: $ PetscErrorCode fun(TS ts,PetscReal t,Vec U,Vec U_t,Vec U_tt,Vec F,ctx);
1376: + t - time at step/stage being solved
1377: . U - state vector
1378: . U_t - time derivative of state vector
1379: . U_tt - second time derivative of state vector
1380: . F - function vector
1381: - ctx - [optional] user-defined context for matrix evaluation routine (may be NULL)
1383: Level: beginner
1385: .seealso: [](chapter_ts), `TS`, `TSSetI2Jacobian()`, `TSSetIFunction()`, `TSCreate()`, `TSSetRHSFunction()`
1386: @*/
1387: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2Function fun, void *ctx)
1388: {
1389: DM dm;
1393: TSSetIFunction(ts, F, NULL, NULL);
1394: TSGetDM(ts, &dm);
1395: DMTSSetI2Function(dm, fun, ctx);
1396: return 0;
1397: }
1399: /*@C
1400: TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.
1402: Not Collective
1404: Input Parameter:
1405: . ts - the `TS` context
1407: Output Parameters:
1408: + r - vector to hold residual (or NULL)
1409: . fun - the function to compute residual (or NULL)
1410: - ctx - the function context (or NULL)
1412: Level: advanced
1414: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
1415: @*/
1416: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2Function *fun, void **ctx)
1417: {
1418: SNES snes;
1419: DM dm;
1422: TSGetSNES(ts, &snes);
1423: SNESGetFunction(snes, r, NULL, NULL);
1424: TSGetDM(ts, &dm);
1425: DMTSGetI2Function(dm, fun, ctx);
1426: return 0;
1427: }
1429: /*@C
1430: TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt
1431: where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.
1433: Logically Collective
1435: Input Parameters:
1436: + ts - the `TS` context obtained from `TSCreate()`
1437: . J - Jacobian matrix
1438: . P - preconditioning matrix for J (may be same as J)
1439: . jac - the Jacobian evaluation routine
1440: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
1442: Calling sequence of jac:
1443: $ PetscErrorCode jac(TS ts,PetscReal t,Vec U,Vec U_t,Vec U_tt,PetscReal v,PetscReal a,Mat J,Mat P,void *ctx);
1445: + t - time at step/stage being solved
1446: . U - state vector
1447: . U_t - time derivative of state vector
1448: . U_tt - second time derivative of state vector
1449: . v - shift for U_t
1450: . a - shift for U_tt
1451: . J - Jacobian of G(U) = F(t,U,W+v*U,W'+a*U), equivalent to dF/dU + v*dF/dU_t + a*dF/dU_tt
1452: . P - preconditioning matrix for J, may be same as J
1453: - ctx - [optional] user-defined context for matrix evaluation routine
1455: Level: beginner
1457: Notes:
1458: The matrices J and P are exactly the matrices that are used by `SNES` for the nonlinear solve.
1460: The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
1461: 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.
1462: The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift"
1463: parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.
1465: .seealso: [](chapter_ts), `TS`, `TSSetI2Function()`, `TSGetI2Jacobian()`
1466: @*/
1467: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2Jacobian jac, void *ctx)
1468: {
1469: DM dm;
1474: TSSetIJacobian(ts, J, P, NULL, NULL);
1475: TSGetDM(ts, &dm);
1476: DMTSSetI2Jacobian(dm, jac, ctx);
1477: return 0;
1478: }
1480: /*@C
1481: TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.
1483: Not Collective, but parallel objects are returned if `TS` is parallel
1485: Input Parameter:
1486: . ts - The `TS` context obtained from `TSCreate()`
1488: Output Parameters:
1489: + J - The (approximate) Jacobian of F(t,U,U_t,U_tt)
1490: . P - The matrix from which the preconditioner is constructed, often the same as J
1491: . jac - The function to compute the Jacobian matrices
1492: - ctx - User-defined context for Jacobian evaluation routine
1494: Level: advanced
1496: Notes:
1497: You can pass in NULL for any return argument you do not need.
1499: .seealso: [](chapter_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
1500: @*/
1501: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2Jacobian *jac, void **ctx)
1502: {
1503: SNES snes;
1504: DM dm;
1506: TSGetSNES(ts, &snes);
1507: SNESSetUpMatrices(snes);
1508: SNESGetJacobian(snes, J, P, NULL, NULL);
1509: TSGetDM(ts, &dm);
1510: DMTSGetI2Jacobian(dm, jac, ctx);
1511: return 0;
1512: }
1514: /*@
1515: TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0
1517: Collective
1519: Input Parameters:
1520: + ts - the `TS` context
1521: . t - current time
1522: . U - state vector
1523: . V - time derivative of state vector (U_t)
1524: - A - second time derivative of state vector (U_tt)
1526: Output Parameter:
1527: . F - the residual vector
1529: Level: developer
1531: Note:
1532: Most users should not need to explicitly call this routine, as it
1533: is used internally within the nonlinear solvers.
1535: .seealso: [](chapter_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
1536: @*/
1537: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
1538: {
1539: DM dm;
1540: TSI2Function I2Function;
1541: void *ctx;
1542: TSRHSFunction rhsfunction;
1550: TSGetDM(ts, &dm);
1551: DMTSGetI2Function(dm, &I2Function, &ctx);
1552: DMTSGetRHSFunction(dm, &rhsfunction, NULL);
1554: if (!I2Function) {
1555: TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE);
1556: return 0;
1557: }
1559: PetscLogEventBegin(TS_FunctionEval, ts, U, V, F);
1561: PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));
1563: if (rhsfunction) {
1564: Vec Frhs;
1565: TSGetRHSVec_Private(ts, &Frhs);
1566: TSComputeRHSFunction(ts, t, U, Frhs);
1567: VecAXPY(F, -1, Frhs);
1568: }
1570: PetscLogEventEnd(TS_FunctionEval, ts, U, V, F);
1571: return 0;
1572: }
1574: /*@
1575: TSComputeI2Jacobian - Evaluates the Jacobian of the DAE
1577: Collective
1579: Input Parameters:
1580: + ts - the `TS` context
1581: . t - current timestep
1582: . U - state vector
1583: . V - time derivative of state vector
1584: . A - second time derivative of state vector
1585: . shiftV - shift to apply, see note below
1586: - shiftA - shift to apply, see note below
1588: Output Parameters:
1589: + J - Jacobian matrix
1590: - P - optional preconditioning matrix
1592: Level: developer
1594: Notes:
1595: If F(t,U,V,A)=0 is the DAE, the required Jacobian is
1597: dF/dU + shiftV*dF/dV + shiftA*dF/dA
1599: Most users should not need to explicitly call this routine, as it
1600: is used internally within the nonlinear solvers.
1602: .seealso: [](chapter_ts), `TS`, `TSSetI2Jacobian()`
1603: @*/
1604: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
1605: {
1606: DM dm;
1607: TSI2Jacobian I2Jacobian;
1608: void *ctx;
1609: TSRHSJacobian rhsjacobian;
1618: TSGetDM(ts, &dm);
1619: DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx);
1620: DMTSGetRHSJacobian(dm, &rhsjacobian, NULL);
1622: if (!I2Jacobian) {
1623: TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE);
1624: return 0;
1625: }
1627: PetscLogEventBegin(TS_JacobianEval, ts, U, J, P);
1628: PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
1629: if (rhsjacobian) {
1630: Mat Jrhs, Prhs;
1631: TSGetRHSMats_Private(ts, &Jrhs, &Prhs);
1632: TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs);
1633: MatAXPY(J, -1, Jrhs, ts->axpy_pattern);
1634: if (P != J) MatAXPY(P, -1, Prhs, ts->axpy_pattern);
1635: }
1637: PetscLogEventEnd(TS_JacobianEval, ts, U, J, P);
1638: return 0;
1639: }
1641: /*@C
1642: TSSetTransientVariable - sets function to transform from state to transient variables
1644: Logically Collective
1646: Input Parameters:
1647: + ts - time stepping context on which to change the transient variable
1648: . tvar - a function that transforms to transient variables
1649: - ctx - a context for tvar
1651: Calling sequence of tvar:
1652: $ PetscErrorCode tvar(TS ts,Vec p,Vec c,void *ctx);
1654: + ts - timestep context
1655: . p - input vector (primitive form)
1656: . c - output vector, transient variables (conservative form)
1657: - ctx - [optional] user-defined function context
1659: Level: advanced
1661: Notes:
1662: This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
1663: can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to
1664: well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is
1665: C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
1666: evaluated via the chain rule, as in
1668: dF/dP + shift * dF/dCdot dC/dP.
1670: .seealso: [](chapter_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
1671: @*/
1672: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariable tvar, void *ctx)
1673: {
1674: DM dm;
1677: TSGetDM(ts, &dm);
1678: DMTSSetTransientVariable(dm, tvar, ctx);
1679: return 0;
1680: }
1682: /*@
1683: TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables
1685: Logically Collective
1687: Input Parameters:
1688: + ts - TS on which to compute
1689: - U - state vector to be transformed to transient variables
1691: Output Parameters:
1692: . C - transient (conservative) variable
1694: Level: developer
1696: Developer Note:
1697: If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = NULL is allowed.
1698: This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are
1699: being used.
1701: .seealso: [](chapter_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
1702: @*/
1703: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
1704: {
1705: DM dm;
1706: DMTS dmts;
1710: TSGetDM(ts, &dm);
1711: DMGetDMTS(dm, &dmts);
1712: if (dmts->ops->transientvar) {
1714: (*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx);
1715: }
1716: return 0;
1717: }
1719: /*@
1720: TSHasTransientVariable - determine whether transient variables have been set
1722: Logically Collective
1724: Input Parameters:
1725: . ts - TS on which to compute
1727: Output Parameters:
1728: . has - `PETSC_TRUE` if transient variables have been set
1730: Level: developer
1732: .seealso: [](chapter_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
1733: @*/
1734: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
1735: {
1736: DM dm;
1737: DMTS dmts;
1740: TSGetDM(ts, &dm);
1741: DMGetDMTS(dm, &dmts);
1742: *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
1743: return 0;
1744: }
1746: /*@
1747: TS2SetSolution - Sets the initial solution and time derivative vectors
1748: for use by the `TS` routines handling second order equations.
1750: Logically Collective
1752: Input Parameters:
1753: + ts - the `TS` context obtained from `TSCreate()`
1754: . u - the solution vector
1755: - v - the time derivative vector
1757: Level: beginner
1759: .seealso: [](chapter_ts), `TS`
1760: @*/
1761: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
1762: {
1766: TSSetSolution(ts, u);
1767: PetscObjectReference((PetscObject)v);
1768: VecDestroy(&ts->vec_dot);
1769: ts->vec_dot = v;
1770: return 0;
1771: }
1773: /*@
1774: TS2GetSolution - Returns the solution and time derivative at the present timestep
1775: for second order equations. It is valid to call this routine inside the function
1776: that you are evaluating in order to move to the new timestep. This vector not
1777: changed until the solution at the next timestep has been calculated.
1779: Not Collective, but Vec returned is parallel if `TS` is parallel
1781: Input Parameter:
1782: . ts - the `TS` context obtained from `TSCreate()`
1784: Output Parameters:
1785: + u - the vector containing the solution
1786: - v - the vector containing the time derivative
1788: Level: intermediate
1790: .seealso: [](chapter_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
1791: @*/
1792: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
1793: {
1797: if (u) *u = ts->vec_sol;
1798: if (v) *v = ts->vec_dot;
1799: return 0;
1800: }
1802: /*@C
1803: TSLoad - Loads a `TS` that has been stored in binary with `TSView()`.
1805: Collective
1807: Input Parameters:
1808: + newdm - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
1809: some related function before a call to `TSLoad()`.
1810: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
1812: Level: intermediate
1814: Note:
1815: The type is determined by the data in the file, any type set into the `TS` before this call is ignored.
1817: .seealso: [](chapter_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
1818: @*/
1819: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
1820: {
1821: PetscBool isbinary;
1822: PetscInt classid;
1823: char type[256];
1824: DMTS sdm;
1825: DM dm;
1829: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
1832: PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT);
1834: PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR);
1835: TSSetType(ts, type);
1836: PetscTryTypeMethod(ts, load, viewer);
1837: DMCreate(PetscObjectComm((PetscObject)ts), &dm);
1838: DMLoad(dm, viewer);
1839: TSSetDM(ts, dm);
1840: DMCreateGlobalVector(ts->dm, &ts->vec_sol);
1841: VecLoad(ts->vec_sol, viewer);
1842: DMGetDMTS(ts->dm, &sdm);
1843: DMTSLoad(sdm, viewer);
1844: return 0;
1845: }
1847: #include <petscdraw.h>
1848: #if defined(PETSC_HAVE_SAWS)
1849: #include <petscviewersaws.h>
1850: #endif
1852: /*@C
1853: TSViewFromOptions - View a `TS` based on values in the options database
1855: Collective
1857: Input Parameters:
1858: + ts - the `TS` context
1859: . obj - Optional object that provides the prefix for the options database keys
1860: - name - command line option string to be passed by user
1862: Level: intermediate
1864: .seealso: [](chapter_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
1865: @*/
1866: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
1867: {
1869: PetscObjectViewFromOptions((PetscObject)ts, obj, name);
1870: return 0;
1871: }
1873: /*@C
1874: TSView - Prints the `TS` data structure.
1876: Collective
1878: Input Parameters:
1879: + ts - the `TS` context obtained from `TSCreate()`
1880: - viewer - visualization context
1882: Options Database Key:
1883: . -ts_view - calls `TSView()` at end of `TSStep()`
1885: Level: beginner
1887: Notes:
1888: The available visualization contexts include
1889: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
1890: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
1891: output where only the first processor opens
1892: the file. All other processors send their
1893: data to the first processor to print.
1895: The user can open an alternative visualization context with
1896: `PetscViewerASCIIOpen()` - output to a specified file.
1898: In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).
1900: .seealso: [](chapter_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
1901: @*/
1902: PetscErrorCode TSView(TS ts, PetscViewer viewer)
1903: {
1904: TSType type;
1905: PetscBool iascii, isstring, isundials, isbinary, isdraw;
1906: DMTS sdm;
1907: #if defined(PETSC_HAVE_SAWS)
1908: PetscBool issaws;
1909: #endif
1912: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer);
1916: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1917: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
1918: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
1919: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1920: #if defined(PETSC_HAVE_SAWS)
1921: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws);
1922: #endif
1923: if (iascii) {
1924: PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer);
1925: if (ts->ops->view) {
1926: PetscViewerASCIIPushTab(viewer);
1927: PetscUseTypeMethod(ts, view, viewer);
1928: PetscViewerASCIIPopTab(viewer);
1929: }
1930: if (ts->max_steps < PETSC_MAX_INT) PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps);
1931: if (ts->max_time < PETSC_MAX_REAL) PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time);
1932: if (ts->ifuncs) PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs);
1933: if (ts->ijacs) PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs);
1934: if (ts->rhsfuncs) PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs);
1935: if (ts->rhsjacs) PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs);
1936: if (ts->usessnes) {
1937: PetscBool lin;
1938: if (ts->problem_type == TS_NONLINEAR) PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its);
1939: PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its);
1940: PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, "");
1941: PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures);
1942: }
1943: PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject);
1944: if (ts->vrtol) PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, ");
1945: else PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol);
1946: if (ts->vatol) PetscViewerASCIIPrintf(viewer, " using vector of absolute error tolerances\n");
1947: else PetscViewerASCIIPrintf(viewer, " using absolute error tolerance of %g\n", (double)ts->atol);
1948: PetscViewerASCIIPushTab(viewer);
1949: TSAdaptView(ts->adapt, viewer);
1950: PetscViewerASCIIPopTab(viewer);
1951: } else if (isstring) {
1952: TSGetType(ts, &type);
1953: PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type);
1954: PetscTryTypeMethod(ts, view, viewer);
1955: } else if (isbinary) {
1956: PetscInt classid = TS_FILE_CLASSID;
1957: MPI_Comm comm;
1958: PetscMPIInt rank;
1959: char type[256];
1961: PetscObjectGetComm((PetscObject)ts, &comm);
1962: MPI_Comm_rank(comm, &rank);
1963: if (rank == 0) {
1964: PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT);
1965: PetscStrncpy(type, ((PetscObject)ts)->type_name, 256);
1966: PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR);
1967: }
1968: PetscTryTypeMethod(ts, view, viewer);
1969: if (ts->adapt) TSAdaptView(ts->adapt, viewer);
1970: DMView(ts->dm, viewer);
1971: VecView(ts->vec_sol, viewer);
1972: DMGetDMTS(ts->dm, &sdm);
1973: DMTSView(sdm, viewer);
1974: } else if (isdraw) {
1975: PetscDraw draw;
1976: char str[36];
1977: PetscReal x, y, bottom, h;
1979: PetscViewerDrawGetDraw(viewer, 0, &draw);
1980: PetscDrawGetCurrentPoint(draw, &x, &y);
1981: PetscStrcpy(str, "TS: ");
1982: PetscStrcat(str, ((PetscObject)ts)->type_name);
1983: PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h);
1984: bottom = y - h;
1985: PetscDrawPushCurrentPoint(draw, x, bottom);
1986: PetscTryTypeMethod(ts, view, viewer);
1987: if (ts->adapt) TSAdaptView(ts->adapt, viewer);
1988: if (ts->snes) SNESView(ts->snes, viewer);
1989: PetscDrawPopCurrentPoint(draw);
1990: #if defined(PETSC_HAVE_SAWS)
1991: } else if (issaws) {
1992: PetscMPIInt rank;
1993: const char *name;
1995: PetscObjectGetName((PetscObject)ts, &name);
1996: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1997: if (!((PetscObject)ts)->amsmem && rank == 0) {
1998: char dir[1024];
2000: PetscObjectViewSAWs((PetscObject)ts, viewer);
2001: PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name);
2002: SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT);
2003: PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name);
2004: SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE);
2005: }
2006: PetscTryTypeMethod(ts, view, viewer);
2007: #endif
2008: }
2009: if (ts->snes && ts->usessnes) {
2010: PetscViewerASCIIPushTab(viewer);
2011: SNESView(ts->snes, viewer);
2012: PetscViewerASCIIPopTab(viewer);
2013: }
2014: DMGetDMTS(ts->dm, &sdm);
2015: DMTSView(sdm, viewer);
2017: PetscViewerASCIIPushTab(viewer);
2018: PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials);
2019: PetscViewerASCIIPopTab(viewer);
2020: return 0;
2021: }
2023: /*@
2024: TSSetApplicationContext - Sets an optional user-defined context for
2025: the timesteppers.
2027: Logically Collective
2029: Input Parameters:
2030: + ts - the `TS` context obtained from `TSCreate()`
2031: - usrP - user context
2033: Level: intermediate
2035: Fortran Note:
2036: To use this from Fortran you must write a Fortran interface definition for this
2037: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
2039: .seealso: [](chapter_ts), `TS`, `TSGetApplicationContext()`
2040: @*/
2041: PetscErrorCode TSSetApplicationContext(TS ts, void *usrP)
2042: {
2044: ts->user = usrP;
2045: return 0;
2046: }
2048: /*@
2049: TSGetApplicationContext - Gets the user-defined context for the
2050: timestepper that was set with `TSSetApplicationContext()`
2052: Not Collective
2054: Input Parameter:
2055: . ts - the `TS` context obtained from `TSCreate()`
2057: Output Parameter:
2058: . usrP - user context
2060: Level: intermediate
2062: Fortran Note:
2063: To use this from Fortran you must write a Fortran interface definition for this
2064: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
2066: .seealso: [](chapter_ts), `TS`, `TSSetApplicationContext()`
2067: @*/
2068: PetscErrorCode TSGetApplicationContext(TS ts, void *usrP)
2069: {
2071: *(void **)usrP = ts->user;
2072: return 0;
2073: }
2075: /*@
2076: TSGetStepNumber - Gets the number of time steps completed.
2078: Not Collective
2080: Input Parameter:
2081: . ts - the `TS` context obtained from `TSCreate()`
2083: Output Parameter:
2084: . steps - number of steps completed so far
2086: Level: intermediate
2088: .seealso: [](chapter_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
2089: @*/
2090: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
2091: {
2094: *steps = ts->steps;
2095: return 0;
2096: }
2098: /*@
2099: TSSetStepNumber - Sets the number of steps completed.
2101: Logically Collective
2103: Input Parameters:
2104: + ts - the `TS` context
2105: - steps - number of steps completed so far
2107: Level: developer
2109: Note:
2110: For most uses of the `TS` solvers the user need not explicitly call
2111: `TSSetStepNumber()`, as the step counter is appropriately updated in
2112: `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
2113: reinitialize timestepping by setting the step counter to zero (and time
2114: to the initial time) to solve a similar problem with different initial
2115: conditions or parameters. Other possible use case is to continue
2116: timestepping from a previously interrupted run in such a way that TS
2117: monitors will be called with a initial nonzero step counter.
2119: .seealso: [](chapter_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
2120: @*/
2121: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
2122: {
2126: ts->steps = steps;
2127: return 0;
2128: }
2130: /*@
2131: TSSetTimeStep - Allows one to reset the timestep at any time,
2132: useful for simple pseudo-timestepping codes.
2134: Logically Collective
2136: Input Parameters:
2137: + ts - the `TS` context obtained from `TSCreate()`
2138: - time_step - the size of the timestep
2140: Level: intermediate
2142: .seealso: [](chapter_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
2143: @*/
2144: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
2145: {
2148: ts->time_step = time_step;
2149: return 0;
2150: }
2152: /*@
2153: TSSetExactFinalTime - Determines whether to adapt the final time step to
2154: match the exact final time, interpolate solution to the exact final time,
2155: or just return at the final time `TS` computed.
2157: Logically Collective
2159: Input Parameters:
2160: + ts - the time-step context
2161: - eftopt - exact final time option
2162: .vb
2163: TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded
2164: TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time
2165: TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time
2166: .ve
2168: Options Database Key:
2169: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step at runtime
2171: Level: beginner
2173: Note:
2174: If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
2175: then the final time you selected.
2177: .seealso: [](chapter_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
2178: @*/
2179: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
2180: {
2183: ts->exact_final_time = eftopt;
2184: return 0;
2185: }
2187: /*@
2188: TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`
2190: Not Collective
2192: Input Parameter:
2193: . ts - the `TS` context
2195: Output Parameter:
2196: . eftopt - exact final time option
2198: Level: beginner
2200: .seealso: [](chapter_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
2201: @*/
2202: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
2203: {
2206: *eftopt = ts->exact_final_time;
2207: return 0;
2208: }
2210: /*@
2211: TSGetTimeStep - Gets the current timestep size.
2213: Not Collective
2215: Input Parameter:
2216: . ts - the `TS` context obtained from `TSCreate()`
2218: Output Parameter:
2219: . dt - the current timestep size
2221: Level: intermediate
2223: .seealso: [](chapter_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
2224: @*/
2225: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
2226: {
2229: *dt = ts->time_step;
2230: return 0;
2231: }
2233: /*@
2234: TSGetSolution - Returns the solution at the present timestep. It
2235: is valid to call this routine inside the function that you are evaluating
2236: in order to move to the new timestep. This vector not changed until
2237: the solution at the next timestep has been calculated.
2239: Not Collective, but v returned is parallel if ts is parallel
2241: Input Parameter:
2242: . ts - the `TS` context obtained from `TSCreate()`
2244: Output Parameter:
2245: . v - the vector containing the solution
2247: Level: intermediate
2249: Note:
2250: If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
2251: final time. It returns the solution at the next timestep.
2253: .seealso: [](chapter_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
2254: @*/
2255: PetscErrorCode TSGetSolution(TS ts, Vec *v)
2256: {
2259: *v = ts->vec_sol;
2260: return 0;
2261: }
2263: /*@
2264: TSGetSolutionComponents - Returns any solution components at the present
2265: timestep, if available for the time integration method being used.
2266: Solution components are quantities that share the same size and
2267: structure as the solution vector.
2269: Not Collective, but v returned is parallel if ts is parallel
2271: Parameters :
2272: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2273: . n - If v is PETSC_NULL, then the number of solution components is
2274: returned through n, else the n-th solution component is
2275: returned in v.
2276: - v - the vector containing the n-th solution component
2277: (may be PETSC_NULL to use this function to find out
2278: the number of solutions components).
2280: Level: advanced
2282: .seealso: [](chapter_ts), `TS`, `TSGetSolution()`
2283: @*/
2284: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
2285: {
2287: if (!ts->ops->getsolutioncomponents) *n = 0;
2288: else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
2289: return 0;
2290: }
2292: /*@
2293: TSGetAuxSolution - Returns an auxiliary solution at the present
2294: timestep, if available for the time integration method being used.
2296: Not Collective, but v returned is parallel if ts is parallel
2298: Parameters :
2299: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2300: - v - the vector containing the auxiliary solution
2302: Level: intermediate
2304: .seealso: [](chapter_ts), `TS`, `TSGetSolution()`
2305: @*/
2306: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
2307: {
2309: if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
2310: else VecZeroEntries(*v);
2311: return 0;
2312: }
2314: /*@
2315: TSGetTimeError - Returns the estimated error vector, if the chosen
2316: `TSType` has an error estimation functionality and `TSSetTimeError()` was called
2318: Not Collective, but v returned is parallel if ts is parallel
2320: Parameters :
2321: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2322: . n - current estimate (n=0) or previous one (n=-1)
2323: - v - the vector containing the error (same size as the solution).
2325: Level: intermediate
2327: Note:
2328: MUST call after `TSSetUp()`
2330: .seealso: [](chapter_ts), `TSGetSolution()`, `TSSetTimeError()`
2331: @*/
2332: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
2333: {
2335: if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
2336: else VecZeroEntries(*v);
2337: return 0;
2338: }
2340: /*@
2341: TSSetTimeError - Sets the estimated error vector, if the chosen
2342: `TSType` has an error estimation functionality. This can be used
2343: to restart such a time integrator with a given error vector.
2345: Not Collective, but v returned is parallel if ts is parallel
2347: Parameters :
2348: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
2349: - v - the vector containing the error (same size as the solution).
2351: Level: intermediate
2353: .seealso: [](chapter_ts), `TS`, `TSSetSolution()`, `TSGetTimeError)`
2354: @*/
2355: PetscErrorCode TSSetTimeError(TS ts, Vec v)
2356: {
2359: PetscTryTypeMethod(ts, settimeerror, v);
2360: return 0;
2361: }
2363: /* ----- Routines to initialize and destroy a timestepper ---- */
2364: /*@
2365: TSSetProblemType - Sets the type of problem to be solved.
2367: Not collective
2369: Input Parameters:
2370: + ts - The `TS`
2371: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2372: .vb
2373: U_t - A U = 0 (linear)
2374: U_t - A(t) U = 0 (linear)
2375: F(t,U,U_t) = 0 (nonlinear)
2376: .ve
2378: Level: beginner
2380: .seealso: [](chapter_ts), `TSSetUp()`, `TSProblemType`, `TS`
2381: @*/
2382: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
2383: {
2385: ts->problem_type = type;
2386: if (type == TS_LINEAR) {
2387: SNES snes;
2388: TSGetSNES(ts, &snes);
2389: SNESSetType(snes, SNESKSPONLY);
2390: }
2391: return 0;
2392: }
2394: /*@C
2395: TSGetProblemType - Gets the type of problem to be solved.
2397: Not collective
2399: Input Parameter:
2400: . ts - The `TS`
2402: Output Parameter:
2403: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
2404: .vb
2405: M U_t = A U
2406: M(t) U_t = A(t) U
2407: F(t,U,U_t)
2408: .ve
2410: Level: beginner
2412: .seealso: [](chapter_ts), `TSSetUp()`, `TSProblemType`, `TS`
2413: @*/
2414: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
2415: {
2418: *type = ts->problem_type;
2419: return 0;
2420: }
2422: /*
2423: Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
2424: */
2425: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
2426: {
2427: PetscBool isnone;
2429: TSGetAdapt(ts, &ts->adapt);
2430: TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type);
2432: PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone);
2433: if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
2434: else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
2435: return 0;
2436: }
2438: /*@
2439: TSSetUp - Sets up the internal data structures for the later use of a timestepper.
2441: Collective
2443: Input Parameter:
2444: . ts - the `TS` context obtained from `TSCreate()`
2446: Level: advanced
2448: Note:
2449: For basic use of the `TS` solvers the user need not explicitly call
2450: `TSSetUp()`, since these actions will automatically occur during
2451: the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this
2452: phase separately, `TSSetUp()` should be called after `TSCreate()`
2453: and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.
2455: .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
2456: @*/
2457: PetscErrorCode TSSetUp(TS ts)
2458: {
2459: DM dm;
2460: PetscErrorCode (*func)(SNES, Vec, Vec, void *);
2461: PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
2462: TSIFunction ifun;
2463: TSIJacobian ijac;
2464: TSI2Jacobian i2jac;
2465: TSRHSJacobian rhsjac;
2468: if (ts->setupcalled) return 0;
2470: if (!((PetscObject)ts)->type_name) {
2471: TSGetIFunction(ts, NULL, &ifun, NULL);
2472: TSSetType(ts, ifun ? TSBEULER : TSEULER);
2473: }
2475: if (!ts->vec_sol) {
2477: DMCreateGlobalVector(ts->dm, &ts->vec_sol);
2478: }
2480: if (ts->tspan) {
2481: if (!ts->tspan->vecs_sol) VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol);
2482: }
2483: if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
2484: PetscObjectReference((PetscObject)ts->Jacprhs);
2485: ts->Jacp = ts->Jacprhs;
2486: }
2488: if (ts->quadraturets) {
2489: TSSetUp(ts->quadraturets);
2490: VecDestroy(&ts->vec_costintegrand);
2491: VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand);
2492: }
2494: TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL);
2495: if (rhsjac == TSComputeRHSJacobianConstant) {
2496: Mat Amat, Pmat;
2497: SNES snes;
2498: TSGetSNES(ts, &snes);
2499: SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL);
2500: /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
2501: * have displaced the RHS matrix */
2502: if (Amat && Amat == ts->Arhs) {
2503: /* 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 */
2504: MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat);
2505: SNESSetJacobian(snes, Amat, NULL, NULL, NULL);
2506: MatDestroy(&Amat);
2507: }
2508: if (Pmat && Pmat == ts->Brhs) {
2509: MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat);
2510: SNESSetJacobian(snes, NULL, Pmat, NULL, NULL);
2511: MatDestroy(&Pmat);
2512: }
2513: }
2515: TSGetAdapt(ts, &ts->adapt);
2516: TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type);
2518: PetscTryTypeMethod(ts, setup);
2520: TSSetExactFinalTimeDefault(ts);
2522: /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
2523: to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
2524: */
2525: TSGetDM(ts, &dm);
2526: DMSNESGetFunction(dm, &func, NULL);
2527: if (!func) DMSNESSetFunction(dm, SNESTSFormFunction, ts);
2529: /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
2530: Otherwise, the SNES will use coloring internally to form the Jacobian.
2531: */
2532: DMSNESGetJacobian(dm, &jac, NULL);
2533: DMTSGetIJacobian(dm, &ijac, NULL);
2534: DMTSGetI2Jacobian(dm, &i2jac, NULL);
2535: DMTSGetRHSJacobian(dm, &rhsjac, NULL);
2536: if (!jac && (ijac || i2jac || rhsjac)) DMSNESSetJacobian(dm, SNESTSFormJacobian, ts);
2538: /* if time integration scheme has a starting method, call it */
2539: PetscTryTypeMethod(ts, startingmethod);
2541: ts->setupcalled = PETSC_TRUE;
2542: return 0;
2543: }
2545: /*@
2546: TSReset - Resets a `TS` context and removes any allocated `Vec`s and `Mat`s.
2548: Collective
2550: Input Parameter:
2551: . ts - the `TS` context obtained from `TSCreate()`
2553: Level: beginner
2555: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()`
2556: @*/
2557: PetscErrorCode TSReset(TS ts)
2558: {
2559: TS_RHSSplitLink ilink = ts->tsrhssplit, next;
2563: PetscTryTypeMethod(ts, reset);
2564: if (ts->snes) SNESReset(ts->snes);
2565: if (ts->adapt) TSAdaptReset(ts->adapt);
2567: MatDestroy(&ts->Arhs);
2568: MatDestroy(&ts->Brhs);
2569: VecDestroy(&ts->Frhs);
2570: VecDestroy(&ts->vec_sol);
2571: VecDestroy(&ts->vec_dot);
2572: VecDestroy(&ts->vatol);
2573: VecDestroy(&ts->vrtol);
2574: VecDestroyVecs(ts->nwork, &ts->work);
2576: MatDestroy(&ts->Jacprhs);
2577: MatDestroy(&ts->Jacp);
2578: if (ts->forward_solve) TSForwardReset(ts);
2579: if (ts->quadraturets) {
2580: TSReset(ts->quadraturets);
2581: VecDestroy(&ts->vec_costintegrand);
2582: }
2583: while (ilink) {
2584: next = ilink->next;
2585: TSDestroy(&ilink->ts);
2586: PetscFree(ilink->splitname);
2587: ISDestroy(&ilink->is);
2588: PetscFree(ilink);
2589: ilink = next;
2590: }
2591: ts->tsrhssplit = NULL;
2592: ts->num_rhs_splits = 0;
2593: if (ts->tspan) {
2594: PetscFree(ts->tspan->span_times);
2595: VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol);
2596: PetscFree(ts->tspan);
2597: }
2598: ts->setupcalled = PETSC_FALSE;
2599: return 0;
2600: }
2602: /*@C
2603: TSDestroy - Destroys the timestepper context that was created
2604: with `TSCreate()`.
2606: Collective
2608: Input Parameter:
2609: . ts - the `TS` context obtained from `TSCreate()`
2611: Level: beginner
2613: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2614: @*/
2615: PetscErrorCode TSDestroy(TS *ts)
2616: {
2617: if (!*ts) return 0;
2619: if (--((PetscObject)(*ts))->refct > 0) {
2620: *ts = NULL;
2621: return 0;
2622: }
2624: TSReset(*ts);
2625: TSAdjointReset(*ts);
2626: if ((*ts)->forward_solve) TSForwardReset(*ts);
2628: /* if memory was published with SAWs then destroy it */
2629: PetscObjectSAWsViewOff((PetscObject)*ts);
2630: PetscTryTypeMethod((*ts), destroy);
2632: TSTrajectoryDestroy(&(*ts)->trajectory);
2634: TSAdaptDestroy(&(*ts)->adapt);
2635: TSEventDestroy(&(*ts)->event);
2637: SNESDestroy(&(*ts)->snes);
2638: DMDestroy(&(*ts)->dm);
2639: TSMonitorCancel((*ts));
2640: TSAdjointMonitorCancel((*ts));
2642: TSDestroy(&(*ts)->quadraturets);
2643: PetscHeaderDestroy(ts);
2644: return 0;
2645: }
2647: /*@
2648: TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
2649: a `TS` (timestepper) context. Valid only for nonlinear problems.
2651: Not Collective, but snes is parallel if ts is parallel
2653: Input Parameter:
2654: . ts - the `TS` context obtained from `TSCreate()`
2656: Output Parameter:
2657: . snes - the nonlinear solver context
2659: Level: beginner
2661: Notes:
2662: The user can then directly manipulate the `SNES` context to set various
2663: options, etc. Likewise, the user can then extract and manipulate the
2664: `KSP`, and `PC` contexts as well.
2666: `TSGetSNES()` does not work for integrators that do not use `SNES`; in
2667: this case `TSGetSNES()` returns NULL in snes.
2669: .seealso: [](chapter_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
2670: @*/
2671: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
2672: {
2675: if (!ts->snes) {
2676: SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes);
2677: PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options);
2678: SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts);
2679: PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1);
2680: if (ts->dm) SNESSetDM(ts->snes, ts->dm);
2681: if (ts->problem_type == TS_LINEAR) SNESSetType(ts->snes, SNESKSPONLY);
2682: }
2683: *snes = ts->snes;
2684: return 0;
2685: }
2687: /*@
2688: TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the timestepping context
2690: Collective
2692: Input Parameters:
2693: + ts - the `TS` context obtained from `TSCreate()`
2694: - snes - the nonlinear solver context
2696: Level: developer
2698: Note:
2699: Most users should have the `TS` created by calling `TSGetSNES()`
2701: .seealso: [](chapter_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2702: @*/
2703: PetscErrorCode TSSetSNES(TS ts, SNES snes)
2704: {
2705: PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);
2709: PetscObjectReference((PetscObject)snes);
2710: SNESDestroy(&ts->snes);
2712: ts->snes = snes;
2714: SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts);
2715: SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL);
2716: if (func == SNESTSFormJacobian) SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts);
2717: return 0;
2718: }
2720: /*@
2721: TSGetKSP - Returns the `KSP` (linear solver) associated with
2722: a `TS` (timestepper) context.
2724: Not Collective, but ksp is parallel if ts is parallel
2726: Input Parameter:
2727: . ts - the `TS` context obtained from `TSCreate()`
2729: Output Parameter:
2730: . ksp - the nonlinear solver context
2732: Level: beginner
2734: Notes:
2735: The user can then directly manipulate the `KSP` context to set various
2736: options, etc. Likewise, the user can then extract and manipulate the
2737: `PC` context as well.
2739: `TSGetKSP()` does not work for integrators that do not use `KSP`;
2740: in this case `TSGetKSP()` returns NULL in ksp.
2742: .seealso: [](chapter_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
2743: @*/
2744: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
2745: {
2746: SNES snes;
2752: TSGetSNES(ts, &snes);
2753: SNESGetKSP(snes, ksp);
2754: return 0;
2755: }
2757: /* ----------- Routines to set solver parameters ---------- */
2759: /*@
2760: TSSetMaxSteps - Sets the maximum number of steps to use.
2762: Logically Collective
2764: Input Parameters:
2765: + ts - the `TS` context obtained from `TSCreate()`
2766: - maxsteps - maximum number of steps to use
2768: Options Database Key:
2769: . -ts_max_steps <maxsteps> - Sets maxsteps
2771: Level: intermediate
2773: Note:
2774: The default maximum number of steps is 5000
2776: .seealso: [](chapter_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
2777: @*/
2778: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
2779: {
2783: ts->max_steps = maxsteps;
2784: return 0;
2785: }
2787: /*@
2788: TSGetMaxSteps - Gets the maximum number of steps to use.
2790: Not Collective
2792: Input Parameters:
2793: . ts - the `TS` context obtained from `TSCreate()`
2795: Output Parameter:
2796: . maxsteps - maximum number of steps to use
2798: Level: advanced
2800: .seealso: [](chapter_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
2801: @*/
2802: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
2803: {
2806: *maxsteps = ts->max_steps;
2807: return 0;
2808: }
2810: /*@
2811: TSSetMaxTime - Sets the maximum (or final) time for timestepping.
2813: Logically Collective
2815: Input Parameters:
2816: + ts - the `TS` context obtained from `TSCreate()`
2817: - maxtime - final time to step to
2819: Options Database Key:
2820: . -ts_max_time <maxtime> - Sets maxtime
2822: Level: intermediate
2824: Notes:
2825: The default maximum time is 5.0
2827: .seealso: [](chapter_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
2828: @*/
2829: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
2830: {
2833: ts->max_time = maxtime;
2834: return 0;
2835: }
2837: /*@
2838: TSGetMaxTime - Gets the maximum (or final) time for timestepping.
2840: Not Collective
2842: Input Parameters:
2843: . ts - the `TS` context obtained from `TSCreate()`
2845: Output Parameter:
2846: . maxtime - final time to step to
2848: Level: advanced
2850: .seealso: [](chapter_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
2851: @*/
2852: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
2853: {
2856: *maxtime = ts->max_time;
2857: return 0;
2858: }
2860: /*@
2861: TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.
2863: Level: deprecated
2865: @*/
2866: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
2867: {
2869: TSSetTime(ts, initial_time);
2870: TSSetTimeStep(ts, time_step);
2871: return 0;
2872: }
2874: /*@
2875: TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.
2877: Level: deprecated
2879: @*/
2880: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
2881: {
2883: if (maxsteps) {
2885: *maxsteps = ts->max_steps;
2886: }
2887: if (maxtime) {
2889: *maxtime = ts->max_time;
2890: }
2891: return 0;
2892: }
2894: /*@
2895: TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.
2897: Level: deprecated
2899: @*/
2900: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
2901: {
2905: if (maxsteps >= 0) ts->max_steps = maxsteps;
2906: if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime;
2907: return 0;
2908: }
2910: /*@
2911: TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.
2913: Level: deprecated
2915: @*/
2916: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
2917: {
2918: return TSGetStepNumber(ts, steps);
2919: }
2921: /*@
2922: TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.
2924: Level: deprecated
2926: @*/
2927: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
2928: {
2929: return TSGetStepNumber(ts, steps);
2930: }
2932: /*@
2933: TSSetSolution - Sets the initial solution vector
2934: for use by the `TS` routines.
2936: Logically Collective
2938: Input Parameters:
2939: + ts - the `TS` context obtained from `TSCreate()`
2940: - u - the solution vector
2942: Level: beginner
2944: .seealso: [](chapter_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
2945: @*/
2946: PetscErrorCode TSSetSolution(TS ts, Vec u)
2947: {
2948: DM dm;
2952: PetscObjectReference((PetscObject)u);
2953: VecDestroy(&ts->vec_sol);
2954: ts->vec_sol = u;
2956: TSGetDM(ts, &dm);
2957: DMShellSetGlobalVector(dm, u);
2958: return 0;
2959: }
2961: /*@C
2962: TSSetPreStep - Sets the general-purpose function
2963: called once at the beginning of each time step.
2965: Logically Collective
2967: Input Parameters:
2968: + ts - The `TS` context obtained from `TSCreate()`
2969: - func - The function
2971: Calling sequence of func:
2972: .vb
2973: PetscErrorCode func (TS ts);
2974: .ve
2976: Level: intermediate
2978: .seealso: [](chapter_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
2979: @*/
2980: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS))
2981: {
2983: ts->prestep = func;
2984: return 0;
2985: }
2987: /*@
2988: TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`
2990: Collective
2992: Input Parameters:
2993: . ts - The `TS` context obtained from `TSCreate()`
2995: Level: developer
2997: Note:
2998: `TSPreStep()` is typically used within time stepping implementations,
2999: so most users would not generally call this routine themselves.
3001: .seealso: [](chapter_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
3002: @*/
3003: PetscErrorCode TSPreStep(TS ts)
3004: {
3006: if (ts->prestep) {
3007: Vec U;
3008: PetscObjectId idprev;
3009: PetscBool sameObject;
3010: PetscObjectState sprev, spost;
3012: TSGetSolution(ts, &U);
3013: PetscObjectGetId((PetscObject)U, &idprev);
3014: PetscObjectStateGet((PetscObject)U, &sprev);
3015: PetscCallBack("TS callback preset", (*ts->prestep)(ts));
3016: TSGetSolution(ts, &U);
3017: PetscObjectCompareId((PetscObject)U, idprev, &sameObject);
3018: PetscObjectStateGet((PetscObject)U, &spost);
3019: if (!sameObject || sprev != spost) TSRestartStep(ts);
3020: }
3021: return 0;
3022: }
3024: /*@C
3025: TSSetPreStage - Sets the general-purpose function
3026: called once at the beginning of each stage.
3028: Logically Collective
3030: Input Parameters:
3031: + ts - The `TS` context obtained from `TSCreate()`
3032: - func - The function
3034: Calling sequence of func:
3035: .vb
3036: PetscErrorCode func(TS ts, PetscReal stagetime);
3037: .ve
3039: Level: intermediate
3041: Note:
3042: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3043: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3044: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3046: .seealso: [](chapter_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3047: @*/
3048: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS, PetscReal))
3049: {
3051: ts->prestage = func;
3052: return 0;
3053: }
3055: /*@C
3056: TSSetPostStage - Sets the general-purpose function, provided with `TSSetPostStep()`,
3057: called once at the end of each stage.
3059: Logically Collective
3061: Input Parameters:
3062: + ts - The `TS` context obtained from `TSCreate()`
3063: - func - The function
3065: Calling sequence of func:
3066: .vb
3067: PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y);
3068: .ve
3070: Level: intermediate
3072: Note:
3073: There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
3074: The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
3075: attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.
3077: .seealso: [](chapter_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3078: @*/
3079: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS, PetscReal, PetscInt, Vec *))
3080: {
3082: ts->poststage = func;
3083: return 0;
3084: }
3086: /*@C
3087: TSSetPostEvaluate - Sets the general-purpose function
3088: called once at the end of each step evaluation.
3090: Logically Collective
3092: Input Parameters:
3093: + ts - The `TS` context obtained from `TSCreate()`
3094: - func - The function
3096: Calling sequence of func:
3097: .vb
3098: PetscErrorCode func(TS ts);
3099: .ve
3101: Level: intermediate
3103: Note:
3104: Semantically, `TSSetPostEvaluate()` differs from `TSSetPostStep()` since the function it sets is called before event-handling
3105: thus guaranteeing the same solution (computed by the time-stepper) will be passed to it. On the other hand, `TSPostStep()`
3106: may be passed a different solution, possibly changed by the event handler. `TSPostEvaluate()` is called after the next step
3107: solution is evaluated allowing to modify it, if need be. The solution can be obtained with `TSGetSolution()`, the time step
3108: with `TSGetTimeStep()`, and the time at the start of the step is available via `TSGetTime()`
3110: .seealso: [](chapter_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
3111: @*/
3112: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS))
3113: {
3115: ts->postevaluate = func;
3116: return 0;
3117: }
3119: /*@
3120: TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`
3122: Collective
3124: Input Parameters:
3125: . ts - The `TS` context obtained from `TSCreate()`
3126: stagetime - The absolute time of the current stage
3128: Level: developer
3130: Note:
3131: `TSPreStage()` is typically used within time stepping implementations,
3132: most users would not generally call this routine themselves.
3134: .seealso: [](chapter_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3135: @*/
3136: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
3137: {
3139: if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
3140: return 0;
3141: }
3143: /*@
3144: TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`
3146: Collective
3148: Input Parameters:
3149: . ts - The `TS` context obtained from `TSCreate()`
3150: stagetime - The absolute time of the current stage
3151: stageindex - Stage number
3152: Y - Array of vectors (of size = total number
3153: of stages) with the stage solutions
3155: Level: developer
3157: Note:
3158: `TSPostStage()` is typically used within time stepping implementations,
3159: most users would not generally call this routine themselves.
3161: .seealso: [](chapter_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3162: @*/
3163: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
3164: {
3166: if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
3167: return 0;
3168: }
3170: /*@
3171: TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`
3173: Collective
3175: Input Parameters:
3176: . ts - The `TS` context obtained from `TSCreate()`
3178: Level: developer
3180: Note:
3181: `TSPostEvaluate()` is typically used within time stepping implementations,
3182: most users would not generally call this routine themselves.
3184: .seealso: [](chapter_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
3185: @*/
3186: PetscErrorCode TSPostEvaluate(TS ts)
3187: {
3189: if (ts->postevaluate) {
3190: Vec U;
3191: PetscObjectState sprev, spost;
3193: TSGetSolution(ts, &U);
3194: PetscObjectStateGet((PetscObject)U, &sprev);
3195: PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
3196: PetscObjectStateGet((PetscObject)U, &spost);
3197: if (sprev != spost) TSRestartStep(ts);
3198: }
3199: return 0;
3200: }
3202: /*@C
3203: TSSetPostStep - Sets the general-purpose function
3204: called once at the end of each time step.
3206: Logically Collective
3208: Input Parameters:
3209: + ts - The `TS` context obtained from `TSCreate()`
3210: - func - The function
3212: Calling sequence of func:
3213: $ func (TS ts);
3215: Level: intermediate
3217: Note:
3218: The function set by `TSSetPostStep()` is called after each successful step. The solution vector X
3219: obtained by `TSGetSolution()` may be different than that computed at the step end if the event handler
3220: locates an event and `TSPostEvent()` modifies it. Use `TSSetPostEvaluate()` if an unmodified solution is needed instead.
3222: .seealso: [](chapter_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
3223: @*/
3224: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS))
3225: {
3227: ts->poststep = func;
3228: return 0;
3229: }
3231: /*@
3232: TSPostStep - Runs the user-defined post-step function that was set with `TSSetPotsStep()`
3234: Collective
3236: Input Parameters:
3237: . ts - The `TS` context obtained from `TSCreate()`
3239: Note:
3240: `TSPostStep()` is typically used within time stepping implementations,
3241: so most users would not generally call this routine themselves.
3243: Level: developer
3245: .seealso: [](chapter_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()`
3246: @*/
3247: PetscErrorCode TSPostStep(TS ts)
3248: {
3250: if (ts->poststep) {
3251: Vec U;
3252: PetscObjectId idprev;
3253: PetscBool sameObject;
3254: PetscObjectState sprev, spost;
3256: TSGetSolution(ts, &U);
3257: PetscObjectGetId((PetscObject)U, &idprev);
3258: PetscObjectStateGet((PetscObject)U, &sprev);
3259: PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
3260: TSGetSolution(ts, &U);
3261: PetscObjectCompareId((PetscObject)U, idprev, &sameObject);
3262: PetscObjectStateGet((PetscObject)U, &spost);
3263: if (!sameObject || sprev != spost) TSRestartStep(ts);
3264: }
3265: return 0;
3266: }
3268: /*@
3269: TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval
3271: Collective
3273: Input Parameters:
3274: + ts - time stepping context
3275: - t - time to interpolate to
3277: Output Parameter:
3278: . U - state at given time
3280: Level: intermediate
3282: Developer Note:
3283: `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.
3285: .seealso: [](chapter_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
3286: @*/
3287: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
3288: {
3292: PetscUseTypeMethod(ts, interpolate, t, U);
3293: return 0;
3294: }
3296: /*@
3297: TSStep - Steps one time step
3299: Collective
3301: Input Parameter:
3302: . ts - the `TS` context obtained from `TSCreate()`
3304: Level: developer
3306: Notes:
3307: The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.
3309: The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
3310: be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.
3312: This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
3313: time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.
3315: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
3316: @*/
3317: PetscErrorCode TSStep(TS ts)
3318: {
3319: static PetscBool cite = PETSC_FALSE;
3320: PetscReal ptime;
3323: PetscCall(PetscCitationsRegister("@article{tspaper,\n"
3324: " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
3325: " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
3326: " journal = {arXiv e-preprints},\n"
3327: " eprint = {1806.01437},\n"
3328: " archivePrefix = {arXiv},\n"
3329: " year = {2018}\n}\n",
3330: &cite));
3331: TSSetUp(ts);
3332: TSTrajectorySetUp(ts->trajectory, ts);
3338: if (!ts->steps) ts->ptime_prev = ts->ptime;
3339: ptime = ts->ptime;
3340: ts->ptime_prev_rollback = ts->ptime_prev;
3341: ts->reason = TS_CONVERGED_ITERATING;
3343: PetscLogEventBegin(TS_Step, ts, 0, 0, 0);
3344: PetscUseTypeMethod(ts, step);
3345: PetscLogEventEnd(TS_Step, ts, 0, 0, 0);
3347: if (ts->tspan && PetscIsCloseAtTol(ts->ptime, ts->tspan->span_times[ts->tspan->spanctr], ts->tspan->reltol * ts->time_step + ts->tspan->abstol, 0) && ts->tspan->spanctr < ts->tspan->num_span_times)
3348: VecCopy(ts->vec_sol, ts->tspan->vecs_sol[ts->tspan->spanctr++]);
3349: if (ts->reason >= 0) {
3350: ts->ptime_prev = ptime;
3351: ts->steps++;
3352: ts->steprollback = PETSC_FALSE;
3353: ts->steprestart = PETSC_FALSE;
3354: }
3355: if (!ts->reason) {
3356: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3357: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3358: }
3362: return 0;
3363: }
3365: /*@
3366: TSEvaluateWLTE - Evaluate the weighted local truncation error norm
3367: at the end of a time step with a given order of accuracy.
3369: Collective
3371: Input Parameters:
3372: + ts - time stepping context
3373: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
3375: Input/Output Parameter:
3376: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
3377: on output, the actual order of the error evaluation
3379: Output Parameter:
3380: . wlte - the weighted local truncation error norm
3382: Level: advanced
3384: Note:
3385: If the timestepper cannot evaluate the error in a particular step
3386: (eg. in the first step or restart steps after event handling),
3387: this routine returns wlte=-1.0 .
3389: .seealso: [](chapter_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
3390: @*/
3391: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
3392: {
3400: PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
3401: return 0;
3402: }
3404: /*@
3405: TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.
3407: Collective
3409: Input Parameters:
3410: + ts - time stepping context
3411: . order - desired order of accuracy
3412: - done - whether the step was evaluated at this order (pass NULL to generate an error if not available)
3414: Output Parameter:
3415: . U - state at the end of the current step
3417: Level: advanced
3419: Notes:
3420: This function cannot be called until all stages have been evaluated.
3422: 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.
3424: .seealso: [](chapter_ts), `TS`, `TSStep()`, `TSAdapt`
3425: @*/
3426: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
3427: {
3431: PetscUseTypeMethod(ts, evaluatestep, order, U, done);
3432: return 0;
3433: }
3435: /*@C
3436: TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.
3438: Not collective
3440: Input Parameter:
3441: . ts - time stepping context
3443: Output Parameter:
3444: . initConditions - The function which computes an initial condition
3446: The calling sequence for the function is
3447: .vb
3448: initCondition(TS ts, Vec u)
3449: ts - The timestepping context
3450: u - The input vector in which the initial condition is stored
3451: .ve
3453: Level: advanced
3455: .seealso: [](chapter_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
3456: @*/
3457: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS, Vec))
3458: {
3461: *initCondition = ts->ops->initcondition;
3462: return 0;
3463: }
3465: /*@C
3466: TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.
3468: Logically collective
3470: Input Parameters:
3471: + ts - time stepping context
3472: - initCondition - The function which computes an initial condition
3474: Calling sequence for initCondition:
3475: $ PetscErrorCode initCondition(TS ts, Vec u)
3476: + ts - The timestepping context
3477: - u - The input vector in which the initial condition is to be stored
3479: Level: advanced
3481: .seealso: [](chapter_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
3482: @*/
3483: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS, Vec))
3484: {
3487: ts->ops->initcondition = initCondition;
3488: return 0;
3489: }
3491: /*@
3492: TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`
3494: Collective
3496: Input Parameters:
3497: + ts - time stepping context
3498: - u - The `Vec` to store the condition in which will be used in `TSSolve()`
3500: Level: advanced
3502: .seealso: [](chapter_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3503: @*/
3504: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
3505: {
3508: PetscTryTypeMethod(ts, initcondition, u);
3509: return 0;
3510: }
3512: /*@C
3513: TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.
3515: Not collective
3517: Input Parameter:
3518: . ts - time stepping context
3520: Output Parameter:
3521: . exactError - The function which computes the solution error
3523: Calling sequence for exactError:
3524: $ PetscErrorCode exactError(TS ts, Vec u)
3525: + ts - The timestepping context
3526: . u - The approximate solution vector
3527: - e - The input vector in which the error is stored
3529: Level: advanced
3531: .seealso: [](chapter_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3532: @*/
3533: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS, Vec, Vec))
3534: {
3537: *exactError = ts->ops->exacterror;
3538: return 0;
3539: }
3541: /*@C
3542: TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.
3544: Logically collective
3546: Input Parameters:
3547: + ts - time stepping context
3548: - exactError - The function which computes the solution error
3550: Calling sequence for exactError:
3551: $ PetscErrorCode exactError(TS ts, Vec u)
3552: + ts - The timestepping context
3553: . u - The approximate solution vector
3554: - e - The input vector in which the error is stored
3556: Level: advanced
3558: .seealso: [](chapter_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
3559: @*/
3560: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS, Vec, Vec))
3561: {
3564: ts->ops->exacterror = exactError;
3565: return 0;
3566: }
3568: /*@
3569: TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`
3571: Collective
3573: Input Parameters:
3574: + ts - time stepping context
3575: . u - The approximate solution
3576: - e - The `Vec` used to store the error
3578: Level: advanced
3580: .seealso: [](chapter_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
3581: @*/
3582: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
3583: {
3587: PetscTryTypeMethod(ts, exacterror, u, e);
3588: return 0;
3589: }
3591: /*@
3592: TSSolve - Steps the requested number of timesteps.
3594: Collective
3596: Input Parameters:
3597: + ts - the `TS` context obtained from `TSCreate()`
3598: - u - the solution vector (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
3599: otherwise must contain the initial conditions and will contain the solution at the final requested time
3601: Level: beginner
3603: Notes:
3604: The final time returned by this function may be different from the time of the internally
3605: held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
3606: stepped over the final time.
3608: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
3609: @*/
3610: PetscErrorCode TSSolve(TS ts, Vec u)
3611: {
3612: Vec solution;
3617: TSSetExactFinalTimeDefault(ts);
3618: 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 */
3619: if (!ts->vec_sol || u == ts->vec_sol) {
3620: VecDuplicate(u, &solution);
3621: TSSetSolution(ts, solution);
3622: VecDestroy(&solution); /* grant ownership */
3623: }
3624: VecCopy(u, ts->vec_sol);
3626: } else if (u) TSSetSolution(ts, u);
3627: TSSetUp(ts);
3628: TSTrajectorySetUp(ts->trajectory, ts);
3635: if (ts->tspan && PetscIsCloseAtTol(ts->ptime, ts->tspan->span_times[0], ts->tspan->reltol * ts->time_step + ts->tspan->abstol, 0)) { /* starting point in time span */
3636: VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]);
3637: ts->tspan->spanctr = 1;
3638: }
3640: if (ts->forward_solve) TSForwardSetUp(ts);
3642: /* reset number of steps only when the step is not restarted. ARKIMEX
3643: restarts the step after an event. Resetting these counters in such case causes
3644: TSTrajectory to incorrectly save the output files
3645: */
3646: /* reset time step and iteration counters */
3647: if (!ts->steps) {
3648: ts->ksp_its = 0;
3649: ts->snes_its = 0;
3650: ts->num_snes_failures = 0;
3651: ts->reject = 0;
3652: ts->steprestart = PETSC_TRUE;
3653: ts->steprollback = PETSC_FALSE;
3654: ts->rhsjacobian.time = PETSC_MIN_REAL;
3655: }
3657: /* make sure initial time step does not overshoot final time or the next point in tspan */
3658: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
3659: PetscReal maxdt;
3660: PetscReal dt = ts->time_step;
3662: if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
3663: else maxdt = ts->max_time - ts->ptime;
3664: ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
3665: }
3666: ts->reason = TS_CONVERGED_ITERATING;
3668: {
3669: PetscViewer viewer;
3670: PetscViewerFormat format;
3671: PetscBool flg;
3672: static PetscBool incall = PETSC_FALSE;
3674: if (!incall) {
3675: /* Estimate the convergence rate of the time discretization */
3676: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg);
3677: if (flg) {
3678: PetscConvEst conv;
3679: DM dm;
3680: PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
3681: PetscInt Nf;
3682: PetscBool checkTemporal = PETSC_TRUE;
3684: incall = PETSC_TRUE;
3685: PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg);
3686: TSGetDM(ts, &dm);
3687: DMGetNumFields(dm, &Nf);
3688: PetscCalloc1(PetscMax(Nf, 1), &alpha);
3689: PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv);
3690: PetscConvEstUseTS(conv, checkTemporal);
3691: PetscConvEstSetSolver(conv, (PetscObject)ts);
3692: PetscConvEstSetFromOptions(conv);
3693: PetscConvEstSetUp(conv);
3694: PetscConvEstGetConvRate(conv, alpha);
3695: PetscViewerPushFormat(viewer, format);
3696: PetscConvEstRateView(conv, alpha, viewer);
3697: PetscViewerPopFormat(viewer);
3698: PetscViewerDestroy(&viewer);
3699: PetscConvEstDestroy(&conv);
3700: PetscFree(alpha);
3701: incall = PETSC_FALSE;
3702: }
3703: }
3704: }
3706: TSViewFromOptions(ts, NULL, "-ts_view_pre");
3708: if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
3709: PetscUseTypeMethod(ts, solve);
3710: if (u) VecCopy(ts->vec_sol, u);
3711: ts->solvetime = ts->ptime;
3712: solution = ts->vec_sol;
3713: } else { /* Step the requested number of timesteps. */
3714: if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
3715: else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
3717: if (!ts->steps) {
3718: TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol);
3719: TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol);
3720: }
3722: while (!ts->reason) {
3723: TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);
3724: if (!ts->steprollback) TSPreStep(ts);
3725: TSStep(ts);
3726: if (ts->testjacobian) TSRHSJacobianTest(ts, NULL);
3727: if (ts->testjacobiantranspose) TSRHSJacobianTestTranspose(ts, NULL);
3728: if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
3729: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
3730: TSForwardCostIntegral(ts);
3731: if (ts->reason >= 0) ts->steps++;
3732: }
3733: if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
3734: if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
3735: TSForwardStep(ts);
3736: if (ts->reason >= 0) ts->steps++;
3737: }
3738: TSPostEvaluate(ts);
3739: TSEventHandler(ts); /* The right-hand side may be changed due to event. Be careful with Any computation using the RHS information after this point. */
3740: if (ts->steprollback) TSPostEvaluate(ts);
3741: if (!ts->steprollback) {
3742: TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol);
3743: TSPostStep(ts);
3744: }
3745: }
3746: TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);
3748: if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
3749: TSInterpolate(ts, ts->max_time, u);
3750: ts->solvetime = ts->max_time;
3751: solution = u;
3752: TSMonitor(ts, -1, ts->solvetime, solution);
3753: } else {
3754: if (u) VecCopy(ts->vec_sol, u);
3755: ts->solvetime = ts->ptime;
3756: solution = ts->vec_sol;
3757: }
3758: }
3760: TSViewFromOptions(ts, NULL, "-ts_view");
3761: VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution");
3762: PetscObjectSAWsBlock((PetscObject)ts);
3763: if (ts->adjoint_solve) TSAdjointSolve(ts);
3764: return 0;
3765: }
3767: /*@
3768: TSGetTime - Gets the time of the most recently completed step.
3770: Not Collective
3772: Input Parameter:
3773: . ts - the `TS` context obtained from `TSCreate()`
3775: Output Parameter:
3776: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.
3778: Level: beginner
3780: Note:
3781: When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
3782: `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.
3784: .seealso: [](chapter_ts), TS`, ``TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
3785: @*/
3786: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
3787: {
3790: *t = ts->ptime;
3791: return 0;
3792: }
3794: /*@
3795: TSGetPrevTime - Gets the starting time of the previously completed step.
3797: Not Collective
3799: Input Parameter:
3800: . ts - the `TS` context obtained from `TSCreate()`
3802: Output Parameter:
3803: . t - the previous time
3805: Level: beginner
3807: .seealso: [](chapter_ts), TS`, ``TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
3808: @*/
3809: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
3810: {
3813: *t = ts->ptime_prev;
3814: return 0;
3815: }
3817: /*@
3818: TSSetTime - Allows one to reset the time.
3820: Logically Collective
3822: Input Parameters:
3823: + ts - the `TS` context obtained from `TSCreate()`
3824: - time - the time
3826: Level: intermediate
3828: .seealso: [](chapter_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
3829: @*/
3830: PetscErrorCode TSSetTime(TS ts, PetscReal t)
3831: {
3834: ts->ptime = t;
3835: return 0;
3836: }
3838: /*@C
3839: TSSetOptionsPrefix - Sets the prefix used for searching for all
3840: TS options in the database.
3842: Logically Collective
3844: Input Parameters:
3845: + ts - The `TS` context
3846: - prefix - The prefix to prepend to all option names
3848: Level: advanced
3850: Note:
3851: A hyphen (-) must NOT be given at the beginning of the prefix name.
3852: The first character of all runtime options is AUTOMATICALLY the
3853: hyphen.
3855: .seealso: [](chapter_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
3856: @*/
3857: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
3858: {
3859: SNES snes;
3862: PetscObjectSetOptionsPrefix((PetscObject)ts, prefix);
3863: TSGetSNES(ts, &snes);
3864: SNESSetOptionsPrefix(snes, prefix);
3865: return 0;
3866: }
3868: /*@C
3869: TSAppendOptionsPrefix - Appends to the prefix used for searching for all
3870: TS options in the database.
3872: Logically Collective
3874: Input Parameters:
3875: + ts - The `TS` context
3876: - prefix - The prefix to prepend to all option names
3878: Level: advanced
3880: Note:
3881: A hyphen (-) must NOT be given at the beginning of the prefix name.
3882: The first character of all runtime options is AUTOMATICALLY the
3883: hyphen.
3885: .seealso: [](chapter_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
3886: @*/
3887: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
3888: {
3889: SNES snes;
3892: PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix);
3893: TSGetSNES(ts, &snes);
3894: SNESAppendOptionsPrefix(snes, prefix);
3895: return 0;
3896: }
3898: /*@C
3899: TSGetOptionsPrefix - Sets the prefix used for searching for all
3900: `TS` options in the database.
3902: Not Collective
3904: Input Parameter:
3905: . ts - The `TS` context
3907: Output Parameter:
3908: . prefix - A pointer to the prefix string used
3910: Level: intermediate
3912: Fortran Note:
3913: On the fortran side, the user should pass in a string 'prefix' of
3914: sufficient length to hold the prefix.
3916: .seealso: [](chapter_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
3917: @*/
3918: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
3919: {
3922: PetscObjectGetOptionsPrefix((PetscObject)ts, prefix);
3923: return 0;
3924: }
3926: /*@C
3927: TSGetRHSJacobian - Returns the Jacobian J at the present timestep.
3929: Not Collective, but parallel objects are returned if ts is parallel
3931: Input Parameter:
3932: . ts - The `TS` context obtained from `TSCreate()`
3934: Output Parameters:
3935: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL)
3936: . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL)
3937: . func - Function to compute the Jacobian of the RHS (or NULL)
3938: - ctx - User-defined context for Jacobian evaluation routine (or NULL)
3940: Level: intermediate
3942: Note:
3943: You can pass in NULL for any return argument you do not need.
3945: .seealso: [](chapter_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
3947: @*/
3948: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobian *func, void **ctx)
3949: {
3950: DM dm;
3952: if (Amat || Pmat) {
3953: SNES snes;
3954: TSGetSNES(ts, &snes);
3955: SNESSetUpMatrices(snes);
3956: SNESGetJacobian(snes, Amat, Pmat, NULL, NULL);
3957: }
3958: TSGetDM(ts, &dm);
3959: DMTSGetRHSJacobian(dm, func, ctx);
3960: return 0;
3961: }
3963: /*@C
3964: TSGetIJacobian - Returns the implicit Jacobian at the present timestep.
3966: Not Collective, but parallel objects are returned if ts is parallel
3968: Input Parameter:
3969: . ts - The `TS` context obtained from `TSCreate()`
3971: Output Parameters:
3972: + Amat - The (approximate) Jacobian of F(t,U,U_t)
3973: . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat
3974: . f - The function to compute the matrices
3975: - ctx - User-defined context for Jacobian evaluation routine
3977: Level: advanced
3979: Note:
3980: You can pass in NULL for any return argument you do not need.
3982: .seealso: [](chapter_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
3984: @*/
3985: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobian *f, void **ctx)
3986: {
3987: DM dm;
3989: if (Amat || Pmat) {
3990: SNES snes;
3991: TSGetSNES(ts, &snes);
3992: SNESSetUpMatrices(snes);
3993: SNESGetJacobian(snes, Amat, Pmat, NULL, NULL);
3994: }
3995: TSGetDM(ts, &dm);
3996: DMTSGetIJacobian(dm, f, ctx);
3997: return 0;
3998: }
4000: #include <petsc/private/dmimpl.h>
4001: /*@
4002: TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`
4004: Logically Collective
4006: Input Parameters:
4007: + ts - the `TS` integrator object
4008: - dm - the dm, cannot be NULL
4010: Level: intermediate
4012: Notes:
4013: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
4014: even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving
4015: different problems using the same function space.
4017: .seealso: [](chapter_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
4018: @*/
4019: PetscErrorCode TSSetDM(TS ts, DM dm)
4020: {
4021: SNES snes;
4022: DMTS tsdm;
4026: PetscObjectReference((PetscObject)dm);
4027: if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
4028: if (ts->dm->dmts && !dm->dmts) {
4029: DMCopyDMTS(ts->dm, dm);
4030: DMGetDMTS(ts->dm, &tsdm);
4031: /* Grant write privileges to the replacement DM */
4032: if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
4033: }
4034: DMDestroy(&ts->dm);
4035: }
4036: ts->dm = dm;
4038: TSGetSNES(ts, &snes);
4039: SNESSetDM(snes, dm);
4040: return 0;
4041: }
4043: /*@
4044: TSGetDM - Gets the `DM` that may be used by some preconditioners
4046: Not Collective
4048: Input Parameter:
4049: . ts - the `TS`
4051: Output Parameter:
4052: . dm - the dm
4054: Level: intermediate
4056: .seealso: [](chapter_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
4057: @*/
4058: PetscErrorCode TSGetDM(TS ts, DM *dm)
4059: {
4061: if (!ts->dm) {
4062: DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm);
4063: if (ts->snes) SNESSetDM(ts->snes, ts->dm);
4064: }
4065: *dm = ts->dm;
4066: return 0;
4067: }
4069: /*@
4070: SNESTSFormFunction - Function to evaluate nonlinear residual
4072: Logically Collective
4074: Input Parameters:
4075: + snes - nonlinear solver
4076: . U - the current state at which to evaluate the residual
4077: - ctx - user context, must be a TS
4079: Output Parameter:
4080: . F - the nonlinear residual
4082: Level: advanced
4084: Note:
4085: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4086: It is most frequently passed to `MatFDColoringSetFunction()`.
4088: .seealso: [](chapter_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
4089: @*/
4090: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
4091: {
4092: TS ts = (TS)ctx;
4098: (ts->ops->snesfunction)(snes, U, F, ts);
4099: return 0;
4100: }
4102: /*@
4103: SNESTSFormJacobian - Function to evaluate the Jacobian
4105: Collective
4107: Input Parameters:
4108: + snes - nonlinear solver
4109: . U - the current state at which to evaluate the residual
4110: - ctx - user context, must be a `TS`
4112: Output Parameters:
4113: + A - the Jacobian
4114: - B - the preconditioning matrix (may be the same as A)
4116: Level: developer
4118: Note:
4119: This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
4121: .seealso: [](chapter_ts), `SNESSetJacobian()`
4122: @*/
4123: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
4124: {
4125: TS ts = (TS)ctx;
4134: (ts->ops->snesjacobian)(snes, U, A, B, ts);
4135: return 0;
4136: }
4138: /*@C
4139: TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems Udot = A U only
4141: Collective
4143: Input Parameters:
4144: + ts - time stepping context
4145: . t - time at which to evaluate
4146: . U - state at which to evaluate
4147: - ctx - context
4149: Output Parameter:
4150: . F - right hand side
4152: Level: intermediate
4154: Note:
4155: This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right hand side for linear problems.
4156: The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.
4158: .seealso: [](chapter_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
4159: @*/
4160: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
4161: {
4162: Mat Arhs, Brhs;
4164: TSGetRHSMats_Private(ts, &Arhs, &Brhs);
4165: /* undo the damage caused by shifting */
4166: TSRecoverRHSJacobian(ts, Arhs, Brhs);
4167: TSComputeRHSJacobian(ts, t, U, Arhs, Brhs);
4168: MatMult(Arhs, U, F);
4169: return 0;
4170: }
4172: /*@C
4173: TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.
4175: Collective
4177: Input Parameters:
4178: + ts - time stepping context
4179: . t - time at which to evaluate
4180: . U - state at which to evaluate
4181: - ctx - context
4183: Output Parameters:
4184: + A - pointer to operator
4185: - B - pointer to preconditioning matrix
4187: Level: intermediate
4189: Note:
4190: This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.
4192: .seealso: [](chapter_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
4193: @*/
4194: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
4195: {
4196: return 0;
4197: }
4199: /*@C
4200: TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only
4202: Collective
4204: Input Parameters:
4205: + ts - time stepping context
4206: . t - time at which to evaluate
4207: . U - state at which to evaluate
4208: . Udot - time derivative of state vector
4209: - ctx - context
4211: Output Parameter:
4212: . F - left hand side
4214: Level: intermediate
4216: Notes:
4217: 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
4218: user is required to write their own `TSComputeIFunction()`.
4219: This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
4220: The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.
4222: Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U
4224: .seealso: [](chapter_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
4225: @*/
4226: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
4227: {
4228: Mat A, B;
4230: TSGetIJacobian(ts, &A, &B, NULL, NULL);
4231: TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE);
4232: MatMult(A, Udot, F);
4233: return 0;
4234: }
4236: /*@C
4237: TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE
4239: Collective
4241: Input Parameters:
4242: + ts - time stepping context
4243: . t - time at which to evaluate
4244: . U - state at which to evaluate
4245: . Udot - time derivative of state vector
4246: . shift - shift to apply
4247: - ctx - context
4249: Output Parameters:
4250: + A - pointer to operator
4251: - B - pointer to preconditioning matrix
4253: Level: advanced
4255: Notes:
4256: This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.
4258: It is only appropriate for problems of the form
4260: $ M Udot = F(U,t)
4262: where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only
4263: works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
4264: an implicit operator of the form
4266: $ shift*M + J
4268: 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
4269: a copy of M or reassemble it when requested.
4271: .seealso: [](chapter_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
4272: @*/
4273: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
4274: {
4275: MatScale(A, shift / ts->ijacobian.shift);
4276: ts->ijacobian.shift = shift;
4277: return 0;
4278: }
4280: /*@
4281: TSGetEquationType - Gets the type of the equation that `TS` is solving.
4283: Not Collective
4285: Input Parameter:
4286: . ts - the `TS` context
4288: Output Parameter:
4289: . equation_type - see `TSEquationType`
4291: Level: beginner
4293: .seealso: [](chapter_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
4294: @*/
4295: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
4296: {
4299: *equation_type = ts->equation_type;
4300: return 0;
4301: }
4303: /*@
4304: TSSetEquationType - Sets the type of the equation that `TS` is solving.
4306: Not Collective
4308: Input Parameters:
4309: + ts - the `TS` context
4310: - equation_type - see `TSEquationType`
4312: Level: advanced
4314: .seealso: [](chapter_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
4315: @*/
4316: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
4317: {
4319: ts->equation_type = equation_type;
4320: return 0;
4321: }
4323: /*@
4324: TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.
4326: Not Collective
4328: Input Parameter:
4329: . ts - the `TS` context
4331: Output Parameter:
4332: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4333: manual pages for the individual convergence tests for complete lists
4335: Level: beginner
4337: Note:
4338: Can only be called after the call to `TSSolve()` is complete.
4340: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4341: @*/
4342: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
4343: {
4346: *reason = ts->reason;
4347: return 0;
4348: }
4350: /*@
4351: TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.
4353: Logically Collective; reason must contain common value
4355: Input Parameters:
4356: + ts - the `TS` context
4357: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
4358: manual pages for the individual convergence tests for complete lists
4360: Level: advanced
4362: Note:
4363: Can only be called while `TSSolve()` is active.
4365: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSConvergedReason`
4366: @*/
4367: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
4368: {
4370: ts->reason = reason;
4371: return 0;
4372: }
4374: /*@
4375: TSGetSolveTime - Gets the time after a call to `TSSolve()`
4377: Not Collective
4379: Input Parameter:
4380: . ts - the `TS` context
4382: Output Parameter:
4383: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`
4385: Level: beginner
4387: Note:
4388: Can only be called after the call to `TSSolve()` is complete.
4390: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason`
4391: @*/
4392: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
4393: {
4396: *ftime = ts->solvetime;
4397: return 0;
4398: }
4400: /*@
4401: TSGetSNESIterations - Gets the total number of nonlinear iterations
4402: used by the time integrator.
4404: Not Collective
4406: Input Parameter:
4407: . ts - `TS` context
4409: Output Parameter:
4410: . nits - number of nonlinear iterations
4412: Level: intermediate
4414: Notes:
4415: This counter is reset to zero for each successive call to `TSSolve()`.
4417: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
4418: @*/
4419: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
4420: {
4423: *nits = ts->snes_its;
4424: return 0;
4425: }
4427: /*@
4428: TSGetKSPIterations - Gets the total number of linear iterations
4429: used by the time integrator.
4431: Not Collective
4433: Input Parameter:
4434: . ts - `TS` context
4436: Output Parameter:
4437: . lits - number of linear iterations
4439: Level: intermediate
4441: Note:
4442: This counter is reset to zero for each successive call to `TSSolve()`.
4444: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()`
4445: @*/
4446: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
4447: {
4450: *lits = ts->ksp_its;
4451: return 0;
4452: }
4454: /*@
4455: TSGetStepRejections - Gets the total number of rejected steps.
4457: Not Collective
4459: Input Parameter:
4460: . ts - `TS` context
4462: Output Parameter:
4463: . rejects - number of steps rejected
4465: Level: intermediate
4467: Note:
4468: This counter is reset to zero for each successive call to `TSSolve()`.
4470: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
4471: @*/
4472: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
4473: {
4476: *rejects = ts->reject;
4477: return 0;
4478: }
4480: /*@
4481: TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`
4483: Not Collective
4485: Input Parameter:
4486: . ts - `TS` context
4488: Output Parameter:
4489: . fails - number of failed nonlinear solves
4491: Level: intermediate
4493: Note:
4494: This counter is reset to zero for each successive call to `TSSolve()`.
4496: .seealso: [](chapter_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
4497: @*/
4498: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
4499: {
4502: *fails = ts->num_snes_failures;
4503: return 0;
4504: }
4506: /*@
4507: TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails
4509: Not Collective
4511: Input Parameters:
4512: + ts - `TS` context
4513: - rejects - maximum number of rejected steps, pass -1 for unlimited
4515: Options Database Key:
4516: . -ts_max_reject - Maximum number of step rejections before a step fails
4518: Level: intermediate
4520: .seealso: [](chapter_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4521: @*/
4522: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
4523: {
4525: ts->max_reject = rejects;
4526: return 0;
4527: }
4529: /*@
4530: TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves
4532: Not Collective
4534: Input Parameters:
4535: + ts - `TS` context
4536: - fails - maximum number of failed nonlinear solves, pass -1 for unlimited
4538: Options Database Key:
4539: . -ts_max_snes_failures - Maximum number of nonlinear solve failures
4541: Level: intermediate
4543: .seealso: [](chapter_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
4544: @*/
4545: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
4546: {
4548: ts->max_snes_failures = fails;
4549: return 0;
4550: }
4552: /*@
4553: TSSetErrorIfStepFails - Immediately error if no step succeeds
4555: Not Collective
4557: Input Parameters:
4558: + ts - `TS` context
4559: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure
4561: Options Database Key:
4562: . -ts_error_if_step_fails - Error if no step succeeds
4564: Level: intermediate
4566: .seealso: [](chapter_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
4567: @*/
4568: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
4569: {
4571: ts->errorifstepfailed = err;
4572: return 0;
4573: }
4575: /*@
4576: TSGetAdapt - Get the adaptive controller context for the current method
4578: Collective on ts if controller has not been created yet
4580: Input Parameter:
4581: . ts - time stepping context
4583: Output Parameter:
4584: . adapt - adaptive controller
4586: Level: intermediate
4588: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
4589: @*/
4590: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
4591: {
4594: if (!ts->adapt) {
4595: TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt);
4596: PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1);
4597: }
4598: *adapt = ts->adapt;
4599: return 0;
4600: }
4602: /*@
4603: TSSetTolerances - Set tolerances for local truncation error when using adaptive controller
4605: Logically Collective
4607: Input Parameters:
4608: + ts - time integration context
4609: . atol - scalar absolute tolerances, `PETSC_DECIDE` to leave current value
4610: . vatol - vector of absolute tolerances or NULL, used in preference to atol if present
4611: . rtol - scalar relative tolerances, `PETSC_DECIDE` to leave current value
4612: - vrtol - vector of relative tolerances or NULL, used in preference to atol if present
4614: Options Database keys:
4615: + -ts_rtol <rtol> - relative tolerance for local truncation error
4616: - -ts_atol <atol> - Absolute tolerance for local truncation error
4618: Level: beginner
4620: Notes:
4621: With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
4622: (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
4623: computed only for the differential or the algebraic part then this can be done using the vector of
4624: tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
4625: differential part and infinity for the algebraic part, the LTE calculation will include only the
4626: differential variables.
4628: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
4629: @*/
4630: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
4631: {
4632: if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol;
4633: if (vatol) {
4634: PetscObjectReference((PetscObject)vatol);
4635: VecDestroy(&ts->vatol);
4636: ts->vatol = vatol;
4637: }
4638: if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol;
4639: if (vrtol) {
4640: PetscObjectReference((PetscObject)vrtol);
4641: VecDestroy(&ts->vrtol);
4642: ts->vrtol = vrtol;
4643: }
4644: return 0;
4645: }
4647: /*@
4648: TSGetTolerances - Get tolerances for local truncation error when using adaptive controller
4650: Logically Collective
4652: Input Parameter:
4653: . ts - time integration context
4655: Output Parameters:
4656: + atol - scalar absolute tolerances, NULL to ignore
4657: . vatol - vector of absolute tolerances, NULL to ignore
4658: . rtol - scalar relative tolerances, NULL to ignore
4659: - vrtol - vector of relative tolerances, NULL to ignore
4661: Level: beginner
4663: .seealso: [](chapter_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
4664: @*/
4665: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
4666: {
4667: if (atol) *atol = ts->atol;
4668: if (vatol) *vatol = ts->vatol;
4669: if (rtol) *rtol = ts->rtol;
4670: if (vrtol) *vrtol = ts->vrtol;
4671: return 0;
4672: }
4674: /*@
4675: TSErrorWeightedNorm2 - compute a weighted 2-norm of the difference between two state vectors
4677: Collective
4679: Input Parameters:
4680: + ts - time stepping context
4681: . U - state vector, usually ts->vec_sol
4682: - Y - state vector to be compared to U
4684: Output Parameters:
4685: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
4686: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
4687: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
4689: Level: developer
4691: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedNorm()`, `TSErrorWeightedNormInfinity()`
4692: @*/
4693: PetscErrorCode TSErrorWeightedNorm2(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4694: {
4695: PetscInt i, n, N, rstart;
4696: PetscInt n_loc, na_loc, nr_loc;
4697: PetscReal n_glb, na_glb, nr_glb;
4698: const PetscScalar *u, *y;
4699: PetscReal sum, suma, sumr, gsum, gsuma, gsumr, diff;
4700: PetscReal tol, tola, tolr;
4701: PetscReal err_loc[6], err_glb[6];
4714: VecGetSize(U, &N);
4715: VecGetLocalSize(U, &n);
4716: VecGetOwnershipRange(U, &rstart, NULL);
4717: VecGetArrayRead(U, &u);
4718: VecGetArrayRead(Y, &y);
4719: sum = 0.;
4720: n_loc = 0;
4721: suma = 0.;
4722: na_loc = 0;
4723: sumr = 0.;
4724: nr_loc = 0;
4725: if (ts->vatol && ts->vrtol) {
4726: const PetscScalar *atol, *rtol;
4727: VecGetArrayRead(ts->vatol, &atol);
4728: VecGetArrayRead(ts->vrtol, &rtol);
4729: for (i = 0; i < n; i++) {
4730: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4731: diff = PetscAbsScalar(y[i] - u[i]);
4732: tola = PetscRealPart(atol[i]);
4733: if (tola > 0.) {
4734: suma += PetscSqr(diff / tola);
4735: na_loc++;
4736: }
4737: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4738: if (tolr > 0.) {
4739: sumr += PetscSqr(diff / tolr);
4740: nr_loc++;
4741: }
4742: tol = tola + tolr;
4743: if (tol > 0.) {
4744: sum += PetscSqr(diff / tol);
4745: n_loc++;
4746: }
4747: }
4748: VecRestoreArrayRead(ts->vatol, &atol);
4749: VecRestoreArrayRead(ts->vrtol, &rtol);
4750: } else if (ts->vatol) { /* vector atol, scalar rtol */
4751: const PetscScalar *atol;
4752: VecGetArrayRead(ts->vatol, &atol);
4753: for (i = 0; i < n; i++) {
4754: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4755: diff = PetscAbsScalar(y[i] - u[i]);
4756: tola = PetscRealPart(atol[i]);
4757: if (tola > 0.) {
4758: suma += PetscSqr(diff / tola);
4759: na_loc++;
4760: }
4761: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4762: if (tolr > 0.) {
4763: sumr += PetscSqr(diff / tolr);
4764: nr_loc++;
4765: }
4766: tol = tola + tolr;
4767: if (tol > 0.) {
4768: sum += PetscSqr(diff / tol);
4769: n_loc++;
4770: }
4771: }
4772: VecRestoreArrayRead(ts->vatol, &atol);
4773: } else if (ts->vrtol) { /* scalar atol, vector rtol */
4774: const PetscScalar *rtol;
4775: VecGetArrayRead(ts->vrtol, &rtol);
4776: for (i = 0; i < n; i++) {
4777: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4778: diff = PetscAbsScalar(y[i] - u[i]);
4779: tola = ts->atol;
4780: if (tola > 0.) {
4781: suma += PetscSqr(diff / tola);
4782: na_loc++;
4783: }
4784: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4785: if (tolr > 0.) {
4786: sumr += PetscSqr(diff / tolr);
4787: nr_loc++;
4788: }
4789: tol = tola + tolr;
4790: if (tol > 0.) {
4791: sum += PetscSqr(diff / tol);
4792: n_loc++;
4793: }
4794: }
4795: VecRestoreArrayRead(ts->vrtol, &rtol);
4796: } else { /* scalar atol, scalar rtol */
4797: for (i = 0; i < n; i++) {
4798: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4799: diff = PetscAbsScalar(y[i] - u[i]);
4800: tola = ts->atol;
4801: if (tola > 0.) {
4802: suma += PetscSqr(diff / tola);
4803: na_loc++;
4804: }
4805: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4806: if (tolr > 0.) {
4807: sumr += PetscSqr(diff / tolr);
4808: nr_loc++;
4809: }
4810: tol = tola + tolr;
4811: if (tol > 0.) {
4812: sum += PetscSqr(diff / tol);
4813: n_loc++;
4814: }
4815: }
4816: }
4817: VecRestoreArrayRead(U, &u);
4818: VecRestoreArrayRead(Y, &y);
4820: err_loc[0] = sum;
4821: err_loc[1] = suma;
4822: err_loc[2] = sumr;
4823: err_loc[3] = (PetscReal)n_loc;
4824: err_loc[4] = (PetscReal)na_loc;
4825: err_loc[5] = (PetscReal)nr_loc;
4827: MPIU_Allreduce(err_loc, err_glb, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts));
4829: gsum = err_glb[0];
4830: gsuma = err_glb[1];
4831: gsumr = err_glb[2];
4832: n_glb = err_glb[3];
4833: na_glb = err_glb[4];
4834: nr_glb = err_glb[5];
4836: *norm = 0.;
4837: if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
4838: *norma = 0.;
4839: if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
4840: *normr = 0.;
4841: if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);
4846: return 0;
4847: }
4849: /*@
4850: TSErrorWeightedNormInfinity - compute a weighted infinity-norm of the difference between two state vectors
4852: Collective
4854: Input Parameters:
4855: + ts - time stepping context
4856: . U - state vector, usually ts->vec_sol
4857: - Y - state vector to be compared to U
4859: Output Parameters:
4860: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
4861: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
4862: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
4864: Level: developer
4866: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedNorm()`, `TSErrorWeightedNorm2()`
4867: @*/
4868: PetscErrorCode TSErrorWeightedNormInfinity(TS ts, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4869: {
4870: PetscInt i, n, N, rstart;
4871: const PetscScalar *u, *y;
4872: PetscReal max, gmax, maxa, gmaxa, maxr, gmaxr;
4873: PetscReal tol, tola, tolr, diff;
4874: PetscReal err_loc[3], err_glb[3];
4887: VecGetSize(U, &N);
4888: VecGetLocalSize(U, &n);
4889: VecGetOwnershipRange(U, &rstart, NULL);
4890: VecGetArrayRead(U, &u);
4891: VecGetArrayRead(Y, &y);
4893: max = 0.;
4894: maxa = 0.;
4895: maxr = 0.;
4897: if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
4898: const PetscScalar *atol, *rtol;
4899: VecGetArrayRead(ts->vatol, &atol);
4900: VecGetArrayRead(ts->vrtol, &rtol);
4902: for (i = 0; i < n; i++) {
4903: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4904: diff = PetscAbsScalar(y[i] - u[i]);
4905: tola = PetscRealPart(atol[i]);
4906: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4907: tol = tola + tolr;
4908: if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4909: if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4910: if (tol > 0.) max = PetscMax(max, diff / tol);
4911: }
4912: VecRestoreArrayRead(ts->vatol, &atol);
4913: VecRestoreArrayRead(ts->vrtol, &rtol);
4914: } else if (ts->vatol) { /* vector atol, scalar rtol */
4915: const PetscScalar *atol;
4916: VecGetArrayRead(ts->vatol, &atol);
4917: for (i = 0; i < n; i++) {
4918: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4919: diff = PetscAbsScalar(y[i] - u[i]);
4920: tola = PetscRealPart(atol[i]);
4921: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4922: tol = tola + tolr;
4923: if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4924: if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4925: if (tol > 0.) max = PetscMax(max, diff / tol);
4926: }
4927: VecRestoreArrayRead(ts->vatol, &atol);
4928: } else if (ts->vrtol) { /* scalar atol, vector rtol */
4929: const PetscScalar *rtol;
4930: VecGetArrayRead(ts->vrtol, &rtol);
4932: for (i = 0; i < n; i++) {
4933: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4934: diff = PetscAbsScalar(y[i] - u[i]);
4935: tola = ts->atol;
4936: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4937: tol = tola + tolr;
4938: if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4939: if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4940: if (tol > 0.) max = PetscMax(max, diff / tol);
4941: }
4942: VecRestoreArrayRead(ts->vrtol, &rtol);
4943: } else { /* scalar atol, scalar rtol */
4945: for (i = 0; i < n; i++) {
4946: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
4947: diff = PetscAbsScalar(y[i] - u[i]);
4948: tola = ts->atol;
4949: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
4950: tol = tola + tolr;
4951: if (tola > 0.) maxa = PetscMax(maxa, diff / tola);
4952: if (tolr > 0.) maxr = PetscMax(maxr, diff / tolr);
4953: if (tol > 0.) max = PetscMax(max, diff / tol);
4954: }
4955: }
4956: VecRestoreArrayRead(U, &u);
4957: VecRestoreArrayRead(Y, &y);
4958: err_loc[0] = max;
4959: err_loc[1] = maxa;
4960: err_loc[2] = maxr;
4961: MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
4962: gmax = err_glb[0];
4963: gmaxa = err_glb[1];
4964: gmaxr = err_glb[2];
4966: *norm = gmax;
4967: *norma = gmaxa;
4968: *normr = gmaxr;
4972: return 0;
4973: }
4975: /*@
4976: TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances
4978: Collective
4980: Input Parameters:
4981: + ts - time stepping context
4982: . U - state vector, usually ts->vec_sol
4983: . Y - state vector to be compared to U
4984: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
4986: Output Parameters:
4987: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
4988: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
4989: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
4991: Options Database Key:
4992: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
4994: Level: developer
4996: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`, `TSErrorWeightedENorm`
4997: @*/
4998: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
4999: {
5000: if (wnormtype == NORM_2) TSErrorWeightedNorm2(ts, U, Y, norm, norma, normr);
5001: else if (wnormtype == NORM_INFINITY) TSErrorWeightedNormInfinity(ts, U, Y, norm, norma, normr);
5002: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5003: return 0;
5004: }
5006: /*@
5007: TSErrorWeightedENorm2 - compute a weighted 2 error norm based on supplied absolute and relative tolerances
5009: Collective
5011: Input Parameters:
5012: + ts - time stepping context
5013: . E - error vector
5014: . U - state vector, usually ts->vec_sol
5015: - Y - state vector, previous time step
5017: Output Parameters:
5018: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5019: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5020: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
5022: Level: developer
5024: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedENorm()`, `TSErrorWeightedENormInfinity()`
5025: @*/
5026: PetscErrorCode TSErrorWeightedENorm2(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5027: {
5028: PetscInt i, n, N, rstart;
5029: PetscInt n_loc, na_loc, nr_loc;
5030: PetscReal n_glb, na_glb, nr_glb;
5031: const PetscScalar *e, *u, *y;
5032: PetscReal err, sum, suma, sumr, gsum, gsuma, gsumr;
5033: PetscReal tol, tola, tolr;
5034: PetscReal err_loc[6], err_glb[6];
5049: VecGetSize(E, &N);
5050: VecGetLocalSize(E, &n);
5051: VecGetOwnershipRange(E, &rstart, NULL);
5052: VecGetArrayRead(E, &e);
5053: VecGetArrayRead(U, &u);
5054: VecGetArrayRead(Y, &y);
5055: sum = 0.;
5056: n_loc = 0;
5057: suma = 0.;
5058: na_loc = 0;
5059: sumr = 0.;
5060: nr_loc = 0;
5061: if (ts->vatol && ts->vrtol) {
5062: const PetscScalar *atol, *rtol;
5063: VecGetArrayRead(ts->vatol, &atol);
5064: VecGetArrayRead(ts->vrtol, &rtol);
5065: for (i = 0; i < n; i++) {
5066: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5067: err = PetscAbsScalar(e[i]);
5068: tola = PetscRealPart(atol[i]);
5069: if (tola > 0.) {
5070: suma += PetscSqr(err / tola);
5071: na_loc++;
5072: }
5073: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5074: if (tolr > 0.) {
5075: sumr += PetscSqr(err / tolr);
5076: nr_loc++;
5077: }
5078: tol = tola + tolr;
5079: if (tol > 0.) {
5080: sum += PetscSqr(err / tol);
5081: n_loc++;
5082: }
5083: }
5084: VecRestoreArrayRead(ts->vatol, &atol);
5085: VecRestoreArrayRead(ts->vrtol, &rtol);
5086: } else if (ts->vatol) { /* vector atol, scalar rtol */
5087: const PetscScalar *atol;
5088: VecGetArrayRead(ts->vatol, &atol);
5089: for (i = 0; i < n; i++) {
5090: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5091: err = PetscAbsScalar(e[i]);
5092: tola = PetscRealPart(atol[i]);
5093: if (tola > 0.) {
5094: suma += PetscSqr(err / tola);
5095: na_loc++;
5096: }
5097: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5098: if (tolr > 0.) {
5099: sumr += PetscSqr(err / tolr);
5100: nr_loc++;
5101: }
5102: tol = tola + tolr;
5103: if (tol > 0.) {
5104: sum += PetscSqr(err / tol);
5105: n_loc++;
5106: }
5107: }
5108: VecRestoreArrayRead(ts->vatol, &atol);
5109: } else if (ts->vrtol) { /* scalar atol, vector rtol */
5110: const PetscScalar *rtol;
5111: VecGetArrayRead(ts->vrtol, &rtol);
5112: for (i = 0; i < n; i++) {
5113: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5114: err = PetscAbsScalar(e[i]);
5115: tola = ts->atol;
5116: if (tola > 0.) {
5117: suma += PetscSqr(err / tola);
5118: na_loc++;
5119: }
5120: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5121: if (tolr > 0.) {
5122: sumr += PetscSqr(err / tolr);
5123: nr_loc++;
5124: }
5125: tol = tola + tolr;
5126: if (tol > 0.) {
5127: sum += PetscSqr(err / tol);
5128: n_loc++;
5129: }
5130: }
5131: VecRestoreArrayRead(ts->vrtol, &rtol);
5132: } else { /* scalar atol, scalar rtol */
5133: for (i = 0; i < n; i++) {
5134: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5135: err = PetscAbsScalar(e[i]);
5136: tola = ts->atol;
5137: if (tola > 0.) {
5138: suma += PetscSqr(err / tola);
5139: na_loc++;
5140: }
5141: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5142: if (tolr > 0.) {
5143: sumr += PetscSqr(err / tolr);
5144: nr_loc++;
5145: }
5146: tol = tola + tolr;
5147: if (tol > 0.) {
5148: sum += PetscSqr(err / tol);
5149: n_loc++;
5150: }
5151: }
5152: }
5153: VecRestoreArrayRead(E, &e);
5154: VecRestoreArrayRead(U, &u);
5155: VecRestoreArrayRead(Y, &y);
5157: err_loc[0] = sum;
5158: err_loc[1] = suma;
5159: err_loc[2] = sumr;
5160: err_loc[3] = (PetscReal)n_loc;
5161: err_loc[4] = (PetscReal)na_loc;
5162: err_loc[5] = (PetscReal)nr_loc;
5164: MPIU_Allreduce(err_loc, err_glb, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts));
5166: gsum = err_glb[0];
5167: gsuma = err_glb[1];
5168: gsumr = err_glb[2];
5169: n_glb = err_glb[3];
5170: na_glb = err_glb[4];
5171: nr_glb = err_glb[5];
5173: *norm = 0.;
5174: if (n_glb > 0.) *norm = PetscSqrtReal(gsum / n_glb);
5175: *norma = 0.;
5176: if (na_glb > 0.) *norma = PetscSqrtReal(gsuma / na_glb);
5177: *normr = 0.;
5178: if (nr_glb > 0.) *normr = PetscSqrtReal(gsumr / nr_glb);
5183: return 0;
5184: }
5186: /*@
5187: TSErrorWeightedENormInfinity - compute a weighted infinity error norm based on supplied absolute and relative tolerances
5189: Collective
5191: Input Parameters:
5192: + ts - time stepping context
5193: . E - error vector
5194: . U - state vector, usually ts->vec_sol
5195: - Y - state vector, previous time step
5197: Output Parameters:
5198: + norm - weighted norm, a value of 1.0 means that the error matches the tolerances
5199: . norma - weighted norm based on the absolute tolerance, a value of 1.0 means that the error matches the tolerances
5200: - normr - weighted norm based on the relative tolerance, a value of 1.0 means that the error matches the tolerances
5202: Level: developer
5204: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedENorm()`, `TSErrorWeightedENorm2()`
5205: @*/
5206: PetscErrorCode TSErrorWeightedENormInfinity(TS ts, Vec E, Vec U, Vec Y, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5207: {
5208: PetscInt i, n, N, rstart;
5209: const PetscScalar *e, *u, *y;
5210: PetscReal err, max, gmax, maxa, gmaxa, maxr, gmaxr;
5211: PetscReal tol, tola, tolr;
5212: PetscReal err_loc[3], err_glb[3];
5227: VecGetSize(E, &N);
5228: VecGetLocalSize(E, &n);
5229: VecGetOwnershipRange(E, &rstart, NULL);
5230: VecGetArrayRead(E, &e);
5231: VecGetArrayRead(U, &u);
5232: VecGetArrayRead(Y, &y);
5234: max = 0.;
5235: maxa = 0.;
5236: maxr = 0.;
5238: if (ts->vatol && ts->vrtol) { /* vector atol, vector rtol */
5239: const PetscScalar *atol, *rtol;
5240: VecGetArrayRead(ts->vatol, &atol);
5241: VecGetArrayRead(ts->vrtol, &rtol);
5243: for (i = 0; i < n; i++) {
5244: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5245: err = PetscAbsScalar(e[i]);
5246: tola = PetscRealPart(atol[i]);
5247: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5248: tol = tola + tolr;
5249: if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5250: if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5251: if (tol > 0.) max = PetscMax(max, err / tol);
5252: }
5253: VecRestoreArrayRead(ts->vatol, &atol);
5254: VecRestoreArrayRead(ts->vrtol, &rtol);
5255: } else if (ts->vatol) { /* vector atol, scalar rtol */
5256: const PetscScalar *atol;
5257: VecGetArrayRead(ts->vatol, &atol);
5258: for (i = 0; i < n; i++) {
5259: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5260: err = PetscAbsScalar(e[i]);
5261: tola = PetscRealPart(atol[i]);
5262: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5263: tol = tola + tolr;
5264: if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5265: if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5266: if (tol > 0.) max = PetscMax(max, err / tol);
5267: }
5268: VecRestoreArrayRead(ts->vatol, &atol);
5269: } else if (ts->vrtol) { /* scalar atol, vector rtol */
5270: const PetscScalar *rtol;
5271: VecGetArrayRead(ts->vrtol, &rtol);
5273: for (i = 0; i < n; i++) {
5274: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5275: err = PetscAbsScalar(e[i]);
5276: tola = ts->atol;
5277: tolr = PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5278: tol = tola + tolr;
5279: if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5280: if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5281: if (tol > 0.) max = PetscMax(max, err / tol);
5282: }
5283: VecRestoreArrayRead(ts->vrtol, &rtol);
5284: } else { /* scalar atol, scalar rtol */
5286: for (i = 0; i < n; i++) {
5287: SkipSmallValue(y[i], u[i], ts->adapt->ignore_max);
5288: err = PetscAbsScalar(e[i]);
5289: tola = ts->atol;
5290: tolr = ts->rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
5291: tol = tola + tolr;
5292: if (tola > 0.) maxa = PetscMax(maxa, err / tola);
5293: if (tolr > 0.) maxr = PetscMax(maxr, err / tolr);
5294: if (tol > 0.) max = PetscMax(max, err / tol);
5295: }
5296: }
5297: VecRestoreArrayRead(E, &e);
5298: VecRestoreArrayRead(U, &u);
5299: VecRestoreArrayRead(Y, &y);
5300: err_loc[0] = max;
5301: err_loc[1] = maxa;
5302: err_loc[2] = maxr;
5303: MPIU_Allreduce(err_loc, err_glb, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
5304: gmax = err_glb[0];
5305: gmaxa = err_glb[1];
5306: gmaxr = err_glb[2];
5308: *norm = gmax;
5309: *norma = gmaxa;
5310: *normr = gmaxr;
5314: return 0;
5315: }
5317: /*@
5318: TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances
5320: Collective
5322: Input Parameters:
5323: + ts - time stepping context
5324: . E - error vector
5325: . U - state vector, usually ts->vec_sol
5326: . Y - state vector, previous time step
5327: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`
5329: Output Parameters:
5330: + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
5331: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
5332: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user
5334: Options Database Key:
5335: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY
5337: Level: developer
5339: .seealso: [](chapter_ts), `TS`, `TSErrorWeightedENormInfinity()`, `TSErrorWeightedENorm2()`, `TSErrorWeightedNormInfinity()`, `TSErrorWeightedNorm2()`
5340: @*/
5341: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
5342: {
5343: if (wnormtype == NORM_2) TSErrorWeightedENorm2(ts, E, U, Y, norm, norma, normr);
5344: else if (wnormtype == NORM_INFINITY) TSErrorWeightedENormInfinity(ts, E, U, Y, norm, norma, normr);
5345: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
5346: return 0;
5347: }
5349: /*@
5350: TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler
5352: Logically Collective
5354: Input Parameters:
5355: + ts - time stepping context
5356: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)
5358: Note:
5359: After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()
5361: Level: intermediate
5363: .seealso: [](chapter_ts), `TSGetCFLTime()`, `TSADAPTCFL`
5364: @*/
5365: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
5366: {
5368: ts->cfltime_local = cfltime;
5369: ts->cfltime = -1.;
5370: return 0;
5371: }
5373: /*@
5374: TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler
5376: Collective
5378: Input Parameter:
5379: . ts - time stepping context
5381: Output Parameter:
5382: . cfltime - maximum stable time step for forward Euler
5384: Level: advanced
5386: .seealso: [](chapter_ts), `TSSetCFLTimeLocal()`
5387: @*/
5388: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
5389: {
5390: if (ts->cfltime < 0) MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts));
5391: *cfltime = ts->cfltime;
5392: return 0;
5393: }
5395: /*@
5396: TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu
5398: Input Parameters:
5399: + ts - the `TS` context.
5400: . xl - lower bound.
5401: - xu - upper bound.
5403: Level: advanced
5405: Note:
5406: If this routine is not called then the lower and upper bounds are set to
5407: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
5409: .seealso: [](chapter_ts), `TS`
5410: @*/
5411: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
5412: {
5413: SNES snes;
5415: TSGetSNES(ts, &snes);
5416: SNESVISetVariableBounds(snes, xl, xu);
5417: return 0;
5418: }
5420: /*@
5421: TSComputeLinearStability - computes the linear stability function at a point
5423: Collective
5425: Input Parameters:
5426: + ts - the `TS` context
5427: - xr,xi - real and imaginary part of input arguments
5429: Output Parameters:
5430: . yr,yi - real and imaginary part of function value
5432: Level: developer
5434: .seealso: [](chapter_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
5435: @*/
5436: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
5437: {
5439: PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
5440: return 0;
5441: }
5443: /*@
5444: TSRestartStep - Flags the solver to restart the next step
5446: Collective
5448: Input Parameter:
5449: . ts - the `TS` context obtained from `TSCreate()`
5451: Level: advanced
5453: Notes:
5454: Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
5455: discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
5456: vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
5457: the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
5458: discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
5459: discontinuous source terms).
5461: .seealso: [](chapter_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
5462: @*/
5463: PetscErrorCode TSRestartStep(TS ts)
5464: {
5466: ts->steprestart = PETSC_TRUE;
5467: return 0;
5468: }
5470: /*@
5471: TSRollBack - Rolls back one time step
5473: Collective
5475: Input Parameter:
5476: . ts - the `TS` context obtained from `TSCreate()`
5478: Level: advanced
5480: .seealso: [](chapter_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
5481: @*/
5482: PetscErrorCode TSRollBack(TS ts)
5483: {
5486: PetscUseTypeMethod(ts, rollback);
5487: ts->time_step = ts->ptime - ts->ptime_prev;
5488: ts->ptime = ts->ptime_prev;
5489: ts->ptime_prev = ts->ptime_prev_rollback;
5490: ts->steps--;
5491: ts->steprollback = PETSC_TRUE;
5492: return 0;
5493: }
5495: /*@
5496: TSGetStages - Get the number of stages and stage values
5498: Input Parameter:
5499: . ts - the `TS` context obtained from `TSCreate()`
5501: Output Parameters:
5502: + ns - the number of stages
5503: - Y - the current stage vectors
5505: Level: advanced
5507: Note:
5508: Both ns and Y can be NULL.
5510: .seealso: [](chapter_ts), `TS`, `TSCreate()`
5511: @*/
5512: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
5513: {
5517: if (!ts->ops->getstages) {
5518: if (ns) *ns = 0;
5519: if (Y) *Y = NULL;
5520: } else PetscUseTypeMethod(ts, getstages, ns, Y);
5521: return 0;
5522: }
5524: /*@C
5525: TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.
5527: Collective
5529: Input Parameters:
5530: + ts - the `TS` context
5531: . t - current timestep
5532: . U - state vector
5533: . Udot - time derivative of state vector
5534: . shift - shift to apply, see note below
5535: - ctx - an optional user context
5537: Output Parameters:
5538: + J - Jacobian matrix (not altered in this routine)
5539: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)
5541: Level: intermediate
5543: Notes:
5544: If F(t,U,Udot)=0 is the DAE, the required Jacobian is
5546: dF/dU + shift*dF/dUdot
5548: Most users should not need to explicitly call this routine, as it
5549: is used internally within the nonlinear solvers.
5551: This will first try to get the coloring from the `DM`. If the `DM` type has no coloring
5552: routine, then it will try to get the coloring from the matrix. This requires that the
5553: matrix have nonzero entries precomputed.
5555: .seealso: [](chapter_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
5556: @*/
5557: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
5558: {
5559: SNES snes;
5560: MatFDColoring color;
5561: PetscBool hascolor, matcolor = PETSC_FALSE;
5563: PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL);
5564: PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color);
5565: if (!color) {
5566: DM dm;
5567: ISColoring iscoloring;
5569: TSGetDM(ts, &dm);
5570: DMHasColoring(dm, &hascolor);
5571: if (hascolor && !matcolor) {
5572: DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring);
5573: MatFDColoringCreate(B, iscoloring, &color);
5574: MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts);
5575: MatFDColoringSetFromOptions(color);
5576: MatFDColoringSetUp(B, iscoloring, color);
5577: ISColoringDestroy(&iscoloring);
5578: } else {
5579: MatColoring mc;
5581: MatColoringCreate(B, &mc);
5582: MatColoringSetDistance(mc, 2);
5583: MatColoringSetType(mc, MATCOLORINGSL);
5584: MatColoringSetFromOptions(mc);
5585: MatColoringApply(mc, &iscoloring);
5586: MatColoringDestroy(&mc);
5587: MatFDColoringCreate(B, iscoloring, &color);
5588: MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts);
5589: MatFDColoringSetFromOptions(color);
5590: MatFDColoringSetUp(B, iscoloring, color);
5591: ISColoringDestroy(&iscoloring);
5592: }
5593: PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color);
5594: PetscObjectDereference((PetscObject)color);
5595: }
5596: TSGetSNES(ts, &snes);
5597: MatFDColoringApply(B, color, U, snes);
5598: if (J != B) {
5599: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
5600: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
5601: }
5602: return 0;
5603: }
5605: /*@
5606: TSSetFunctionDomainError - Set a function that tests if the current state vector is valid
5608: Input Parameters:
5609: + ts - the `TS` context
5610: - func - function called within `TSFunctionDomainError()`
5612: Calling sequence of func:
5613: $ PetscErrorCode func(TS ts,PetscReal time,Vec state,PetscBool reject)
5615: + ts - the TS context
5616: . time - the current time (of the stage)
5617: . state - the state to check if it is valid
5618: - reject - (output parameter) PETSC_FALSE if the state is acceptable, PETSC_TRUE if not acceptable
5620: Level: intermediate
5622: Notes:
5623: If an implicit ODE solver is being used then, in addition to providing this routine, the
5624: user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
5625: function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
5626: Use `TSGetSNES()` to obtain the `SNES` object
5628: Developer Note:
5629: The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
5630: since one takes a function pointer and the other does not.
5632: .seealso: [](chapter_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
5633: @*/
5635: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS, PetscReal, Vec, PetscBool *))
5636: {
5638: ts->functiondomainerror = func;
5639: return 0;
5640: }
5642: /*@
5643: TSFunctionDomainError - Checks if the current state is valid
5645: Input Parameters:
5646: + ts - the `TS` context
5647: . stagetime - time of the simulation
5648: - Y - state vector to check.
5650: Output Parameter:
5651: . accept - Set to `PETSC_FALSE` if the current state vector is valid.
5653: Level: developer
5655: Note:
5656: This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
5657: to check if the current state is valid.
5659: .seealso: [](chapter_ts), `TS`, `TSSetFunctionDomainError()`
5660: @*/
5661: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
5662: {
5664: *accept = PETSC_TRUE;
5665: if (ts->functiondomainerror) (*ts->functiondomainerror)(ts, stagetime, Y, accept);
5666: return 0;
5667: }
5669: /*@C
5670: TSClone - This function clones a time step object.
5672: Collective
5674: Input Parameter:
5675: . tsin - The input `TS`
5677: Output Parameter:
5678: . tsout - The output `TS` (cloned)
5680: Level: developer
5682: Notes:
5683: 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.
5684: It will likely be replaced in the future with a mechanism of switching methods on the fly.
5686: 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 `SNES` snes_dup=NULL; `TSGetSNES`(ts,&snes_dup); `TSSetSNES`(ts,snes_dup);
5688: .seealso: [](chapter_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
5689: @*/
5690: PetscErrorCode TSClone(TS tsin, TS *tsout)
5691: {
5692: TS t;
5693: SNES snes_start;
5694: DM dm;
5695: TSType type;
5698: *tsout = NULL;
5700: PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView);
5702: /* General TS description */
5703: t->numbermonitors = 0;
5704: t->monitorFrequency = 1;
5705: t->setupcalled = 0;
5706: t->ksp_its = 0;
5707: t->snes_its = 0;
5708: t->nwork = 0;
5709: t->rhsjacobian.time = PETSC_MIN_REAL;
5710: t->rhsjacobian.scale = 1.;
5711: t->ijacobian.shift = 1.;
5713: TSGetSNES(tsin, &snes_start);
5714: TSSetSNES(t, snes_start);
5716: TSGetDM(tsin, &dm);
5717: TSSetDM(t, dm);
5719: t->adapt = tsin->adapt;
5720: PetscObjectReference((PetscObject)t->adapt);
5722: t->trajectory = tsin->trajectory;
5723: PetscObjectReference((PetscObject)t->trajectory);
5725: t->event = tsin->event;
5726: if (t->event) t->event->refct++;
5728: t->problem_type = tsin->problem_type;
5729: t->ptime = tsin->ptime;
5730: t->ptime_prev = tsin->ptime_prev;
5731: t->time_step = tsin->time_step;
5732: t->max_time = tsin->max_time;
5733: t->steps = tsin->steps;
5734: t->max_steps = tsin->max_steps;
5735: t->equation_type = tsin->equation_type;
5736: t->atol = tsin->atol;
5737: t->rtol = tsin->rtol;
5738: t->max_snes_failures = tsin->max_snes_failures;
5739: t->max_reject = tsin->max_reject;
5740: t->errorifstepfailed = tsin->errorifstepfailed;
5742: TSGetType(tsin, &type);
5743: TSSetType(t, type);
5745: t->vec_sol = NULL;
5747: t->cfltime = tsin->cfltime;
5748: t->cfltime_local = tsin->cfltime_local;
5749: t->exact_final_time = tsin->exact_final_time;
5751: PetscMemcpy(t->ops, tsin->ops, sizeof(struct _TSOps));
5753: if (((PetscObject)tsin)->fortran_func_pointers) {
5754: PetscInt i;
5755: PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers);
5756: for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
5757: }
5758: *tsout = t;
5759: return 0;
5760: }
5762: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
5763: {
5764: TS ts = (TS)ctx;
5766: TSComputeRHSFunction(ts, 0, x, y);
5767: return 0;
5768: }
5770: /*@
5771: TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5773: Logically Collective
5775: Input Parameters:
5776: TS - the time stepping routine
5778: Output Parameter:
5779: . flg - `PETSC_TRUE` if the multiply is likely correct
5781: Options Database Key:
5782: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator
5784: Level: advanced
5786: Note:
5787: This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian
5789: .seealso: [](chapter_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
5790: @*/
5791: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
5792: {
5793: Mat J, B;
5794: TSRHSJacobian func;
5795: void *ctx;
5797: TSGetRHSJacobian(ts, &J, &B, &func, &ctx);
5798: (*func)(ts, 0.0, ts->vec_sol, J, B, ctx);
5799: MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg);
5800: return 0;
5801: }
5803: /*@C
5804: TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.
5806: Logically Collective
5808: Input Parameters:
5809: TS - the time stepping routine
5811: Output Parameter:
5812: . flg - `PETSC_TRUE` if the multiply is likely correct
5814: Options Database Key:
5815: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator
5817: Level: advanced
5819: Notes:
5820: This only works for problems defined only the RHS function and Jacobian NOT IFunction and IJacobian
5822: .seealso: [](chapter_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
5823: @*/
5824: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
5825: {
5826: Mat J, B;
5827: void *ctx;
5828: TSRHSJacobian func;
5830: TSGetRHSJacobian(ts, &J, &B, &func, &ctx);
5831: (*func)(ts, 0.0, ts->vec_sol, J, B, ctx);
5832: MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg);
5833: return 0;
5834: }
5836: /*@
5837: TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.
5839: Logically collective
5841: Input Parameters:
5842: + ts - timestepping context
5843: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5845: Options Database Key:
5846: . -ts_use_splitrhsfunction - <true,false>
5848: Level: intermediate
5850: Note:
5851: This is only useful for multirate methods
5853: .seealso: [](chapter_ts), `TS`, `TSGetUseSplitRHSFunction()`
5854: @*/
5855: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
5856: {
5858: ts->use_splitrhsfunction = use_splitrhsfunction;
5859: return 0;
5860: }
5862: /*@
5863: TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.
5865: Not collective
5867: Input Parameter:
5868: . ts - timestepping context
5870: Output Parameter:
5871: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used
5873: Level: intermediate
5875: .seealso: [](chapter_ts), `TS`, `TSSetUseSplitRHSFunction()`
5876: @*/
5877: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
5878: {
5880: *use_splitrhsfunction = ts->use_splitrhsfunction;
5881: return 0;
5882: }
5884: /*@
5885: TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.
5887: Logically Collective
5889: Input Parameters:
5890: + ts - the time-stepper
5891: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)
5893: Level: intermediate
5895: Note:
5896: When the relationship between the nonzero structures is known and supplied the solution process can be much faster
5898: .seealso: [](chapter_ts), `TS`, `MatAXPY()`, `MatStructure`
5899: @*/
5900: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
5901: {
5903: ts->axpy_pattern = str;
5904: return 0;
5905: }
5907: /*@
5908: TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span
5910: Collective
5912: Input Parameters:
5913: + ts - the time-stepper
5914: . n - number of the time points (>=2)
5915: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5917: Options Database Key:
5918: . -ts_time_span <t0,...tf> - Sets the time span
5920: Level: intermediate
5922: Notes:
5923: The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
5924: `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
5925: The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may
5926: pressure the memory system when using a large number of span points.
5928: .seealso: [](chapter_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()`
5929: @*/
5930: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
5931: {
5934: if (ts->tspan && n != ts->tspan->num_span_times) {
5935: PetscFree(ts->tspan->span_times);
5936: VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol);
5937: PetscMalloc1(n, &ts->tspan->span_times);
5938: }
5939: if (!ts->tspan) {
5940: TSTimeSpan tspan;
5941: PetscNew(&tspan);
5942: PetscMalloc1(n, &tspan->span_times);
5943: tspan->reltol = 1e-6;
5944: tspan->abstol = 10 * PETSC_MACHINE_EPSILON;
5945: ts->tspan = tspan;
5946: }
5947: ts->tspan->num_span_times = n;
5948: PetscArraycpy(ts->tspan->span_times, span_times, n);
5949: TSSetTime(ts, ts->tspan->span_times[0]);
5950: TSSetMaxTime(ts, ts->tspan->span_times[n - 1]);
5951: return 0;
5952: }
5954: /*@C
5955: TSGetTimeSpan - gets the time span.
5957: Not Collective
5959: Input Parameter:
5960: . ts - the time-stepper
5962: Output Parameters:
5963: + n - number of the time points (>=2)
5964: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.
5965: The values are valid until the `TS` object is destroyed.
5967: Level: beginner
5969: Note:
5970: Both n and span_times can be NULL.
5972: .seealso: [](chapter_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()`
5973: @*/
5974: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal **span_times)
5975: {
5979: if (!ts->tspan) {
5980: if (n) *n = 0;
5981: if (span_times) *span_times = NULL;
5982: } else {
5983: if (n) *n = ts->tspan->num_span_times;
5984: if (span_times) *span_times = ts->tspan->span_times;
5985: }
5986: return 0;
5987: }
5989: /*@
5990: TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.
5992: Input Parameter:
5993: . ts - the `TS` context obtained from `TSCreate()`
5995: Output Parameters:
5996: + nsol - the number of solutions
5997: - Sols - the solution vectors
5999: Level: intermediate
6001: Notes:
6002: Both nsol and Sols can be NULL.
6004: Some time points in the time span may be skipped by TS so that nsol is less than the number of points specified by `TSSetTimeSpan()`.
6005: For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span.
6007: .seealso: [](chapter_ts), `TS`, `TSSetTimeSpan()`
6008: @*/
6009: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
6010: {
6014: if (!ts->tspan) {
6015: if (nsol) *nsol = 0;
6016: if (Sols) *Sols = NULL;
6017: } else {
6018: if (nsol) *nsol = ts->tspan->spanctr;
6019: if (Sols) *Sols = ts->tspan->vecs_sol;
6020: }
6021: return 0;
6022: }