1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscdmda.h> 3 #include <petscdmshell.h> 4 #include <petscdmplex.h> // For TSSetFromOptions() 5 #include <petscdmswarm.h> // For TSSetFromOptions() 6 #include <petscviewer.h> 7 #include <petscdraw.h> 8 #include <petscconvest.h> 9 10 /* Logging support */ 11 PetscClassId TS_CLASSID, DMTS_CLASSID; 12 PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval; 13 14 const char *const TSExactFinalTimeOptions[] = {"UNSPECIFIED", "STEPOVER", "INTERPOLATE", "MATCHSTEP", "TSExactFinalTimeOption", "TS_EXACTFINALTIME_", NULL}; 15 16 static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt, TSAdaptType default_type) 17 { 18 PetscFunctionBegin; 19 PetscValidHeaderSpecific(adapt, TSADAPT_CLASSID, 1); 20 PetscValidCharPointer(default_type, 2); 21 if (!((PetscObject)adapt)->type_name) PetscCall(TSAdaptSetType(adapt, default_type)); 22 PetscFunctionReturn(PETSC_SUCCESS); 23 } 24 25 /*@ 26 TSSetFromOptions - Sets various `TS` parameters from the options database 27 28 Collective 29 30 Input Parameter: 31 . ts - the `TS` context obtained from `TSCreate()` 32 33 Options Database Keys: 34 + -ts_type <type> - EULER, BEULER, SUNDIALS, PSEUDO, CN, RK, THETA, ALPHA, GLLE, SSP, GLEE, BSYMP, IRK, see `TSType` 35 . -ts_save_trajectory - checkpoint the solution at each time-step 36 . -ts_max_time <time> - maximum time to compute to 37 . -ts_time_span <t0,...tf> - sets the time span, solutions are computed and stored for each indicated time 38 . -ts_max_steps <steps> - maximum number of time-steps to take 39 . -ts_init_time <time> - initial time to start computation 40 . -ts_final_time <time> - final time to compute to (deprecated: use `-ts_max_time`) 41 . -ts_dt <dt> - initial time step 42 . -ts_exact_final_time <stepover,interpolate,matchstep> - whether to stop at the exact given final time and how to compute the solution at that time 43 . -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed 44 . -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails 45 . -ts_error_if_step_fails <true,false> - Error if no step succeeds 46 . -ts_rtol <rtol> - relative tolerance for local truncation error 47 . -ts_atol <atol> - Absolute tolerance for local truncation error 48 . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - test the Jacobian at each iteration against finite difference with RHS function 49 . -ts_rhs_jacobian_test_mult_transpose -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_interval <interval> - output once every interval (default=1) time steps; used with -ts_monitor_solution 67 . -ts_monitor_solution_vtk <filename.vts,filename.vtu> - Save each time step to a binary file, use filename-%%03" PetscInt_FMT ".vts (filename-%%03" PetscInt_FMT ".vtu) 68 - -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time 69 70 Level: beginner 71 72 Notes: 73 See `SNESSetFromOptions()` and `KSPSetFromOptions()` for how to control the nonlinear and linear solves used by the time-stepper. 74 75 Certain `SNES` options get reset for each new nonlinear solver, for example `-snes_lag_jacobian its` and `-snes_lag_preconditioner its`, in order 76 to retain them over the multiple nonlinear solves that `TS` uses you mush also provide `-snes_lag_jacobian_persists true` and 77 `-snes_lag_preconditioner_persists true` 78 79 Developer Note: 80 We should unify all the -ts_monitor options in the way that -xxx_view has been unified 81 82 .seealso: [](ch_ts), `TS`, `TSGetType()` 83 @*/ 84 PetscErrorCode TSSetFromOptions(TS ts) 85 { 86 PetscBool opt, flg, tflg; 87 char monfilename[PETSC_MAX_PATH_LEN]; 88 PetscReal time_step, tspan[100]; 89 PetscInt nt = PETSC_STATIC_ARRAY_LENGTH(tspan); 90 TSExactFinalTimeOption eftopt; 91 char dir[16]; 92 TSIFunction ifun; 93 const char *defaultType; 94 char typeName[256]; 95 96 PetscFunctionBegin; 97 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 98 99 PetscCall(TSRegisterAll()); 100 PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL)); 101 102 PetscObjectOptionsBegin((PetscObject)ts); 103 if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name; 104 else defaultType = ifun ? TSBEULER : TSEULER; 105 PetscCall(PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, 256, &opt)); 106 if (opt) PetscCall(TSSetType(ts, typeName)); 107 else PetscCall(TSSetType(ts, defaultType)); 108 109 /* Handle generic TS options */ 110 PetscCall(PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL)); 111 PetscCall(PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL)); 112 PetscCall(PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", tspan, &nt, &flg)); 113 if (flg) PetscCall(TSSetTimeSpan(ts, nt, tspan)); 114 PetscCall(PetscOptionsInt("-ts_max_steps", "Maximum number of time steps", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL)); 115 PetscCall(PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL)); 116 PetscCall(PetscOptionsReal("-ts_dt", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg)); 117 if (flg) PetscCall(TSSetTimeStep(ts, time_step)); 118 PetscCall(PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg)); 119 if (flg) PetscCall(TSSetExactFinalTime(ts, eftopt)); 120 PetscCall(PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, NULL)); 121 PetscCall(PetscOptionsInt("-ts_max_reject", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, NULL)); 122 PetscCall(PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL)); 123 PetscCall(PetscOptionsReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL)); 124 PetscCall(PetscOptionsReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL)); 125 126 PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL)); 127 PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult_transpose", "Test the RHS Jacobian transpose for consistency with RHS at each solve ", "None", ts->testjacobiantranspose, &ts->testjacobiantranspose, NULL)); 128 PetscCall(PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL)); 129 #if defined(PETSC_HAVE_SAWS) 130 { 131 PetscBool set; 132 flg = PETSC_FALSE; 133 PetscCall(PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set)); 134 if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)ts, flg)); 135 } 136 #endif 137 138 /* Monitor options */ 139 PetscCall(PetscOptionsInt("-ts_monitor_frequency", "Number of time steps between monitor output", "TSMonitorSetFrequency", ts->monitorFrequency, &ts->monitorFrequency, NULL)); 140 PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL)); 141 PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL)); 142 PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, NULL)); 143 PetscCall(TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL)); 144 145 PetscCall(PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg)); 146 if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)ts, monfilename)); 147 148 PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt)); 149 if (opt) { 150 PetscInt howoften = 1; 151 DM dm; 152 PetscBool net; 153 154 PetscCall(PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL)); 155 PetscCall(TSGetDM(ts, &dm)); 156 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net)); 157 if (net) { 158 TSMonitorLGCtxNetwork ctx; 159 PetscCall(TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx)); 160 PetscCall(TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxNetworkDestroy)); 161 PetscCall(PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL)); 162 } else { 163 TSMonitorLGCtx ctx; 164 PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx)); 165 PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy)); 166 } 167 } 168 169 PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt)); 170 if (opt) { 171 TSMonitorLGCtx ctx; 172 PetscInt howoften = 1; 173 174 PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL)); 175 PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx)); 176 PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy)); 177 } 178 PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL)); 179 180 PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt)); 181 if (opt) { 182 TSMonitorLGCtx ctx; 183 PetscInt howoften = 1; 184 185 PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL)); 186 PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx)); 187 PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy)); 188 } 189 PetscCall(PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt)); 190 if (opt) { 191 TSMonitorLGCtx ctx; 192 PetscInt howoften = 1; 193 194 PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL)); 195 PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx)); 196 PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy)); 197 ctx->semilogy = PETSC_TRUE; 198 } 199 200 PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt)); 201 if (opt) { 202 TSMonitorLGCtx ctx; 203 PetscInt howoften = 1; 204 205 PetscCall(PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL)); 206 PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx)); 207 PetscCall(TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy)); 208 } 209 PetscCall(PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt)); 210 if (opt) { 211 TSMonitorLGCtx ctx; 212 PetscInt howoften = 1; 213 214 PetscCall(PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL)); 215 PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx)); 216 PetscCall(TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscErrorCode(*)(void **))TSMonitorLGCtxDestroy)); 217 } 218 PetscCall(PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt)); 219 if (opt) { 220 TSMonitorSPEigCtx ctx; 221 PetscInt howoften = 1; 222 223 PetscCall(PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL)); 224 PetscCall(TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx)); 225 PetscCall(TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscErrorCode(*)(void **))TSMonitorSPEigCtxDestroy)); 226 } 227 PetscCall(PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase space from the DMSwarm", "TSMonitorSPSwarm", &opt)); 228 if (opt) { 229 TSMonitorSPCtx ctx; 230 PetscInt howoften = 1, retain = 0; 231 PetscBool phase = PETSC_TRUE, create = PETSC_TRUE, multispecies = PETSC_FALSE; 232 233 for (PetscInt i = 0; i < ts->numbermonitors; ++i) 234 if (ts->monitor[i] == TSMonitorSPSwarmSolution) { 235 create = PETSC_FALSE; 236 break; 237 } 238 if (create) { 239 PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase space from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL)); 240 PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL)); 241 PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL)); 242 PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_multi_species", "Color particles by particle species", "TSMonitorSPSwarm", multispecies, &multispecies, NULL)); 243 PetscCall(TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, multispecies, &ctx)); 244 PetscCall(TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorSPCtxDestroy)); 245 } 246 } 247 PetscCall(PetscOptionsName("-ts_monitor_hg_swarm", "Display particle histogram from the DMSwarm", "TSMonitorHGSwarm", &opt)); 248 if (opt) { 249 TSMonitorHGCtx ctx; 250 PetscInt howoften = 1, Ns = 1; 251 PetscBool velocity = PETSC_FALSE, create = PETSC_TRUE; 252 253 for (PetscInt i = 0; i < ts->numbermonitors; ++i) 254 if (ts->monitor[i] == TSMonitorHGSwarmSolution) { 255 create = PETSC_FALSE; 256 break; 257 } 258 if (create) { 259 DM sw, dm; 260 PetscInt Nc, Nb; 261 262 PetscCall(TSGetDM(ts, &sw)); 263 PetscCall(DMSwarmGetCellDM(sw, &dm)); 264 PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc)); 265 Nb = PetscMin(20, PetscMax(10, Nc)); 266 PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL)); 267 PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL)); 268 PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL)); 269 PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL)); 270 PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx)); 271 PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorHGCtxDestroy)); 272 } 273 } 274 opt = PETSC_FALSE; 275 PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt)); 276 if (opt) { 277 TSMonitorDrawCtx ctx; 278 PetscInt howoften = 1; 279 280 PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL)); 281 PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx)); 282 PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy)); 283 } 284 opt = PETSC_FALSE; 285 PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt)); 286 if (opt) { 287 TSMonitorDrawCtx ctx; 288 PetscReal bounds[4]; 289 PetscInt n = 4; 290 PetscDraw draw; 291 PetscDrawAxis axis; 292 293 PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL)); 294 PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field"); 295 PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx)); 296 PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw)); 297 PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis)); 298 PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3])); 299 PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2")); 300 PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy)); 301 } 302 opt = PETSC_FALSE; 303 PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt)); 304 if (opt) { 305 TSMonitorDrawCtx ctx; 306 PetscInt howoften = 1; 307 308 PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL)); 309 PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx)); 310 PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy)); 311 } 312 opt = PETSC_FALSE; 313 PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt)); 314 if (opt) { 315 TSMonitorDrawCtx ctx; 316 PetscInt howoften = 1; 317 318 PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL)); 319 PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx)); 320 PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscErrorCode(*)(void **))TSMonitorDrawCtxDestroy)); 321 } 322 323 opt = PETSC_FALSE; 324 PetscCall(PetscOptionsString("-ts_monitor_solution_vtk", "Save each time step to a binary file, use filename-%%03" PetscInt_FMT ".vts", "TSMonitorSolutionVTK", NULL, monfilename, sizeof(monfilename), &flg)); 325 if (flg) { 326 const char *ptr = NULL, *ptr2 = NULL; 327 char *filetemplate; 328 PetscCheck(monfilename[0], PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03" PetscInt_FMT ".vts"); 329 /* Do some cursory validation of the input. */ 330 PetscCall(PetscStrstr(monfilename, "%", (char **)&ptr)); 331 PetscCheck(ptr, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03" PetscInt_FMT ".vts"); 332 for (ptr++; ptr && *ptr; ptr++) { 333 PetscCall(PetscStrchr("DdiouxX", *ptr, (char **)&ptr2)); 334 PetscCheck(ptr2 || (*ptr >= '0' && *ptr <= '9'), PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Invalid file template argument to -ts_monitor_solution_vtk, should look like filename-%%03" PetscInt_FMT ".vts"); 335 if (ptr2) break; 336 } 337 PetscCall(PetscStrallocpy(monfilename, &filetemplate)); 338 PetscCall(TSMonitorSet(ts, TSMonitorSolutionVTK, filetemplate, (PetscErrorCode(*)(void **))TSMonitorSolutionVTKDestroy)); 339 } 340 341 PetscCall(PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg)); 342 if (flg) { 343 TSMonitorDMDARayCtx *rayctx; 344 int ray = 0; 345 DMDirection ddir; 346 DM da; 347 PetscMPIInt rank; 348 349 PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir); 350 if (dir[0] == 'x') ddir = DM_X; 351 else if (dir[0] == 'y') ddir = DM_Y; 352 else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir); 353 sscanf(dir + 2, "%d", &ray); 354 355 PetscCall(PetscInfo(((PetscObject)ts), "Displaying DMDA ray %c = %d\n", dir[0], ray)); 356 PetscCall(PetscNew(&rayctx)); 357 PetscCall(TSGetDM(ts, &da)); 358 PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter)); 359 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank)); 360 if (rank == 0) PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer)); 361 rayctx->lgctx = NULL; 362 PetscCall(TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy)); 363 } 364 PetscCall(PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg)); 365 if (flg) { 366 TSMonitorDMDARayCtx *rayctx; 367 int ray = 0; 368 DMDirection ddir; 369 DM da; 370 PetscInt howoften = 1; 371 372 PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir); 373 if (dir[0] == 'x') ddir = DM_X; 374 else if (dir[0] == 'y') ddir = DM_Y; 375 else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir); 376 sscanf(dir + 2, "%d", &ray); 377 378 PetscCall(PetscInfo(((PetscObject)ts), "Displaying LG DMDA ray %c = %d\n", dir[0], ray)); 379 PetscCall(PetscNew(&rayctx)); 380 PetscCall(TSGetDM(ts, &da)); 381 PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter)); 382 PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx)); 383 PetscCall(TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy)); 384 } 385 386 PetscCall(PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt)); 387 if (opt) { 388 TSMonitorEnvelopeCtx ctx; 389 390 PetscCall(TSMonitorEnvelopeCtxCreate(ts, &ctx)); 391 PetscCall(TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscErrorCode(*)(void **))TSMonitorEnvelopeCtxDestroy)); 392 } 393 flg = PETSC_FALSE; 394 PetscCall(PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt)); 395 if (opt && flg) PetscCall(TSMonitorCancel(ts)); 396 397 flg = PETSC_FALSE; 398 PetscCall(PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeIJacobianDefaultColor", flg, &flg, NULL)); 399 if (flg) { 400 DM dm; 401 402 PetscCall(TSGetDM(ts, &dm)); 403 PetscCall(DMTSUnsetIJacobianContext_Internal(dm)); 404 PetscCall(TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL)); 405 PetscCall(PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n")); 406 } 407 408 /* Handle specific TS options */ 409 PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject); 410 411 /* Handle TSAdapt options */ 412 PetscCall(TSGetAdapt(ts, &ts->adapt)); 413 PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type)); 414 PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject)); 415 416 /* TS trajectory must be set after TS, since it may use some TS options above */ 417 tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE; 418 PetscCall(PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL)); 419 if (tflg) PetscCall(TSSetSaveTrajectory(ts)); 420 421 PetscCall(TSAdjointSetFromOptions(ts, PetscOptionsObject)); 422 423 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 424 PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject)); 425 PetscOptionsEnd(); 426 427 if (ts->trajectory) PetscCall(TSTrajectorySetFromOptions(ts->trajectory, ts)); 428 429 /* why do we have to do this here and not during TSSetUp? */ 430 PetscCall(TSGetSNES(ts, &ts->snes)); 431 if (ts->problem_type == TS_LINEAR) { 432 PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, "")); 433 if (!flg) PetscCall(SNESSetType(ts->snes, SNESKSPONLY)); 434 } 435 PetscCall(SNESSetFromOptions(ts->snes)); 436 PetscFunctionReturn(PETSC_SUCCESS); 437 } 438 439 /*@ 440 TSGetTrajectory - Gets the trajectory from a `TS` if it exists 441 442 Collective 443 444 Input Parameter: 445 . ts - the `TS` context obtained from `TSCreate()` 446 447 Output Parameter: 448 . tr - the `TSTrajectory` object, if it exists 449 450 Level: advanced 451 452 Note: 453 This routine should be called after all `TS` options have been set 454 455 .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectory`, `TSTrajectoryCreate()` 456 @*/ 457 PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr) 458 { 459 PetscFunctionBegin; 460 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 461 *tr = ts->trajectory; 462 PetscFunctionReturn(PETSC_SUCCESS); 463 } 464 465 /*@ 466 TSSetSaveTrajectory - Causes the `TS` to save its solutions as it iterates forward in time in a `TSTrajectory` object 467 468 Collective 469 470 Input Parameter: 471 . ts - the `TS` context obtained from `TSCreate()` 472 473 Options Database Keys: 474 + -ts_save_trajectory - saves the trajectory to a file 475 - -ts_trajectory_type type - set trajectory type 476 477 Level: intermediate 478 479 Notes: 480 This routine should be called after all `TS` options have been set 481 482 The `TSTRAJECTORYVISUALIZATION` files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and 483 MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m 484 485 .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()` 486 @*/ 487 PetscErrorCode TSSetSaveTrajectory(TS ts) 488 { 489 PetscFunctionBegin; 490 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 491 if (!ts->trajectory) PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory)); 492 PetscFunctionReturn(PETSC_SUCCESS); 493 } 494 495 /*@ 496 TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object 497 498 Collective 499 500 Input Parameter: 501 . ts - the `TS` context obtained from `TSCreate()` 502 503 Level: intermediate 504 505 .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()` 506 @*/ 507 PetscErrorCode TSResetTrajectory(TS ts) 508 { 509 PetscFunctionBegin; 510 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 511 if (ts->trajectory) { 512 PetscCall(TSTrajectoryDestroy(&ts->trajectory)); 513 PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory)); 514 } 515 PetscFunctionReturn(PETSC_SUCCESS); 516 } 517 518 /*@ 519 TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from a `TS` 520 521 Collective 522 523 Input Parameter: 524 . ts - the `TS` context obtained from `TSCreate()` 525 526 Level: intermediate 527 528 .seealso: [](ch_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()` 529 @*/ 530 PetscErrorCode TSRemoveTrajectory(TS ts) 531 { 532 PetscFunctionBegin; 533 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 534 if (ts->trajectory) PetscCall(TSTrajectoryDestroy(&ts->trajectory)); 535 PetscFunctionReturn(PETSC_SUCCESS); 536 } 537 538 /*@ 539 TSComputeRHSJacobian - Computes the Jacobian matrix that has been 540 set with `TSSetRHSJacobian()`. 541 542 Collective 543 544 Input Parameters: 545 + ts - the `TS` context 546 . t - current timestep 547 - U - input vector 548 549 Output Parameters: 550 + A - Jacobian matrix 551 - B - optional preconditioning matrix 552 553 Level: developer 554 555 Note: 556 Most users should not need to explicitly call this routine, as it 557 is used internally within the nonlinear solvers. 558 559 .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()` 560 @*/ 561 PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B) 562 { 563 PetscObjectState Ustate; 564 PetscObjectId Uid; 565 DM dm; 566 DMTS tsdm; 567 TSRHSJacobian rhsjacobianfunc; 568 void *ctx; 569 TSRHSFunction rhsfunction; 570 571 PetscFunctionBegin; 572 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 573 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 574 PetscCheckSameComm(ts, 1, U, 3); 575 PetscCall(TSGetDM(ts, &dm)); 576 PetscCall(DMGetDMTS(dm, &tsdm)); 577 PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL)); 578 PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx)); 579 PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate)); 580 PetscCall(PetscObjectGetId((PetscObject)U, &Uid)); 581 582 if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) PetscFunctionReturn(PETSC_SUCCESS); 583 584 PetscCheck(ts->rhsjacobian.shift == 0.0 || !ts->rhsjacobian.reuse, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Should not call TSComputeRHSJacobian() on a shifted matrix (shift=%lf) when RHSJacobian is reusable.", (double)ts->rhsjacobian.shift); 585 if (rhsjacobianfunc) { 586 PetscCall(PetscLogEventBegin(TS_JacobianEval, ts, U, A, B)); 587 PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx)); 588 ts->rhsjacs++; 589 PetscCall(PetscLogEventEnd(TS_JacobianEval, ts, U, A, B)); 590 } else { 591 PetscCall(MatZeroEntries(A)); 592 if (B && A != B) PetscCall(MatZeroEntries(B)); 593 } 594 ts->rhsjacobian.time = t; 595 ts->rhsjacobian.shift = 0; 596 ts->rhsjacobian.scale = 1.; 597 PetscCall(PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid)); 598 PetscCall(PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate)); 599 PetscFunctionReturn(PETSC_SUCCESS); 600 } 601 602 /*@ 603 TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS` 604 605 Collective 606 607 Input Parameters: 608 + ts - the `TS` context 609 . t - current time 610 - U - state vector 611 612 Output Parameter: 613 . y - right hand side 614 615 Level: developer 616 617 Note: 618 Most users should not need to explicitly call this routine, as it 619 is used internally within the nonlinear solvers. 620 621 .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()` 622 @*/ 623 PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y) 624 { 625 TSRHSFunction rhsfunction; 626 TSIFunction ifunction; 627 void *ctx; 628 DM dm; 629 630 PetscFunctionBegin; 631 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 632 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 633 PetscValidHeaderSpecific(y, VEC_CLASSID, 4); 634 PetscCall(TSGetDM(ts, &dm)); 635 PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx)); 636 PetscCall(DMTSGetIFunction(dm, &ifunction, NULL)); 637 638 PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()"); 639 640 if (rhsfunction) { 641 PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, y, 0)); 642 PetscCall(VecLockReadPush(U)); 643 PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx)); 644 PetscCall(VecLockReadPop(U)); 645 ts->rhsfuncs++; 646 PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, y, 0)); 647 } else PetscCall(VecZeroEntries(y)); 648 PetscFunctionReturn(PETSC_SUCCESS); 649 } 650 651 /*@ 652 TSComputeSolutionFunction - Evaluates the solution function. 653 654 Collective 655 656 Input Parameters: 657 + ts - the `TS` context 658 - t - current time 659 660 Output Parameter: 661 . U - the solution 662 663 Level: developer 664 665 .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()` 666 @*/ 667 PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U) 668 { 669 TSSolutionFunction solutionfunction; 670 void *ctx; 671 DM dm; 672 673 PetscFunctionBegin; 674 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 675 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 676 PetscCall(TSGetDM(ts, &dm)); 677 PetscCall(DMTSGetSolutionFunction(dm, &solutionfunction, &ctx)); 678 if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx)); 679 PetscFunctionReturn(PETSC_SUCCESS); 680 } 681 /*@ 682 TSComputeForcingFunction - Evaluates the forcing function. 683 684 Collective 685 686 Input Parameters: 687 + ts - the `TS` context 688 - t - current time 689 690 Output Parameter: 691 . U - the function value 692 693 Level: developer 694 695 .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()` 696 @*/ 697 PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U) 698 { 699 void *ctx; 700 DM dm; 701 TSForcingFunction forcing; 702 703 PetscFunctionBegin; 704 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 705 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 706 PetscCall(TSGetDM(ts, &dm)); 707 PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx)); 708 709 if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx)); 710 PetscFunctionReturn(PETSC_SUCCESS); 711 } 712 713 static PetscErrorCode TSGetRHSVec_Private(TS ts, Vec *Frhs) 714 { 715 Vec F; 716 717 PetscFunctionBegin; 718 *Frhs = NULL; 719 PetscCall(TSGetIFunction(ts, &F, NULL, NULL)); 720 if (!ts->Frhs) PetscCall(VecDuplicate(F, &ts->Frhs)); 721 *Frhs = ts->Frhs; 722 PetscFunctionReturn(PETSC_SUCCESS); 723 } 724 725 PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs) 726 { 727 Mat A, B; 728 TSIJacobian ijacobian; 729 730 PetscFunctionBegin; 731 if (Arhs) *Arhs = NULL; 732 if (Brhs) *Brhs = NULL; 733 PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL)); 734 if (Arhs) { 735 if (!ts->Arhs) { 736 if (ijacobian) { 737 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs)); 738 PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN)); 739 } else { 740 ts->Arhs = A; 741 PetscCall(PetscObjectReference((PetscObject)A)); 742 } 743 } else { 744 PetscBool flg; 745 PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg)); 746 /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */ 747 if (flg && !ijacobian && ts->Arhs == ts->Brhs) { 748 PetscCall(PetscObjectDereference((PetscObject)ts->Arhs)); 749 ts->Arhs = A; 750 PetscCall(PetscObjectReference((PetscObject)A)); 751 } 752 } 753 *Arhs = ts->Arhs; 754 } 755 if (Brhs) { 756 if (!ts->Brhs) { 757 if (A != B) { 758 if (ijacobian) { 759 PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs)); 760 } else { 761 ts->Brhs = B; 762 PetscCall(PetscObjectReference((PetscObject)B)); 763 } 764 } else { 765 PetscCall(PetscObjectReference((PetscObject)ts->Arhs)); 766 ts->Brhs = ts->Arhs; 767 } 768 } 769 *Brhs = ts->Brhs; 770 } 771 PetscFunctionReturn(PETSC_SUCCESS); 772 } 773 774 /*@ 775 TSComputeIFunction - Evaluates the DAE residual written in the implicit form F(t,U,Udot)=0 776 777 Collective 778 779 Input Parameters: 780 + ts - the `TS` context 781 . t - current time 782 . U - state vector 783 . Udot - time derivative of state vector 784 - imex - flag indicates if the method is `TSIMEX` so that the RHSFunction should be kept separate 785 786 Output Parameter: 787 . Y - right hand side 788 789 Level: developer 790 791 Note: 792 Most users should not need to explicitly call this routine, as it 793 is used internally within the nonlinear solvers. 794 795 If the user did did not write their equations in implicit form, this 796 function recasts them in implicit form. 797 798 .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()` 799 @*/ 800 PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex) 801 { 802 TSIFunction ifunction; 803 TSRHSFunction rhsfunction; 804 void *ctx; 805 DM dm; 806 807 PetscFunctionBegin; 808 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 809 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 810 PetscValidHeaderSpecific(Udot, VEC_CLASSID, 4); 811 PetscValidHeaderSpecific(Y, VEC_CLASSID, 5); 812 813 PetscCall(TSGetDM(ts, &dm)); 814 PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx)); 815 PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL)); 816 817 PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()"); 818 819 PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, Udot, Y)); 820 if (ifunction) { 821 PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx)); 822 ts->ifuncs++; 823 } 824 if (imex) { 825 if (!ifunction) PetscCall(VecCopy(Udot, Y)); 826 } else if (rhsfunction) { 827 if (ifunction) { 828 Vec Frhs; 829 PetscCall(TSGetRHSVec_Private(ts, &Frhs)); 830 PetscCall(TSComputeRHSFunction(ts, t, U, Frhs)); 831 PetscCall(VecAXPY(Y, -1, Frhs)); 832 } else { 833 PetscCall(TSComputeRHSFunction(ts, t, U, Y)); 834 PetscCall(VecAYPX(Y, -1, Udot)); 835 } 836 } 837 PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, Udot, Y)); 838 PetscFunctionReturn(PETSC_SUCCESS); 839 } 840 841 /* 842 TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call `TSComputeRHSJacobian()` on it. 843 844 Note: 845 This routine is needed when one switches from `TSComputeIJacobian()` to `TSComputeRHSJacobian()` because the Jacobian matrix may be shifted or scaled in `TSComputeIJacobian()`. 846 847 */ 848 static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B) 849 { 850 PetscFunctionBegin; 851 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 852 PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat"); 853 PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat"); 854 855 if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift)); 856 if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1)); 857 if (B && B == ts->Brhs && A != B) { 858 if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift)); 859 if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1)); 860 } 861 ts->rhsjacobian.shift = 0; 862 ts->rhsjacobian.scale = 1.; 863 PetscFunctionReturn(PETSC_SUCCESS); 864 } 865 866 /*@ 867 TSComputeIJacobian - Evaluates the Jacobian of the DAE 868 869 Collective 870 871 Input 872 Input Parameters: 873 + ts - the `TS` context 874 . t - current timestep 875 . U - state vector 876 . Udot - time derivative of state vector 877 . shift - shift to apply, see note below 878 - imex - flag indicates if the method is `TSIMEX` so that the RHSJacobian should be kept separate 879 880 Output Parameters: 881 + A - Jacobian matrix 882 - B - matrix from which the preconditioner is constructed; often the same as `A` 883 884 Level: developer 885 886 Notes: 887 If F(t,U,Udot)=0 is the DAE, the required Jacobian is 888 .vb 889 dF/dU + shift*dF/dUdot 890 .ve 891 Most users should not need to explicitly call this routine, as it 892 is used internally within the nonlinear solvers. 893 894 .seealso: [](ch_ts), `TS`, `TSSetIJacobian()` 895 @*/ 896 PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex) 897 { 898 TSIJacobian ijacobian; 899 TSRHSJacobian rhsjacobian; 900 DM dm; 901 void *ctx; 902 903 PetscFunctionBegin; 904 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 905 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 906 PetscValidHeaderSpecific(Udot, VEC_CLASSID, 4); 907 PetscValidPointer(A, 6); 908 PetscValidHeaderSpecific(A, MAT_CLASSID, 6); 909 PetscValidPointer(B, 7); 910 PetscValidHeaderSpecific(B, MAT_CLASSID, 7); 911 912 PetscCall(TSGetDM(ts, &dm)); 913 PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx)); 914 PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL)); 915 916 PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()"); 917 918 PetscCall(PetscLogEventBegin(TS_JacobianEval, ts, U, A, B)); 919 if (ijacobian) { 920 PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx)); 921 ts->ijacs++; 922 } 923 if (imex) { 924 if (!ijacobian) { /* system was written as Udot = G(t,U) */ 925 PetscBool assembled; 926 if (rhsjacobian) { 927 Mat Arhs = NULL; 928 PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL)); 929 if (A == Arhs) { 930 PetscCheck(rhsjacobian != TSComputeRHSJacobianConstant, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Unsupported operation! cannot use TSComputeRHSJacobianConstant"); /* there is no way to reconstruct shift*M-J since J cannot be reevaluated */ 931 ts->rhsjacobian.time = PETSC_MIN_REAL; 932 } 933 } 934 PetscCall(MatZeroEntries(A)); 935 PetscCall(MatAssembled(A, &assembled)); 936 if (!assembled) { 937 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 938 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 939 } 940 PetscCall(MatShift(A, shift)); 941 if (A != B) { 942 PetscCall(MatZeroEntries(B)); 943 PetscCall(MatAssembled(B, &assembled)); 944 if (!assembled) { 945 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 946 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 947 } 948 PetscCall(MatShift(B, shift)); 949 } 950 } 951 } else { 952 Mat Arhs = NULL, Brhs = NULL; 953 954 /* RHSJacobian needs to be converted to part of IJacobian if exists */ 955 if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs)); 956 if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */ 957 PetscObjectState Ustate; 958 PetscObjectId Uid; 959 TSRHSFunction rhsfunction; 960 961 PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL)); 962 PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate)); 963 PetscCall(PetscObjectGetId((PetscObject)U, &Uid)); 964 if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) && 965 ts->rhsjacobian.scale == -1.) { /* No need to recompute RHSJacobian */ 966 PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */ 967 if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift)); 968 } else { 969 PetscBool flg; 970 971 if (ts->rhsjacobian.reuse) { /* Undo the damage */ 972 /* MatScale has a short path for this case. 973 However, this code path is taken the first time TSComputeRHSJacobian is called 974 and the matrices have not been assembled yet */ 975 PetscCall(TSRecoverRHSJacobian(ts, A, B)); 976 } 977 PetscCall(TSComputeRHSJacobian(ts, t, U, A, B)); 978 PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg)); 979 /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */ 980 if (!flg) { 981 PetscCall(MatScale(A, -1)); 982 PetscCall(MatShift(A, shift)); 983 } 984 if (A != B) { 985 PetscCall(MatScale(B, -1)); 986 PetscCall(MatShift(B, shift)); 987 } 988 } 989 ts->rhsjacobian.scale = -1; 990 ts->rhsjacobian.shift = shift; 991 } else if (Arhs) { /* Both IJacobian and RHSJacobian */ 992 if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */ 993 PetscCall(MatZeroEntries(A)); 994 PetscCall(MatShift(A, shift)); 995 if (A != B) { 996 PetscCall(MatZeroEntries(B)); 997 PetscCall(MatShift(B, shift)); 998 } 999 } 1000 PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs)); 1001 PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern)); 1002 if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern)); 1003 } 1004 } 1005 PetscCall(PetscLogEventEnd(TS_JacobianEval, ts, U, A, B)); 1006 PetscFunctionReturn(PETSC_SUCCESS); 1007 } 1008 1009 /*@C 1010 TSSetRHSFunction - Sets the routine for evaluating the function, 1011 where U_t = G(t,u). 1012 1013 Logically Collective 1014 1015 Input Parameters: 1016 + ts - the `TS` context obtained from `TSCreate()` 1017 . r - vector to put the computed right hand side (or `NULL` to have it created) 1018 . f - routine for evaluating the right-hand-side function 1019 - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`) 1020 1021 Calling sequence of `f`: 1022 $ PetscErrorCode f(TS ts, PetscReal t, Vec u, Vec F, void *ctx) 1023 + ts - timestep context 1024 . t - current timestep 1025 . u - input vector 1026 . F - function vector 1027 - ctx - [optional] user-defined function context 1028 1029 Level: beginner 1030 1031 Note: 1032 You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE. 1033 1034 .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()` 1035 @*/ 1036 PetscErrorCode TSSetRHSFunction(TS ts, Vec r, PetscErrorCode (*f)(TS, PetscReal, Vec, Vec, void *), void *ctx) 1037 { 1038 SNES snes; 1039 Vec ralloc = NULL; 1040 DM dm; 1041 1042 PetscFunctionBegin; 1043 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1044 if (r) PetscValidHeaderSpecific(r, VEC_CLASSID, 2); 1045 1046 PetscCall(TSGetDM(ts, &dm)); 1047 PetscCall(DMTSSetRHSFunction(dm, f, ctx)); 1048 PetscCall(TSGetSNES(ts, &snes)); 1049 if (!r && !ts->dm && ts->vec_sol) { 1050 PetscCall(VecDuplicate(ts->vec_sol, &ralloc)); 1051 r = ralloc; 1052 } 1053 PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts)); 1054 PetscCall(VecDestroy(&ralloc)); 1055 PetscFunctionReturn(PETSC_SUCCESS); 1056 } 1057 1058 /*@C 1059 TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE 1060 1061 Logically Collective 1062 1063 Input Parameters: 1064 + ts - the `TS` context obtained from `TSCreate()` 1065 . f - routine for evaluating the solution 1066 - ctx - [optional] user-defined context for private data for the 1067 function evaluation routine (may be `NULL`) 1068 1069 Calling sequence of `f`: 1070 $ PetscErrorCode f(TS ts, PetscReal t, Vec u, void *ctx) 1071 + ts - the integrator context 1072 . t - current timestep 1073 . u - output vector 1074 - ctx - [optional] user-defined function context 1075 1076 Options Database Keys: 1077 + -ts_monitor_lg_error - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()` 1078 - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()` 1079 1080 Level: intermediate 1081 1082 Notes: 1083 This routine is used for testing accuracy of time integration schemes when you already know the solution. 1084 If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to 1085 create closed-form solutions with non-physical forcing terms. 1086 1087 For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history. 1088 1089 .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()` 1090 @*/ 1091 PetscErrorCode TSSetSolutionFunction(TS ts, PetscErrorCode (*f)(TS, PetscReal, Vec, void *), void *ctx) 1092 { 1093 DM dm; 1094 1095 PetscFunctionBegin; 1096 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1097 PetscCall(TSGetDM(ts, &dm)); 1098 PetscCall(DMTSSetSolutionFunction(dm, f, ctx)); 1099 PetscFunctionReturn(PETSC_SUCCESS); 1100 } 1101 1102 /*@C 1103 TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE 1104 1105 Logically Collective 1106 1107 Input Parameters: 1108 + ts - the `TS` context obtained from `TSCreate()` 1109 . func - routine for evaluating the forcing function 1110 - ctx - [optional] user-defined context for private data for the 1111 function evaluation routine (may be `NULL`) 1112 1113 Calling sequence of func: 1114 $ PetscErrorCode func(TS ts, PetscReal t, Vec f, void *ctx) 1115 + ts - the integrator context 1116 . t - current timestep 1117 . f - output vector 1118 - ctx - [optional] user-defined function context 1119 1120 Level: intermediate 1121 1122 Notes: 1123 This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to 1124 create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the 1125 definition of the problem you are solving and hence possibly introducing bugs. 1126 1127 This replaces the ODE F(u,u_t,t) = 0 the TS is solving with F(u,u_t,t) - func(t) = 0 1128 1129 This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the 1130 parameters can be passed in the ctx variable. 1131 1132 For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history. 1133 1134 .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetSolutionFunction()` 1135 @*/ 1136 PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFunction func, void *ctx) 1137 { 1138 DM dm; 1139 1140 PetscFunctionBegin; 1141 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1142 PetscCall(TSGetDM(ts, &dm)); 1143 PetscCall(DMTSSetForcingFunction(dm, func, ctx)); 1144 PetscFunctionReturn(PETSC_SUCCESS); 1145 } 1146 1147 /*@C 1148 TSSetRHSJacobian - Sets the function to compute the Jacobian of G, 1149 where U_t = G(U,t), as well as the location to store the matrix. 1150 1151 Logically Collective 1152 1153 Input Parameters: 1154 + ts - the `TS` context obtained from `TSCreate()` 1155 . Amat - (approximate) location to store Jacobian matrix entries computed by `f` 1156 . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`) 1157 . f - the Jacobian evaluation routine 1158 - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`) 1159 1160 Calling sequence of `f`: 1161 $ PetscErrorCode f(TS ts, PetscReal t, Vec u, Mat A, Mat B, void *ctx) 1162 + ts - the `TS` context obtained from `TSCreate()` 1163 . t - current timestep 1164 . u - input vector 1165 . Amat - (approximate) Jacobian matrix 1166 . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`) 1167 - ctx - [optional] user-defined context for matrix evaluation routine 1168 1169 Level: beginner 1170 1171 Notes: 1172 You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value 1173 1174 The `TS` solver may modify the nonzero structure and the entries of the matrices Amat and Pmat between the calls to f() 1175 You should not assume the values are the same in the next call to f() as you set them in the previous call. 1176 1177 .seealso: [](ch_ts), `TS`, `SNESComputeJacobianDefaultColor()`, `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()` 1178 @*/ 1179 PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobian f, void *ctx) 1180 { 1181 SNES snes; 1182 DM dm; 1183 TSIJacobian ijacobian; 1184 1185 PetscFunctionBegin; 1186 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1187 if (Amat) PetscValidHeaderSpecific(Amat, MAT_CLASSID, 2); 1188 if (Pmat) PetscValidHeaderSpecific(Pmat, MAT_CLASSID, 3); 1189 if (Amat) PetscCheckSameComm(ts, 1, Amat, 2); 1190 if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3); 1191 1192 PetscCall(TSGetDM(ts, &dm)); 1193 PetscCall(DMTSSetRHSJacobian(dm, f, ctx)); 1194 PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL)); 1195 PetscCall(TSGetSNES(ts, &snes)); 1196 if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts)); 1197 if (Amat) { 1198 PetscCall(PetscObjectReference((PetscObject)Amat)); 1199 PetscCall(MatDestroy(&ts->Arhs)); 1200 ts->Arhs = Amat; 1201 } 1202 if (Pmat) { 1203 PetscCall(PetscObjectReference((PetscObject)Pmat)); 1204 PetscCall(MatDestroy(&ts->Brhs)); 1205 ts->Brhs = Pmat; 1206 } 1207 PetscFunctionReturn(PETSC_SUCCESS); 1208 } 1209 1210 /*@C 1211 TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved. 1212 1213 Logically Collective 1214 1215 Input Parameters: 1216 + ts - the `TS` context obtained from `TSCreate()` 1217 . r - vector to hold the residual (or `NULL` to have it created internally) 1218 . f - the function evaluation routine 1219 - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`) 1220 1221 Calling sequence of `f`: 1222 $ PetscErrorCode f(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx) 1223 + ts - the `TS` context obtained from `TSCreate()` 1224 . t - time at step/stage being solved 1225 . u - state vector 1226 . u_t - time derivative of state vector 1227 . F - function vector 1228 - ctx - [optional] user-defined context for matrix evaluation routine 1229 1230 Level: beginner 1231 1232 Note: 1233 The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE. When solving DAEs you must use this function. 1234 1235 .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`, `TSSetIJacobian()` 1236 @*/ 1237 PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunction f, void *ctx) 1238 { 1239 SNES snes; 1240 Vec ralloc = NULL; 1241 DM dm; 1242 1243 PetscFunctionBegin; 1244 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1245 if (r) PetscValidHeaderSpecific(r, VEC_CLASSID, 2); 1246 1247 PetscCall(TSGetDM(ts, &dm)); 1248 PetscCall(DMTSSetIFunction(dm, f, ctx)); 1249 1250 PetscCall(TSGetSNES(ts, &snes)); 1251 if (!r && !ts->dm && ts->vec_sol) { 1252 PetscCall(VecDuplicate(ts->vec_sol, &ralloc)); 1253 r = ralloc; 1254 } 1255 PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts)); 1256 PetscCall(VecDestroy(&ralloc)); 1257 PetscFunctionReturn(PETSC_SUCCESS); 1258 } 1259 1260 /*@C 1261 TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it. 1262 1263 Not Collective 1264 1265 Input Parameter: 1266 . ts - the `TS` context 1267 1268 Output Parameters: 1269 + r - vector to hold residual (or `NULL`) 1270 . func - the function to compute residual (or `NULL`) 1271 - ctx - the function context (or `NULL`) 1272 1273 Level: advanced 1274 1275 .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()` 1276 @*/ 1277 PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunction *func, void **ctx) 1278 { 1279 SNES snes; 1280 DM dm; 1281 1282 PetscFunctionBegin; 1283 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1284 PetscCall(TSGetSNES(ts, &snes)); 1285 PetscCall(SNESGetFunction(snes, r, NULL, NULL)); 1286 PetscCall(TSGetDM(ts, &dm)); 1287 PetscCall(DMTSGetIFunction(dm, func, ctx)); 1288 PetscFunctionReturn(PETSC_SUCCESS); 1289 } 1290 1291 /*@C 1292 TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it. 1293 1294 Not Collective 1295 1296 Input Parameter: 1297 . ts - the `TS` context 1298 1299 Output Parameters: 1300 + r - vector to hold computed right hand side (or `NULL`) 1301 . func - the function to compute right hand side (or `NULL`) 1302 - ctx - the function context (or `NULL`) 1303 1304 Level: advanced 1305 1306 .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()` 1307 @*/ 1308 PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunction *func, void **ctx) 1309 { 1310 SNES snes; 1311 DM dm; 1312 1313 PetscFunctionBegin; 1314 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1315 PetscCall(TSGetSNES(ts, &snes)); 1316 PetscCall(SNESGetFunction(snes, r, NULL, NULL)); 1317 PetscCall(TSGetDM(ts, &dm)); 1318 PetscCall(DMTSGetRHSFunction(dm, func, ctx)); 1319 PetscFunctionReturn(PETSC_SUCCESS); 1320 } 1321 1322 /*@C 1323 TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function 1324 provided with `TSSetIFunction()`. 1325 1326 Logically Collective 1327 1328 Input Parameters: 1329 + ts - the `TS` context obtained from `TSCreate()` 1330 . Amat - (approximate) matrix to store Jacobian entries computed by `f` 1331 . Pmat - matrix used to compute preconditioner (usually the same as `Amat`) 1332 . f - the Jacobian evaluation routine 1333 - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`) 1334 1335 Calling sequence of `f`: 1336 $ PetscErrorCode f(TS ts, PetscReal t, Vec U, Vec U_t, PetscReal a, Mat Amat, Mat Pmat, void *ctx) 1337 + ts - the `TS` context obtained from `TSCreate()` 1338 . t - time at step/stage being solved 1339 . U - state vector 1340 . U_t - time derivative of state vector 1341 . a - shift 1342 . Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t 1343 . Pmat - matrix used for constructing preconditioner, usually the same as `Amat` 1344 - ctx - [optional] user-defined context for matrix evaluation routine 1345 1346 Level: beginner 1347 1348 Notes: 1349 The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve. 1350 1351 If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null 1352 space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process. 1353 1354 The matrix dF/dU + a*dF/dU_t you provide turns out to be 1355 the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved. 1356 The time integrator internally approximates U_t by W+a*U where the positive "shift" 1357 a and vector W depend on the integration method, step size, and past states. For example with 1358 the backward Euler method a = 1/dt and W = -a*U(previous timestep) so 1359 W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt 1360 1361 You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value 1362 1363 The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f` 1364 You should not assume the values are the same in the next call to `f` as you set them in the previous call. 1365 1366 .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetRHSJacobian()`, `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()` 1367 @*/ 1368 PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobian f, void *ctx) 1369 { 1370 SNES snes; 1371 DM dm; 1372 1373 PetscFunctionBegin; 1374 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1375 if (Amat) PetscValidHeaderSpecific(Amat, MAT_CLASSID, 2); 1376 if (Pmat) PetscValidHeaderSpecific(Pmat, MAT_CLASSID, 3); 1377 if (Amat) PetscCheckSameComm(ts, 1, Amat, 2); 1378 if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3); 1379 1380 PetscCall(TSGetDM(ts, &dm)); 1381 PetscCall(DMTSSetIJacobian(dm, f, ctx)); 1382 1383 PetscCall(TSGetSNES(ts, &snes)); 1384 PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts)); 1385 PetscFunctionReturn(PETSC_SUCCESS); 1386 } 1387 1388 /*@ 1389 TSRHSJacobianSetReuse - restore RHS Jacobian before re-evaluating. Without this flag, `TS` will change the sign and 1390 shift the RHS Jacobian for a finite-time-step implicit solve, in which case the user function will need to recompute 1391 the entire Jacobian. The reuse flag must be set if the evaluation function will assume that the matrix entries have 1392 not been changed by the `TS`. 1393 1394 Logically Collective 1395 1396 Input Parameters: 1397 + ts - `TS` context obtained from `TSCreate()` 1398 - reuse - `PETSC_TRUE` if the RHS Jacobian 1399 1400 Level: intermediate 1401 1402 .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()` 1403 @*/ 1404 PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse) 1405 { 1406 PetscFunctionBegin; 1407 ts->rhsjacobian.reuse = reuse; 1408 PetscFunctionReturn(PETSC_SUCCESS); 1409 } 1410 1411 /*@C 1412 TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved. 1413 1414 Logically Collective 1415 1416 Input Parameters: 1417 + ts - the `TS` context obtained from `TSCreate()` 1418 . F - vector to hold the residual (or `NULL` to have it created internally) 1419 . fun - the function evaluation routine 1420 - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`) 1421 1422 Calling sequence of `fun`: 1423 $ PetscErrorCode fun(TS ts, PetscReal t, Vec U, Vec U_t, Vec U_tt, Vec F,void *ctx); 1424 + ts - the `TS` context obtained from `TSCreate()` 1425 . t - time at step/stage being solved 1426 . U - state vector 1427 . U_t - time derivative of state vector 1428 . U_tt - second time derivative of state vector 1429 . F - function vector 1430 - ctx - [optional] user-defined context for matrix evaluation routine (may be `NULL`) 1431 1432 Level: beginner 1433 1434 .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`, `TSSetIFunction()`, `TSCreate()`, `TSSetRHSFunction()` 1435 @*/ 1436 PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2Function fun, void *ctx) 1437 { 1438 DM dm; 1439 1440 PetscFunctionBegin; 1441 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1442 if (F) PetscValidHeaderSpecific(F, VEC_CLASSID, 2); 1443 PetscCall(TSSetIFunction(ts, F, NULL, NULL)); 1444 PetscCall(TSGetDM(ts, &dm)); 1445 PetscCall(DMTSSetI2Function(dm, fun, ctx)); 1446 PetscFunctionReturn(PETSC_SUCCESS); 1447 } 1448 1449 /*@C 1450 TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it. 1451 1452 Not Collective 1453 1454 Input Parameter: 1455 . ts - the `TS` context 1456 1457 Output Parameters: 1458 + r - vector to hold residual (or `NULL`) 1459 . fun - the function to compute residual (or `NULL`) 1460 - ctx - the function context (or `NULL`) 1461 1462 Level: advanced 1463 1464 .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()` 1465 @*/ 1466 PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2Function *fun, void **ctx) 1467 { 1468 SNES snes; 1469 DM dm; 1470 1471 PetscFunctionBegin; 1472 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1473 PetscCall(TSGetSNES(ts, &snes)); 1474 PetscCall(SNESGetFunction(snes, r, NULL, NULL)); 1475 PetscCall(TSGetDM(ts, &dm)); 1476 PetscCall(DMTSGetI2Function(dm, fun, ctx)); 1477 PetscFunctionReturn(PETSC_SUCCESS); 1478 } 1479 1480 /*@C 1481 TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t + a*dF/dU_tt 1482 where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`. 1483 1484 Logically Collective 1485 1486 Input Parameters: 1487 + ts - the `TS` context obtained from `TSCreate()` 1488 . J - matrix to hold the Jacobian values 1489 . P - matrix for constructing the preconditioner (may be same as `J`) 1490 . jac - the Jacobian evaluation routine 1491 - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`) 1492 1493 Calling sequence of `jac`: 1494 $ PetscErrorCode jac(TS ts, PetscReal t, Vec U, Vec U_t, Vec U_tt, PetscReal v, PetscReal a, Mat J, Mat P, void *ctx) 1495 + ts - the `TS` context obtained from `TSCreate()` 1496 . t - time at step/stage being solved 1497 . U - state vector 1498 . U_t - time derivative of state vector 1499 . U_tt - second time derivative of state vector 1500 . v - shift for U_t 1501 . a - shift for U_tt 1502 . 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 1503 . P - preconditioning matrix for J, may be same as J 1504 - ctx - [optional] user-defined context for matrix evaluation routine 1505 1506 Level: beginner 1507 1508 Notes: 1509 The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve. 1510 1511 The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be 1512 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. 1513 The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U where the positive "shift" 1514 parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states. 1515 1516 .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Jacobian()` 1517 @*/ 1518 PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2Jacobian jac, void *ctx) 1519 { 1520 DM dm; 1521 1522 PetscFunctionBegin; 1523 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1524 if (J) PetscValidHeaderSpecific(J, MAT_CLASSID, 2); 1525 if (P) PetscValidHeaderSpecific(P, MAT_CLASSID, 3); 1526 PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL)); 1527 PetscCall(TSGetDM(ts, &dm)); 1528 PetscCall(DMTSSetI2Jacobian(dm, jac, ctx)); 1529 PetscFunctionReturn(PETSC_SUCCESS); 1530 } 1531 1532 /*@C 1533 TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep. 1534 1535 Not Collective, but parallel objects are returned if `TS` is parallel 1536 1537 Input Parameter: 1538 . ts - The `TS` context obtained from `TSCreate()` 1539 1540 Output Parameters: 1541 + J - The (approximate) Jacobian of F(t,U,U_t,U_tt) 1542 . P - The matrix from which the preconditioner is constructed, often the same as `J` 1543 . jac - The function to compute the Jacobian matrices 1544 - ctx - User-defined context for Jacobian evaluation routine 1545 1546 Level: advanced 1547 1548 Note: 1549 You can pass in `NULL` for any return argument you do not need. 1550 1551 .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()` 1552 @*/ 1553 PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2Jacobian *jac, void **ctx) 1554 { 1555 SNES snes; 1556 DM dm; 1557 1558 PetscFunctionBegin; 1559 PetscCall(TSGetSNES(ts, &snes)); 1560 PetscCall(SNESSetUpMatrices(snes)); 1561 PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL)); 1562 PetscCall(TSGetDM(ts, &dm)); 1563 PetscCall(DMTSGetI2Jacobian(dm, jac, ctx)); 1564 PetscFunctionReturn(PETSC_SUCCESS); 1565 } 1566 1567 /*@ 1568 TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0 1569 1570 Collective 1571 1572 Input Parameters: 1573 + ts - the `TS` context 1574 . t - current time 1575 . U - state vector 1576 . V - time derivative of state vector (U_t) 1577 - A - second time derivative of state vector (U_tt) 1578 1579 Output Parameter: 1580 . F - the residual vector 1581 1582 Level: developer 1583 1584 Note: 1585 Most users should not need to explicitly call this routine, as it 1586 is used internally within the nonlinear solvers. 1587 1588 .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()` 1589 @*/ 1590 PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F) 1591 { 1592 DM dm; 1593 TSI2Function I2Function; 1594 void *ctx; 1595 TSRHSFunction rhsfunction; 1596 1597 PetscFunctionBegin; 1598 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1599 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 1600 PetscValidHeaderSpecific(V, VEC_CLASSID, 4); 1601 PetscValidHeaderSpecific(A, VEC_CLASSID, 5); 1602 PetscValidHeaderSpecific(F, VEC_CLASSID, 6); 1603 1604 PetscCall(TSGetDM(ts, &dm)); 1605 PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx)); 1606 PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL)); 1607 1608 if (!I2Function) { 1609 PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE)); 1610 PetscFunctionReturn(PETSC_SUCCESS); 1611 } 1612 1613 PetscCall(PetscLogEventBegin(TS_FunctionEval, ts, U, V, F)); 1614 1615 PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx)); 1616 1617 if (rhsfunction) { 1618 Vec Frhs; 1619 PetscCall(TSGetRHSVec_Private(ts, &Frhs)); 1620 PetscCall(TSComputeRHSFunction(ts, t, U, Frhs)); 1621 PetscCall(VecAXPY(F, -1, Frhs)); 1622 } 1623 1624 PetscCall(PetscLogEventEnd(TS_FunctionEval, ts, U, V, F)); 1625 PetscFunctionReturn(PETSC_SUCCESS); 1626 } 1627 1628 /*@ 1629 TSComputeI2Jacobian - Evaluates the Jacobian of the DAE 1630 1631 Collective 1632 1633 Input Parameters: 1634 + ts - the `TS` context 1635 . t - current timestep 1636 . U - state vector 1637 . V - time derivative of state vector 1638 . A - second time derivative of state vector 1639 . shiftV - shift to apply, see note below 1640 - shiftA - shift to apply, see note below 1641 1642 Output Parameters: 1643 + J - Jacobian matrix 1644 - P - optional preconditioning matrix 1645 1646 Level: developer 1647 1648 Notes: 1649 If F(t,U,V,A)=0 is the DAE, the required Jacobian is 1650 1651 dF/dU + shiftV*dF/dV + shiftA*dF/dA 1652 1653 Most users should not need to explicitly call this routine, as it 1654 is used internally within the nonlinear solvers. 1655 1656 .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()` 1657 @*/ 1658 PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P) 1659 { 1660 DM dm; 1661 TSI2Jacobian I2Jacobian; 1662 void *ctx; 1663 TSRHSJacobian rhsjacobian; 1664 1665 PetscFunctionBegin; 1666 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1667 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 1668 PetscValidHeaderSpecific(V, VEC_CLASSID, 4); 1669 PetscValidHeaderSpecific(A, VEC_CLASSID, 5); 1670 PetscValidHeaderSpecific(J, MAT_CLASSID, 8); 1671 PetscValidHeaderSpecific(P, MAT_CLASSID, 9); 1672 1673 PetscCall(TSGetDM(ts, &dm)); 1674 PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx)); 1675 PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL)); 1676 1677 if (!I2Jacobian) { 1678 PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE)); 1679 PetscFunctionReturn(PETSC_SUCCESS); 1680 } 1681 1682 PetscCall(PetscLogEventBegin(TS_JacobianEval, ts, U, J, P)); 1683 PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx)); 1684 if (rhsjacobian) { 1685 Mat Jrhs, Prhs; 1686 PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs)); 1687 PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs)); 1688 PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern)); 1689 if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern)); 1690 } 1691 1692 PetscCall(PetscLogEventEnd(TS_JacobianEval, ts, U, J, P)); 1693 PetscFunctionReturn(PETSC_SUCCESS); 1694 } 1695 1696 /*@C 1697 TSSetTransientVariable - sets function to transform from state to transient variables 1698 1699 Logically Collective 1700 1701 Input Parameters: 1702 + ts - time stepping context on which to change the transient variable 1703 . tvar - a function that transforms to transient variables 1704 - ctx - a context for tvar 1705 1706 Calling sequence of `tvar`: 1707 $ PetscErrorCode tvar(TS ts, Vec p, Vec c, void *ctx) 1708 + ts - timestep context 1709 . p - input vector (primitive form) 1710 . c - output vector, transient variables (conservative form) 1711 - ctx - [optional] user-defined function context 1712 1713 Level: advanced 1714 1715 Notes: 1716 This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`) 1717 can be conservative. In this context, primitive variables P are used to model the state (e.g., because they lead to 1718 well-conditioned formulations even in limiting cases such as low-Mach or zero porosity). The transient variable is 1719 C(P), specified by calling this function. An IFunction thus receives arguments (P, Cdot) and the IJacobian must be 1720 evaluated via the chain rule, as in 1721 .vb 1722 dF/dP + shift * dF/dCdot dC/dP. 1723 .ve 1724 1725 .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()` 1726 @*/ 1727 PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariable tvar, void *ctx) 1728 { 1729 DM dm; 1730 1731 PetscFunctionBegin; 1732 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1733 PetscCall(TSGetDM(ts, &dm)); 1734 PetscCall(DMTSSetTransientVariable(dm, tvar, ctx)); 1735 PetscFunctionReturn(PETSC_SUCCESS); 1736 } 1737 1738 /*@ 1739 TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables 1740 1741 Logically Collective 1742 1743 Input Parameters: 1744 + ts - TS on which to compute 1745 - U - state vector to be transformed to transient variables 1746 1747 Output Parameter: 1748 . C - transient (conservative) variable 1749 1750 Level: developer 1751 1752 Developer Note: 1753 If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed. 1754 This makes it safe to call without a guard. One can use `TSHasTransientVariable()` to check if transient variables are 1755 being used. 1756 1757 .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()` 1758 @*/ 1759 PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C) 1760 { 1761 DM dm; 1762 DMTS dmts; 1763 1764 PetscFunctionBegin; 1765 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1766 PetscValidHeaderSpecific(U, VEC_CLASSID, 2); 1767 PetscCall(TSGetDM(ts, &dm)); 1768 PetscCall(DMGetDMTS(dm, &dmts)); 1769 if (dmts->ops->transientvar) { 1770 PetscValidHeaderSpecific(C, VEC_CLASSID, 3); 1771 PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx)); 1772 } 1773 PetscFunctionReturn(PETSC_SUCCESS); 1774 } 1775 1776 /*@ 1777 TSHasTransientVariable - determine whether transient variables have been set 1778 1779 Logically Collective 1780 1781 Input Parameter: 1782 . ts - `TS` on which to compute 1783 1784 Output Parameter: 1785 . has - `PETSC_TRUE` if transient variables have been set 1786 1787 Level: developer 1788 1789 .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()` 1790 @*/ 1791 PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has) 1792 { 1793 DM dm; 1794 DMTS dmts; 1795 1796 PetscFunctionBegin; 1797 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1798 PetscCall(TSGetDM(ts, &dm)); 1799 PetscCall(DMGetDMTS(dm, &dmts)); 1800 *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE; 1801 PetscFunctionReturn(PETSC_SUCCESS); 1802 } 1803 1804 /*@ 1805 TS2SetSolution - Sets the initial solution and time derivative vectors 1806 for use by the `TS` routines handling second order equations. 1807 1808 Logically Collective 1809 1810 Input Parameters: 1811 + ts - the `TS` context obtained from `TSCreate()` 1812 . u - the solution vector 1813 - v - the time derivative vector 1814 1815 Level: beginner 1816 1817 .seealso: [](ch_ts), `TS` 1818 @*/ 1819 PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v) 1820 { 1821 PetscFunctionBegin; 1822 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1823 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 1824 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 1825 PetscCall(TSSetSolution(ts, u)); 1826 PetscCall(PetscObjectReference((PetscObject)v)); 1827 PetscCall(VecDestroy(&ts->vec_dot)); 1828 ts->vec_dot = v; 1829 PetscFunctionReturn(PETSC_SUCCESS); 1830 } 1831 1832 /*@ 1833 TS2GetSolution - Returns the solution and time derivative at the present timestep 1834 for second order equations. It is valid to call this routine inside the function 1835 that you are evaluating in order to move to the new timestep. This vector not 1836 changed until the solution at the next timestep has been calculated. 1837 1838 Not Collective, but `u` returned is parallel if `TS` is parallel 1839 1840 Input Parameter: 1841 . ts - the `TS` context obtained from `TSCreate()` 1842 1843 Output Parameters: 1844 + u - the vector containing the solution 1845 - v - the vector containing the time derivative 1846 1847 Level: intermediate 1848 1849 .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()` 1850 @*/ 1851 PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v) 1852 { 1853 PetscFunctionBegin; 1854 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1855 if (u) PetscValidPointer(u, 2); 1856 if (v) PetscValidPointer(v, 3); 1857 if (u) *u = ts->vec_sol; 1858 if (v) *v = ts->vec_dot; 1859 PetscFunctionReturn(PETSC_SUCCESS); 1860 } 1861 1862 /*@C 1863 TSLoad - Loads a `TS` that has been stored in binary with `TSView()`. 1864 1865 Collective 1866 1867 Input Parameters: 1868 + newdm - the newly loaded `TS`, this needs to have been created with `TSCreate()` or 1869 some related function before a call to `TSLoad()`. 1870 - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` 1871 1872 Level: intermediate 1873 1874 Note: 1875 The type is determined by the data in the file, any type set into the `TS` before this call is ignored. 1876 1877 .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()` 1878 @*/ 1879 PetscErrorCode TSLoad(TS ts, PetscViewer viewer) 1880 { 1881 PetscBool isbinary; 1882 PetscInt classid; 1883 char type[256]; 1884 DMTS sdm; 1885 DM dm; 1886 1887 PetscFunctionBegin; 1888 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1889 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 1890 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 1891 PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 1892 1893 PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT)); 1894 PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file"); 1895 PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR)); 1896 PetscCall(TSSetType(ts, type)); 1897 PetscTryTypeMethod(ts, load, viewer); 1898 PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm)); 1899 PetscCall(DMLoad(dm, viewer)); 1900 PetscCall(TSSetDM(ts, dm)); 1901 PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol)); 1902 PetscCall(VecLoad(ts->vec_sol, viewer)); 1903 PetscCall(DMGetDMTS(ts->dm, &sdm)); 1904 PetscCall(DMTSLoad(sdm, viewer)); 1905 PetscFunctionReturn(PETSC_SUCCESS); 1906 } 1907 1908 #include <petscdraw.h> 1909 #if defined(PETSC_HAVE_SAWS) 1910 #include <petscviewersaws.h> 1911 #endif 1912 1913 /*@C 1914 TSViewFromOptions - View a `TS` based on values in the options database 1915 1916 Collective 1917 1918 Input Parameters: 1919 + ts - the `TS` context 1920 . obj - Optional object that provides the prefix for the options database keys 1921 - name - command line option string to be passed by user 1922 1923 Level: intermediate 1924 1925 .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()` 1926 @*/ 1927 PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[]) 1928 { 1929 PetscFunctionBegin; 1930 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1931 PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name)); 1932 PetscFunctionReturn(PETSC_SUCCESS); 1933 } 1934 1935 /*@C 1936 TSView - Prints the `TS` data structure. 1937 1938 Collective 1939 1940 Input Parameters: 1941 + ts - the `TS` context obtained from `TSCreate()` 1942 - viewer - visualization context 1943 1944 Options Database Key: 1945 . -ts_view - calls `TSView()` at end of `TSStep()` 1946 1947 Level: beginner 1948 1949 Notes: 1950 The available visualization contexts include 1951 + `PETSC_VIEWER_STDOUT_SELF` - standard output (default) 1952 - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard 1953 output where only the first processor opens 1954 the file. All other processors send their 1955 data to the first processor to print. 1956 1957 The user can open an alternative visualization context with 1958 `PetscViewerASCIIOpen()` - output to a specified file. 1959 1960 In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer). 1961 1962 .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()` 1963 @*/ 1964 PetscErrorCode TSView(TS ts, PetscViewer viewer) 1965 { 1966 TSType type; 1967 PetscBool iascii, isstring, isundials, isbinary, isdraw; 1968 DMTS sdm; 1969 #if defined(PETSC_HAVE_SAWS) 1970 PetscBool issaws; 1971 #endif 1972 1973 PetscFunctionBegin; 1974 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1975 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer)); 1976 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 1977 PetscCheckSameComm(ts, 1, viewer, 2); 1978 1979 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1980 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring)); 1981 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 1982 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 1983 #if defined(PETSC_HAVE_SAWS) 1984 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws)); 1985 #endif 1986 if (iascii) { 1987 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer)); 1988 if (ts->ops->view) { 1989 PetscCall(PetscViewerASCIIPushTab(viewer)); 1990 PetscUseTypeMethod(ts, view, viewer); 1991 PetscCall(PetscViewerASCIIPopTab(viewer)); 1992 } 1993 if (ts->max_steps < PETSC_MAX_INT) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum steps=%" PetscInt_FMT "\n", ts->max_steps)); 1994 if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, " maximum time=%g\n", (double)ts->max_time)); 1995 if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs)); 1996 if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs)); 1997 if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs)); 1998 if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs)); 1999 if (ts->usessnes) { 2000 PetscBool lin; 2001 if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its)); 2002 PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its)); 2003 PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, "")); 2004 PetscCall(PetscViewerASCIIPrintf(viewer, " total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures)); 2005 } 2006 PetscCall(PetscViewerASCIIPrintf(viewer, " total number of rejected steps=%" PetscInt_FMT "\n", ts->reject)); 2007 if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of relative error tolerances, ")); 2008 else PetscCall(PetscViewerASCIIPrintf(viewer, " using relative error tolerance of %g, ", (double)ts->rtol)); 2009 if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, " using vector of absolute error tolerances\n")); 2010 else PetscCall(PetscViewerASCIIPrintf(viewer, " using absolute error tolerance of %g\n", (double)ts->atol)); 2011 PetscCall(PetscViewerASCIIPushTab(viewer)); 2012 PetscCall(TSAdaptView(ts->adapt, viewer)); 2013 PetscCall(PetscViewerASCIIPopTab(viewer)); 2014 } else if (isstring) { 2015 PetscCall(TSGetType(ts, &type)); 2016 PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type)); 2017 PetscTryTypeMethod(ts, view, viewer); 2018 } else if (isbinary) { 2019 PetscInt classid = TS_FILE_CLASSID; 2020 MPI_Comm comm; 2021 PetscMPIInt rank; 2022 char type[256]; 2023 2024 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm)); 2025 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2026 if (rank == 0) { 2027 PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT)); 2028 PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256)); 2029 PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR)); 2030 } 2031 PetscTryTypeMethod(ts, view, viewer); 2032 if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer)); 2033 PetscCall(DMView(ts->dm, viewer)); 2034 PetscCall(VecView(ts->vec_sol, viewer)); 2035 PetscCall(DMGetDMTS(ts->dm, &sdm)); 2036 PetscCall(DMTSView(sdm, viewer)); 2037 } else if (isdraw) { 2038 PetscDraw draw; 2039 char str[36]; 2040 PetscReal x, y, bottom, h; 2041 2042 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 2043 PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); 2044 PetscCall(PetscStrncpy(str, "TS: ", sizeof(str))); 2045 PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str))); 2046 PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h)); 2047 bottom = y - h; 2048 PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom)); 2049 PetscTryTypeMethod(ts, view, viewer); 2050 if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer)); 2051 if (ts->snes) PetscCall(SNESView(ts->snes, viewer)); 2052 PetscCall(PetscDrawPopCurrentPoint(draw)); 2053 #if defined(PETSC_HAVE_SAWS) 2054 } else if (issaws) { 2055 PetscMPIInt rank; 2056 const char *name; 2057 2058 PetscCall(PetscObjectGetName((PetscObject)ts, &name)); 2059 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 2060 if (!((PetscObject)ts)->amsmem && rank == 0) { 2061 char dir[1024]; 2062 2063 PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer)); 2064 PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name)); 2065 PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT)); 2066 PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name)); 2067 PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE)); 2068 } 2069 PetscTryTypeMethod(ts, view, viewer); 2070 #endif 2071 } 2072 if (ts->snes && ts->usessnes) { 2073 PetscCall(PetscViewerASCIIPushTab(viewer)); 2074 PetscCall(SNESView(ts->snes, viewer)); 2075 PetscCall(PetscViewerASCIIPopTab(viewer)); 2076 } 2077 PetscCall(DMGetDMTS(ts->dm, &sdm)); 2078 PetscCall(DMTSView(sdm, viewer)); 2079 2080 PetscCall(PetscViewerASCIIPushTab(viewer)); 2081 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials)); 2082 PetscCall(PetscViewerASCIIPopTab(viewer)); 2083 PetscFunctionReturn(PETSC_SUCCESS); 2084 } 2085 2086 /*@ 2087 TSSetApplicationContext - Sets an optional user-defined context for 2088 the timesteppers. 2089 2090 Logically Collective 2091 2092 Input Parameters: 2093 + ts - the `TS` context obtained from `TSCreate()` 2094 - usrP - user context 2095 2096 Level: intermediate 2097 2098 Fortran Note: 2099 You must write a Fortran interface definition for this 2100 function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument. 2101 2102 .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()` 2103 @*/ 2104 PetscErrorCode TSSetApplicationContext(TS ts, void *usrP) 2105 { 2106 PetscFunctionBegin; 2107 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2108 ts->user = usrP; 2109 PetscFunctionReturn(PETSC_SUCCESS); 2110 } 2111 2112 /*@ 2113 TSGetApplicationContext - Gets the user-defined context for the 2114 timestepper that was set with `TSSetApplicationContext()` 2115 2116 Not Collective 2117 2118 Input Parameter: 2119 . ts - the `TS` context obtained from `TSCreate()` 2120 2121 Output Parameter: 2122 . usrP - user context 2123 2124 Level: intermediate 2125 2126 Fortran Note: 2127 You must write a Fortran interface definition for this 2128 function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument. 2129 2130 .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()` 2131 @*/ 2132 PetscErrorCode TSGetApplicationContext(TS ts, void *usrP) 2133 { 2134 PetscFunctionBegin; 2135 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2136 *(void **)usrP = ts->user; 2137 PetscFunctionReturn(PETSC_SUCCESS); 2138 } 2139 2140 /*@ 2141 TSGetStepNumber - Gets the number of time steps completed. 2142 2143 Not Collective 2144 2145 Input Parameter: 2146 . ts - the `TS` context obtained from `TSCreate()` 2147 2148 Output Parameter: 2149 . steps - number of steps completed so far 2150 2151 Level: intermediate 2152 2153 .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()` 2154 @*/ 2155 PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps) 2156 { 2157 PetscFunctionBegin; 2158 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2159 PetscValidIntPointer(steps, 2); 2160 *steps = ts->steps; 2161 PetscFunctionReturn(PETSC_SUCCESS); 2162 } 2163 2164 /*@ 2165 TSSetStepNumber - Sets the number of steps completed. 2166 2167 Logically Collective 2168 2169 Input Parameters: 2170 + ts - the `TS` context 2171 - steps - number of steps completed so far 2172 2173 Level: developer 2174 2175 Note: 2176 For most uses of the `TS` solvers the user need not explicitly call 2177 `TSSetStepNumber()`, as the step counter is appropriately updated in 2178 `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to 2179 reinitialize timestepping by setting the step counter to zero (and time 2180 to the initial time) to solve a similar problem with different initial 2181 conditions or parameters. Other possible use case is to continue 2182 timestepping from a previously interrupted run in such a way that `TS` 2183 monitors will be called with a initial nonzero step counter. 2184 2185 .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()` 2186 @*/ 2187 PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps) 2188 { 2189 PetscFunctionBegin; 2190 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2191 PetscValidLogicalCollectiveInt(ts, steps, 2); 2192 PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative"); 2193 ts->steps = steps; 2194 PetscFunctionReturn(PETSC_SUCCESS); 2195 } 2196 2197 /*@ 2198 TSSetTimeStep - Allows one to reset the timestep at any time, 2199 useful for simple pseudo-timestepping codes. 2200 2201 Logically Collective 2202 2203 Input Parameters: 2204 + ts - the `TS` context obtained from `TSCreate()` 2205 - time_step - the size of the timestep 2206 2207 Level: intermediate 2208 2209 .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()` 2210 @*/ 2211 PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step) 2212 { 2213 PetscFunctionBegin; 2214 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2215 PetscValidLogicalCollectiveReal(ts, time_step, 2); 2216 ts->time_step = time_step; 2217 PetscFunctionReturn(PETSC_SUCCESS); 2218 } 2219 2220 /*@ 2221 TSSetExactFinalTime - Determines whether to adapt the final time step to 2222 match the exact final time, interpolate solution to the exact final time, 2223 or just return at the final time `TS` computed. 2224 2225 Logically Collective 2226 2227 Input Parameters: 2228 + ts - the time-step context 2229 - eftopt - exact final time option 2230 .vb 2231 TS_EXACTFINALTIME_STEPOVER - Don't do anything if final time is exceeded 2232 TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time 2233 TS_EXACTFINALTIME_MATCHSTEP - Adapt final time step to match the final time 2234 .ve 2235 2236 Options Database Key: 2237 . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step at runtime 2238 2239 Level: beginner 2240 2241 Note: 2242 If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time 2243 then the final time you selected. 2244 2245 .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()` 2246 @*/ 2247 PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt) 2248 { 2249 PetscFunctionBegin; 2250 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2251 PetscValidLogicalCollectiveEnum(ts, eftopt, 2); 2252 ts->exact_final_time = eftopt; 2253 PetscFunctionReturn(PETSC_SUCCESS); 2254 } 2255 2256 /*@ 2257 TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()` 2258 2259 Not Collective 2260 2261 Input Parameter: 2262 . ts - the `TS` context 2263 2264 Output Parameter: 2265 . eftopt - exact final time option 2266 2267 Level: beginner 2268 2269 .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()` 2270 @*/ 2271 PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt) 2272 { 2273 PetscFunctionBegin; 2274 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2275 PetscValidPointer(eftopt, 2); 2276 *eftopt = ts->exact_final_time; 2277 PetscFunctionReturn(PETSC_SUCCESS); 2278 } 2279 2280 /*@ 2281 TSGetTimeStep - Gets the current timestep size. 2282 2283 Not Collective 2284 2285 Input Parameter: 2286 . ts - the `TS` context obtained from `TSCreate()` 2287 2288 Output Parameter: 2289 . dt - the current timestep size 2290 2291 Level: intermediate 2292 2293 .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()` 2294 @*/ 2295 PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt) 2296 { 2297 PetscFunctionBegin; 2298 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2299 PetscValidRealPointer(dt, 2); 2300 *dt = ts->time_step; 2301 PetscFunctionReturn(PETSC_SUCCESS); 2302 } 2303 2304 /*@ 2305 TSGetSolution - Returns the solution at the present timestep. It 2306 is valid to call this routine inside the function that you are evaluating 2307 in order to move to the new timestep. This vector not changed until 2308 the solution at the next timestep has been calculated. 2309 2310 Not Collective, but v returned is parallel if ts is parallel 2311 2312 Input Parameter: 2313 . ts - the `TS` context obtained from `TSCreate()` 2314 2315 Output Parameter: 2316 . v - the vector containing the solution 2317 2318 Level: intermediate 2319 2320 Note: 2321 If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested 2322 final time. It returns the solution at the next timestep. 2323 2324 .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()` 2325 @*/ 2326 PetscErrorCode TSGetSolution(TS ts, Vec *v) 2327 { 2328 PetscFunctionBegin; 2329 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2330 PetscValidPointer(v, 2); 2331 *v = ts->vec_sol; 2332 PetscFunctionReturn(PETSC_SUCCESS); 2333 } 2334 2335 /*@ 2336 TSGetSolutionComponents - Returns any solution components at the present 2337 timestep, if available for the time integration method being used. 2338 Solution components are quantities that share the same size and 2339 structure as the solution vector. 2340 2341 Not Collective, but v returned is parallel if ts is parallel 2342 2343 Parameters : 2344 + ts - the `TS` context obtained from `TSCreate()` (input parameter). 2345 . n - If v is `NULL`, then the number of solution components is 2346 returned through n, else the n-th solution component is 2347 returned in v. 2348 - v - the vector containing the n-th solution component 2349 (may be `NULL` to use this function to find out 2350 the number of solutions components). 2351 2352 Level: advanced 2353 2354 .seealso: [](ch_ts), `TS`, `TSGetSolution()` 2355 @*/ 2356 PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v) 2357 { 2358 PetscFunctionBegin; 2359 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2360 if (!ts->ops->getsolutioncomponents) *n = 0; 2361 else PetscUseTypeMethod(ts, getsolutioncomponents, n, v); 2362 PetscFunctionReturn(PETSC_SUCCESS); 2363 } 2364 2365 /*@ 2366 TSGetAuxSolution - Returns an auxiliary solution at the present 2367 timestep, if available for the time integration method being used. 2368 2369 Not Collective, but v returned is parallel if ts is parallel 2370 2371 Parameters : 2372 + ts - the `TS` context obtained from `TSCreate()` (input parameter). 2373 - v - the vector containing the auxiliary solution 2374 2375 Level: intermediate 2376 2377 .seealso: [](ch_ts), `TS`, `TSGetSolution()` 2378 @*/ 2379 PetscErrorCode TSGetAuxSolution(TS ts, Vec *v) 2380 { 2381 PetscFunctionBegin; 2382 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2383 if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v); 2384 else PetscCall(VecZeroEntries(*v)); 2385 PetscFunctionReturn(PETSC_SUCCESS); 2386 } 2387 2388 /*@ 2389 TSGetTimeError - Returns the estimated error vector, if the chosen 2390 `TSType` has an error estimation functionality and `TSSetTimeError()` was called 2391 2392 Not Collective, but v returned is parallel if ts is parallel 2393 2394 Parameters : 2395 + ts - the `TS` context obtained from `TSCreate()` (input parameter). 2396 . n - current estimate (n=0) or previous one (n=-1) 2397 - v - the vector containing the error (same size as the solution). 2398 2399 Level: intermediate 2400 2401 Note: 2402 MUST call after `TSSetUp()` 2403 2404 .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()` 2405 @*/ 2406 PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v) 2407 { 2408 PetscFunctionBegin; 2409 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2410 if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v); 2411 else PetscCall(VecZeroEntries(*v)); 2412 PetscFunctionReturn(PETSC_SUCCESS); 2413 } 2414 2415 /*@ 2416 TSSetTimeError - Sets the estimated error vector, if the chosen 2417 `TSType` has an error estimation functionality. This can be used 2418 to restart such a time integrator with a given error vector. 2419 2420 Not Collective, but v returned is parallel if ts is parallel 2421 2422 Parameters : 2423 + ts - the `TS` context obtained from `TSCreate()` (input parameter). 2424 - v - the vector containing the error (same size as the solution). 2425 2426 Level: intermediate 2427 2428 .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError)` 2429 @*/ 2430 PetscErrorCode TSSetTimeError(TS ts, Vec v) 2431 { 2432 PetscFunctionBegin; 2433 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2434 PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first"); 2435 PetscTryTypeMethod(ts, settimeerror, v); 2436 PetscFunctionReturn(PETSC_SUCCESS); 2437 } 2438 2439 /* ----- Routines to initialize and destroy a timestepper ---- */ 2440 /*@ 2441 TSSetProblemType - Sets the type of problem to be solved. 2442 2443 Not collective 2444 2445 Input Parameters: 2446 + ts - The `TS` 2447 - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms 2448 .vb 2449 U_t - A U = 0 (linear) 2450 U_t - A(t) U = 0 (linear) 2451 F(t,U,U_t) = 0 (nonlinear) 2452 .ve 2453 2454 Level: beginner 2455 2456 .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS` 2457 @*/ 2458 PetscErrorCode TSSetProblemType(TS ts, TSProblemType type) 2459 { 2460 PetscFunctionBegin; 2461 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2462 ts->problem_type = type; 2463 if (type == TS_LINEAR) { 2464 SNES snes; 2465 PetscCall(TSGetSNES(ts, &snes)); 2466 PetscCall(SNESSetType(snes, SNESKSPONLY)); 2467 } 2468 PetscFunctionReturn(PETSC_SUCCESS); 2469 } 2470 2471 /*@C 2472 TSGetProblemType - Gets the type of problem to be solved. 2473 2474 Not collective 2475 2476 Input Parameter: 2477 . ts - The `TS` 2478 2479 Output Parameter: 2480 . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms 2481 .vb 2482 M U_t = A U 2483 M(t) U_t = A(t) U 2484 F(t,U,U_t) 2485 .ve 2486 2487 Level: beginner 2488 2489 .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS` 2490 @*/ 2491 PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type) 2492 { 2493 PetscFunctionBegin; 2494 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2495 PetscValidIntPointer(type, 2); 2496 *type = ts->problem_type; 2497 PetscFunctionReturn(PETSC_SUCCESS); 2498 } 2499 2500 /* 2501 Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp() 2502 */ 2503 static PetscErrorCode TSSetExactFinalTimeDefault(TS ts) 2504 { 2505 PetscBool isnone; 2506 2507 PetscFunctionBegin; 2508 PetscCall(TSGetAdapt(ts, &ts->adapt)); 2509 PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type)); 2510 2511 PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone)); 2512 if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP; 2513 else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE; 2514 PetscFunctionReturn(PETSC_SUCCESS); 2515 } 2516 2517 /*@ 2518 TSSetUp - Sets up the internal data structures for the later use of a timestepper. 2519 2520 Collective 2521 2522 Input Parameter: 2523 . ts - the `TS` context obtained from `TSCreate()` 2524 2525 Level: advanced 2526 2527 Note: 2528 For basic use of the `TS` solvers the user need not explicitly call 2529 `TSSetUp()`, since these actions will automatically occur during 2530 the call to `TSStep()` or `TSSolve()`. However, if one wishes to control this 2531 phase separately, `TSSetUp()` should be called after `TSCreate()` 2532 and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`. 2533 2534 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()` 2535 @*/ 2536 PetscErrorCode TSSetUp(TS ts) 2537 { 2538 DM dm; 2539 PetscErrorCode (*func)(SNES, Vec, Vec, void *); 2540 PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *); 2541 TSIFunction ifun; 2542 TSIJacobian ijac; 2543 TSI2Jacobian i2jac; 2544 TSRHSJacobian rhsjac; 2545 2546 PetscFunctionBegin; 2547 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2548 if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS); 2549 2550 if (!((PetscObject)ts)->type_name) { 2551 PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL)); 2552 PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER)); 2553 } 2554 2555 if (!ts->vec_sol) { 2556 PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first"); 2557 PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol)); 2558 } 2559 2560 if (ts->tspan) { 2561 if (!ts->tspan->vecs_sol) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol)); 2562 } 2563 if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */ 2564 PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs)); 2565 ts->Jacp = ts->Jacprhs; 2566 } 2567 2568 if (ts->quadraturets) { 2569 PetscCall(TSSetUp(ts->quadraturets)); 2570 PetscCall(VecDestroy(&ts->vec_costintegrand)); 2571 PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand)); 2572 } 2573 2574 PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL)); 2575 if (rhsjac == TSComputeRHSJacobianConstant) { 2576 Mat Amat, Pmat; 2577 SNES snes; 2578 PetscCall(TSGetSNES(ts, &snes)); 2579 PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL)); 2580 /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would 2581 * have displaced the RHS matrix */ 2582 if (Amat && Amat == ts->Arhs) { 2583 /* 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 */ 2584 PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat)); 2585 PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL)); 2586 PetscCall(MatDestroy(&Amat)); 2587 } 2588 if (Pmat && Pmat == ts->Brhs) { 2589 PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat)); 2590 PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL)); 2591 PetscCall(MatDestroy(&Pmat)); 2592 } 2593 } 2594 2595 PetscCall(TSGetAdapt(ts, &ts->adapt)); 2596 PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type)); 2597 2598 PetscTryTypeMethod(ts, setup); 2599 2600 PetscCall(TSSetExactFinalTimeDefault(ts)); 2601 2602 /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction 2603 to be set right but can't do it elsewhere due to the overreliance on ctx=ts. 2604 */ 2605 PetscCall(TSGetDM(ts, &dm)); 2606 PetscCall(DMSNESGetFunction(dm, &func, NULL)); 2607 if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts)); 2608 2609 /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it. 2610 Otherwise, the SNES will use coloring internally to form the Jacobian. 2611 */ 2612 PetscCall(DMSNESGetJacobian(dm, &jac, NULL)); 2613 PetscCall(DMTSGetIJacobian(dm, &ijac, NULL)); 2614 PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL)); 2615 PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL)); 2616 if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts)); 2617 2618 /* if time integration scheme has a starting method, call it */ 2619 PetscTryTypeMethod(ts, startingmethod); 2620 2621 ts->setupcalled = PETSC_TRUE; 2622 PetscFunctionReturn(PETSC_SUCCESS); 2623 } 2624 2625 /*@ 2626 TSReset - Resets a `TS` context and removes any allocated `Vec`s and `Mat`s. 2627 2628 Collective 2629 2630 Input Parameter: 2631 . ts - the `TS` context obtained from `TSCreate()` 2632 2633 Level: beginner 2634 2635 .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()` 2636 @*/ 2637 PetscErrorCode TSReset(TS ts) 2638 { 2639 TS_RHSSplitLink ilink = ts->tsrhssplit, next; 2640 2641 PetscFunctionBegin; 2642 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2643 2644 PetscTryTypeMethod(ts, reset); 2645 if (ts->snes) PetscCall(SNESReset(ts->snes)); 2646 if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt)); 2647 2648 PetscCall(MatDestroy(&ts->Arhs)); 2649 PetscCall(MatDestroy(&ts->Brhs)); 2650 PetscCall(VecDestroy(&ts->Frhs)); 2651 PetscCall(VecDestroy(&ts->vec_sol)); 2652 PetscCall(VecDestroy(&ts->vec_dot)); 2653 PetscCall(VecDestroy(&ts->vatol)); 2654 PetscCall(VecDestroy(&ts->vrtol)); 2655 PetscCall(VecDestroyVecs(ts->nwork, &ts->work)); 2656 2657 PetscCall(MatDestroy(&ts->Jacprhs)); 2658 PetscCall(MatDestroy(&ts->Jacp)); 2659 if (ts->forward_solve) PetscCall(TSForwardReset(ts)); 2660 if (ts->quadraturets) { 2661 PetscCall(TSReset(ts->quadraturets)); 2662 PetscCall(VecDestroy(&ts->vec_costintegrand)); 2663 } 2664 while (ilink) { 2665 next = ilink->next; 2666 PetscCall(TSDestroy(&ilink->ts)); 2667 PetscCall(PetscFree(ilink->splitname)); 2668 PetscCall(ISDestroy(&ilink->is)); 2669 PetscCall(PetscFree(ilink)); 2670 ilink = next; 2671 } 2672 ts->tsrhssplit = NULL; 2673 ts->num_rhs_splits = 0; 2674 if (ts->tspan) { 2675 PetscCall(PetscFree(ts->tspan->span_times)); 2676 PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol)); 2677 PetscCall(PetscFree(ts->tspan)); 2678 } 2679 ts->setupcalled = PETSC_FALSE; 2680 PetscFunctionReturn(PETSC_SUCCESS); 2681 } 2682 2683 /*@C 2684 TSDestroy - Destroys the timestepper context that was created 2685 with `TSCreate()`. 2686 2687 Collective 2688 2689 Input Parameter: 2690 . ts - the `TS` context obtained from `TSCreate()` 2691 2692 Level: beginner 2693 2694 .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()` 2695 @*/ 2696 PetscErrorCode TSDestroy(TS *ts) 2697 { 2698 PetscFunctionBegin; 2699 if (!*ts) PetscFunctionReturn(PETSC_SUCCESS); 2700 PetscValidHeaderSpecific(*ts, TS_CLASSID, 1); 2701 if (--((PetscObject)(*ts))->refct > 0) { 2702 *ts = NULL; 2703 PetscFunctionReturn(PETSC_SUCCESS); 2704 } 2705 2706 PetscCall(TSReset(*ts)); 2707 PetscCall(TSAdjointReset(*ts)); 2708 if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts)); 2709 2710 /* if memory was published with SAWs then destroy it */ 2711 PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts)); 2712 PetscTryTypeMethod((*ts), destroy); 2713 2714 PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory)); 2715 2716 PetscCall(TSAdaptDestroy(&(*ts)->adapt)); 2717 PetscCall(TSEventDestroy(&(*ts)->event)); 2718 2719 PetscCall(SNESDestroy(&(*ts)->snes)); 2720 PetscCall(DMDestroy(&(*ts)->dm)); 2721 PetscCall(TSMonitorCancel((*ts))); 2722 PetscCall(TSAdjointMonitorCancel((*ts))); 2723 2724 PetscCall(TSDestroy(&(*ts)->quadraturets)); 2725 PetscCall(PetscHeaderDestroy(ts)); 2726 PetscFunctionReturn(PETSC_SUCCESS); 2727 } 2728 2729 /*@ 2730 TSGetSNES - Returns the `SNES` (nonlinear solver) associated with 2731 a `TS` (timestepper) context. Valid only for nonlinear problems. 2732 2733 Not Collective, but snes is parallel if ts is parallel 2734 2735 Input Parameter: 2736 . ts - the `TS` context obtained from `TSCreate()` 2737 2738 Output Parameter: 2739 . snes - the nonlinear solver context 2740 2741 Level: beginner 2742 2743 Notes: 2744 The user can then directly manipulate the `SNES` context to set various 2745 options, etc. Likewise, the user can then extract and manipulate the 2746 `KSP`, and `PC` contexts as well. 2747 2748 `TSGetSNES()` does not work for integrators that do not use `SNES`; in 2749 this case `TSGetSNES()` returns `NULL` in `snes`. 2750 2751 .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()` 2752 @*/ 2753 PetscErrorCode TSGetSNES(TS ts, SNES *snes) 2754 { 2755 PetscFunctionBegin; 2756 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2757 PetscValidPointer(snes, 2); 2758 if (!ts->snes) { 2759 PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes)); 2760 PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options)); 2761 PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts)); 2762 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1)); 2763 if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm)); 2764 if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY)); 2765 } 2766 *snes = ts->snes; 2767 PetscFunctionReturn(PETSC_SUCCESS); 2768 } 2769 2770 /*@ 2771 TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the timestepping context 2772 2773 Collective 2774 2775 Input Parameters: 2776 + ts - the `TS` context obtained from `TSCreate()` 2777 - snes - the nonlinear solver context 2778 2779 Level: developer 2780 2781 Note: 2782 Most users should have the `TS` created by calling `TSGetSNES()` 2783 2784 .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()` 2785 @*/ 2786 PetscErrorCode TSSetSNES(TS ts, SNES snes) 2787 { 2788 PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *); 2789 2790 PetscFunctionBegin; 2791 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2792 PetscValidHeaderSpecific(snes, SNES_CLASSID, 2); 2793 PetscCall(PetscObjectReference((PetscObject)snes)); 2794 PetscCall(SNESDestroy(&ts->snes)); 2795 2796 ts->snes = snes; 2797 2798 PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts)); 2799 PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL)); 2800 if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts)); 2801 PetscFunctionReturn(PETSC_SUCCESS); 2802 } 2803 2804 /*@ 2805 TSGetKSP - Returns the `KSP` (linear solver) associated with 2806 a `TS` (timestepper) context. 2807 2808 Not Collective, but `ksp` is parallel if `ts` is parallel 2809 2810 Input Parameter: 2811 . ts - the `TS` context obtained from `TSCreate()` 2812 2813 Output Parameter: 2814 . ksp - the nonlinear solver context 2815 2816 Level: beginner 2817 2818 Notes: 2819 The user can then directly manipulate the `KSP` context to set various 2820 options, etc. Likewise, the user can then extract and manipulate the 2821 `PC` context as well. 2822 2823 `TSGetKSP()` does not work for integrators that do not use `KSP`; 2824 in this case `TSGetKSP()` returns `NULL` in `ksp`. 2825 2826 .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()` 2827 @*/ 2828 PetscErrorCode TSGetKSP(TS ts, KSP *ksp) 2829 { 2830 SNES snes; 2831 2832 PetscFunctionBegin; 2833 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2834 PetscValidPointer(ksp, 2); 2835 PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first"); 2836 PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()"); 2837 PetscCall(TSGetSNES(ts, &snes)); 2838 PetscCall(SNESGetKSP(snes, ksp)); 2839 PetscFunctionReturn(PETSC_SUCCESS); 2840 } 2841 2842 /* ----------- Routines to set solver parameters ---------- */ 2843 2844 /*@ 2845 TSSetMaxSteps - Sets the maximum number of steps to use. 2846 2847 Logically Collective 2848 2849 Input Parameters: 2850 + ts - the `TS` context obtained from `TSCreate()` 2851 - maxsteps - maximum number of steps to use 2852 2853 Options Database Key: 2854 . -ts_max_steps <maxsteps> - Sets maxsteps 2855 2856 Level: intermediate 2857 2858 Note: 2859 The default maximum number of steps is 5000 2860 2861 .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()` 2862 @*/ 2863 PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps) 2864 { 2865 PetscFunctionBegin; 2866 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2867 PetscValidLogicalCollectiveInt(ts, maxsteps, 2); 2868 PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative"); 2869 ts->max_steps = maxsteps; 2870 PetscFunctionReturn(PETSC_SUCCESS); 2871 } 2872 2873 /*@ 2874 TSGetMaxSteps - Gets the maximum number of steps to use. 2875 2876 Not Collective 2877 2878 Input Parameter: 2879 . ts - the `TS` context obtained from `TSCreate()` 2880 2881 Output Parameter: 2882 . maxsteps - maximum number of steps to use 2883 2884 Level: advanced 2885 2886 .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()` 2887 @*/ 2888 PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps) 2889 { 2890 PetscFunctionBegin; 2891 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2892 PetscValidIntPointer(maxsteps, 2); 2893 *maxsteps = ts->max_steps; 2894 PetscFunctionReturn(PETSC_SUCCESS); 2895 } 2896 2897 /*@ 2898 TSSetMaxTime - Sets the maximum (or final) time for timestepping. 2899 2900 Logically Collective 2901 2902 Input Parameters: 2903 + ts - the `TS` context obtained from `TSCreate()` 2904 - maxtime - final time to step to 2905 2906 Options Database Key: 2907 . -ts_max_time <maxtime> - Sets maxtime 2908 2909 Level: intermediate 2910 2911 Notes: 2912 The default maximum time is 5.0 2913 2914 .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()` 2915 @*/ 2916 PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime) 2917 { 2918 PetscFunctionBegin; 2919 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2920 PetscValidLogicalCollectiveReal(ts, maxtime, 2); 2921 ts->max_time = maxtime; 2922 PetscFunctionReturn(PETSC_SUCCESS); 2923 } 2924 2925 /*@ 2926 TSGetMaxTime - Gets the maximum (or final) time for timestepping. 2927 2928 Not Collective 2929 2930 Input Parameter: 2931 . ts - the `TS` context obtained from `TSCreate()` 2932 2933 Output Parameter: 2934 . maxtime - final time to step to 2935 2936 Level: advanced 2937 2938 .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()` 2939 @*/ 2940 PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime) 2941 { 2942 PetscFunctionBegin; 2943 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2944 PetscValidRealPointer(maxtime, 2); 2945 *maxtime = ts->max_time; 2946 PetscFunctionReturn(PETSC_SUCCESS); 2947 } 2948 2949 /*@ 2950 TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`. 2951 2952 Level: deprecated 2953 2954 @*/ 2955 PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step) 2956 { 2957 PetscFunctionBegin; 2958 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2959 PetscCall(TSSetTime(ts, initial_time)); 2960 PetscCall(TSSetTimeStep(ts, time_step)); 2961 PetscFunctionReturn(PETSC_SUCCESS); 2962 } 2963 2964 /*@ 2965 TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`. 2966 2967 Level: deprecated 2968 2969 @*/ 2970 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 2971 { 2972 PetscFunctionBegin; 2973 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2974 if (maxsteps) { 2975 PetscValidIntPointer(maxsteps, 2); 2976 *maxsteps = ts->max_steps; 2977 } 2978 if (maxtime) { 2979 PetscValidRealPointer(maxtime, 3); 2980 *maxtime = ts->max_time; 2981 } 2982 PetscFunctionReturn(PETSC_SUCCESS); 2983 } 2984 2985 /*@ 2986 TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`. 2987 2988 Level: deprecated 2989 2990 @*/ 2991 PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime) 2992 { 2993 PetscFunctionBegin; 2994 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 2995 PetscValidLogicalCollectiveInt(ts, maxsteps, 2); 2996 PetscValidLogicalCollectiveReal(ts, maxtime, 3); 2997 if (maxsteps >= 0) ts->max_steps = maxsteps; 2998 if (maxtime != (PetscReal)PETSC_DEFAULT) ts->max_time = maxtime; 2999 PetscFunctionReturn(PETSC_SUCCESS); 3000 } 3001 3002 /*@ 3003 TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`. 3004 3005 Level: deprecated 3006 3007 @*/ 3008 PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps) 3009 { 3010 return TSGetStepNumber(ts, steps); 3011 } 3012 3013 /*@ 3014 TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`. 3015 3016 Level: deprecated 3017 3018 @*/ 3019 PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps) 3020 { 3021 return TSGetStepNumber(ts, steps); 3022 } 3023 3024 /*@ 3025 TSSetSolution - Sets the initial solution vector 3026 for use by the `TS` routines. 3027 3028 Logically Collective 3029 3030 Input Parameters: 3031 + ts - the `TS` context obtained from `TSCreate()` 3032 - u - the solution vector 3033 3034 Level: beginner 3035 3036 .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()` 3037 @*/ 3038 PetscErrorCode TSSetSolution(TS ts, Vec u) 3039 { 3040 DM dm; 3041 3042 PetscFunctionBegin; 3043 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3044 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 3045 PetscCall(PetscObjectReference((PetscObject)u)); 3046 PetscCall(VecDestroy(&ts->vec_sol)); 3047 ts->vec_sol = u; 3048 3049 PetscCall(TSGetDM(ts, &dm)); 3050 PetscCall(DMShellSetGlobalVector(dm, u)); 3051 PetscFunctionReturn(PETSC_SUCCESS); 3052 } 3053 3054 /*@C 3055 TSSetPreStep - Sets the general-purpose function 3056 called once at the beginning of each time step. 3057 3058 Logically Collective 3059 3060 Input Parameters: 3061 + ts - The `TS` context obtained from `TSCreate()` 3062 - func - The function 3063 3064 Calling sequence of `func`: 3065 .vb 3066 PetscErrorCode func (TS ts) 3067 .ve 3068 3069 Level: intermediate 3070 3071 .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()` 3072 @*/ 3073 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 3074 { 3075 PetscFunctionBegin; 3076 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3077 ts->prestep = func; 3078 PetscFunctionReturn(PETSC_SUCCESS); 3079 } 3080 3081 /*@ 3082 TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()` 3083 3084 Collective 3085 3086 Input Parameter: 3087 . ts - The `TS` context obtained from `TSCreate()` 3088 3089 Level: developer 3090 3091 Note: 3092 `TSPreStep()` is typically used within time stepping implementations, 3093 so most users would not generally call this routine themselves. 3094 3095 .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()` 3096 @*/ 3097 PetscErrorCode TSPreStep(TS ts) 3098 { 3099 PetscFunctionBegin; 3100 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3101 if (ts->prestep) { 3102 Vec U; 3103 PetscObjectId idprev; 3104 PetscBool sameObject; 3105 PetscObjectState sprev, spost; 3106 3107 PetscCall(TSGetSolution(ts, &U)); 3108 PetscCall(PetscObjectGetId((PetscObject)U, &idprev)); 3109 PetscCall(PetscObjectStateGet((PetscObject)U, &sprev)); 3110 PetscCallBack("TS callback preset", (*ts->prestep)(ts)); 3111 PetscCall(TSGetSolution(ts, &U)); 3112 PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject)); 3113 PetscCall(PetscObjectStateGet((PetscObject)U, &spost)); 3114 if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts)); 3115 } 3116 PetscFunctionReturn(PETSC_SUCCESS); 3117 } 3118 3119 /*@C 3120 TSSetPreStage - Sets the general-purpose function 3121 called once at the beginning of each stage. 3122 3123 Logically Collective 3124 3125 Input Parameters: 3126 + ts - The `TS` context obtained from `TSCreate()` 3127 - func - The function 3128 3129 Calling sequence of `func`: 3130 .vb 3131 PetscErrorCode func(TS ts, PetscReal stagetime) 3132 .ve 3133 3134 Level: intermediate 3135 3136 Note: 3137 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 3138 The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being 3139 attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`. 3140 3141 .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()` 3142 @*/ 3143 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS, PetscReal)) 3144 { 3145 PetscFunctionBegin; 3146 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3147 ts->prestage = func; 3148 PetscFunctionReturn(PETSC_SUCCESS); 3149 } 3150 3151 /*@C 3152 TSSetPostStage - Sets the general-purpose function, provided with `TSSetPostStep()`, 3153 called once at the end of each stage. 3154 3155 Logically Collective 3156 3157 Input Parameters: 3158 + ts - The `TS` context obtained from `TSCreate()` 3159 - func - The function 3160 3161 Calling sequence of `func`: 3162 .vb 3163 PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y) 3164 .ve 3165 3166 Level: intermediate 3167 3168 Note: 3169 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 3170 The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being 3171 attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`. 3172 3173 .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()` 3174 @*/ 3175 PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS, PetscReal, PetscInt, Vec *)) 3176 { 3177 PetscFunctionBegin; 3178 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3179 ts->poststage = func; 3180 PetscFunctionReturn(PETSC_SUCCESS); 3181 } 3182 3183 /*@C 3184 TSSetPostEvaluate - Sets the general-purpose function 3185 called once at the end of each step evaluation. 3186 3187 Logically Collective 3188 3189 Input Parameters: 3190 + ts - The `TS` context obtained from `TSCreate()` 3191 - func - The function 3192 3193 Calling sequence of `func`: 3194 .vb 3195 PetscErrorCode func(TS ts) 3196 .ve 3197 3198 Level: intermediate 3199 3200 Note: 3201 Semantically, `TSSetPostEvaluate()` differs from `TSSetPostStep()` since the function it sets is called before event-handling 3202 thus guaranteeing the same solution (computed by the time-stepper) will be passed to it. On the other hand, `TSPostStep()` 3203 may be passed a different solution, possibly changed by the event handler. `TSPostEvaluate()` is called after the next step 3204 solution is evaluated allowing to modify it, if need be. The solution can be obtained with `TSGetSolution()`, the time step 3205 with `TSGetTimeStep()`, and the time at the start of the step is available via `TSGetTime()` 3206 3207 .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()` 3208 @*/ 3209 PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS)) 3210 { 3211 PetscFunctionBegin; 3212 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3213 ts->postevaluate = func; 3214 PetscFunctionReturn(PETSC_SUCCESS); 3215 } 3216 3217 /*@ 3218 TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()` 3219 3220 Collective 3221 3222 Input Parameters: 3223 . ts - The `TS` context obtained from `TSCreate()` 3224 stagetime - The absolute time of the current stage 3225 3226 Level: developer 3227 3228 Note: 3229 `TSPreStage()` is typically used within time stepping implementations, 3230 most users would not generally call this routine themselves. 3231 3232 .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()` 3233 @*/ 3234 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 3235 { 3236 PetscFunctionBegin; 3237 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3238 if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime)); 3239 PetscFunctionReturn(PETSC_SUCCESS); 3240 } 3241 3242 /*@ 3243 TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()` 3244 3245 Collective 3246 3247 Input Parameters: 3248 . ts - The `TS` context obtained from `TSCreate()` 3249 stagetime - The absolute time of the current stage 3250 stageindex - Stage number 3251 Y - Array of vectors (of size = total number 3252 of stages) with the stage solutions 3253 3254 Level: developer 3255 3256 Note: 3257 `TSPostStage()` is typically used within time stepping implementations, 3258 most users would not generally call this routine themselves. 3259 3260 .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()` 3261 @*/ 3262 PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y) 3263 { 3264 PetscFunctionBegin; 3265 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3266 if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y)); 3267 PetscFunctionReturn(PETSC_SUCCESS); 3268 } 3269 3270 /*@ 3271 TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()` 3272 3273 Collective 3274 3275 Input Parameter: 3276 . ts - The `TS` context obtained from `TSCreate()` 3277 3278 Level: developer 3279 3280 Note: 3281 `TSPostEvaluate()` is typically used within time stepping implementations, 3282 most users would not generally call this routine themselves. 3283 3284 .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()` 3285 @*/ 3286 PetscErrorCode TSPostEvaluate(TS ts) 3287 { 3288 PetscFunctionBegin; 3289 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3290 if (ts->postevaluate) { 3291 Vec U; 3292 PetscObjectState sprev, spost; 3293 3294 PetscCall(TSGetSolution(ts, &U)); 3295 PetscCall(PetscObjectStateGet((PetscObject)U, &sprev)); 3296 PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts)); 3297 PetscCall(PetscObjectStateGet((PetscObject)U, &spost)); 3298 if (sprev != spost) PetscCall(TSRestartStep(ts)); 3299 } 3300 PetscFunctionReturn(PETSC_SUCCESS); 3301 } 3302 3303 /*@C 3304 TSSetPostStep - Sets the general-purpose function 3305 called once at the end of each time step. 3306 3307 Logically Collective 3308 3309 Input Parameters: 3310 + ts - The `TS` context obtained from `TSCreate()` 3311 - func - The function 3312 3313 Calling sequence of `func`: 3314 $ PetscErrorCode func(TS ts) 3315 3316 Level: intermediate 3317 3318 Note: 3319 The function set by `TSSetPostStep()` is called after each successful step. The solution vector 3320 obtained by `TSGetSolution()` may be different than that computed at the step end if the event handler 3321 locates an event and `TSPostEvent()` modifies it. Use `TSSetPostEvaluate()` if an unmodified solution is needed instead. 3322 3323 .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()` 3324 @*/ 3325 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 3326 { 3327 PetscFunctionBegin; 3328 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3329 ts->poststep = func; 3330 PetscFunctionReturn(PETSC_SUCCESS); 3331 } 3332 3333 /*@ 3334 TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()` 3335 3336 Collective 3337 3338 Input Parameter: 3339 . ts - The `TS` context obtained from `TSCreate()` 3340 3341 Note: 3342 `TSPostStep()` is typically used within time stepping implementations, 3343 so most users would not generally call this routine themselves. 3344 3345 Level: developer 3346 3347 .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()` 3348 @*/ 3349 PetscErrorCode TSPostStep(TS ts) 3350 { 3351 PetscFunctionBegin; 3352 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3353 if (ts->poststep) { 3354 Vec U; 3355 PetscObjectId idprev; 3356 PetscBool sameObject; 3357 PetscObjectState sprev, spost; 3358 3359 PetscCall(TSGetSolution(ts, &U)); 3360 PetscCall(PetscObjectGetId((PetscObject)U, &idprev)); 3361 PetscCall(PetscObjectStateGet((PetscObject)U, &sprev)); 3362 PetscCallBack("TS callback poststep", (*ts->poststep)(ts)); 3363 PetscCall(TSGetSolution(ts, &U)); 3364 PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject)); 3365 PetscCall(PetscObjectStateGet((PetscObject)U, &spost)); 3366 if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts)); 3367 } 3368 PetscFunctionReturn(PETSC_SUCCESS); 3369 } 3370 3371 /*@ 3372 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 3373 3374 Collective 3375 3376 Input Parameters: 3377 + ts - time stepping context 3378 - t - time to interpolate to 3379 3380 Output Parameter: 3381 . U - state at given time 3382 3383 Level: intermediate 3384 3385 Developer Note: 3386 `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 3387 3388 .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()` 3389 @*/ 3390 PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U) 3391 { 3392 PetscFunctionBegin; 3393 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3394 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 3395 PetscCheck(t >= ts->ptime_prev && t <= ts->ptime, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Requested time %g not in last time steps [%g,%g]", (double)t, (double)ts->ptime_prev, (double)ts->ptime); 3396 PetscUseTypeMethod(ts, interpolate, t, U); 3397 PetscFunctionReturn(PETSC_SUCCESS); 3398 } 3399 3400 /*@ 3401 TSStep - Steps one time step 3402 3403 Collective 3404 3405 Input Parameter: 3406 . ts - the `TS` context obtained from `TSCreate()` 3407 3408 Level: developer 3409 3410 Notes: 3411 The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine. 3412 3413 The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may 3414 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 3415 3416 This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the 3417 time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep. 3418 3419 .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()` 3420 @*/ 3421 PetscErrorCode TSStep(TS ts) 3422 { 3423 static PetscBool cite = PETSC_FALSE; 3424 PetscReal ptime; 3425 3426 PetscFunctionBegin; 3427 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3428 PetscCall(PetscCitationsRegister("@article{tspaper,\n" 3429 " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n" 3430 " author = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n" 3431 " journal = {arXiv e-preprints},\n" 3432 " eprint = {1806.01437},\n" 3433 " archivePrefix = {arXiv},\n" 3434 " year = {2018}\n}\n", 3435 &cite)); 3436 PetscCall(TSSetUp(ts)); 3437 PetscCall(TSTrajectorySetUp(ts->trajectory, ts)); 3438 3439 PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_MAX_INT, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>"); 3440 PetscCheck(ts->exact_final_time != TS_EXACTFINALTIME_UNSPECIFIED, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetExactFinalTime() or use -ts_exact_final_time <stepover,interpolate,matchstep> before calling TSStep()"); 3441 PetscCheck(ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP || ts->adapt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Since TS is not adaptive you cannot use TS_EXACTFINALTIME_MATCHSTEP, suggest TS_EXACTFINALTIME_INTERPOLATE"); 3442 3443 if (!ts->steps) ts->ptime_prev = ts->ptime; 3444 ptime = ts->ptime; 3445 ts->ptime_prev_rollback = ts->ptime_prev; 3446 ts->reason = TS_CONVERGED_ITERATING; 3447 3448 PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0)); 3449 PetscUseTypeMethod(ts, step); 3450 PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0)); 3451 3452 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) 3453 PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[ts->tspan->spanctr++])); 3454 if (ts->reason >= 0) { 3455 ts->ptime_prev = ptime; 3456 ts->steps++; 3457 ts->steprollback = PETSC_FALSE; 3458 ts->steprestart = PETSC_FALSE; 3459 } 3460 if (!ts->reason) { 3461 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 3462 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3463 } 3464 3465 if (ts->reason < 0 && ts->errorifstepfailed) { 3466 PetscCall(TSMonitorCancel(ts)); 3467 PetscCheck(ts->reason != TS_DIVERGED_NONLINEAR_SOLVE, PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s, increase -ts_max_snes_failures or make negative to attempt recovery", TSConvergedReasons[ts->reason]); 3468 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]); 3469 } 3470 PetscFunctionReturn(PETSC_SUCCESS); 3471 } 3472 3473 /*@ 3474 TSEvaluateWLTE - Evaluate the weighted local truncation error norm 3475 at the end of a time step with a given order of accuracy. 3476 3477 Collective 3478 3479 Input Parameters: 3480 + ts - time stepping context 3481 - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY` 3482 3483 Input/Output Parameter: 3484 . order - optional, desired order for the error evaluation or `PETSC_DECIDE`; 3485 on output, the actual order of the error evaluation 3486 3487 Output Parameter: 3488 . wlte - the weighted local truncation error norm 3489 3490 Level: advanced 3491 3492 Note: 3493 If the timestepper cannot evaluate the error in a particular step 3494 (eg. in the first step or restart steps after event handling), 3495 this routine returns wlte=-1.0 . 3496 3497 .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()` 3498 @*/ 3499 PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte) 3500 { 3501 PetscFunctionBegin; 3502 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3503 PetscValidType(ts, 1); 3504 PetscValidLogicalCollectiveEnum(ts, wnormtype, 2); 3505 if (order) PetscValidIntPointer(order, 3); 3506 if (order) PetscValidLogicalCollectiveInt(ts, *order, 3); 3507 PetscValidRealPointer(wlte, 4); 3508 PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]); 3509 PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte); 3510 PetscFunctionReturn(PETSC_SUCCESS); 3511 } 3512 3513 /*@ 3514 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 3515 3516 Collective 3517 3518 Input Parameters: 3519 + ts - time stepping context 3520 . order - desired order of accuracy 3521 - done - whether the step was evaluated at this order (pass `NULL` to generate an error if not available) 3522 3523 Output Parameter: 3524 . U - state at the end of the current step 3525 3526 Level: advanced 3527 3528 Notes: 3529 This function cannot be called until all stages have been evaluated. 3530 3531 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. 3532 3533 .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt` 3534 @*/ 3535 PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done) 3536 { 3537 PetscFunctionBegin; 3538 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3539 PetscValidType(ts, 1); 3540 PetscValidHeaderSpecific(U, VEC_CLASSID, 3); 3541 PetscUseTypeMethod(ts, evaluatestep, order, U, done); 3542 PetscFunctionReturn(PETSC_SUCCESS); 3543 } 3544 3545 /*@C 3546 TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping. 3547 3548 Not collective 3549 3550 Input Parameter: 3551 . ts - time stepping context 3552 3553 Output Parameter: 3554 . initConditions - The function which computes an initial condition 3555 3556 Calling sequence of `initCondition`: 3557 .vb 3558 PetscErrorCode initCondition(TS ts, Vec u) 3559 .ve 3560 + ts - The timestepping context 3561 - u - The input vector in which the initial condition is stored 3562 3563 Level: advanced 3564 3565 .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()` 3566 @*/ 3567 PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS, Vec)) 3568 { 3569 PetscFunctionBegin; 3570 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3571 PetscValidPointer(initCondition, 2); 3572 *initCondition = ts->ops->initcondition; 3573 PetscFunctionReturn(PETSC_SUCCESS); 3574 } 3575 3576 /*@C 3577 TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping. 3578 3579 Logically collective 3580 3581 Input Parameters: 3582 + ts - time stepping context 3583 - initCondition - The function which computes an initial condition 3584 3585 Calling sequence of `initCondition`: 3586 $ PetscErrorCode initCondition(TS ts, Vec u) 3587 + ts - The timestepping context 3588 - u - The input vector in which the initial condition is to be stored 3589 3590 Level: advanced 3591 3592 .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()` 3593 @*/ 3594 PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS, Vec)) 3595 { 3596 PetscFunctionBegin; 3597 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3598 PetscValidFunction(initCondition, 2); 3599 ts->ops->initcondition = initCondition; 3600 PetscFunctionReturn(PETSC_SUCCESS); 3601 } 3602 3603 /*@ 3604 TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()` 3605 3606 Collective 3607 3608 Input Parameters: 3609 + ts - time stepping context 3610 - u - The `Vec` to store the condition in which will be used in `TSSolve()` 3611 3612 Level: advanced 3613 3614 .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()` 3615 @*/ 3616 PetscErrorCode TSComputeInitialCondition(TS ts, Vec u) 3617 { 3618 PetscFunctionBegin; 3619 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3620 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 3621 PetscTryTypeMethod(ts, initcondition, u); 3622 PetscFunctionReturn(PETSC_SUCCESS); 3623 } 3624 3625 /*@C 3626 TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping. 3627 3628 Not collective 3629 3630 Input Parameter: 3631 . ts - time stepping context 3632 3633 Output Parameter: 3634 . exactError - The function which computes the solution error 3635 3636 Calling sequence of `exactError`: 3637 $ PetscErrorCode exactError(TS ts, Vec u, Vec e) 3638 + ts - The timestepping context 3639 . u - The approximate solution vector 3640 - e - The vector in which the error is stored 3641 3642 Level: advanced 3643 3644 .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()` 3645 @*/ 3646 PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS, Vec, Vec)) 3647 { 3648 PetscFunctionBegin; 3649 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3650 PetscValidPointer(exactError, 2); 3651 *exactError = ts->ops->exacterror; 3652 PetscFunctionReturn(PETSC_SUCCESS); 3653 } 3654 3655 /*@C 3656 TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping. 3657 3658 Logically collective 3659 3660 Input Parameters: 3661 + ts - time stepping context 3662 - exactError - The function which computes the solution error 3663 3664 Calling sequence of `exactError`: 3665 $ PetscErrorCode exactError(TS ts, Vec u) 3666 + ts - The timestepping context 3667 . u - The approximate solution vector 3668 - e - The vector in which the error is stored 3669 3670 Level: advanced 3671 3672 .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()` 3673 @*/ 3674 PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS, Vec, Vec)) 3675 { 3676 PetscFunctionBegin; 3677 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3678 PetscValidFunction(exactError, 2); 3679 ts->ops->exacterror = exactError; 3680 PetscFunctionReturn(PETSC_SUCCESS); 3681 } 3682 3683 /*@ 3684 TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()` 3685 3686 Collective 3687 3688 Input Parameters: 3689 + ts - time stepping context 3690 . u - The approximate solution 3691 - e - The `Vec` used to store the error 3692 3693 Level: advanced 3694 3695 .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()` 3696 @*/ 3697 PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e) 3698 { 3699 PetscFunctionBegin; 3700 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3701 PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 3702 PetscValidHeaderSpecific(e, VEC_CLASSID, 3); 3703 PetscTryTypeMethod(ts, exacterror, u, e); 3704 PetscFunctionReturn(PETSC_SUCCESS); 3705 } 3706 3707 /*@ 3708 TSSolve - Steps the requested number of timesteps. 3709 3710 Collective 3711 3712 Input Parameters: 3713 + ts - the `TS` context obtained from `TSCreate()` 3714 - u - the solution vector (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used, 3715 otherwise must contain the initial conditions and will contain the solution at the final requested time 3716 3717 Level: beginner 3718 3719 Notes: 3720 The final time returned by this function may be different from the time of the internally 3721 held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have 3722 stepped over the final time. 3723 3724 .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()` 3725 @*/ 3726 PetscErrorCode TSSolve(TS ts, Vec u) 3727 { 3728 Vec solution; 3729 3730 PetscFunctionBegin; 3731 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3732 if (u) PetscValidHeaderSpecific(u, VEC_CLASSID, 2); 3733 3734 PetscCall(TSSetExactFinalTimeDefault(ts)); 3735 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 */ 3736 if (!ts->vec_sol || u == ts->vec_sol) { 3737 PetscCall(VecDuplicate(u, &solution)); 3738 PetscCall(TSSetSolution(ts, solution)); 3739 PetscCall(VecDestroy(&solution)); /* grant ownership */ 3740 } 3741 PetscCall(VecCopy(u, ts->vec_sol)); 3742 PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE"); 3743 } else if (u) PetscCall(TSSetSolution(ts, u)); 3744 PetscCall(TSSetUp(ts)); 3745 PetscCall(TSTrajectorySetUp(ts->trajectory, ts)); 3746 3747 PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_MAX_INT, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>"); 3748 PetscCheck(ts->exact_final_time != TS_EXACTFINALTIME_UNSPECIFIED, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetExactFinalTime() or use -ts_exact_final_time <stepover,interpolate,matchstep> before calling TSSolve()"); 3749 PetscCheck(ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP || ts->adapt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Since TS is not adaptive you cannot use TS_EXACTFINALTIME_MATCHSTEP, suggest TS_EXACTFINALTIME_INTERPOLATE"); 3750 PetscCheck(!(ts->tspan && ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP), PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "You must use TS_EXACTFINALTIME_MATCHSTEP when using time span"); 3751 3752 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 */ 3753 PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0])); 3754 ts->tspan->spanctr = 1; 3755 } 3756 3757 if (ts->forward_solve) PetscCall(TSForwardSetUp(ts)); 3758 3759 /* reset number of steps only when the step is not restarted. ARKIMEX 3760 restarts the step after an event. Resetting these counters in such case causes 3761 TSTrajectory to incorrectly save the output files 3762 */ 3763 /* reset time step and iteration counters */ 3764 if (!ts->steps) { 3765 ts->ksp_its = 0; 3766 ts->snes_its = 0; 3767 ts->num_snes_failures = 0; 3768 ts->reject = 0; 3769 ts->steprestart = PETSC_TRUE; 3770 ts->steprollback = PETSC_FALSE; 3771 ts->rhsjacobian.time = PETSC_MIN_REAL; 3772 } 3773 3774 /* make sure initial time step does not overshoot final time or the next point in tspan */ 3775 if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) { 3776 PetscReal maxdt; 3777 PetscReal dt = ts->time_step; 3778 3779 if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime; 3780 else maxdt = ts->max_time - ts->ptime; 3781 ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt); 3782 } 3783 ts->reason = TS_CONVERGED_ITERATING; 3784 3785 { 3786 PetscViewer viewer; 3787 PetscViewerFormat format; 3788 PetscBool flg; 3789 static PetscBool incall = PETSC_FALSE; 3790 3791 if (!incall) { 3792 /* Estimate the convergence rate of the time discretization */ 3793 PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg)); 3794 if (flg) { 3795 PetscConvEst conv; 3796 DM dm; 3797 PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */ 3798 PetscInt Nf; 3799 PetscBool checkTemporal = PETSC_TRUE; 3800 3801 incall = PETSC_TRUE; 3802 PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg)); 3803 PetscCall(TSGetDM(ts, &dm)); 3804 PetscCall(DMGetNumFields(dm, &Nf)); 3805 PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha)); 3806 PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv)); 3807 PetscCall(PetscConvEstUseTS(conv, checkTemporal)); 3808 PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts)); 3809 PetscCall(PetscConvEstSetFromOptions(conv)); 3810 PetscCall(PetscConvEstSetUp(conv)); 3811 PetscCall(PetscConvEstGetConvRate(conv, alpha)); 3812 PetscCall(PetscViewerPushFormat(viewer, format)); 3813 PetscCall(PetscConvEstRateView(conv, alpha, viewer)); 3814 PetscCall(PetscViewerPopFormat(viewer)); 3815 PetscCall(PetscViewerDestroy(&viewer)); 3816 PetscCall(PetscConvEstDestroy(&conv)); 3817 PetscCall(PetscFree(alpha)); 3818 incall = PETSC_FALSE; 3819 } 3820 } 3821 } 3822 3823 PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre")); 3824 3825 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 3826 PetscUseTypeMethod(ts, solve); 3827 if (u) PetscCall(VecCopy(ts->vec_sol, u)); 3828 ts->solvetime = ts->ptime; 3829 solution = ts->vec_sol; 3830 } else { /* Step the requested number of timesteps. */ 3831 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 3832 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3833 3834 if (!ts->steps) { 3835 PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol)); 3836 PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol)); 3837 } 3838 3839 while (!ts->reason) { 3840 PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol)); 3841 if (!ts->steprollback) PetscCall(TSPreStep(ts)); 3842 PetscCall(TSStep(ts)); 3843 if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL)); 3844 if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL)); 3845 if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */ 3846 if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */ 3847 PetscCall(TSForwardCostIntegral(ts)); 3848 if (ts->reason >= 0) ts->steps++; 3849 } 3850 if (ts->forward_solve) { /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */ 3851 if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */ 3852 PetscCall(TSForwardStep(ts)); 3853 if (ts->reason >= 0) ts->steps++; 3854 } 3855 PetscCall(TSPostEvaluate(ts)); 3856 PetscCall(TSEventHandler(ts)); /* The right-hand side may be changed due to event. Be careful with Any computation using the RHS information after this point. */ 3857 if (ts->steprollback) PetscCall(TSPostEvaluate(ts)); 3858 if (!ts->steprollback) { 3859 PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol)); 3860 PetscCall(TSPostStep(ts)); 3861 } 3862 } 3863 PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol)); 3864 3865 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 3866 PetscCall(TSInterpolate(ts, ts->max_time, u)); 3867 ts->solvetime = ts->max_time; 3868 solution = u; 3869 PetscCall(TSMonitor(ts, -1, ts->solvetime, solution)); 3870 } else { 3871 if (u) PetscCall(VecCopy(ts->vec_sol, u)); 3872 ts->solvetime = ts->ptime; 3873 solution = ts->vec_sol; 3874 } 3875 } 3876 3877 PetscCall(TSViewFromOptions(ts, NULL, "-ts_view")); 3878 PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution")); 3879 PetscCall(PetscObjectSAWsBlock((PetscObject)ts)); 3880 if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts)); 3881 PetscFunctionReturn(PETSC_SUCCESS); 3882 } 3883 3884 /*@ 3885 TSGetTime - Gets the time of the most recently completed step. 3886 3887 Not Collective 3888 3889 Input Parameter: 3890 . ts - the `TS` context obtained from `TSCreate()` 3891 3892 Output Parameter: 3893 . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`. 3894 3895 Level: beginner 3896 3897 Note: 3898 When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`, 3899 `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated. 3900 3901 .seealso: [](ch_ts), TS`, ``TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()` 3902 @*/ 3903 PetscErrorCode TSGetTime(TS ts, PetscReal *t) 3904 { 3905 PetscFunctionBegin; 3906 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3907 PetscValidRealPointer(t, 2); 3908 *t = ts->ptime; 3909 PetscFunctionReturn(PETSC_SUCCESS); 3910 } 3911 3912 /*@ 3913 TSGetPrevTime - Gets the starting time of the previously completed step. 3914 3915 Not Collective 3916 3917 Input Parameter: 3918 . ts - the `TS` context obtained from `TSCreate()` 3919 3920 Output Parameter: 3921 . t - the previous time 3922 3923 Level: beginner 3924 3925 .seealso: [](ch_ts), TS`, ``TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()` 3926 @*/ 3927 PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t) 3928 { 3929 PetscFunctionBegin; 3930 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3931 PetscValidRealPointer(t, 2); 3932 *t = ts->ptime_prev; 3933 PetscFunctionReturn(PETSC_SUCCESS); 3934 } 3935 3936 /*@ 3937 TSSetTime - Allows one to reset the time. 3938 3939 Logically Collective 3940 3941 Input Parameters: 3942 + ts - the `TS` context obtained from `TSCreate()` 3943 - time - the time 3944 3945 Level: intermediate 3946 3947 .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()` 3948 @*/ 3949 PetscErrorCode TSSetTime(TS ts, PetscReal t) 3950 { 3951 PetscFunctionBegin; 3952 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3953 PetscValidLogicalCollectiveReal(ts, t, 2); 3954 ts->ptime = t; 3955 PetscFunctionReturn(PETSC_SUCCESS); 3956 } 3957 3958 /*@C 3959 TSSetOptionsPrefix - Sets the prefix used for searching for all 3960 TS options in the database. 3961 3962 Logically Collective 3963 3964 Input Parameters: 3965 + ts - The `TS` context 3966 - prefix - The prefix to prepend to all option names 3967 3968 Level: advanced 3969 3970 Note: 3971 A hyphen (-) must NOT be given at the beginning of the prefix name. 3972 The first character of all runtime options is AUTOMATICALLY the 3973 hyphen. 3974 3975 .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()` 3976 @*/ 3977 PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[]) 3978 { 3979 SNES snes; 3980 3981 PetscFunctionBegin; 3982 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 3983 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix)); 3984 PetscCall(TSGetSNES(ts, &snes)); 3985 PetscCall(SNESSetOptionsPrefix(snes, prefix)); 3986 PetscFunctionReturn(PETSC_SUCCESS); 3987 } 3988 3989 /*@C 3990 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 3991 TS options in the database. 3992 3993 Logically Collective 3994 3995 Input Parameters: 3996 + ts - The `TS` context 3997 - prefix - The prefix to prepend to all option names 3998 3999 Level: advanced 4000 4001 Note: 4002 A hyphen (-) must NOT be given at the beginning of the prefix name. 4003 The first character of all runtime options is AUTOMATICALLY the 4004 hyphen. 4005 4006 .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()` 4007 @*/ 4008 PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[]) 4009 { 4010 SNES snes; 4011 4012 PetscFunctionBegin; 4013 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4014 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix)); 4015 PetscCall(TSGetSNES(ts, &snes)); 4016 PetscCall(SNESAppendOptionsPrefix(snes, prefix)); 4017 PetscFunctionReturn(PETSC_SUCCESS); 4018 } 4019 4020 /*@C 4021 TSGetOptionsPrefix - Sets the prefix used for searching for all 4022 `TS` options in the database. 4023 4024 Not Collective 4025 4026 Input Parameter: 4027 . ts - The `TS` context 4028 4029 Output Parameter: 4030 . prefix - A pointer to the prefix string used 4031 4032 Level: intermediate 4033 4034 Fortran Note: 4035 The user should pass in a string 'prefix' of 4036 sufficient length to hold the prefix. 4037 4038 .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()` 4039 @*/ 4040 PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[]) 4041 { 4042 PetscFunctionBegin; 4043 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4044 PetscValidPointer(prefix, 2); 4045 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix)); 4046 PetscFunctionReturn(PETSC_SUCCESS); 4047 } 4048 4049 /*@C 4050 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 4051 4052 Not Collective, but parallel objects are returned if ts is parallel 4053 4054 Input Parameter: 4055 . ts - The `TS` context obtained from `TSCreate()` 4056 4057 Output Parameters: 4058 + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or `NULL`) 4059 . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat` (or `NULL`) 4060 . func - Function to compute the Jacobian of the RHS (or `NULL`) 4061 - ctx - User-defined context for Jacobian evaluation routine (or `NULL`) 4062 4063 Level: intermediate 4064 4065 Note: 4066 You can pass in `NULL` for any return argument you do not need. 4067 4068 .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()` 4069 4070 @*/ 4071 PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobian *func, void **ctx) 4072 { 4073 DM dm; 4074 4075 PetscFunctionBegin; 4076 if (Amat || Pmat) { 4077 SNES snes; 4078 PetscCall(TSGetSNES(ts, &snes)); 4079 PetscCall(SNESSetUpMatrices(snes)); 4080 PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL)); 4081 } 4082 PetscCall(TSGetDM(ts, &dm)); 4083 PetscCall(DMTSGetRHSJacobian(dm, func, ctx)); 4084 PetscFunctionReturn(PETSC_SUCCESS); 4085 } 4086 4087 /*@C 4088 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 4089 4090 Not Collective, but parallel objects are returned if ts is parallel 4091 4092 Input Parameter: 4093 . ts - The `TS` context obtained from `TSCreate()` 4094 4095 Output Parameters: 4096 + Amat - The (approximate) Jacobian of F(t,U,U_t) 4097 . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat` 4098 . f - The function to compute the matrices 4099 - ctx - User-defined context for Jacobian evaluation routine 4100 4101 Level: advanced 4102 4103 Note: 4104 You can pass in `NULL` for any return argument you do not need. 4105 4106 .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()` 4107 @*/ 4108 PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobian *f, void **ctx) 4109 { 4110 DM dm; 4111 4112 PetscFunctionBegin; 4113 if (Amat || Pmat) { 4114 SNES snes; 4115 PetscCall(TSGetSNES(ts, &snes)); 4116 PetscCall(SNESSetUpMatrices(snes)); 4117 PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL)); 4118 } 4119 PetscCall(TSGetDM(ts, &dm)); 4120 PetscCall(DMTSGetIJacobian(dm, f, ctx)); 4121 PetscFunctionReturn(PETSC_SUCCESS); 4122 } 4123 4124 #include <petsc/private/dmimpl.h> 4125 /*@ 4126 TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS` 4127 4128 Logically Collective 4129 4130 Input Parameters: 4131 + ts - the `TS` integrator object 4132 - dm - the dm, cannot be `NULL` 4133 4134 Level: intermediate 4135 4136 Notes: 4137 A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`, 4138 even when not using interfaces like `DMTSSetIFunction()`. Use `DMClone()` to get a distinct `DM` when solving 4139 different problems using the same function space. 4140 4141 .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()` 4142 @*/ 4143 PetscErrorCode TSSetDM(TS ts, DM dm) 4144 { 4145 SNES snes; 4146 DMTS tsdm; 4147 4148 PetscFunctionBegin; 4149 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4150 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 4151 PetscCall(PetscObjectReference((PetscObject)dm)); 4152 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 4153 if (ts->dm->dmts && !dm->dmts) { 4154 PetscCall(DMCopyDMTS(ts->dm, dm)); 4155 PetscCall(DMGetDMTS(ts->dm, &tsdm)); 4156 /* Grant write privileges to the replacement DM */ 4157 if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm; 4158 } 4159 PetscCall(DMDestroy(&ts->dm)); 4160 } 4161 ts->dm = dm; 4162 4163 PetscCall(TSGetSNES(ts, &snes)); 4164 PetscCall(SNESSetDM(snes, dm)); 4165 PetscFunctionReturn(PETSC_SUCCESS); 4166 } 4167 4168 /*@ 4169 TSGetDM - Gets the `DM` that may be used by some preconditioners 4170 4171 Not Collective 4172 4173 Input Parameter: 4174 . ts - the `TS` 4175 4176 Output Parameter: 4177 . dm - the `DM` 4178 4179 Level: intermediate 4180 4181 .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()` 4182 @*/ 4183 PetscErrorCode TSGetDM(TS ts, DM *dm) 4184 { 4185 PetscFunctionBegin; 4186 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4187 if (!ts->dm) { 4188 PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm)); 4189 if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm)); 4190 } 4191 *dm = ts->dm; 4192 PetscFunctionReturn(PETSC_SUCCESS); 4193 } 4194 4195 /*@ 4196 SNESTSFormFunction - Function to evaluate nonlinear residual 4197 4198 Logically Collective 4199 4200 Input Parameters: 4201 + snes - nonlinear solver 4202 . U - the current state at which to evaluate the residual 4203 - ctx - user context, must be a TS 4204 4205 Output Parameter: 4206 . F - the nonlinear residual 4207 4208 Level: advanced 4209 4210 Note: 4211 This function is not normally called by users and is automatically registered with the `SNES` used by `TS`. 4212 It is most frequently passed to `MatFDColoringSetFunction()`. 4213 4214 .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()` 4215 @*/ 4216 PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx) 4217 { 4218 TS ts = (TS)ctx; 4219 4220 PetscFunctionBegin; 4221 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4222 PetscValidHeaderSpecific(U, VEC_CLASSID, 2); 4223 PetscValidHeaderSpecific(F, VEC_CLASSID, 3); 4224 PetscValidHeaderSpecific(ts, TS_CLASSID, 4); 4225 PetscCall((ts->ops->snesfunction)(snes, U, F, ts)); 4226 PetscFunctionReturn(PETSC_SUCCESS); 4227 } 4228 4229 /*@ 4230 SNESTSFormJacobian - Function to evaluate the Jacobian 4231 4232 Collective 4233 4234 Input Parameters: 4235 + snes - nonlinear solver 4236 . U - the current state at which to evaluate the residual 4237 - ctx - user context, must be a `TS` 4238 4239 Output Parameters: 4240 + A - the Jacobian 4241 - B - the preconditioning matrix (may be the same as A) 4242 4243 Level: developer 4244 4245 Note: 4246 This function is not normally called by users and is automatically registered with the `SNES` used by `TS`. 4247 4248 .seealso: [](ch_ts), `SNESSetJacobian()` 4249 @*/ 4250 PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx) 4251 { 4252 TS ts = (TS)ctx; 4253 4254 PetscFunctionBegin; 4255 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4256 PetscValidHeaderSpecific(U, VEC_CLASSID, 2); 4257 PetscValidPointer(A, 3); 4258 PetscValidHeaderSpecific(A, MAT_CLASSID, 3); 4259 PetscValidPointer(B, 4); 4260 PetscValidHeaderSpecific(B, MAT_CLASSID, 4); 4261 PetscValidHeaderSpecific(ts, TS_CLASSID, 5); 4262 PetscCall((ts->ops->snesjacobian)(snes, U, A, B, ts)); 4263 PetscFunctionReturn(PETSC_SUCCESS); 4264 } 4265 4266 /*@C 4267 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems Udot = A U only 4268 4269 Collective 4270 4271 Input Parameters: 4272 + ts - time stepping context 4273 . t - time at which to evaluate 4274 . U - state at which to evaluate 4275 - ctx - context 4276 4277 Output Parameter: 4278 . F - right hand side 4279 4280 Level: intermediate 4281 4282 Note: 4283 This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right hand side for linear problems. 4284 The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`. 4285 4286 .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()` 4287 @*/ 4288 PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx) 4289 { 4290 Mat Arhs, Brhs; 4291 4292 PetscFunctionBegin; 4293 PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs)); 4294 /* undo the damage caused by shifting */ 4295 PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs)); 4296 PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs)); 4297 PetscCall(MatMult(Arhs, U, F)); 4298 PetscFunctionReturn(PETSC_SUCCESS); 4299 } 4300 4301 /*@C 4302 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 4303 4304 Collective 4305 4306 Input Parameters: 4307 + ts - time stepping context 4308 . t - time at which to evaluate 4309 . U - state at which to evaluate 4310 - ctx - context 4311 4312 Output Parameters: 4313 + A - pointer to operator 4314 - B - pointer to preconditioning matrix 4315 4316 Level: intermediate 4317 4318 Note: 4319 This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems. 4320 4321 .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()` 4322 @*/ 4323 PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx) 4324 { 4325 PetscFunctionBegin; 4326 PetscFunctionReturn(PETSC_SUCCESS); 4327 } 4328 4329 /*@C 4330 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 4331 4332 Collective 4333 4334 Input Parameters: 4335 + ts - time stepping context 4336 . t - time at which to evaluate 4337 . U - state at which to evaluate 4338 . Udot - time derivative of state vector 4339 - ctx - context 4340 4341 Output Parameter: 4342 . F - left hand side 4343 4344 Level: intermediate 4345 4346 Notes: 4347 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 4348 user is required to write their own `TSComputeIFunction()`. 4349 This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems. 4350 The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`. 4351 4352 Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U 4353 4354 .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()` 4355 @*/ 4356 PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx) 4357 { 4358 Mat A, B; 4359 4360 PetscFunctionBegin; 4361 PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL)); 4362 PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE)); 4363 PetscCall(MatMult(A, Udot, F)); 4364 PetscFunctionReturn(PETSC_SUCCESS); 4365 } 4366 4367 /*@C 4368 TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE 4369 4370 Collective 4371 4372 Input Parameters: 4373 + ts - time stepping context 4374 . t - time at which to evaluate 4375 . U - state at which to evaluate 4376 . Udot - time derivative of state vector 4377 . shift - shift to apply 4378 - ctx - context 4379 4380 Output Parameters: 4381 + A - pointer to operator 4382 - B - pointer to preconditioning matrix 4383 4384 Level: advanced 4385 4386 Notes: 4387 This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems. 4388 4389 It is only appropriate for problems of the form 4390 4391 $ M Udot = F(U,t) 4392 4393 where M is constant and F is non-stiff. The user must pass M to `TSSetIJacobian()`. The current implementation only 4394 works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing 4395 an implicit operator of the form 4396 4397 $ shift*M + J 4398 4399 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 4400 a copy of M or reassemble it when requested. 4401 4402 .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()` 4403 @*/ 4404 PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx) 4405 { 4406 PetscFunctionBegin; 4407 PetscCall(MatScale(A, shift / ts->ijacobian.shift)); 4408 ts->ijacobian.shift = shift; 4409 PetscFunctionReturn(PETSC_SUCCESS); 4410 } 4411 4412 /*@ 4413 TSGetEquationType - Gets the type of the equation that `TS` is solving. 4414 4415 Not Collective 4416 4417 Input Parameter: 4418 . ts - the `TS` context 4419 4420 Output Parameter: 4421 . equation_type - see `TSEquationType` 4422 4423 Level: beginner 4424 4425 .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType` 4426 @*/ 4427 PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type) 4428 { 4429 PetscFunctionBegin; 4430 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4431 PetscValidPointer(equation_type, 2); 4432 *equation_type = ts->equation_type; 4433 PetscFunctionReturn(PETSC_SUCCESS); 4434 } 4435 4436 /*@ 4437 TSSetEquationType - Sets the type of the equation that `TS` is solving. 4438 4439 Not Collective 4440 4441 Input Parameters: 4442 + ts - the `TS` context 4443 - equation_type - see `TSEquationType` 4444 4445 Level: advanced 4446 4447 .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType` 4448 @*/ 4449 PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type) 4450 { 4451 PetscFunctionBegin; 4452 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4453 ts->equation_type = equation_type; 4454 PetscFunctionReturn(PETSC_SUCCESS); 4455 } 4456 4457 /*@ 4458 TSGetConvergedReason - Gets the reason the `TS` iteration was stopped. 4459 4460 Not Collective 4461 4462 Input Parameter: 4463 . ts - the `TS` context 4464 4465 Output Parameter: 4466 . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the 4467 manual pages for the individual convergence tests for complete lists 4468 4469 Level: beginner 4470 4471 Note: 4472 Can only be called after the call to `TSSolve()` is complete. 4473 4474 .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason` 4475 @*/ 4476 PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason) 4477 { 4478 PetscFunctionBegin; 4479 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4480 PetscValidPointer(reason, 2); 4481 *reason = ts->reason; 4482 PetscFunctionReturn(PETSC_SUCCESS); 4483 } 4484 4485 /*@ 4486 TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`. 4487 4488 Logically Collective; reason must contain common value 4489 4490 Input Parameters: 4491 + ts - the `TS` context 4492 - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the 4493 manual pages for the individual convergence tests for complete lists 4494 4495 Level: advanced 4496 4497 Note: 4498 Can only be called while `TSSolve()` is active. 4499 4500 .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason` 4501 @*/ 4502 PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason) 4503 { 4504 PetscFunctionBegin; 4505 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4506 ts->reason = reason; 4507 PetscFunctionReturn(PETSC_SUCCESS); 4508 } 4509 4510 /*@ 4511 TSGetSolveTime - Gets the time after a call to `TSSolve()` 4512 4513 Not Collective 4514 4515 Input Parameter: 4516 . ts - the `TS` context 4517 4518 Output Parameter: 4519 . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()` 4520 4521 Level: beginner 4522 4523 Note: 4524 Can only be called after the call to `TSSolve()` is complete. 4525 4526 .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSSetConvergenceTest()`, `TSConvergedReason` 4527 @*/ 4528 PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime) 4529 { 4530 PetscFunctionBegin; 4531 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4532 PetscValidRealPointer(ftime, 2); 4533 *ftime = ts->solvetime; 4534 PetscFunctionReturn(PETSC_SUCCESS); 4535 } 4536 4537 /*@ 4538 TSGetSNESIterations - Gets the total number of nonlinear iterations 4539 used by the time integrator. 4540 4541 Not Collective 4542 4543 Input Parameter: 4544 . ts - `TS` context 4545 4546 Output Parameter: 4547 . nits - number of nonlinear iterations 4548 4549 Level: intermediate 4550 4551 Note: 4552 This counter is reset to zero for each successive call to `TSSolve()`. 4553 4554 .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()` 4555 @*/ 4556 PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits) 4557 { 4558 PetscFunctionBegin; 4559 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4560 PetscValidIntPointer(nits, 2); 4561 *nits = ts->snes_its; 4562 PetscFunctionReturn(PETSC_SUCCESS); 4563 } 4564 4565 /*@ 4566 TSGetKSPIterations - Gets the total number of linear iterations 4567 used by the time integrator. 4568 4569 Not Collective 4570 4571 Input Parameter: 4572 . ts - `TS` context 4573 4574 Output Parameter: 4575 . lits - number of linear iterations 4576 4577 Level: intermediate 4578 4579 Note: 4580 This counter is reset to zero for each successive call to `TSSolve()`. 4581 4582 .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()` 4583 @*/ 4584 PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits) 4585 { 4586 PetscFunctionBegin; 4587 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4588 PetscValidIntPointer(lits, 2); 4589 *lits = ts->ksp_its; 4590 PetscFunctionReturn(PETSC_SUCCESS); 4591 } 4592 4593 /*@ 4594 TSGetStepRejections - Gets the total number of rejected steps. 4595 4596 Not Collective 4597 4598 Input Parameter: 4599 . ts - `TS` context 4600 4601 Output Parameter: 4602 . rejects - number of steps rejected 4603 4604 Level: intermediate 4605 4606 Note: 4607 This counter is reset to zero for each successive call to `TSSolve()`. 4608 4609 .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()` 4610 @*/ 4611 PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects) 4612 { 4613 PetscFunctionBegin; 4614 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4615 PetscValidIntPointer(rejects, 2); 4616 *rejects = ts->reject; 4617 PetscFunctionReturn(PETSC_SUCCESS); 4618 } 4619 4620 /*@ 4621 TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS` 4622 4623 Not Collective 4624 4625 Input Parameter: 4626 . ts - `TS` context 4627 4628 Output Parameter: 4629 . fails - number of failed nonlinear solves 4630 4631 Level: intermediate 4632 4633 Note: 4634 This counter is reset to zero for each successive call to `TSSolve()`. 4635 4636 .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()` 4637 @*/ 4638 PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails) 4639 { 4640 PetscFunctionBegin; 4641 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4642 PetscValidIntPointer(fails, 2); 4643 *fails = ts->num_snes_failures; 4644 PetscFunctionReturn(PETSC_SUCCESS); 4645 } 4646 4647 /*@ 4648 TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails 4649 4650 Not Collective 4651 4652 Input Parameters: 4653 + ts - `TS` context 4654 - rejects - maximum number of rejected steps, pass -1 for unlimited 4655 4656 Options Database Key: 4657 . -ts_max_reject - Maximum number of step rejections before a step fails 4658 4659 Level: intermediate 4660 4661 .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()` 4662 @*/ 4663 PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects) 4664 { 4665 PetscFunctionBegin; 4666 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4667 ts->max_reject = rejects; 4668 PetscFunctionReturn(PETSC_SUCCESS); 4669 } 4670 4671 /*@ 4672 TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves 4673 4674 Not Collective 4675 4676 Input Parameters: 4677 + ts - `TS` context 4678 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 4679 4680 Options Database Key: 4681 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 4682 4683 Level: intermediate 4684 4685 .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()` 4686 @*/ 4687 PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails) 4688 { 4689 PetscFunctionBegin; 4690 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4691 ts->max_snes_failures = fails; 4692 PetscFunctionReturn(PETSC_SUCCESS); 4693 } 4694 4695 /*@ 4696 TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()` 4697 4698 Not Collective 4699 4700 Input Parameters: 4701 + ts - `TS` context 4702 - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure 4703 4704 Options Database Key: 4705 . -ts_error_if_step_fails - Error if no step succeeds 4706 4707 Level: intermediate 4708 4709 .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()` 4710 @*/ 4711 PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err) 4712 { 4713 PetscFunctionBegin; 4714 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4715 ts->errorifstepfailed = err; 4716 PetscFunctionReturn(PETSC_SUCCESS); 4717 } 4718 4719 /*@ 4720 TSGetAdapt - Get the adaptive controller context for the current method 4721 4722 Collective on `ts` if controller has not been created yet 4723 4724 Input Parameter: 4725 . ts - time stepping context 4726 4727 Output Parameter: 4728 . adapt - adaptive controller 4729 4730 Level: intermediate 4731 4732 .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()` 4733 @*/ 4734 PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt) 4735 { 4736 PetscFunctionBegin; 4737 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4738 PetscValidPointer(adapt, 2); 4739 if (!ts->adapt) { 4740 PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt)); 4741 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1)); 4742 } 4743 *adapt = ts->adapt; 4744 PetscFunctionReturn(PETSC_SUCCESS); 4745 } 4746 4747 /*@ 4748 TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller 4749 4750 Logically Collective 4751 4752 Input Parameters: 4753 + ts - time integration context 4754 . atol - scalar absolute tolerances, `PETSC_DECIDE` to leave current value 4755 . vatol - vector of absolute tolerances or `NULL`, used in preference to atol if present 4756 . rtol - scalar relative tolerances, `PETSC_DECIDE` to leave current value 4757 - vrtol - vector of relative tolerances or `NULL`, used in preference to atol if present 4758 4759 Options Database Keys: 4760 + -ts_rtol <rtol> - relative tolerance for local truncation error 4761 - -ts_atol <atol> - Absolute tolerance for local truncation error 4762 4763 Level: beginner 4764 4765 Notes: 4766 With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error 4767 (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be 4768 computed only for the differential or the algebraic part then this can be done using the vector of 4769 tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the 4770 differential part and infinity for the algebraic part, the LTE calculation will include only the 4771 differential variables. 4772 4773 .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()` 4774 @*/ 4775 PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol) 4776 { 4777 PetscFunctionBegin; 4778 if (atol != (PetscReal)PETSC_DECIDE && atol != (PetscReal)PETSC_DEFAULT) ts->atol = atol; 4779 if (vatol) { 4780 PetscCall(PetscObjectReference((PetscObject)vatol)); 4781 PetscCall(VecDestroy(&ts->vatol)); 4782 ts->vatol = vatol; 4783 } 4784 if (rtol != (PetscReal)PETSC_DECIDE && rtol != (PetscReal)PETSC_DEFAULT) ts->rtol = rtol; 4785 if (vrtol) { 4786 PetscCall(PetscObjectReference((PetscObject)vrtol)); 4787 PetscCall(VecDestroy(&ts->vrtol)); 4788 ts->vrtol = vrtol; 4789 } 4790 PetscFunctionReturn(PETSC_SUCCESS); 4791 } 4792 4793 /*@ 4794 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 4795 4796 Logically Collective 4797 4798 Input Parameter: 4799 . ts - time integration context 4800 4801 Output Parameters: 4802 + atol - scalar absolute tolerances, `NULL` to ignore 4803 . vatol - vector of absolute tolerances, `NULL` to ignore 4804 . rtol - scalar relative tolerances, `NULL` to ignore 4805 - vrtol - vector of relative tolerances, `NULL` to ignore 4806 4807 Level: beginner 4808 4809 .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()` 4810 @*/ 4811 PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol) 4812 { 4813 PetscFunctionBegin; 4814 if (atol) *atol = ts->atol; 4815 if (vatol) *vatol = ts->vatol; 4816 if (rtol) *rtol = ts->rtol; 4817 if (vrtol) *vrtol = ts->vrtol; 4818 PetscFunctionReturn(PETSC_SUCCESS); 4819 } 4820 4821 /*@ 4822 TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances 4823 4824 Collective 4825 4826 Input Parameters: 4827 + ts - time stepping context 4828 . U - state vector, usually ts->vec_sol 4829 . Y - state vector to be compared to U 4830 - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY` 4831 4832 Output Parameters: 4833 + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances 4834 . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user 4835 - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user 4836 4837 Options Database Key: 4838 . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY 4839 4840 Level: developer 4841 4842 .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()` 4843 @*/ 4844 PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr) 4845 { 4846 PetscInt norma_loc, norm_loc, normr_loc; 4847 4848 PetscFunctionBegin; 4849 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4850 PetscCall(VecErrorWeightedNorms(U, Y, NULL, wnormtype, ts->atol, ts->vatol, ts->rtol, ts->vrtol, ts->adapt->ignore_max, norm, &norm_loc, norma, &norma_loc, normr, &normr_loc)); 4851 if (wnormtype == NORM_2) { 4852 if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc); 4853 if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc); 4854 if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc); 4855 } 4856 PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm"); 4857 PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma"); 4858 PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr"); 4859 PetscFunctionReturn(PETSC_SUCCESS); 4860 } 4861 4862 /*@ 4863 TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances 4864 4865 Collective 4866 4867 Input Parameters: 4868 + ts - time stepping context 4869 . E - error vector 4870 . U - state vector, usually ts->vec_sol 4871 . Y - state vector, previous time step 4872 - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY` 4873 4874 Output Parameters: 4875 + norm - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances 4876 . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user 4877 - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user 4878 4879 Options Database Key: 4880 . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY 4881 4882 Level: developer 4883 4884 .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()` 4885 @*/ 4886 PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr) 4887 { 4888 PetscInt norma_loc, norm_loc, normr_loc; 4889 4890 PetscFunctionBegin; 4891 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4892 PetscCall(VecErrorWeightedNorms(U, Y, E, wnormtype, ts->atol, ts->vatol, ts->rtol, ts->vrtol, ts->adapt->ignore_max, norm, &norm_loc, norma, &norma_loc, normr, &normr_loc)); 4893 if (wnormtype == NORM_2) { 4894 if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc); 4895 if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc); 4896 if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc); 4897 } 4898 PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm"); 4899 PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma"); 4900 PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr"); 4901 PetscFunctionReturn(PETSC_SUCCESS); 4902 } 4903 4904 /*@ 4905 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 4906 4907 Logically Collective 4908 4909 Input Parameters: 4910 + ts - time stepping context 4911 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 4912 4913 Note: 4914 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 4915 4916 Level: intermediate 4917 4918 .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL` 4919 @*/ 4920 PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime) 4921 { 4922 PetscFunctionBegin; 4923 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 4924 ts->cfltime_local = cfltime; 4925 ts->cfltime = -1.; 4926 PetscFunctionReturn(PETSC_SUCCESS); 4927 } 4928 4929 /*@ 4930 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 4931 4932 Collective 4933 4934 Input Parameter: 4935 . ts - time stepping context 4936 4937 Output Parameter: 4938 . cfltime - maximum stable time step for forward Euler 4939 4940 Level: advanced 4941 4942 .seealso: [](ch_ts), `TSSetCFLTimeLocal()` 4943 @*/ 4944 PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime) 4945 { 4946 PetscFunctionBegin; 4947 if (ts->cfltime < 0) PetscCall(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts))); 4948 *cfltime = ts->cfltime; 4949 PetscFunctionReturn(PETSC_SUCCESS); 4950 } 4951 4952 /*@ 4953 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 4954 4955 Input Parameters: 4956 + ts - the `TS` context. 4957 . xl - lower bound. 4958 - xu - upper bound. 4959 4960 Level: advanced 4961 4962 Note: 4963 If this routine is not called then the lower and upper bounds are set to 4964 `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`. 4965 4966 .seealso: [](ch_ts), `TS` 4967 @*/ 4968 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 4969 { 4970 SNES snes; 4971 4972 PetscFunctionBegin; 4973 PetscCall(TSGetSNES(ts, &snes)); 4974 PetscCall(SNESVISetVariableBounds(snes, xl, xu)); 4975 PetscFunctionReturn(PETSC_SUCCESS); 4976 } 4977 4978 /*@ 4979 TSComputeLinearStability - computes the linear stability function at a point 4980 4981 Collective 4982 4983 Input Parameters: 4984 + ts - the `TS` context 4985 . xr - real part of input argument 4986 - xi - imaginary part of input argument 4987 4988 Output Parameters: 4989 + yr - real part of function value 4990 - yi - imaginary part of function value 4991 4992 Level: developer 4993 4994 .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()` 4995 @*/ 4996 PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi) 4997 { 4998 PetscFunctionBegin; 4999 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5000 PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi); 5001 PetscFunctionReturn(PETSC_SUCCESS); 5002 } 5003 5004 /*@ 5005 TSRestartStep - Flags the solver to restart the next step 5006 5007 Collective 5008 5009 Input Parameter: 5010 . ts - the `TS` context obtained from `TSCreate()` 5011 5012 Level: advanced 5013 5014 Notes: 5015 Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of 5016 discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution 5017 vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For 5018 the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce 5019 discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with 5020 discontinuous source terms). 5021 5022 .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()` 5023 @*/ 5024 PetscErrorCode TSRestartStep(TS ts) 5025 { 5026 PetscFunctionBegin; 5027 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5028 ts->steprestart = PETSC_TRUE; 5029 PetscFunctionReturn(PETSC_SUCCESS); 5030 } 5031 5032 /*@ 5033 TSRollBack - Rolls back one time step 5034 5035 Collective 5036 5037 Input Parameter: 5038 . ts - the `TS` context obtained from `TSCreate()` 5039 5040 Level: advanced 5041 5042 .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()` 5043 @*/ 5044 PetscErrorCode TSRollBack(TS ts) 5045 { 5046 PetscFunctionBegin; 5047 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5048 PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called"); 5049 PetscUseTypeMethod(ts, rollback); 5050 ts->time_step = ts->ptime - ts->ptime_prev; 5051 ts->ptime = ts->ptime_prev; 5052 ts->ptime_prev = ts->ptime_prev_rollback; 5053 ts->steps--; 5054 ts->steprollback = PETSC_TRUE; 5055 PetscFunctionReturn(PETSC_SUCCESS); 5056 } 5057 5058 /*@ 5059 TSGetStages - Get the number of stages and stage values 5060 5061 Input Parameter: 5062 . ts - the `TS` context obtained from `TSCreate()` 5063 5064 Output Parameters: 5065 + ns - the number of stages 5066 - Y - the current stage vectors 5067 5068 Level: advanced 5069 5070 Note: 5071 Both `ns` and `Y` can be `NULL`. 5072 5073 .seealso: [](ch_ts), `TS`, `TSCreate()` 5074 @*/ 5075 PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y) 5076 { 5077 PetscFunctionBegin; 5078 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5079 if (ns) PetscValidIntPointer(ns, 2); 5080 if (Y) PetscValidPointer(Y, 3); 5081 if (!ts->ops->getstages) { 5082 if (ns) *ns = 0; 5083 if (Y) *Y = NULL; 5084 } else PetscUseTypeMethod(ts, getstages, ns, Y); 5085 PetscFunctionReturn(PETSC_SUCCESS); 5086 } 5087 5088 /*@C 5089 TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity. 5090 5091 Collective 5092 5093 Input Parameters: 5094 + ts - the `TS` context 5095 . t - current timestep 5096 . U - state vector 5097 . Udot - time derivative of state vector 5098 . shift - shift to apply, see note below 5099 - ctx - an optional user context 5100 5101 Output Parameters: 5102 + J - Jacobian matrix (not altered in this routine) 5103 - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`) 5104 5105 Level: intermediate 5106 5107 Notes: 5108 If F(t,U,Udot)=0 is the DAE, the required Jacobian is 5109 5110 dF/dU + shift*dF/dUdot 5111 5112 Most users should not need to explicitly call this routine, as it 5113 is used internally within the nonlinear solvers. 5114 5115 This will first try to get the coloring from the `DM`. If the `DM` type has no coloring 5116 routine, then it will try to get the coloring from the matrix. This requires that the 5117 matrix have nonzero entries precomputed. 5118 5119 .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()` 5120 @*/ 5121 PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx) 5122 { 5123 SNES snes; 5124 MatFDColoring color; 5125 PetscBool hascolor, matcolor = PETSC_FALSE; 5126 5127 PetscFunctionBegin; 5128 PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL)); 5129 PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color)); 5130 if (!color) { 5131 DM dm; 5132 ISColoring iscoloring; 5133 5134 PetscCall(TSGetDM(ts, &dm)); 5135 PetscCall(DMHasColoring(dm, &hascolor)); 5136 if (hascolor && !matcolor) { 5137 PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring)); 5138 PetscCall(MatFDColoringCreate(B, iscoloring, &color)); 5139 PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts)); 5140 PetscCall(MatFDColoringSetFromOptions(color)); 5141 PetscCall(MatFDColoringSetUp(B, iscoloring, color)); 5142 PetscCall(ISColoringDestroy(&iscoloring)); 5143 } else { 5144 MatColoring mc; 5145 5146 PetscCall(MatColoringCreate(B, &mc)); 5147 PetscCall(MatColoringSetDistance(mc, 2)); 5148 PetscCall(MatColoringSetType(mc, MATCOLORINGSL)); 5149 PetscCall(MatColoringSetFromOptions(mc)); 5150 PetscCall(MatColoringApply(mc, &iscoloring)); 5151 PetscCall(MatColoringDestroy(&mc)); 5152 PetscCall(MatFDColoringCreate(B, iscoloring, &color)); 5153 PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts)); 5154 PetscCall(MatFDColoringSetFromOptions(color)); 5155 PetscCall(MatFDColoringSetUp(B, iscoloring, color)); 5156 PetscCall(ISColoringDestroy(&iscoloring)); 5157 } 5158 PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color)); 5159 PetscCall(PetscObjectDereference((PetscObject)color)); 5160 } 5161 PetscCall(TSGetSNES(ts, &snes)); 5162 PetscCall(MatFDColoringApply(B, color, U, snes)); 5163 if (J != B) { 5164 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 5165 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 5166 } 5167 PetscFunctionReturn(PETSC_SUCCESS); 5168 } 5169 5170 /*@ 5171 TSSetFunctionDomainError - Set a function that tests if the current state vector is valid 5172 5173 Input Parameters: 5174 + ts - the `TS` context 5175 - func - function called within `TSFunctionDomainError()` 5176 5177 Calling sequence of `func`: 5178 $ PetscErrorCode func(TS ts, PetscReal time, Vec state, PetscBool reject) 5179 + ts - the TS context 5180 . time - the current time (of the stage) 5181 . state - the state to check if it is valid 5182 - reject - (output parameter) `PETSC_FALSE` if the state is acceptable, `PETSC_TRUE` if not acceptable 5183 5184 Level: intermediate 5185 5186 Notes: 5187 If an implicit ODE solver is being used then, in addition to providing this routine, the 5188 user's code should call `SNESSetFunctionDomainError()` when domain errors occur during 5189 function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`. 5190 Use `TSGetSNES()` to obtain the `SNES` object 5191 5192 Developer Note: 5193 The naming of this function is inconsistent with the `SNESSetFunctionDomainError()` 5194 since one takes a function pointer and the other does not. 5195 5196 .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()` 5197 @*/ 5198 5199 PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS, PetscReal, Vec, PetscBool *)) 5200 { 5201 PetscFunctionBegin; 5202 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5203 ts->functiondomainerror = func; 5204 PetscFunctionReturn(PETSC_SUCCESS); 5205 } 5206 5207 /*@ 5208 TSFunctionDomainError - Checks if the current state is valid 5209 5210 Input Parameters: 5211 + ts - the `TS` context 5212 . stagetime - time of the simulation 5213 - Y - state vector to check. 5214 5215 Output Parameter: 5216 . accept - Set to `PETSC_FALSE` if the current state vector is valid. 5217 5218 Level: developer 5219 5220 Note: 5221 This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`) 5222 to check if the current state is valid. 5223 5224 .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()` 5225 @*/ 5226 PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept) 5227 { 5228 PetscFunctionBegin; 5229 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5230 *accept = PETSC_TRUE; 5231 if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept)); 5232 PetscFunctionReturn(PETSC_SUCCESS); 5233 } 5234 5235 /*@C 5236 TSClone - This function clones a time step `TS` object. 5237 5238 Collective 5239 5240 Input Parameter: 5241 . tsin - The input `TS` 5242 5243 Output Parameter: 5244 . tsout - The output `TS` (cloned) 5245 5246 Level: developer 5247 5248 Notes: 5249 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. 5250 It will likely be replaced in the future with a mechanism of switching methods on the fly. 5251 5252 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 5253 .vb 5254 SNES snes_dup = NULL; 5255 TSGetSNES(ts,&snes_dup); 5256 TSSetSNES(ts,snes_dup); 5257 .ve 5258 5259 .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()` 5260 @*/ 5261 PetscErrorCode TSClone(TS tsin, TS *tsout) 5262 { 5263 TS t; 5264 SNES snes_start; 5265 DM dm; 5266 TSType type; 5267 5268 PetscFunctionBegin; 5269 PetscValidPointer(tsin, 1); 5270 *tsout = NULL; 5271 5272 PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView)); 5273 5274 /* General TS description */ 5275 t->numbermonitors = 0; 5276 t->monitorFrequency = 1; 5277 t->setupcalled = 0; 5278 t->ksp_its = 0; 5279 t->snes_its = 0; 5280 t->nwork = 0; 5281 t->rhsjacobian.time = PETSC_MIN_REAL; 5282 t->rhsjacobian.scale = 1.; 5283 t->ijacobian.shift = 1.; 5284 5285 PetscCall(TSGetSNES(tsin, &snes_start)); 5286 PetscCall(TSSetSNES(t, snes_start)); 5287 5288 PetscCall(TSGetDM(tsin, &dm)); 5289 PetscCall(TSSetDM(t, dm)); 5290 5291 t->adapt = tsin->adapt; 5292 PetscCall(PetscObjectReference((PetscObject)t->adapt)); 5293 5294 t->trajectory = tsin->trajectory; 5295 PetscCall(PetscObjectReference((PetscObject)t->trajectory)); 5296 5297 t->event = tsin->event; 5298 if (t->event) t->event->refct++; 5299 5300 t->problem_type = tsin->problem_type; 5301 t->ptime = tsin->ptime; 5302 t->ptime_prev = tsin->ptime_prev; 5303 t->time_step = tsin->time_step; 5304 t->max_time = tsin->max_time; 5305 t->steps = tsin->steps; 5306 t->max_steps = tsin->max_steps; 5307 t->equation_type = tsin->equation_type; 5308 t->atol = tsin->atol; 5309 t->rtol = tsin->rtol; 5310 t->max_snes_failures = tsin->max_snes_failures; 5311 t->max_reject = tsin->max_reject; 5312 t->errorifstepfailed = tsin->errorifstepfailed; 5313 5314 PetscCall(TSGetType(tsin, &type)); 5315 PetscCall(TSSetType(t, type)); 5316 5317 t->vec_sol = NULL; 5318 5319 t->cfltime = tsin->cfltime; 5320 t->cfltime_local = tsin->cfltime_local; 5321 t->exact_final_time = tsin->exact_final_time; 5322 5323 t->ops[0] = tsin->ops[0]; 5324 5325 if (((PetscObject)tsin)->fortran_func_pointers) { 5326 PetscInt i; 5327 PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers)); 5328 for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i]; 5329 } 5330 *tsout = t; 5331 PetscFunctionReturn(PETSC_SUCCESS); 5332 } 5333 5334 static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y) 5335 { 5336 TS ts = (TS)ctx; 5337 5338 PetscFunctionBegin; 5339 PetscCall(TSComputeRHSFunction(ts, 0, x, y)); 5340 PetscFunctionReturn(PETSC_SUCCESS); 5341 } 5342 5343 /*@ 5344 TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function. 5345 5346 Logically Collective 5347 5348 Input Parameter: 5349 . TS - the time stepping routine 5350 5351 Output Parameter: 5352 . flg - `PETSC_TRUE` if the multiply is likely correct 5353 5354 Options Database Key: 5355 . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator 5356 5357 Level: advanced 5358 5359 Note: 5360 This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian 5361 5362 .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()` 5363 @*/ 5364 PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg) 5365 { 5366 Mat J, B; 5367 TSRHSJacobian func; 5368 void *ctx; 5369 5370 PetscFunctionBegin; 5371 PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx)); 5372 PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx)); 5373 PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg)); 5374 PetscFunctionReturn(PETSC_SUCCESS); 5375 } 5376 5377 /*@C 5378 TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function. 5379 5380 Logically Collective 5381 5382 Input Parameter: 5383 . TS - the time stepping routine 5384 5385 Output Parameter: 5386 . flg - `PETSC_TRUE` if the multiply is likely correct 5387 5388 Options Database Key: 5389 . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator 5390 5391 Level: advanced 5392 5393 Notes: 5394 This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian 5395 5396 .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()` 5397 @*/ 5398 PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg) 5399 { 5400 Mat J, B; 5401 void *ctx; 5402 TSRHSJacobian func; 5403 5404 PetscFunctionBegin; 5405 PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx)); 5406 PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx)); 5407 PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg)); 5408 PetscFunctionReturn(PETSC_SUCCESS); 5409 } 5410 5411 /*@ 5412 TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used. 5413 5414 Logically Collective 5415 5416 Input Parameters: 5417 + ts - timestepping context 5418 - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used 5419 5420 Options Database Key: 5421 . -ts_use_splitrhsfunction - <true,false> 5422 5423 Level: intermediate 5424 5425 Note: 5426 This is only for multirate methods 5427 5428 .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()` 5429 @*/ 5430 PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction) 5431 { 5432 PetscFunctionBegin; 5433 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5434 ts->use_splitrhsfunction = use_splitrhsfunction; 5435 PetscFunctionReturn(PETSC_SUCCESS); 5436 } 5437 5438 /*@ 5439 TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used. 5440 5441 Not Collective 5442 5443 Input Parameter: 5444 . ts - timestepping context 5445 5446 Output Parameter: 5447 . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used 5448 5449 Level: intermediate 5450 5451 .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()` 5452 @*/ 5453 PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction) 5454 { 5455 PetscFunctionBegin; 5456 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5457 *use_splitrhsfunction = ts->use_splitrhsfunction; 5458 PetscFunctionReturn(PETSC_SUCCESS); 5459 } 5460 5461 /*@ 5462 TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix. 5463 5464 Logically Collective 5465 5466 Input Parameters: 5467 + ts - the time-stepper 5468 - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`) 5469 5470 Level: intermediate 5471 5472 Note: 5473 When the relationship between the nonzero structures is known and supplied the solution process can be much faster 5474 5475 .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure` 5476 @*/ 5477 PetscErrorCode TSSetMatStructure(TS ts, MatStructure str) 5478 { 5479 PetscFunctionBegin; 5480 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5481 ts->axpy_pattern = str; 5482 PetscFunctionReturn(PETSC_SUCCESS); 5483 } 5484 5485 /*@ 5486 TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span 5487 5488 Collective 5489 5490 Input Parameters: 5491 + ts - the time-stepper 5492 . n - number of the time points (>=2) 5493 - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively. 5494 5495 Options Database Key: 5496 . -ts_time_span <t0,...tf> - Sets the time span 5497 5498 Level: intermediate 5499 5500 Notes: 5501 The elements in tspan must be all increasing. They correspond to the intermediate points for time integration. 5502 `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified. 5503 The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may 5504 pressure the memory system when using a large number of span points. 5505 5506 .seealso: [](ch_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()` 5507 @*/ 5508 PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times) 5509 { 5510 PetscFunctionBegin; 5511 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5512 PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n); 5513 if (ts->tspan && n != ts->tspan->num_span_times) { 5514 PetscCall(PetscFree(ts->tspan->span_times)); 5515 PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol)); 5516 PetscCall(PetscMalloc1(n, &ts->tspan->span_times)); 5517 } 5518 if (!ts->tspan) { 5519 TSTimeSpan tspan; 5520 PetscCall(PetscNew(&tspan)); 5521 PetscCall(PetscMalloc1(n, &tspan->span_times)); 5522 tspan->reltol = 1e-6; 5523 tspan->abstol = 10 * PETSC_MACHINE_EPSILON; 5524 ts->tspan = tspan; 5525 } 5526 ts->tspan->num_span_times = n; 5527 PetscCall(PetscArraycpy(ts->tspan->span_times, span_times, n)); 5528 PetscCall(TSSetTime(ts, ts->tspan->span_times[0])); 5529 PetscCall(TSSetMaxTime(ts, ts->tspan->span_times[n - 1])); 5530 PetscFunctionReturn(PETSC_SUCCESS); 5531 } 5532 5533 /*@C 5534 TSGetTimeSpan - gets the time span set with `TSSetTimeSpan()` 5535 5536 Not Collective 5537 5538 Input Parameter: 5539 . ts - the time-stepper 5540 5541 Output Parameters: 5542 + n - number of the time points (>=2) 5543 - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively. 5544 5545 Level: beginner 5546 5547 Note: 5548 The values obtained are valid until the `TS` object is destroyed. 5549 5550 Both `n` and `span_times` can be `NULL`. 5551 5552 .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()` 5553 @*/ 5554 PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal **span_times) 5555 { 5556 PetscFunctionBegin; 5557 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5558 if (n) PetscValidIntPointer(n, 2); 5559 if (span_times) PetscValidPointer(span_times, 3); 5560 if (!ts->tspan) { 5561 if (n) *n = 0; 5562 if (span_times) *span_times = NULL; 5563 } else { 5564 if (n) *n = ts->tspan->num_span_times; 5565 if (span_times) *span_times = ts->tspan->span_times; 5566 } 5567 PetscFunctionReturn(PETSC_SUCCESS); 5568 } 5569 5570 /*@ 5571 TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span. 5572 5573 Input Parameter: 5574 . ts - the `TS` context obtained from `TSCreate()` 5575 5576 Output Parameters: 5577 + nsol - the number of solutions 5578 - Sols - the solution vectors 5579 5580 Level: intermediate 5581 5582 Notes: 5583 Both `nsol` and `Sols` can be `NULL`. 5584 5585 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()`. 5586 For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span. 5587 5588 .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()` 5589 @*/ 5590 PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols) 5591 { 5592 PetscFunctionBegin; 5593 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 5594 if (nsol) PetscValidIntPointer(nsol, 2); 5595 if (Sols) PetscValidPointer(Sols, 3); 5596 if (!ts->tspan) { 5597 if (nsol) *nsol = 0; 5598 if (Sols) *Sols = NULL; 5599 } else { 5600 if (nsol) *nsol = ts->tspan->spanctr; 5601 if (Sols) *Sols = ts->tspan->vecs_sol; 5602 } 5603 PetscFunctionReturn(PETSC_SUCCESS); 5604 } 5605 5606 /*@C 5607 TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information. 5608 5609 Collective 5610 5611 Input Parameters: 5612 + ts - the `TS` context 5613 . J - Jacobian matrix (not altered in this routine) 5614 - B - newly computed Jacobian matrix to use with preconditioner 5615 5616 Level: intermediate 5617 5618 Notes: 5619 This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains 5620 many constant zeros entries, which is typically the case when the matrix is generated by a `DM` 5621 and multiple fields are involved. 5622 5623 Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity 5624 structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can 5625 usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian. 5626 `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`. 5627 5628 .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()` 5629 @*/ 5630 PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B) 5631 { 5632 MatColoring mc = NULL; 5633 ISColoring iscoloring = NULL; 5634 MatFDColoring matfdcoloring = NULL; 5635 5636 PetscFunctionBegin; 5637 /* Generate new coloring after eliminating zeros in the matrix */ 5638 PetscCall(MatEliminateZeros(B)); 5639 PetscCall(MatColoringCreate(B, &mc)); 5640 PetscCall(MatColoringSetDistance(mc, 2)); 5641 PetscCall(MatColoringSetType(mc, MATCOLORINGSL)); 5642 PetscCall(MatColoringSetFromOptions(mc)); 5643 PetscCall(MatColoringApply(mc, &iscoloring)); 5644 PetscCall(MatColoringDestroy(&mc)); 5645 /* Replace the old coloring with the new one */ 5646 PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring)); 5647 PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts)); 5648 PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 5649 PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring)); 5650 PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring)); 5651 PetscCall(PetscObjectDereference((PetscObject)matfdcoloring)); 5652 PetscCall(ISColoringDestroy(&iscoloring)); 5653 PetscFunctionReturn(PETSC_SUCCESS); 5654 } 5655