Actual source code: ts.c

      1: #include <petsc/private/tsimpl.h>
      2: #include <petscdmda.h>
      3: #include <petscdmshell.h>
      4: #include <petscdmplex.h>
      5: #include <petscdmswarm.h>
      6: #include <petscviewer.h>
      7: #include <petscdraw.h>
      8: #include <petscconvest.h>

     10: /* Logging support */
     11: PetscClassId  TS_CLASSID, DMTS_CLASSID;
     12: PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;

     14: const char *const TSExactFinalTimeOptions[] = {"UNSPECIFIED", "STEPOVER", "INTERPOLATE", "MATCHSTEP", "TSExactFinalTimeOption", "TS_EXACTFINALTIME_", NULL};

     16: static PetscErrorCode TSAdaptSetDefaultType(TSAdapt adapt, TSAdaptType default_type)
     17: {
     18:   PetscFunctionBegin;
     20:   PetscAssertPointer(default_type, 2);
     21:   if (!((PetscObject)adapt)->type_name) PetscCall(TSAdaptSetType(adapt, default_type));
     22:   PetscFunctionReturn(PETSC_SUCCESS);
     23: }

     25: /*@
     26:   TSSetFromOptions - Sets various `TS` parameters from the options database

     28:   Collective

     30:   Input Parameter:
     31: . ts - the `TS` context obtained from `TSCreate()`

     33:   Options Database Keys:
     34: + -ts_type <type>                                                    - EULER, BEULER, SUNDIALS, PSEUDO, CN, RK, THETA, ALPHA, GLLE,  SSP, GLEE, BSYMP, IRK, see `TSType`
     35: . -ts_save_trajectory                                                - checkpoint the solution at each time-step
     36: . -ts_max_time <time>                                                - maximum time to compute to
     37: . -ts_time_span <t0,...tf>                                           - sets the time span, solutions are computed and stored for each indicated time
     38: . -ts_max_steps <steps>                                              - maximum number of time-steps to take
     39: . -ts_init_time <time>                                               - initial time to start computation
     40: . -ts_final_time <time>                                              - final time to compute to (deprecated: use `-ts_max_time`)
     41: . -ts_dt <dt>                                                        - initial time step
     42: . -ts_exact_final_time <stepover,interpolate,matchstep>              - whether to stop at the exact given final time and how to compute the solution at that time
     43: . -ts_max_snes_failures <maxfailures>                                - Maximum number of nonlinear solve failures allowed
     44: . -ts_max_reject <maxrejects>                                        - Maximum number of step rejections before step fails
     45: . -ts_error_if_step_fails <true,false>                               - Error if no step succeeds
     46: . -ts_rtol <rtol>                                                    - relative tolerance for local truncation error
     47: . -ts_atol <atol>                                                    - Absolute tolerance for local truncation error
     48: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view               - test the Jacobian at each iteration against finite difference with RHS function
     49: . -ts_rhs_jacobian_test_mult_transpose                               - test the Jacobian at each iteration against finite difference with RHS function
     50: . -ts_adjoint_solve <yes,no>                                         - After solving the ODE/DAE solve the adjoint problem (requires `-ts_save_trajectory`)
     51: . -ts_fd_color                                                       - Use finite differences with coloring to compute IJacobian
     52: . -ts_monitor                                                        - print information at each timestep
     53: . -ts_monitor_cancel                                                 - Cancel all monitors
     54: . -ts_monitor_lg_solution                                            - Monitor solution graphically
     55: . -ts_monitor_lg_error                                               - Monitor error graphically
     56: . -ts_monitor_error                                                  - Monitors norm of error
     57: . -ts_monitor_lg_timestep                                            - Monitor timestep size graphically
     58: . -ts_monitor_lg_timestep_log                                        - Monitor log timestep size graphically
     59: . -ts_monitor_lg_snes_iterations                                     - Monitor number nonlinear iterations for each timestep graphically
     60: . -ts_monitor_lg_ksp_iterations                                      - Monitor number nonlinear iterations for each timestep graphically
     61: . -ts_monitor_sp_eig                                                 - Monitor eigenvalues of linearized operator graphically
     62: . -ts_monitor_draw_solution                                          - Monitor solution graphically
     63: . -ts_monitor_draw_solution_phase  <xleft,yleft,xright,yright>       - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom
     64: . -ts_monitor_draw_error                                             - Monitor error graphically, requires use to have provided TSSetSolutionFunction()
     65: . -ts_monitor_solution [ascii binary draw][:filename][:viewerformat] - monitors the solution at each timestep
     66: . -ts_monitor_solution_interval <interval>                           - output once every interval (default=1) time steps. Use -1 to only output at the end of the simulation
     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_solution_vtk_interval <interval>                       - output once every interval (default=1) time steps. Use -1 to only output at the end of the simulation
     69: - -ts_monitor_envelope                                               - determine maximum and minimum value of each component of the solution over the solution time

     71:   Level: beginner

     73:   Notes:
     74:   See `SNESSetFromOptions()` and `KSPSetFromOptions()` for how to control the nonlinear and linear solves used by the time-stepper.

     76:   Certain `SNES` options get reset for each new nonlinear solver, for example `-snes_lag_jacobian its` and `-snes_lag_preconditioner its`, in order
     77:   to retain them over the multiple nonlinear solves that `TS` uses you must also provide `-snes_lag_jacobian_persists true` and
     78:   `-snes_lag_preconditioner_persists true`

     80:   Developer Notes:
     81:   We should unify all the -ts_monitor options in the way that -xxx_view has been unified

     83: .seealso: [](ch_ts), `TS`, `TSGetType()`
     84: @*/
     85: PetscErrorCode TSSetFromOptions(TS ts)
     86: {
     87:   PetscBool              opt, flg, tflg;
     88:   char                   monfilename[PETSC_MAX_PATH_LEN];
     89:   PetscReal              time_step, tspan[100];
     90:   PetscInt               nt = PETSC_STATIC_ARRAY_LENGTH(tspan);
     91:   TSExactFinalTimeOption eftopt;
     92:   char                   dir[16];
     93:   TSIFunctionFn         *ifun;
     94:   const char            *defaultType;
     95:   char                   typeName[256];

     97:   PetscFunctionBegin;

    100:   PetscCall(TSRegisterAll());
    101:   PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));

    103:   PetscObjectOptionsBegin((PetscObject)ts);
    104:   if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name;
    105:   else defaultType = ifun ? TSBEULER : TSEULER;
    106:   PetscCall(PetscOptionsFList("-ts_type", "TS method", "TSSetType", TSList, defaultType, typeName, 256, &opt));
    107:   if (opt) PetscCall(TSSetType(ts, typeName));
    108:   else PetscCall(TSSetType(ts, defaultType));

    110:   /* Handle generic TS options */
    111:   PetscCall(PetscOptionsDeprecated("-ts_final_time", "-ts_max_time", "3.10", NULL));
    112:   PetscCall(PetscOptionsReal("-ts_max_time", "Maximum time to run to", "TSSetMaxTime", ts->max_time, &ts->max_time, NULL));
    113:   PetscCall(PetscOptionsRealArray("-ts_time_span", "Time span", "TSSetTimeSpan", tspan, &nt, &flg));
    114:   if (flg) PetscCall(TSSetTimeSpan(ts, nt, tspan));
    115:   PetscCall(PetscOptionsInt("-ts_max_steps", "Maximum number of time steps", "TSSetMaxSteps", ts->max_steps, &ts->max_steps, NULL));
    116:   PetscCall(PetscOptionsReal("-ts_init_time", "Initial time", "TSSetTime", ts->ptime, &ts->ptime, NULL));
    117:   PetscCall(PetscOptionsReal("-ts_dt", "Initial time step", "TSSetTimeStep", ts->time_step, &time_step, &flg));
    118:   if (flg) PetscCall(TSSetTimeStep(ts, time_step));
    119:   PetscCall(PetscOptionsEnum("-ts_exact_final_time", "Option for handling of final time step", "TSSetExactFinalTime", TSExactFinalTimeOptions, (PetscEnum)ts->exact_final_time, (PetscEnum *)&eftopt, &flg));
    120:   if (flg) PetscCall(TSSetExactFinalTime(ts, eftopt));
    121:   PetscCall(PetscOptionsInt("-ts_max_snes_failures", "Maximum number of nonlinear solve failures", "TSSetMaxSNESFailures", ts->max_snes_failures, &ts->max_snes_failures, &flg));
    122:   if (flg) PetscCall(TSSetMaxSNESFailures(ts, ts->max_snes_failures));
    123:   PetscCall(PetscOptionsInt("-ts_max_reject", "Maximum number of step rejections before step fails", "TSSetMaxStepRejections", ts->max_reject, &ts->max_reject, &flg));
    124:   if (flg) PetscCall(TSSetMaxStepRejections(ts, ts->max_reject));
    125:   PetscCall(PetscOptionsBool("-ts_error_if_step_fails", "Error if no step succeeds", "TSSetErrorIfStepFails", ts->errorifstepfailed, &ts->errorifstepfailed, NULL));
    126:   PetscCall(PetscOptionsBoundedReal("-ts_rtol", "Relative tolerance for local truncation error", "TSSetTolerances", ts->rtol, &ts->rtol, NULL, 0));
    127:   PetscCall(PetscOptionsBoundedReal("-ts_atol", "Absolute tolerance for local truncation error", "TSSetTolerances", ts->atol, &ts->atol, NULL, 0));

    129:   PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult", "Test the RHS Jacobian for consistency with RHS at each solve ", "None", ts->testjacobian, &ts->testjacobian, NULL));
    130:   PetscCall(PetscOptionsBool("-ts_rhs_jacobian_test_mult_transpose", "Test the RHS Jacobian transpose for consistency with RHS at each solve ", "None", ts->testjacobiantranspose, &ts->testjacobiantranspose, NULL));
    131:   PetscCall(PetscOptionsBool("-ts_use_splitrhsfunction", "Use the split RHS function for multirate solvers ", "TSSetUseSplitRHSFunction", ts->use_splitrhsfunction, &ts->use_splitrhsfunction, NULL));
    132: #if defined(PETSC_HAVE_SAWS)
    133:   {
    134:     PetscBool set;
    135:     flg = PETSC_FALSE;
    136:     PetscCall(PetscOptionsBool("-ts_saws_block", "Block for SAWs memory snooper at end of TSSolve", "PetscObjectSAWsBlock", ((PetscObject)ts)->amspublishblock, &flg, &set));
    137:     if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)ts, flg));
    138:   }
    139: #endif

    141:   /* Monitor options */
    142:   PetscCall(PetscOptionsInt("-ts_monitor_frequency", "Number of time steps between monitor output", "TSMonitorSetFrequency", ts->monitorFrequency, &ts->monitorFrequency, NULL));
    143:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor", "Monitor time and timestep size", "TSMonitorDefault", TSMonitorDefault, NULL));
    144:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_extreme", "Monitor extreme values of the solution", "TSMonitorExtreme", TSMonitorExtreme, NULL));
    145:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_solution", "View the solution at each timestep", "TSMonitorSolution", TSMonitorSolution, NULL));
    146:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_dmswarm_monitor_moments", "Monitor moments of particle distribution", "TSDMSwarmMonitorMoments", TSDMSwarmMonitorMoments, NULL));

    148:   PetscCall(PetscOptionsString("-ts_monitor_python", "Use Python function", "TSMonitorSet", NULL, monfilename, sizeof(monfilename), &flg));
    149:   if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)ts, monfilename));

    151:   PetscCall(PetscOptionsName("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", &opt));
    152:   if (opt) {
    153:     PetscInt  howoften = 1;
    154:     DM        dm;
    155:     PetscBool net;

    157:     PetscCall(PetscOptionsInt("-ts_monitor_lg_solution", "Monitor solution graphically", "TSMonitorLGSolution", howoften, &howoften, NULL));
    158:     PetscCall(TSGetDM(ts, &dm));
    159:     PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMNETWORK, &net));
    160:     if (net) {
    161:       TSMonitorLGCtxNetwork ctx;
    162:       PetscCall(TSMonitorLGCtxNetworkCreate(ts, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &ctx));
    163:       PetscCall(TSMonitorSet(ts, TSMonitorLGCtxNetworkSolution, ctx, (PetscErrorCode (*)(void **))TSMonitorLGCtxNetworkDestroy));
    164:       PetscCall(PetscOptionsBool("-ts_monitor_lg_solution_semilogy", "Plot the solution with a semi-log axis", "", ctx->semilogy, &ctx->semilogy, NULL));
    165:     } else {
    166:       TSMonitorLGCtx ctx;
    167:       PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
    168:       PetscCall(TSMonitorSet(ts, TSMonitorLGSolution, ctx, (PetscErrorCode (*)(void **))TSMonitorLGCtxDestroy));
    169:     }
    170:   }

    172:   PetscCall(PetscOptionsName("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", &opt));
    173:   if (opt) {
    174:     TSMonitorLGCtx ctx;
    175:     PetscInt       howoften = 1;

    177:     PetscCall(PetscOptionsInt("-ts_monitor_lg_error", "Monitor error graphically", "TSMonitorLGError", howoften, &howoften, NULL));
    178:     PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
    179:     PetscCall(TSMonitorSet(ts, TSMonitorLGError, ctx, (PetscErrorCode (*)(void **))TSMonitorLGCtxDestroy));
    180:   }
    181:   PetscCall(TSMonitorSetFromOptions(ts, "-ts_monitor_error", "View the error at each timestep", "TSMonitorError", TSMonitorError, NULL));

    183:   PetscCall(PetscOptionsName("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", &opt));
    184:   if (opt) {
    185:     TSMonitorLGCtx ctx;
    186:     PetscInt       howoften = 1;

    188:     PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep", "Monitor timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
    189:     PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
    190:     PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode (*)(void **))TSMonitorLGCtxDestroy));
    191:   }
    192:   PetscCall(PetscOptionsName("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", &opt));
    193:   if (opt) {
    194:     TSMonitorLGCtx ctx;
    195:     PetscInt       howoften = 1;

    197:     PetscCall(PetscOptionsInt("-ts_monitor_lg_timestep_log", "Monitor log timestep size graphically", "TSMonitorLGTimeStep", howoften, &howoften, NULL));
    198:     PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
    199:     PetscCall(TSMonitorSet(ts, TSMonitorLGTimeStep, ctx, (PetscErrorCode (*)(void **))TSMonitorLGCtxDestroy));
    200:     ctx->semilogy = PETSC_TRUE;
    201:   }

    203:   PetscCall(PetscOptionsName("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", &opt));
    204:   if (opt) {
    205:     TSMonitorLGCtx ctx;
    206:     PetscInt       howoften = 1;

    208:     PetscCall(PetscOptionsInt("-ts_monitor_lg_snes_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGSNESIterations", howoften, &howoften, NULL));
    209:     PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
    210:     PetscCall(TSMonitorSet(ts, TSMonitorLGSNESIterations, ctx, (PetscErrorCode (*)(void **))TSMonitorLGCtxDestroy));
    211:   }
    212:   PetscCall(PetscOptionsName("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", &opt));
    213:   if (opt) {
    214:     TSMonitorLGCtx ctx;
    215:     PetscInt       howoften = 1;

    217:     PetscCall(PetscOptionsInt("-ts_monitor_lg_ksp_iterations", "Monitor number nonlinear iterations for each timestep graphically", "TSMonitorLGKSPIterations", howoften, &howoften, NULL));
    218:     PetscCall(TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, howoften, &ctx));
    219:     PetscCall(TSMonitorSet(ts, TSMonitorLGKSPIterations, ctx, (PetscErrorCode (*)(void **))TSMonitorLGCtxDestroy));
    220:   }
    221:   PetscCall(PetscOptionsName("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", &opt));
    222:   if (opt) {
    223:     TSMonitorSPEigCtx ctx;
    224:     PetscInt          howoften = 1;

    226:     PetscCall(PetscOptionsInt("-ts_monitor_sp_eig", "Monitor eigenvalues of linearized operator graphically", "TSMonitorSPEig", howoften, &howoften, NULL));
    227:     PetscCall(TSMonitorSPEigCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
    228:     PetscCall(TSMonitorSet(ts, TSMonitorSPEig, ctx, (PetscErrorCode (*)(void **))TSMonitorSPEigCtxDestroy));
    229:   }
    230:   PetscCall(PetscOptionsName("-ts_monitor_sp_swarm", "Display particle phase space from the DMSwarm", "TSMonitorSPSwarm", &opt));
    231:   if (opt) {
    232:     TSMonitorSPCtx ctx;
    233:     PetscInt       howoften = 1, retain = 0;
    234:     PetscBool      phase = PETSC_TRUE, create = PETSC_TRUE, multispecies = PETSC_FALSE;

    236:     for (PetscInt i = 0; i < ts->numbermonitors; ++i)
    237:       if (ts->monitor[i] == TSMonitorSPSwarmSolution) {
    238:         create = PETSC_FALSE;
    239:         break;
    240:       }
    241:     if (create) {
    242:       PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm", "Display particles phase space from the DMSwarm", "TSMonitorSPSwarm", howoften, &howoften, NULL));
    243:       PetscCall(PetscOptionsInt("-ts_monitor_sp_swarm_retain", "Retain n points plotted to show trajectory, -1 for all points", "TSMonitorSPSwarm", retain, &retain, NULL));
    244:       PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_phase", "Plot in phase space rather than coordinate space", "TSMonitorSPSwarm", phase, &phase, NULL));
    245:       PetscCall(PetscOptionsBool("-ts_monitor_sp_swarm_multi_species", "Color particles by particle species", "TSMonitorSPSwarm", multispecies, &multispecies, NULL));
    246:       PetscCall(TSMonitorSPCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, retain, phase, multispecies, &ctx));
    247:       PetscCall(TSMonitorSet(ts, TSMonitorSPSwarmSolution, ctx, (PetscErrorCode (*)(void **))TSMonitorSPCtxDestroy));
    248:     }
    249:   }
    250:   PetscCall(PetscOptionsName("-ts_monitor_hg_swarm", "Display particle histogram from the DMSwarm", "TSMonitorHGSwarm", &opt));
    251:   if (opt) {
    252:     TSMonitorHGCtx ctx;
    253:     PetscInt       howoften = 1, Ns = 1;
    254:     PetscBool      velocity = PETSC_FALSE, create = PETSC_TRUE;

    256:     for (PetscInt i = 0; i < ts->numbermonitors; ++i)
    257:       if (ts->monitor[i] == TSMonitorHGSwarmSolution) {
    258:         create = PETSC_FALSE;
    259:         break;
    260:       }
    261:     if (create) {
    262:       DM       sw, dm;
    263:       PetscInt Nc, Nb;

    265:       PetscCall(TSGetDM(ts, &sw));
    266:       PetscCall(DMSwarmGetCellDM(sw, &dm));
    267:       PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &Nc));
    268:       Nb = PetscMin(20, PetscMax(10, Nc));
    269:       PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm", "Display particles histogram from the DMSwarm", "TSMonitorHGSwarm", howoften, &howoften, NULL));
    270:       PetscCall(PetscOptionsBool("-ts_monitor_hg_swarm_velocity", "Plot in velocity space rather than coordinate space", "TSMonitorHGSwarm", velocity, &velocity, NULL));
    271:       PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_species", "Number of species to histogram", "TSMonitorHGSwarm", Ns, &Ns, NULL));
    272:       PetscCall(PetscOptionsInt("-ts_monitor_hg_swarm_bins", "Number of histogram bins", "TSMonitorHGSwarm", Nb, &Nb, NULL));
    273:       PetscCall(TSMonitorHGCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, Ns, Nb, velocity, &ctx));
    274:       PetscCall(TSMonitorSet(ts, TSMonitorHGSwarmSolution, ctx, (PetscErrorCode (*)(void **))TSMonitorHGCtxDestroy));
    275:     }
    276:   }
    277:   opt = PETSC_FALSE;
    278:   PetscCall(PetscOptionsName("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", &opt));
    279:   if (opt) {
    280:     TSMonitorDrawCtx ctx;
    281:     PetscInt         howoften = 1;

    283:     PetscCall(PetscOptionsInt("-ts_monitor_draw_solution", "Monitor solution graphically", "TSMonitorDrawSolution", howoften, &howoften, NULL));
    284:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Computed Solution", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
    285:     PetscCall(TSMonitorSet(ts, TSMonitorDrawSolution, ctx, (PetscErrorCode (*)(void **))TSMonitorDrawCtxDestroy));
    286:   }
    287:   opt = PETSC_FALSE;
    288:   PetscCall(PetscOptionsName("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", &opt));
    289:   if (opt) {
    290:     TSMonitorDrawCtx ctx;
    291:     PetscReal        bounds[4];
    292:     PetscInt         n = 4;
    293:     PetscDraw        draw;
    294:     PetscDrawAxis    axis;

    296:     PetscCall(PetscOptionsRealArray("-ts_monitor_draw_solution_phase", "Monitor solution graphically", "TSMonitorDrawSolutionPhase", bounds, &n, NULL));
    297:     PetscCheck(n == 4, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Must provide bounding box of phase field");
    298:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 300, 300, 1, &ctx));
    299:     PetscCall(PetscViewerDrawGetDraw(ctx->viewer, 0, &draw));
    300:     PetscCall(PetscViewerDrawGetDrawAxis(ctx->viewer, 0, &axis));
    301:     PetscCall(PetscDrawAxisSetLimits(axis, bounds[0], bounds[2], bounds[1], bounds[3]));
    302:     PetscCall(PetscDrawAxisSetLabels(axis, "Phase Diagram", "Variable 1", "Variable 2"));
    303:     PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionPhase, ctx, (PetscErrorCode (*)(void **))TSMonitorDrawCtxDestroy));
    304:   }
    305:   opt = PETSC_FALSE;
    306:   PetscCall(PetscOptionsName("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", &opt));
    307:   if (opt) {
    308:     TSMonitorDrawCtx ctx;
    309:     PetscInt         howoften = 1;

    311:     PetscCall(PetscOptionsInt("-ts_monitor_draw_error", "Monitor error graphically", "TSMonitorDrawError", howoften, &howoften, NULL));
    312:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Error", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
    313:     PetscCall(TSMonitorSet(ts, TSMonitorDrawError, ctx, (PetscErrorCode (*)(void **))TSMonitorDrawCtxDestroy));
    314:   }
    315:   opt = PETSC_FALSE;
    316:   PetscCall(PetscOptionsName("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", &opt));
    317:   if (opt) {
    318:     TSMonitorDrawCtx ctx;
    319:     PetscInt         howoften = 1;

    321:     PetscCall(PetscOptionsInt("-ts_monitor_draw_solution_function", "Monitor solution provided by TSMonitorSetSolutionFunction() graphically", "TSMonitorDrawSolutionFunction", howoften, &howoften, NULL));
    322:     PetscCall(TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts), NULL, "Solution provided by user function", PETSC_DECIDE, PETSC_DECIDE, 300, 300, howoften, &ctx));
    323:     PetscCall(TSMonitorSet(ts, TSMonitorDrawSolutionFunction, ctx, (PetscErrorCode (*)(void **))TSMonitorDrawCtxDestroy));
    324:   }

    326:   opt = PETSC_FALSE;
    327:   PetscCall(PetscOptionsString("-ts_monitor_solution_vtk", "Save each time step to a binary file, use filename-%%03" PetscInt_FMT ".vts", "TSMonitorSolutionVTK", NULL, monfilename, sizeof(monfilename), &flg));
    328:   if (flg) {
    329:     TSMonitorVTKCtx ctx;

    331:     PetscCall(TSMonitorSolutionVTKCtxCreate(monfilename, &ctx));
    332:     PetscCall(PetscOptionsInt("-ts_monitor_solution_vtk_interval", "Save every interval time step (-1 for last step only)", NULL, ctx->interval, &ctx->interval, NULL));
    333:     PetscCall(TSMonitorSet(ts, (PetscErrorCode (*)(TS, PetscInt, PetscReal, Vec, void *))TSMonitorSolutionVTK, ctx, (PetscErrorCode (*)(void **))TSMonitorSolutionVTKDestroy));
    334:   }

    336:   PetscCall(PetscOptionsString("-ts_monitor_dmda_ray", "Display a ray of the solution", "None", "y=0", dir, sizeof(dir), &flg));
    337:   if (flg) {
    338:     TSMonitorDMDARayCtx *rayctx;
    339:     int                  ray = 0;
    340:     DMDirection          ddir;
    341:     DM                   da;
    342:     PetscMPIInt          rank;

    344:     PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
    345:     if (dir[0] == 'x') ddir = DM_X;
    346:     else if (dir[0] == 'y') ddir = DM_Y;
    347:     else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray %s", dir);
    348:     sscanf(dir + 2, "%d", &ray);

    350:     PetscCall(PetscInfo((PetscObject)ts, "Displaying DMDA ray %c = %d\n", dir[0], ray));
    351:     PetscCall(PetscNew(&rayctx));
    352:     PetscCall(TSGetDM(ts, &da));
    353:     PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
    354:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ts), &rank));
    355:     if (rank == 0) PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, NULL, 0, 0, 600, 300, &rayctx->viewer));
    356:     rayctx->lgctx = NULL;
    357:     PetscCall(TSMonitorSet(ts, TSMonitorDMDARay, rayctx, TSMonitorDMDARayDestroy));
    358:   }
    359:   PetscCall(PetscOptionsString("-ts_monitor_lg_dmda_ray", "Display a ray of the solution", "None", "x=0", dir, sizeof(dir), &flg));
    360:   if (flg) {
    361:     TSMonitorDMDARayCtx *rayctx;
    362:     int                  ray = 0;
    363:     DMDirection          ddir;
    364:     DM                   da;
    365:     PetscInt             howoften = 1;

    367:     PetscCheck(dir[1] == '=', PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir);
    368:     if (dir[0] == 'x') ddir = DM_X;
    369:     else if (dir[0] == 'y') ddir = DM_Y;
    370:     else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir);
    371:     sscanf(dir + 2, "%d", &ray);

    373:     PetscCall(PetscInfo((PetscObject)ts, "Displaying LG DMDA ray %c = %d\n", dir[0], ray));
    374:     PetscCall(PetscNew(&rayctx));
    375:     PetscCall(TSGetDM(ts, &da));
    376:     PetscCall(DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter));
    377:     PetscCall(TSMonitorLGCtxCreate(PETSC_COMM_SELF, NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 600, 400, howoften, &rayctx->lgctx));
    378:     PetscCall(TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy));
    379:   }

    381:   PetscCall(PetscOptionsName("-ts_monitor_envelope", "Monitor maximum and minimum value of each component of the solution", "TSMonitorEnvelope", &opt));
    382:   if (opt) {
    383:     TSMonitorEnvelopeCtx ctx;

    385:     PetscCall(TSMonitorEnvelopeCtxCreate(ts, &ctx));
    386:     PetscCall(TSMonitorSet(ts, TSMonitorEnvelope, ctx, (PetscErrorCode (*)(void **))TSMonitorEnvelopeCtxDestroy));
    387:   }
    388:   flg = PETSC_FALSE;
    389:   PetscCall(PetscOptionsBool("-ts_monitor_cancel", "Remove all monitors", "TSMonitorCancel", flg, &flg, &opt));
    390:   if (opt && flg) PetscCall(TSMonitorCancel(ts));

    392:   flg = PETSC_FALSE;
    393:   PetscCall(PetscOptionsBool("-ts_fd_color", "Use finite differences with coloring to compute IJacobian", "TSComputeIJacobianDefaultColor", flg, &flg, NULL));
    394:   if (flg) {
    395:     DM dm;

    397:     PetscCall(TSGetDM(ts, &dm));
    398:     PetscCall(DMTSUnsetIJacobianContext_Internal(dm));
    399:     PetscCall(TSSetIJacobian(ts, NULL, NULL, TSComputeIJacobianDefaultColor, NULL));
    400:     PetscCall(PetscInfo(ts, "Setting default finite difference coloring Jacobian matrix\n"));
    401:   }

    403:   /* Handle specific TS options */
    404:   PetscTryTypeMethod(ts, setfromoptions, PetscOptionsObject);

    406:   /* Handle TSAdapt options */
    407:   PetscCall(TSGetAdapt(ts, &ts->adapt));
    408:   PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));
    409:   PetscCall(TSAdaptSetFromOptions(ts->adapt, PetscOptionsObject));

    411:   /* TS trajectory must be set after TS, since it may use some TS options above */
    412:   tflg = ts->trajectory ? PETSC_TRUE : PETSC_FALSE;
    413:   PetscCall(PetscOptionsBool("-ts_save_trajectory", "Save the solution at each timestep", "TSSetSaveTrajectory", tflg, &tflg, NULL));
    414:   if (tflg) PetscCall(TSSetSaveTrajectory(ts));

    416:   PetscCall(TSAdjointSetFromOptions(ts, PetscOptionsObject));

    418:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
    419:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)ts, PetscOptionsObject));
    420:   PetscOptionsEnd();

    422:   if (ts->trajectory) PetscCall(TSTrajectorySetFromOptions(ts->trajectory, ts));

    424:   /* why do we have to do this here and not during TSSetUp? */
    425:   PetscCall(TSGetSNES(ts, &ts->snes));
    426:   if (ts->problem_type == TS_LINEAR) {
    427:     PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &flg, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
    428:     if (!flg) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
    429:   }
    430:   PetscCall(SNESSetFromOptions(ts->snes));
    431:   PetscFunctionReturn(PETSC_SUCCESS);
    432: }

    434: /*@
    435:   TSGetTrajectory - Gets the trajectory from a `TS` if it exists

    437:   Collective

    439:   Input Parameter:
    440: . ts - the `TS` context obtained from `TSCreate()`

    442:   Output Parameter:
    443: . tr - the `TSTrajectory` object, if it exists

    445:   Level: advanced

    447:   Note:
    448:   This routine should be called after all `TS` options have been set

    450: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSAdjointSolve()`, `TSTrajectoryCreate()`
    451: @*/
    452: PetscErrorCode TSGetTrajectory(TS ts, TSTrajectory *tr)
    453: {
    454:   PetscFunctionBegin;
    456:   *tr = ts->trajectory;
    457:   PetscFunctionReturn(PETSC_SUCCESS);
    458: }

    460: /*@
    461:   TSSetSaveTrajectory - Causes the `TS` to save its solutions as it iterates forward in time in a `TSTrajectory` object

    463:   Collective

    465:   Input Parameter:
    466: . ts - the `TS` context obtained from `TSCreate()`

    468:   Options Database Keys:
    469: + -ts_save_trajectory      - saves the trajectory to a file
    470: - -ts_trajectory_type type - set trajectory type

    472:   Level: intermediate

    474:   Notes:
    475:   This routine should be called after all `TS` options have been set

    477:   The `TSTRAJECTORYVISUALIZATION` files can be loaded into Python with $PETSC_DIR/lib/petsc/bin/PetscBinaryIOTrajectory.py and
    478:   MATLAB with $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m

    480: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`
    481: @*/
    482: PetscErrorCode TSSetSaveTrajectory(TS ts)
    483: {
    484:   PetscFunctionBegin;
    486:   if (!ts->trajectory) PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
    487:   PetscFunctionReturn(PETSC_SUCCESS);
    488: }

    490: /*@
    491:   TSResetTrajectory - Destroys and recreates the internal `TSTrajectory` object

    493:   Collective

    495:   Input Parameter:
    496: . ts - the `TS` context obtained from `TSCreate()`

    498:   Level: intermediate

    500: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`, `TSAdjointSolve()`, `TSRemoveTrajectory()`
    501: @*/
    502: PetscErrorCode TSResetTrajectory(TS ts)
    503: {
    504:   PetscFunctionBegin;
    506:   if (ts->trajectory) {
    507:     PetscCall(TSTrajectoryDestroy(&ts->trajectory));
    508:     PetscCall(TSTrajectoryCreate(PetscObjectComm((PetscObject)ts), &ts->trajectory));
    509:   }
    510:   PetscFunctionReturn(PETSC_SUCCESS);
    511: }

    513: /*@
    514:   TSRemoveTrajectory - Destroys and removes the internal `TSTrajectory` object from a `TS`

    516:   Collective

    518:   Input Parameter:
    519: . ts - the `TS` context obtained from `TSCreate()`

    521:   Level: intermediate

    523: .seealso: [](ch_ts), `TSTrajectory`, `TSResetTrajectory()`, `TSAdjointSolve()`
    524: @*/
    525: PetscErrorCode TSRemoveTrajectory(TS ts)
    526: {
    527:   PetscFunctionBegin;
    529:   if (ts->trajectory) PetscCall(TSTrajectoryDestroy(&ts->trajectory));
    530:   PetscFunctionReturn(PETSC_SUCCESS);
    531: }

    533: /*@
    534:   TSComputeRHSJacobian - Computes the Jacobian matrix that has been
    535:   set with `TSSetRHSJacobian()`.

    537:   Collective

    539:   Input Parameters:
    540: + ts - the `TS` context
    541: . t  - current timestep
    542: - U  - input vector

    544:   Output Parameters:
    545: + A - Jacobian matrix
    546: - B - optional preconditioning matrix

    548:   Level: developer

    550:   Note:
    551:   Most users should not need to explicitly call this routine, as it
    552:   is used internally within the nonlinear solvers.

    554: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `KSPSetOperators()`
    555: @*/
    556: PetscErrorCode TSComputeRHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B)
    557: {
    558:   PetscObjectState Ustate;
    559:   PetscObjectId    Uid;
    560:   DM               dm;
    561:   DMTS             tsdm;
    562:   TSRHSJacobianFn *rhsjacobianfunc;
    563:   void            *ctx;
    564:   TSRHSFunctionFn *rhsfunction;

    566:   PetscFunctionBegin;
    569:   PetscCheckSameComm(ts, 1, U, 3);
    570:   PetscCall(TSGetDM(ts, &dm));
    571:   PetscCall(DMGetDMTS(dm, &tsdm));
    572:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
    573:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobianfunc, &ctx));
    574:   PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
    575:   PetscCall(PetscObjectGetId((PetscObject)U, &Uid));

    577:   if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) PetscFunctionReturn(PETSC_SUCCESS);

    579:   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);
    580:   if (rhsjacobianfunc) {
    581:     PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
    582:     PetscCallBack("TS callback Jacobian", (*rhsjacobianfunc)(ts, t, U, A, B, ctx));
    583:     ts->rhsjacs++;
    584:     PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
    585:   } else {
    586:     PetscCall(MatZeroEntries(A));
    587:     if (B && A != B) PetscCall(MatZeroEntries(B));
    588:   }
    589:   ts->rhsjacobian.time  = t;
    590:   ts->rhsjacobian.shift = 0;
    591:   ts->rhsjacobian.scale = 1.;
    592:   PetscCall(PetscObjectGetId((PetscObject)U, &ts->rhsjacobian.Xid));
    593:   PetscCall(PetscObjectStateGet((PetscObject)U, &ts->rhsjacobian.Xstate));
    594:   PetscFunctionReturn(PETSC_SUCCESS);
    595: }

    597: /*@
    598:   TSComputeRHSFunction - Evaluates the right-hand-side function for a `TS`

    600:   Collective

    602:   Input Parameters:
    603: + ts - the `TS` context
    604: . t  - current time
    605: - U  - state vector

    607:   Output Parameter:
    608: . y - right-hand side

    610:   Level: developer

    612:   Note:
    613:   Most users should not need to explicitly call this routine, as it
    614:   is used internally within the nonlinear solvers.

    616: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
    617: @*/
    618: PetscErrorCode TSComputeRHSFunction(TS ts, PetscReal t, Vec U, Vec y)
    619: {
    620:   TSRHSFunctionFn *rhsfunction;
    621:   TSIFunctionFn   *ifunction;
    622:   void            *ctx;
    623:   DM               dm;

    625:   PetscFunctionBegin;
    629:   PetscCall(TSGetDM(ts, &dm));
    630:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, &ctx));
    631:   PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));

    633:   PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");

    635:   if (rhsfunction) {
    636:     PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, y, 0));
    637:     PetscCall(VecLockReadPush(U));
    638:     PetscCallBack("TS callback right-hand-side", (*rhsfunction)(ts, t, U, y, ctx));
    639:     PetscCall(VecLockReadPop(U));
    640:     ts->rhsfuncs++;
    641:     PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, y, 0));
    642:   } else PetscCall(VecZeroEntries(y));
    643:   PetscFunctionReturn(PETSC_SUCCESS);
    644: }

    646: /*@
    647:   TSComputeSolutionFunction - Evaluates the solution function.

    649:   Collective

    651:   Input Parameters:
    652: + ts - the `TS` context
    653: - t  - current time

    655:   Output Parameter:
    656: . U - the solution

    658:   Level: developer

    660: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
    661: @*/
    662: PetscErrorCode TSComputeSolutionFunction(TS ts, PetscReal t, Vec U)
    663: {
    664:   TSSolutionFn *solutionfunction;
    665:   void         *ctx;
    666:   DM            dm;

    668:   PetscFunctionBegin;
    671:   PetscCall(TSGetDM(ts, &dm));
    672:   PetscCall(DMTSGetSolutionFunction(dm, &solutionfunction, &ctx));
    673:   if (solutionfunction) PetscCallBack("TS callback solution", (*solutionfunction)(ts, t, U, ctx));
    674:   PetscFunctionReturn(PETSC_SUCCESS);
    675: }
    676: /*@
    677:   TSComputeForcingFunction - Evaluates the forcing function.

    679:   Collective

    681:   Input Parameters:
    682: + ts - the `TS` context
    683: - t  - current time

    685:   Output Parameter:
    686: . U - the function value

    688:   Level: developer

    690: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSSetRHSFunction()`, `TSComputeIFunction()`
    691: @*/
    692: PetscErrorCode TSComputeForcingFunction(TS ts, PetscReal t, Vec U)
    693: {
    694:   void        *ctx;
    695:   DM           dm;
    696:   TSForcingFn *forcing;

    698:   PetscFunctionBegin;
    701:   PetscCall(TSGetDM(ts, &dm));
    702:   PetscCall(DMTSGetForcingFunction(dm, &forcing, &ctx));

    704:   if (forcing) PetscCallBack("TS callback forcing function", (*forcing)(ts, t, U, ctx));
    705:   PetscFunctionReturn(PETSC_SUCCESS);
    706: }

    708: PetscErrorCode TSGetRHSMats_Private(TS ts, Mat *Arhs, Mat *Brhs)
    709: {
    710:   Mat            A, B;
    711:   TSIJacobianFn *ijacobian;

    713:   PetscFunctionBegin;
    714:   if (Arhs) *Arhs = NULL;
    715:   if (Brhs) *Brhs = NULL;
    716:   PetscCall(TSGetIJacobian(ts, &A, &B, &ijacobian, NULL));
    717:   if (Arhs) {
    718:     if (!ts->Arhs) {
    719:       if (ijacobian) {
    720:         PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &ts->Arhs));
    721:         PetscCall(TSSetMatStructure(ts, SAME_NONZERO_PATTERN));
    722:       } else {
    723:         ts->Arhs = A;
    724:         PetscCall(PetscObjectReference((PetscObject)A));
    725:       }
    726:     } else {
    727:       PetscBool flg;
    728:       PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
    729:       /* Handle case where user provided only RHSJacobian and used -snes_mf_operator */
    730:       if (flg && !ijacobian && ts->Arhs == ts->Brhs) {
    731:         PetscCall(PetscObjectDereference((PetscObject)ts->Arhs));
    732:         ts->Arhs = A;
    733:         PetscCall(PetscObjectReference((PetscObject)A));
    734:       }
    735:     }
    736:     *Arhs = ts->Arhs;
    737:   }
    738:   if (Brhs) {
    739:     if (!ts->Brhs) {
    740:       if (A != B) {
    741:         if (ijacobian) {
    742:           PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &ts->Brhs));
    743:         } else {
    744:           ts->Brhs = B;
    745:           PetscCall(PetscObjectReference((PetscObject)B));
    746:         }
    747:       } else {
    748:         PetscCall(PetscObjectReference((PetscObject)ts->Arhs));
    749:         ts->Brhs = ts->Arhs;
    750:       }
    751:     }
    752:     *Brhs = ts->Brhs;
    753:   }
    754:   PetscFunctionReturn(PETSC_SUCCESS);
    755: }

    757: /*@
    758:   TSComputeIFunction - Evaluates the DAE residual written in the implicit form F(t,U,Udot)=0

    760:   Collective

    762:   Input Parameters:
    763: + ts   - the `TS` context
    764: . t    - current time
    765: . U    - state vector
    766: . Udot - time derivative of state vector
    767: - imex - flag indicates if the method is `TSARKIMEX` so that the RHSFunction should be kept separate

    769:   Output Parameter:
    770: . Y - right-hand side

    772:   Level: developer

    774:   Note:
    775:   Most users should not need to explicitly call this routine, as it
    776:   is used internally within the nonlinear solvers.

    778:   If the user did not write their equations in implicit form, this
    779:   function recasts them in implicit form.

    781: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSComputeRHSFunction()`
    782: @*/
    783: PetscErrorCode TSComputeIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec Y, PetscBool imex)
    784: {
    785:   TSIFunctionFn   *ifunction;
    786:   TSRHSFunctionFn *rhsfunction;
    787:   void            *ctx;
    788:   DM               dm;

    790:   PetscFunctionBegin;

    796:   PetscCall(TSGetDM(ts, &dm));
    797:   PetscCall(DMTSGetIFunction(dm, &ifunction, &ctx));
    798:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));

    800:   PetscCheck(rhsfunction || ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSFunction() and / or TSSetIFunction()");

    802:   PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, Udot, Y));
    803:   if (ifunction) {
    804:     PetscCallBack("TS callback implicit function", (*ifunction)(ts, t, U, Udot, Y, ctx));
    805:     ts->ifuncs++;
    806:   }
    807:   if (imex) {
    808:     if (!ifunction) PetscCall(VecCopy(Udot, Y));
    809:   } else if (rhsfunction) {
    810:     if (ifunction) {
    811:       Vec Frhs;

    813:       PetscCall(DMGetGlobalVector(dm, &Frhs));
    814:       PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
    815:       PetscCall(VecAXPY(Y, -1, Frhs));
    816:       PetscCall(DMRestoreGlobalVector(dm, &Frhs));
    817:     } else {
    818:       PetscCall(TSComputeRHSFunction(ts, t, U, Y));
    819:       PetscCall(VecAYPX(Y, -1, Udot));
    820:     }
    821:   }
    822:   PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, Udot, Y));
    823:   PetscFunctionReturn(PETSC_SUCCESS);
    824: }

    826: /*
    827:    TSRecoverRHSJacobian - Recover the Jacobian matrix so that one can call `TSComputeRHSJacobian()` on it.

    829:    Note:
    830:    This routine is needed when one switches from `TSComputeIJacobian()` to `TSComputeRHSJacobian()` because the Jacobian matrix may be shifted or scaled in `TSComputeIJacobian()`.

    832: */
    833: static PetscErrorCode TSRecoverRHSJacobian(TS ts, Mat A, Mat B)
    834: {
    835:   PetscFunctionBegin;
    837:   PetscCheck(A == ts->Arhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Amat");
    838:   PetscCheck(B == ts->Brhs, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Invalid Bmat");

    840:   if (ts->rhsjacobian.shift) PetscCall(MatShift(A, -ts->rhsjacobian.shift));
    841:   if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(A, -1));
    842:   if (B && B == ts->Brhs && A != B) {
    843:     if (ts->rhsjacobian.shift) PetscCall(MatShift(B, -ts->rhsjacobian.shift));
    844:     if (ts->rhsjacobian.scale == -1.) PetscCall(MatScale(B, -1));
    845:   }
    846:   ts->rhsjacobian.shift = 0;
    847:   ts->rhsjacobian.scale = 1.;
    848:   PetscFunctionReturn(PETSC_SUCCESS);
    849: }

    851: /*@
    852:   TSComputeIJacobian - Evaluates the Jacobian of the DAE

    854:   Collective

    856:   Input Parameters:
    857: + ts    - the `TS` context
    858: . t     - current timestep
    859: . U     - state vector
    860: . Udot  - time derivative of state vector
    861: . shift - shift to apply, see note below
    862: - imex  - flag indicates if the method is `TSARKIMEX` so that the RHSJacobian should be kept separate

    864:   Output Parameters:
    865: + A - Jacobian matrix
    866: - B - matrix from which the preconditioner is constructed; often the same as `A`

    868:   Level: developer

    870:   Notes:
    871:   If F(t,U,Udot)=0 is the DAE, the required Jacobian is
    872: .vb
    873:    dF/dU + shift*dF/dUdot
    874: .ve
    875:   Most users should not need to explicitly call this routine, as it
    876:   is used internally within the nonlinear solvers.

    878: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`
    879: @*/
    880: PetscErrorCode TSComputeIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, PetscBool imex)
    881: {
    882:   TSIJacobianFn   *ijacobian;
    883:   TSRHSJacobianFn *rhsjacobian;
    884:   DM               dm;
    885:   void            *ctx;

    887:   PetscFunctionBegin;

    894:   PetscCall(TSGetDM(ts, &dm));
    895:   PetscCall(DMTSGetIJacobian(dm, &ijacobian, &ctx));
    896:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));

    898:   PetscCheck(rhsjacobian || ijacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_USER, "Must call TSSetRHSJacobian() and / or TSSetIJacobian()");

    900:   PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, A, B));
    901:   if (ijacobian) {
    902:     PetscCallBack("TS callback implicit Jacobian", (*ijacobian)(ts, t, U, Udot, shift, A, B, ctx));
    903:     ts->ijacs++;
    904:   }
    905:   if (imex) {
    906:     if (!ijacobian) { /* system was written as Udot = G(t,U) */
    907:       PetscBool assembled;
    908:       if (rhsjacobian) {
    909:         Mat Arhs = NULL;
    910:         PetscCall(TSGetRHSMats_Private(ts, &Arhs, NULL));
    911:         if (A == Arhs) {
    912:           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 */
    913:           ts->rhsjacobian.time = PETSC_MIN_REAL;
    914:         }
    915:       }
    916:       PetscCall(MatZeroEntries(A));
    917:       PetscCall(MatAssembled(A, &assembled));
    918:       if (!assembled) {
    919:         PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    920:         PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    921:       }
    922:       PetscCall(MatShift(A, shift));
    923:       if (A != B) {
    924:         PetscCall(MatZeroEntries(B));
    925:         PetscCall(MatAssembled(B, &assembled));
    926:         if (!assembled) {
    927:           PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    928:           PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    929:         }
    930:         PetscCall(MatShift(B, shift));
    931:       }
    932:     }
    933:   } else {
    934:     Mat Arhs = NULL, Brhs = NULL;

    936:     /* RHSJacobian needs to be converted to part of IJacobian if exists */
    937:     if (rhsjacobian) PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
    938:     if (Arhs == A) { /* No IJacobian matrix, so we only have the RHS matrix */
    939:       PetscObjectState Ustate;
    940:       PetscObjectId    Uid;
    941:       TSRHSFunctionFn *rhsfunction;

    943:       PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));
    944:       PetscCall(PetscObjectStateGet((PetscObject)U, &Ustate));
    945:       PetscCall(PetscObjectGetId((PetscObject)U, &Uid));
    946:       if ((rhsjacobian == TSComputeRHSJacobianConstant || (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && rhsfunction != TSComputeRHSFunctionLinear)) &&
    947:           ts->rhsjacobian.scale == -1.) {                      /* No need to recompute RHSJacobian */
    948:         PetscCall(MatShift(A, shift - ts->rhsjacobian.shift)); /* revert the old shift and add the new shift with a single call to MatShift */
    949:         if (A != B) PetscCall(MatShift(B, shift - ts->rhsjacobian.shift));
    950:       } else {
    951:         PetscBool flg;

    953:         if (ts->rhsjacobian.reuse) { /* Undo the damage */
    954:           /* MatScale has a short path for this case.
    955:              However, this code path is taken the first time TSComputeRHSJacobian is called
    956:              and the matrices have not been assembled yet */
    957:           PetscCall(TSRecoverRHSJacobian(ts, A, B));
    958:         }
    959:         PetscCall(TSComputeRHSJacobian(ts, t, U, A, B));
    960:         PetscCall(SNESGetUseMatrixFree(ts->snes, NULL, &flg));
    961:         /* since -snes_mf_operator uses the full SNES function it does not need to be shifted or scaled here */
    962:         if (!flg) {
    963:           PetscCall(MatScale(A, -1));
    964:           PetscCall(MatShift(A, shift));
    965:         }
    966:         if (A != B) {
    967:           PetscCall(MatScale(B, -1));
    968:           PetscCall(MatShift(B, shift));
    969:         }
    970:       }
    971:       ts->rhsjacobian.scale = -1;
    972:       ts->rhsjacobian.shift = shift;
    973:     } else if (Arhs) {  /* Both IJacobian and RHSJacobian */
    974:       if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
    975:         PetscCall(MatZeroEntries(A));
    976:         PetscCall(MatShift(A, shift));
    977:         if (A != B) {
    978:           PetscCall(MatZeroEntries(B));
    979:           PetscCall(MatShift(B, shift));
    980:         }
    981:       }
    982:       PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
    983:       PetscCall(MatAXPY(A, -1, Arhs, ts->axpy_pattern));
    984:       if (A != B) PetscCall(MatAXPY(B, -1, Brhs, ts->axpy_pattern));
    985:     }
    986:   }
    987:   PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, A, B));
    988:   PetscFunctionReturn(PETSC_SUCCESS);
    989: }

    991: /*@C
    992:   TSSetRHSFunction - Sets the routine for evaluating the function,
    993:   where U_t = G(t,u).

    995:   Logically Collective

    997:   Input Parameters:
    998: + ts  - the `TS` context obtained from `TSCreate()`
    999: . r   - vector to put the computed right-hand side (or `NULL` to have it created)
    1000: . f   - routine for evaluating the right-hand-side function
    1001: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`)

    1003:   Level: beginner

    1005:   Note:
    1006:   You must call this function or `TSSetIFunction()` to define your ODE. You cannot use this function when solving a DAE.

    1008: .seealso: [](ch_ts), `TS`, `TSRHSFunctionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSSetIFunction()`
    1009: @*/
    1010: PetscErrorCode TSSetRHSFunction(TS ts, Vec r, TSRHSFunctionFn *f, void *ctx)
    1011: {
    1012:   SNES snes;
    1013:   Vec  ralloc = NULL;
    1014:   DM   dm;

    1016:   PetscFunctionBegin;

    1020:   PetscCall(TSGetDM(ts, &dm));
    1021:   PetscCall(DMTSSetRHSFunction(dm, f, ctx));
    1022:   PetscCall(TSGetSNES(ts, &snes));
    1023:   if (!r && !ts->dm && ts->vec_sol) {
    1024:     PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
    1025:     r = ralloc;
    1026:   }
    1027:   PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
    1028:   PetscCall(VecDestroy(&ralloc));
    1029:   PetscFunctionReturn(PETSC_SUCCESS);
    1030: }

    1032: /*@C
    1033:   TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE

    1035:   Logically Collective

    1037:   Input Parameters:
    1038: + ts  - the `TS` context obtained from `TSCreate()`
    1039: . f   - routine for evaluating the solution
    1040: - ctx - [optional] user-defined context for private data for the
    1041:           function evaluation routine (may be `NULL`)

    1043:   Options Database Keys:
    1044: + -ts_monitor_lg_error   - create a graphical monitor of error history, requires user to have provided `TSSetSolutionFunction()`
    1045: - -ts_monitor_draw_error - Monitor error graphically, requires user to have provided `TSSetSolutionFunction()`

    1047:   Level: intermediate

    1049:   Notes:
    1050:   This routine is used for testing accuracy of time integration schemes when you already know the solution.
    1051:   If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to
    1052:   create closed-form solutions with non-physical forcing terms.

    1054:   For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.

    1056: .seealso: [](ch_ts), `TS`, `TSSolutionFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`, `TSComputeSolutionFunction()`, `TSSetForcingFunction()`, `TSSetSolution()`, `TSGetSolution()`, `TSMonitorLGError()`, `TSMonitorDrawError()`
    1057: @*/
    1058: PetscErrorCode TSSetSolutionFunction(TS ts, TSSolutionFn *f, void *ctx)
    1059: {
    1060:   DM dm;

    1062:   PetscFunctionBegin;
    1064:   PetscCall(TSGetDM(ts, &dm));
    1065:   PetscCall(DMTSSetSolutionFunction(dm, f, ctx));
    1066:   PetscFunctionReturn(PETSC_SUCCESS);
    1067: }

    1069: /*@C
    1070:   TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE

    1072:   Logically Collective

    1074:   Input Parameters:
    1075: + ts   - the `TS` context obtained from `TSCreate()`
    1076: . func - routine for evaluating the forcing function
    1077: - ctx  - [optional] user-defined context for private data for the function evaluation routine
    1078:          (may be `NULL`)

    1080:   Level: intermediate

    1082:   Notes:
    1083:   This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to
    1084:   create closed-form solutions with a non-physical forcing term. It allows you to use the Method of Manufactored Solution without directly editing the
    1085:   definition of the problem you are solving and hence possibly introducing bugs.

    1087:   This replaces the ODE F(u,u_t,t) = 0 the `TS` is solving with F(u,u_t,t) - func(t) = 0

    1089:   This forcing function does not depend on the solution to the equations, it can only depend on spatial location, time, and possibly parameters, the
    1090:   parameters can be passed in the ctx variable.

    1092:   For low-dimensional problems solved in serial, such as small discrete systems, `TSMonitorLGError()` can be used to monitor the error history.

    1094: .seealso: [](ch_ts), `TS`, `TSForcingFn`, `TSSetRHSJacobian()`, `TSSetIJacobian()`,
    1095: `TSComputeSolutionFunction()`, `TSSetSolutionFunction()`
    1096: @*/
    1097: PetscErrorCode TSSetForcingFunction(TS ts, TSForcingFn *func, void *ctx)
    1098: {
    1099:   DM dm;

    1101:   PetscFunctionBegin;
    1103:   PetscCall(TSGetDM(ts, &dm));
    1104:   PetscCall(DMTSSetForcingFunction(dm, func, ctx));
    1105:   PetscFunctionReturn(PETSC_SUCCESS);
    1106: }

    1108: /*@C
    1109:   TSSetRHSJacobian - Sets the function to compute the Jacobian of G,
    1110:   where U_t = G(U,t), as well as the location to store the matrix.

    1112:   Logically Collective

    1114:   Input Parameters:
    1115: + ts   - the `TS` context obtained from `TSCreate()`
    1116: . Amat - (approximate) location to store Jacobian matrix entries computed by `f`
    1117: . Pmat - matrix from which preconditioner is to be constructed (usually the same as `Amat`)
    1118: . f    - the Jacobian evaluation routine
    1119: - ctx  - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)

    1121:   Level: beginner

    1123:   Notes:
    1124:   You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value

    1126:   The `TS` solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f()`
    1127:   You should not assume the values are the same in the next call to f() as you set them in the previous call.

    1129: .seealso: [](ch_ts), `TS`, `TSRHSJacobianFn`, `SNESComputeJacobianDefaultColor()`,
    1130: `TSSetRHSFunction()`, `TSRHSJacobianSetReuse()`, `TSSetIJacobian()`, `TSRHSFunctionFn`, `TSIFunctionFn`
    1131: @*/
    1132: PetscErrorCode TSSetRHSJacobian(TS ts, Mat Amat, Mat Pmat, TSRHSJacobianFn *f, void *ctx)
    1133: {
    1134:   SNES           snes;
    1135:   DM             dm;
    1136:   TSIJacobianFn *ijacobian;

    1138:   PetscFunctionBegin;
    1142:   if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
    1143:   if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);

    1145:   PetscCall(TSGetDM(ts, &dm));
    1146:   PetscCall(DMTSSetRHSJacobian(dm, f, ctx));
    1147:   PetscCall(DMTSGetIJacobian(dm, &ijacobian, NULL));
    1148:   PetscCall(TSGetSNES(ts, &snes));
    1149:   if (!ijacobian) PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
    1150:   if (Amat) {
    1151:     PetscCall(PetscObjectReference((PetscObject)Amat));
    1152:     PetscCall(MatDestroy(&ts->Arhs));
    1153:     ts->Arhs = Amat;
    1154:   }
    1155:   if (Pmat) {
    1156:     PetscCall(PetscObjectReference((PetscObject)Pmat));
    1157:     PetscCall(MatDestroy(&ts->Brhs));
    1158:     ts->Brhs = Pmat;
    1159:   }
    1160:   PetscFunctionReturn(PETSC_SUCCESS);
    1161: }

    1163: /*@C
    1164:   TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved.

    1166:   Logically Collective

    1168:   Input Parameters:
    1169: + ts  - the `TS` context obtained from `TSCreate()`
    1170: . r   - vector to hold the residual (or `NULL` to have it created internally)
    1171: . f   - the function evaluation routine
    1172: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)

    1174:   Level: beginner

    1176:   Note:
    1177:   The user MUST call either this routine or `TSSetRHSFunction()` to define the ODE.  When solving DAEs you must use this function.

    1179: .seealso: [](ch_ts), `TS`, `TSIFunctionFn`, `TSSetRHSJacobian()`, `TSSetRHSFunction()`,
    1180: `TSSetIJacobian()`
    1181: @*/
    1182: PetscErrorCode TSSetIFunction(TS ts, Vec r, TSIFunctionFn *f, void *ctx)
    1183: {
    1184:   SNES snes;
    1185:   Vec  ralloc = NULL;
    1186:   DM   dm;

    1188:   PetscFunctionBegin;

    1192:   PetscCall(TSGetDM(ts, &dm));
    1193:   PetscCall(DMTSSetIFunction(dm, f, ctx));

    1195:   PetscCall(TSGetSNES(ts, &snes));
    1196:   if (!r && !ts->dm && ts->vec_sol) {
    1197:     PetscCall(VecDuplicate(ts->vec_sol, &ralloc));
    1198:     r = ralloc;
    1199:   }
    1200:   PetscCall(SNESSetFunction(snes, r, SNESTSFormFunction, ts));
    1201:   PetscCall(VecDestroy(&ralloc));
    1202:   PetscFunctionReturn(PETSC_SUCCESS);
    1203: }

    1205: /*@C
    1206:   TSGetIFunction - Returns the vector where the implicit residual is stored and the function/context to compute it.

    1208:   Not Collective

    1210:   Input Parameter:
    1211: . ts - the `TS` context

    1213:   Output Parameters:
    1214: + r    - vector to hold residual (or `NULL`)
    1215: . func - the function to compute residual (or `NULL`)
    1216: - ctx  - the function context (or `NULL`)

    1218:   Level: advanced

    1220: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`
    1221: @*/
    1222: PetscErrorCode TSGetIFunction(TS ts, Vec *r, TSIFunctionFn **func, void **ctx)
    1223: {
    1224:   SNES snes;
    1225:   DM   dm;

    1227:   PetscFunctionBegin;
    1229:   PetscCall(TSGetSNES(ts, &snes));
    1230:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
    1231:   PetscCall(TSGetDM(ts, &dm));
    1232:   PetscCall(DMTSGetIFunction(dm, func, ctx));
    1233:   PetscFunctionReturn(PETSC_SUCCESS);
    1234: }

    1236: /*@C
    1237:   TSGetRHSFunction - Returns the vector where the right-hand side is stored and the function/context to compute it.

    1239:   Not Collective

    1241:   Input Parameter:
    1242: . ts - the `TS` context

    1244:   Output Parameters:
    1245: + r    - vector to hold computed right-hand side (or `NULL`)
    1246: . func - the function to compute right-hand side (or `NULL`)
    1247: - ctx  - the function context (or `NULL`)

    1249:   Level: advanced

    1251: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `SNESGetFunction()`
    1252: @*/
    1253: PetscErrorCode TSGetRHSFunction(TS ts, Vec *r, TSRHSFunctionFn **func, void **ctx)
    1254: {
    1255:   SNES snes;
    1256:   DM   dm;

    1258:   PetscFunctionBegin;
    1260:   PetscCall(TSGetSNES(ts, &snes));
    1261:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
    1262:   PetscCall(TSGetDM(ts, &dm));
    1263:   PetscCall(DMTSGetRHSFunction(dm, func, ctx));
    1264:   PetscFunctionReturn(PETSC_SUCCESS);
    1265: }

    1267: /*@C
    1268:   TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function
    1269:   provided with `TSSetIFunction()`.

    1271:   Logically Collective

    1273:   Input Parameters:
    1274: + ts   - the `TS` context obtained from `TSCreate()`
    1275: . Amat - (approximate) matrix to store Jacobian entries computed by `f`
    1276: . Pmat - matrix used to compute preconditioner (usually the same as `Amat`)
    1277: . f    - the Jacobian evaluation routine
    1278: - ctx  - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)

    1280:   Level: beginner

    1282:   Notes:
    1283:   The matrices `Amat` and `Pmat` are exactly the matrices that are used by `SNES` for the nonlinear solve.

    1285:   If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
    1286:   space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.

    1288:   The matrix dF/dU + a*dF/dU_t you provide turns out to be
    1289:   the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved.
    1290:   The time integrator internally approximates U_t by W+a*U where the positive "shift"
    1291:   a and vector W depend on the integration method, step size, and past states. For example with
    1292:   the backward Euler method a = 1/dt and W = -a*U(previous timestep) so
    1293:   W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt

    1295:   You must set all the diagonal entries of the matrices, if they are zero you must still set them with a zero value

    1297:   The TS solver may modify the nonzero structure and the entries of the matrices `Amat` and `Pmat` between the calls to `f`
    1298:   You should not assume the values are the same in the next call to `f` as you set them in the previous call.

    1300:   In case `TSSetRHSJacobian()` is also used in conjunction with a fully-implicit solver,
    1301:   multilevel linear solvers, e.g. `PCMG`, will likely not work due to the way `TS` handles rhs matrices.

    1303: .seealso: [](ch_ts), `TS`, `TSIJacobianFn`, `TSSetIFunction()`, `TSSetRHSJacobian()`,
    1304: `SNESComputeJacobianDefaultColor()`, `SNESComputeJacobianDefault()`, `TSSetRHSFunction()`
    1305: @*/
    1306: PetscErrorCode TSSetIJacobian(TS ts, Mat Amat, Mat Pmat, TSIJacobianFn *f, void *ctx)
    1307: {
    1308:   SNES snes;
    1309:   DM   dm;

    1311:   PetscFunctionBegin;
    1315:   if (Amat) PetscCheckSameComm(ts, 1, Amat, 2);
    1316:   if (Pmat) PetscCheckSameComm(ts, 1, Pmat, 3);

    1318:   PetscCall(TSGetDM(ts, &dm));
    1319:   PetscCall(DMTSSetIJacobian(dm, f, ctx));

    1321:   PetscCall(TSGetSNES(ts, &snes));
    1322:   PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESTSFormJacobian, ts));
    1323:   PetscFunctionReturn(PETSC_SUCCESS);
    1324: }

    1326: /*@
    1327:   TSRHSJacobianSetReuse - restore the RHS Jacobian before calling the user-provided `TSRHSJacobianFn` function again

    1329:   Logically Collective

    1331:   Input Parameters:
    1332: + ts    - `TS` context obtained from `TSCreate()`
    1333: - reuse - `PETSC_TRUE` if the RHS Jacobian

    1335:   Level: intermediate

    1337:   Notes:
    1338:   Without this flag, `TS` will change the sign and shift the RHS Jacobian for a
    1339:   finite-time-step implicit solve, in which case the user function will need to recompute the
    1340:   entire Jacobian.  The `reuse `flag must be set if the evaluation function assumes that the
    1341:   matrix entries have not been changed by the `TS`.

    1343: .seealso: [](ch_ts), `TS`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
    1344: @*/
    1345: PetscErrorCode TSRHSJacobianSetReuse(TS ts, PetscBool reuse)
    1346: {
    1347:   PetscFunctionBegin;
    1348:   ts->rhsjacobian.reuse = reuse;
    1349:   PetscFunctionReturn(PETSC_SUCCESS);
    1350: }

    1352: /*@C
    1353:   TSSetI2Function - Set the function to compute F(t,U,U_t,U_tt) where F = 0 is the DAE to be solved.

    1355:   Logically Collective

    1357:   Input Parameters:
    1358: + ts  - the `TS` context obtained from `TSCreate()`
    1359: . F   - vector to hold the residual (or `NULL` to have it created internally)
    1360: . fun - the function evaluation routine
    1361: - ctx - user-defined context for private data for the function evaluation routine (may be `NULL`)

    1363:   Level: beginner

    1365: .seealso: [](ch_ts), `TS`, `TSI2FunctionFn`, `TSSetI2Jacobian()`, `TSSetIFunction()`,
    1366: `TSCreate()`, `TSSetRHSFunction()`
    1367: @*/
    1368: PetscErrorCode TSSetI2Function(TS ts, Vec F, TSI2FunctionFn *fun, void *ctx)
    1369: {
    1370:   DM dm;

    1372:   PetscFunctionBegin;
    1375:   PetscCall(TSSetIFunction(ts, F, NULL, NULL));
    1376:   PetscCall(TSGetDM(ts, &dm));
    1377:   PetscCall(DMTSSetI2Function(dm, fun, ctx));
    1378:   PetscFunctionReturn(PETSC_SUCCESS);
    1379: }

    1381: /*@C
    1382:   TSGetI2Function - Returns the vector where the implicit residual is stored and the function/context to compute it.

    1384:   Not Collective

    1386:   Input Parameter:
    1387: . ts - the `TS` context

    1389:   Output Parameters:
    1390: + r   - vector to hold residual (or `NULL`)
    1391: . fun - the function to compute residual (or `NULL`)
    1392: - ctx - the function context (or `NULL`)

    1394:   Level: advanced

    1396: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `SNESGetFunction()`, `TSCreate()`
    1397: @*/
    1398: PetscErrorCode TSGetI2Function(TS ts, Vec *r, TSI2FunctionFn **fun, void **ctx)
    1399: {
    1400:   SNES snes;
    1401:   DM   dm;

    1403:   PetscFunctionBegin;
    1405:   PetscCall(TSGetSNES(ts, &snes));
    1406:   PetscCall(SNESGetFunction(snes, r, NULL, NULL));
    1407:   PetscCall(TSGetDM(ts, &dm));
    1408:   PetscCall(DMTSGetI2Function(dm, fun, ctx));
    1409:   PetscFunctionReturn(PETSC_SUCCESS);
    1410: }

    1412: /*@C
    1413:   TSSetI2Jacobian - Set the function to compute the matrix dF/dU + v*dF/dU_t  + a*dF/dU_tt
    1414:   where F(t,U,U_t,U_tt) is the function you provided with `TSSetI2Function()`.

    1416:   Logically Collective

    1418:   Input Parameters:
    1419: + ts  - the `TS` context obtained from `TSCreate()`
    1420: . J   - matrix to hold the Jacobian values
    1421: . P   - matrix for constructing the preconditioner (may be same as `J`)
    1422: . jac - the Jacobian evaluation routine, see `TSI2JacobianFn` for the calling sequence
    1423: - ctx - user-defined context for private data for the Jacobian evaluation routine (may be `NULL`)

    1425:   Level: beginner

    1427:   Notes:
    1428:   The matrices `J` and `P` are exactly the matrices that are used by `SNES` for the nonlinear solve.

    1430:   The matrix dF/dU + v*dF/dU_t + a*dF/dU_tt you provide turns out to be
    1431:   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.
    1432:   The time integrator internally approximates U_t by W+v*U and U_tt by W'+a*U  where the positive "shift"
    1433:   parameters 'v' and 'a' and vectors W, W' depend on the integration method, step size, and past states.

    1435: .seealso: [](ch_ts), `TS`, `TSI2JacobianFn`, `TSSetI2Function()`, `TSGetI2Jacobian()`
    1436: @*/
    1437: PetscErrorCode TSSetI2Jacobian(TS ts, Mat J, Mat P, TSI2JacobianFn *jac, void *ctx)
    1438: {
    1439:   DM dm;

    1441:   PetscFunctionBegin;
    1445:   PetscCall(TSSetIJacobian(ts, J, P, NULL, NULL));
    1446:   PetscCall(TSGetDM(ts, &dm));
    1447:   PetscCall(DMTSSetI2Jacobian(dm, jac, ctx));
    1448:   PetscFunctionReturn(PETSC_SUCCESS);
    1449: }

    1451: /*@C
    1452:   TSGetI2Jacobian - Returns the implicit Jacobian at the present timestep.

    1454:   Not Collective, but parallel objects are returned if `TS` is parallel

    1456:   Input Parameter:
    1457: . ts - The `TS` context obtained from `TSCreate()`

    1459:   Output Parameters:
    1460: + J   - The (approximate) Jacobian of F(t,U,U_t,U_tt)
    1461: . P   - The matrix from which the preconditioner is constructed, often the same as `J`
    1462: . jac - The function to compute the Jacobian matrices
    1463: - ctx - User-defined context for Jacobian evaluation routine

    1465:   Level: advanced

    1467:   Note:
    1468:   You can pass in `NULL` for any return argument you do not need.

    1470: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`, `TSSetI2Jacobian()`, `TSGetI2Function()`, `TSCreate()`
    1471: @*/
    1472: PetscErrorCode TSGetI2Jacobian(TS ts, Mat *J, Mat *P, TSI2JacobianFn **jac, void **ctx)
    1473: {
    1474:   SNES snes;
    1475:   DM   dm;

    1477:   PetscFunctionBegin;
    1478:   PetscCall(TSGetSNES(ts, &snes));
    1479:   PetscCall(SNESSetUpMatrices(snes));
    1480:   PetscCall(SNESGetJacobian(snes, J, P, NULL, NULL));
    1481:   PetscCall(TSGetDM(ts, &dm));
    1482:   PetscCall(DMTSGetI2Jacobian(dm, jac, ctx));
    1483:   PetscFunctionReturn(PETSC_SUCCESS);
    1484: }

    1486: /*@
    1487:   TSComputeI2Function - Evaluates the DAE residual written in implicit form F(t,U,U_t,U_tt) = 0

    1489:   Collective

    1491:   Input Parameters:
    1492: + ts - the `TS` context
    1493: . t  - current time
    1494: . U  - state vector
    1495: . V  - time derivative of state vector (U_t)
    1496: - A  - second time derivative of state vector (U_tt)

    1498:   Output Parameter:
    1499: . F - the residual vector

    1501:   Level: developer

    1503:   Note:
    1504:   Most users should not need to explicitly call this routine, as it
    1505:   is used internally within the nonlinear solvers.

    1507: .seealso: [](ch_ts), `TS`, `TSSetI2Function()`, `TSGetI2Function()`
    1508: @*/
    1509: PetscErrorCode TSComputeI2Function(TS ts, PetscReal t, Vec U, Vec V, Vec A, Vec F)
    1510: {
    1511:   DM               dm;
    1512:   TSI2FunctionFn  *I2Function;
    1513:   void            *ctx;
    1514:   TSRHSFunctionFn *rhsfunction;

    1516:   PetscFunctionBegin;

    1523:   PetscCall(TSGetDM(ts, &dm));
    1524:   PetscCall(DMTSGetI2Function(dm, &I2Function, &ctx));
    1525:   PetscCall(DMTSGetRHSFunction(dm, &rhsfunction, NULL));

    1527:   if (!I2Function) {
    1528:     PetscCall(TSComputeIFunction(ts, t, U, A, F, PETSC_FALSE));
    1529:     PetscFunctionReturn(PETSC_SUCCESS);
    1530:   }

    1532:   PetscCall(PetscLogEventBegin(TS_FunctionEval, U, ts, V, F));

    1534:   PetscCallBack("TS callback implicit function", I2Function(ts, t, U, V, A, F, ctx));

    1536:   if (rhsfunction) {
    1537:     Vec Frhs;

    1539:     PetscCall(DMGetGlobalVector(dm, &Frhs));
    1540:     PetscCall(TSComputeRHSFunction(ts, t, U, Frhs));
    1541:     PetscCall(VecAXPY(F, -1, Frhs));
    1542:     PetscCall(DMRestoreGlobalVector(dm, &Frhs));
    1543:   }

    1545:   PetscCall(PetscLogEventEnd(TS_FunctionEval, U, ts, V, F));
    1546:   PetscFunctionReturn(PETSC_SUCCESS);
    1547: }

    1549: /*@
    1550:   TSComputeI2Jacobian - Evaluates the Jacobian of the DAE

    1552:   Collective

    1554:   Input Parameters:
    1555: + ts     - the `TS` context
    1556: . t      - current timestep
    1557: . U      - state vector
    1558: . V      - time derivative of state vector
    1559: . A      - second time derivative of state vector
    1560: . shiftV - shift to apply, see note below
    1561: - shiftA - shift to apply, see note below

    1563:   Output Parameters:
    1564: + J - Jacobian matrix
    1565: - P - optional preconditioning matrix

    1567:   Level: developer

    1569:   Notes:
    1570:   If F(t,U,V,A)=0 is the DAE, the required Jacobian is

    1572:   dF/dU + shiftV*dF/dV + shiftA*dF/dA

    1574:   Most users should not need to explicitly call this routine, as it
    1575:   is used internally within the nonlinear solvers.

    1577: .seealso: [](ch_ts), `TS`, `TSSetI2Jacobian()`
    1578: @*/
    1579: PetscErrorCode TSComputeI2Jacobian(TS ts, PetscReal t, Vec U, Vec V, Vec A, PetscReal shiftV, PetscReal shiftA, Mat J, Mat P)
    1580: {
    1581:   DM               dm;
    1582:   TSI2JacobianFn  *I2Jacobian;
    1583:   void            *ctx;
    1584:   TSRHSJacobianFn *rhsjacobian;

    1586:   PetscFunctionBegin;

    1594:   PetscCall(TSGetDM(ts, &dm));
    1595:   PetscCall(DMTSGetI2Jacobian(dm, &I2Jacobian, &ctx));
    1596:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL));

    1598:   if (!I2Jacobian) {
    1599:     PetscCall(TSComputeIJacobian(ts, t, U, A, shiftA, J, P, PETSC_FALSE));
    1600:     PetscFunctionReturn(PETSC_SUCCESS);
    1601:   }

    1603:   PetscCall(PetscLogEventBegin(TS_JacobianEval, U, ts, J, P));
    1604:   PetscCallBack("TS callback implicit Jacobian", I2Jacobian(ts, t, U, V, A, shiftV, shiftA, J, P, ctx));
    1605:   if (rhsjacobian) {
    1606:     Mat Jrhs, Prhs;
    1607:     PetscCall(TSGetRHSMats_Private(ts, &Jrhs, &Prhs));
    1608:     PetscCall(TSComputeRHSJacobian(ts, t, U, Jrhs, Prhs));
    1609:     PetscCall(MatAXPY(J, -1, Jrhs, ts->axpy_pattern));
    1610:     if (P != J) PetscCall(MatAXPY(P, -1, Prhs, ts->axpy_pattern));
    1611:   }

    1613:   PetscCall(PetscLogEventEnd(TS_JacobianEval, U, ts, J, P));
    1614:   PetscFunctionReturn(PETSC_SUCCESS);
    1615: }

    1617: /*@C
    1618:   TSSetTransientVariable - sets function to transform from state to transient variables

    1620:   Logically Collective

    1622:   Input Parameters:
    1623: + ts   - time stepping context on which to change the transient variable
    1624: . tvar - a function that transforms to transient variables, see `TSTransientVariableFn` for the calling sequence
    1625: - ctx  - a context for tvar

    1627:   Level: advanced

    1629:   Notes:
    1630:   This is typically used to transform from primitive to conservative variables so that a time integrator (e.g., `TSBDF`)
    1631:   can be conservative.  In this context, primitive variables P are used to model the state (e.g., because they lead to
    1632:   well-conditioned formulations even in limiting cases such as low-Mach or zero porosity).  The transient variable is
    1633:   C(P), specified by calling this function.  An IFunction thus receives arguments (P, Cdot) and the IJacobian must be
    1634:   evaluated via the chain rule, as in
    1635: .vb
    1636:      dF/dP + shift * dF/dCdot dC/dP.
    1637: .ve

    1639: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSTransientVariableFn`, `DMTSSetTransientVariable()`, `DMTSGetTransientVariable()`, `TSSetIFunction()`, `TSSetIJacobian()`
    1640: @*/
    1641: PetscErrorCode TSSetTransientVariable(TS ts, TSTransientVariableFn *tvar, void *ctx)
    1642: {
    1643:   DM dm;

    1645:   PetscFunctionBegin;
    1647:   PetscCall(TSGetDM(ts, &dm));
    1648:   PetscCall(DMTSSetTransientVariable(dm, tvar, ctx));
    1649:   PetscFunctionReturn(PETSC_SUCCESS);
    1650: }

    1652: /*@
    1653:   TSComputeTransientVariable - transforms state (primitive) variables to transient (conservative) variables

    1655:   Logically Collective

    1657:   Input Parameters:
    1658: + ts - TS on which to compute
    1659: - U  - state vector to be transformed to transient variables

    1661:   Output Parameter:
    1662: . C - transient (conservative) variable

    1664:   Level: developer

    1666:   Developer Notes:
    1667:   If `DMTSSetTransientVariable()` has not been called, then C is not modified in this routine and C = `NULL` is allowed.
    1668:   This makes it safe to call without a guard.  One can use `TSHasTransientVariable()` to check if transient variables are
    1669:   being used.

    1671: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeIFunction()`, `TSComputeIJacobian()`
    1672: @*/
    1673: PetscErrorCode TSComputeTransientVariable(TS ts, Vec U, Vec C)
    1674: {
    1675:   DM   dm;
    1676:   DMTS dmts;

    1678:   PetscFunctionBegin;
    1681:   PetscCall(TSGetDM(ts, &dm));
    1682:   PetscCall(DMGetDMTS(dm, &dmts));
    1683:   if (dmts->ops->transientvar) {
    1685:     PetscCall((*dmts->ops->transientvar)(ts, U, C, dmts->transientvarctx));
    1686:   }
    1687:   PetscFunctionReturn(PETSC_SUCCESS);
    1688: }

    1690: /*@
    1691:   TSHasTransientVariable - determine whether transient variables have been set

    1693:   Logically Collective

    1695:   Input Parameter:
    1696: . ts - `TS` on which to compute

    1698:   Output Parameter:
    1699: . has - `PETSC_TRUE` if transient variables have been set

    1701:   Level: developer

    1703: .seealso: [](ch_ts), `TS`, `TSBDF`, `DMTSSetTransientVariable()`, `TSComputeTransientVariable()`
    1704: @*/
    1705: PetscErrorCode TSHasTransientVariable(TS ts, PetscBool *has)
    1706: {
    1707:   DM   dm;
    1708:   DMTS dmts;

    1710:   PetscFunctionBegin;
    1712:   PetscCall(TSGetDM(ts, &dm));
    1713:   PetscCall(DMGetDMTS(dm, &dmts));
    1714:   *has = dmts->ops->transientvar ? PETSC_TRUE : PETSC_FALSE;
    1715:   PetscFunctionReturn(PETSC_SUCCESS);
    1716: }

    1718: /*@
    1719:   TS2SetSolution - Sets the initial solution and time derivative vectors
    1720:   for use by the `TS` routines handling second order equations.

    1722:   Logically Collective

    1724:   Input Parameters:
    1725: + ts - the `TS` context obtained from `TSCreate()`
    1726: . u  - the solution vector
    1727: - v  - the time derivative vector

    1729:   Level: beginner

    1731: .seealso: [](ch_ts), `TS`
    1732: @*/
    1733: PetscErrorCode TS2SetSolution(TS ts, Vec u, Vec v)
    1734: {
    1735:   PetscFunctionBegin;
    1739:   PetscCall(TSSetSolution(ts, u));
    1740:   PetscCall(PetscObjectReference((PetscObject)v));
    1741:   PetscCall(VecDestroy(&ts->vec_dot));
    1742:   ts->vec_dot = v;
    1743:   PetscFunctionReturn(PETSC_SUCCESS);
    1744: }

    1746: /*@
    1747:   TS2GetSolution - Returns the solution and time derivative at the present timestep
    1748:   for second order equations.

    1750:   Not Collective

    1752:   Input Parameter:
    1753: . ts - the `TS` context obtained from `TSCreate()`

    1755:   Output Parameters:
    1756: + u - the vector containing the solution
    1757: - v - the vector containing the time derivative

    1759:   Level: intermediate

    1761:   Notes:
    1762:   It is valid to call this routine inside the function
    1763:   that you are evaluating in order to move to the new timestep. This vector not
    1764:   changed until the solution at the next timestep has been calculated.

    1766: .seealso: [](ch_ts), `TS`, `TS2SetSolution()`, `TSGetTimeStep()`, `TSGetTime()`
    1767: @*/
    1768: PetscErrorCode TS2GetSolution(TS ts, Vec *u, Vec *v)
    1769: {
    1770:   PetscFunctionBegin;
    1772:   if (u) PetscAssertPointer(u, 2);
    1773:   if (v) PetscAssertPointer(v, 3);
    1774:   if (u) *u = ts->vec_sol;
    1775:   if (v) *v = ts->vec_dot;
    1776:   PetscFunctionReturn(PETSC_SUCCESS);
    1777: }

    1779: /*@
    1780:   TSLoad - Loads a `TS` that has been stored in binary  with `TSView()`.

    1782:   Collective

    1784:   Input Parameters:
    1785: + ts     - the newly loaded `TS`, this needs to have been created with `TSCreate()` or
    1786:            some related function before a call to `TSLoad()`.
    1787: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`

    1789:   Level: intermediate

    1791:   Note:
    1792:   The type is determined by the data in the file, any type set into the `TS` before this call is ignored.

    1794: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerBinaryOpen()`, `TSView()`, `MatLoad()`, `VecLoad()`
    1795: @*/
    1796: PetscErrorCode TSLoad(TS ts, PetscViewer viewer)
    1797: {
    1798:   PetscBool isbinary;
    1799:   PetscInt  classid;
    1800:   char      type[256];
    1801:   DMTS      sdm;
    1802:   DM        dm;

    1804:   PetscFunctionBegin;
    1807:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
    1808:   PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");

    1810:   PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
    1811:   PetscCheck(classid == TS_FILE_CLASSID, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Not TS next in file");
    1812:   PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
    1813:   PetscCall(TSSetType(ts, type));
    1814:   PetscTryTypeMethod(ts, load, viewer);
    1815:   PetscCall(DMCreate(PetscObjectComm((PetscObject)ts), &dm));
    1816:   PetscCall(DMLoad(dm, viewer));
    1817:   PetscCall(TSSetDM(ts, dm));
    1818:   PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
    1819:   PetscCall(VecLoad(ts->vec_sol, viewer));
    1820:   PetscCall(DMGetDMTS(ts->dm, &sdm));
    1821:   PetscCall(DMTSLoad(sdm, viewer));
    1822:   PetscFunctionReturn(PETSC_SUCCESS);
    1823: }

    1825: #include <petscdraw.h>
    1826: #if defined(PETSC_HAVE_SAWS)
    1827: #include <petscviewersaws.h>
    1828: #endif

    1830: /*@
    1831:   TSViewFromOptions - View a `TS` based on values in the options database

    1833:   Collective

    1835:   Input Parameters:
    1836: + ts   - the `TS` context
    1837: . obj  - Optional object that provides the prefix for the options database keys
    1838: - name - command line option string to be passed by user

    1840:   Level: intermediate

    1842: .seealso: [](ch_ts), `TS`, `TSView`, `PetscObjectViewFromOptions()`, `TSCreate()`
    1843: @*/
    1844: PetscErrorCode TSViewFromOptions(TS ts, PetscObject obj, const char name[])
    1845: {
    1846:   PetscFunctionBegin;
    1848:   PetscCall(PetscObjectViewFromOptions((PetscObject)ts, obj, name));
    1849:   PetscFunctionReturn(PETSC_SUCCESS);
    1850: }

    1852: /*@
    1853:   TSView - Prints the `TS` data structure.

    1855:   Collective

    1857:   Input Parameters:
    1858: + ts     - the `TS` context obtained from `TSCreate()`
    1859: - viewer - visualization context

    1861:   Options Database Key:
    1862: . -ts_view - calls `TSView()` at end of `TSStep()`

    1864:   Level: beginner

    1866:   Notes:
    1867:   The available visualization contexts include
    1868: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
    1869: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
    1870:   output where only the first processor opens
    1871:   the file.  All other processors send their
    1872:   data to the first processor to print.

    1874:   The user can open an alternative visualization context with
    1875:   `PetscViewerASCIIOpen()` - output to a specified file.

    1877:   In the debugger you can do call `TSView`(ts,0) to display the `TS` solver. (The same holds for any PETSc object viewer).

    1879: .seealso: [](ch_ts), `TS`, `PetscViewer`, `PetscViewerASCIIOpen()`
    1880: @*/
    1881: PetscErrorCode TSView(TS ts, PetscViewer viewer)
    1882: {
    1883:   TSType    type;
    1884:   PetscBool iascii, isstring, isundials, isbinary, isdraw;
    1885:   DMTS      sdm;
    1886: #if defined(PETSC_HAVE_SAWS)
    1887:   PetscBool issaws;
    1888: #endif

    1890:   PetscFunctionBegin;
    1892:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts), &viewer));
    1894:   PetscCheckSameComm(ts, 1, viewer, 2);

    1896:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    1897:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
    1898:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
    1899:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    1900: #if defined(PETSC_HAVE_SAWS)
    1901:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
    1902: #endif
    1903:   if (iascii) {
    1904:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ts, viewer));
    1905:     if (ts->ops->view) {
    1906:       PetscCall(PetscViewerASCIIPushTab(viewer));
    1907:       PetscUseTypeMethod(ts, view, viewer);
    1908:       PetscCall(PetscViewerASCIIPopTab(viewer));
    1909:     }
    1910:     if (ts->max_steps < PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, "  maximum steps=%" PetscInt_FMT "\n", ts->max_steps));
    1911:     if (ts->max_time < PETSC_MAX_REAL) PetscCall(PetscViewerASCIIPrintf(viewer, "  maximum time=%g\n", (double)ts->max_time));
    1912:     if (ts->ifuncs) PetscCall(PetscViewerASCIIPrintf(viewer, "  total number of I function evaluations=%" PetscInt_FMT "\n", ts->ifuncs));
    1913:     if (ts->ijacs) PetscCall(PetscViewerASCIIPrintf(viewer, "  total number of I Jacobian evaluations=%" PetscInt_FMT "\n", ts->ijacs));
    1914:     if (ts->rhsfuncs) PetscCall(PetscViewerASCIIPrintf(viewer, "  total number of RHS function evaluations=%" PetscInt_FMT "\n", ts->rhsfuncs));
    1915:     if (ts->rhsjacs) PetscCall(PetscViewerASCIIPrintf(viewer, "  total number of RHS Jacobian evaluations=%" PetscInt_FMT "\n", ts->rhsjacs));
    1916:     if (ts->usessnes) {
    1917:       PetscBool lin;
    1918:       if (ts->problem_type == TS_NONLINEAR) PetscCall(PetscViewerASCIIPrintf(viewer, "  total number of nonlinear solver iterations=%" PetscInt_FMT "\n", ts->snes_its));
    1919:       PetscCall(PetscViewerASCIIPrintf(viewer, "  total number of linear solver iterations=%" PetscInt_FMT "\n", ts->ksp_its));
    1920:       PetscCall(PetscObjectTypeCompareAny((PetscObject)ts->snes, &lin, SNESKSPONLY, SNESKSPTRANSPOSEONLY, ""));
    1921:       PetscCall(PetscViewerASCIIPrintf(viewer, "  total number of %slinear solve failures=%" PetscInt_FMT "\n", lin ? "" : "non", ts->num_snes_failures));
    1922:     }
    1923:     PetscCall(PetscViewerASCIIPrintf(viewer, "  total number of rejected steps=%" PetscInt_FMT "\n", ts->reject));
    1924:     if (ts->vrtol) PetscCall(PetscViewerASCIIPrintf(viewer, "  using vector of relative error tolerances, "));
    1925:     else PetscCall(PetscViewerASCIIPrintf(viewer, "  using relative error tolerance of %g, ", (double)ts->rtol));
    1926:     if (ts->vatol) PetscCall(PetscViewerASCIIPrintf(viewer, "  using vector of absolute error tolerances\n"));
    1927:     else PetscCall(PetscViewerASCIIPrintf(viewer, "  using absolute error tolerance of %g\n", (double)ts->atol));
    1928:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1929:     PetscCall(TSAdaptView(ts->adapt, viewer));
    1930:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1931:   } else if (isstring) {
    1932:     PetscCall(TSGetType(ts, &type));
    1933:     PetscCall(PetscViewerStringSPrintf(viewer, " TSType: %-7.7s", type));
    1934:     PetscTryTypeMethod(ts, view, viewer);
    1935:   } else if (isbinary) {
    1936:     PetscInt    classid = TS_FILE_CLASSID;
    1937:     MPI_Comm    comm;
    1938:     PetscMPIInt rank;
    1939:     char        type[256];

    1941:     PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
    1942:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1943:     if (rank == 0) {
    1944:       PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
    1945:       PetscCall(PetscStrncpy(type, ((PetscObject)ts)->type_name, 256));
    1946:       PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
    1947:     }
    1948:     PetscTryTypeMethod(ts, view, viewer);
    1949:     if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
    1950:     PetscCall(DMView(ts->dm, viewer));
    1951:     PetscCall(VecView(ts->vec_sol, viewer));
    1952:     PetscCall(DMGetDMTS(ts->dm, &sdm));
    1953:     PetscCall(DMTSView(sdm, viewer));
    1954:   } else if (isdraw) {
    1955:     PetscDraw draw;
    1956:     char      str[36];
    1957:     PetscReal x, y, bottom, h;

    1959:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    1960:     PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
    1961:     PetscCall(PetscStrncpy(str, "TS: ", sizeof(str)));
    1962:     PetscCall(PetscStrlcat(str, ((PetscObject)ts)->type_name, sizeof(str)));
    1963:     PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLACK, PETSC_DRAW_BLACK, str, NULL, &h));
    1964:     bottom = y - h;
    1965:     PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
    1966:     PetscTryTypeMethod(ts, view, viewer);
    1967:     if (ts->adapt) PetscCall(TSAdaptView(ts->adapt, viewer));
    1968:     if (ts->snes) PetscCall(SNESView(ts->snes, viewer));
    1969:     PetscCall(PetscDrawPopCurrentPoint(draw));
    1970: #if defined(PETSC_HAVE_SAWS)
    1971:   } else if (issaws) {
    1972:     PetscMPIInt rank;
    1973:     const char *name;

    1975:     PetscCall(PetscObjectGetName((PetscObject)ts, &name));
    1976:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    1977:     if (!((PetscObject)ts)->amsmem && rank == 0) {
    1978:       char dir[1024];

    1980:       PetscCall(PetscObjectViewSAWs((PetscObject)ts, viewer));
    1981:       PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time_step", name));
    1982:       PetscCallSAWs(SAWs_Register, (dir, &ts->steps, 1, SAWs_READ, SAWs_INT));
    1983:       PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/time", name));
    1984:       PetscCallSAWs(SAWs_Register, (dir, &ts->ptime, 1, SAWs_READ, SAWs_DOUBLE));
    1985:     }
    1986:     PetscTryTypeMethod(ts, view, viewer);
    1987: #endif
    1988:   }
    1989:   if (ts->snes && ts->usessnes) {
    1990:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1991:     PetscCall(SNESView(ts->snes, viewer));
    1992:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1993:   }
    1994:   PetscCall(DMGetDMTS(ts->dm, &sdm));
    1995:   PetscCall(DMTSView(sdm, viewer));

    1997:   PetscCall(PetscViewerASCIIPushTab(viewer));
    1998:   PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSSUNDIALS, &isundials));
    1999:   PetscCall(PetscViewerASCIIPopTab(viewer));
    2000:   PetscFunctionReturn(PETSC_SUCCESS);
    2001: }

    2003: /*@
    2004:   TSSetApplicationContext - Sets an optional user-defined context for
    2005:   the timesteppers.

    2007:   Logically Collective

    2009:   Input Parameters:
    2010: + ts   - the `TS` context obtained from `TSCreate()`
    2011: - usrP - user context

    2013:   Level: intermediate

    2015:   Fortran Notes:
    2016:   You must write a Fortran interface definition for this
    2017:   function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.

    2019: .seealso: [](ch_ts), `TS`, `TSGetApplicationContext()`
    2020: @*/
    2021: PetscErrorCode TSSetApplicationContext(TS ts, void *usrP)
    2022: {
    2023:   PetscFunctionBegin;
    2025:   ts->user = usrP;
    2026:   PetscFunctionReturn(PETSC_SUCCESS);
    2027: }

    2029: /*@
    2030:   TSGetApplicationContext - Gets the user-defined context for the
    2031:   timestepper that was set with `TSSetApplicationContext()`

    2033:   Not Collective

    2035:   Input Parameter:
    2036: . ts - the `TS` context obtained from `TSCreate()`

    2038:   Output Parameter:
    2039: . usrP - user context

    2041:   Level: intermediate

    2043:   Fortran Notes:
    2044:   You must write a Fortran interface definition for this
    2045:   function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.

    2047: .seealso: [](ch_ts), `TS`, `TSSetApplicationContext()`
    2048: @*/
    2049: PetscErrorCode TSGetApplicationContext(TS ts, void *usrP)
    2050: {
    2051:   PetscFunctionBegin;
    2053:   *(void **)usrP = ts->user;
    2054:   PetscFunctionReturn(PETSC_SUCCESS);
    2055: }

    2057: /*@
    2058:   TSGetStepNumber - Gets the number of time steps completed.

    2060:   Not Collective

    2062:   Input Parameter:
    2063: . ts - the `TS` context obtained from `TSCreate()`

    2065:   Output Parameter:
    2066: . steps - number of steps completed so far

    2068:   Level: intermediate

    2070: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetTimeStep()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`
    2071: @*/
    2072: PetscErrorCode TSGetStepNumber(TS ts, PetscInt *steps)
    2073: {
    2074:   PetscFunctionBegin;
    2076:   PetscAssertPointer(steps, 2);
    2077:   *steps = ts->steps;
    2078:   PetscFunctionReturn(PETSC_SUCCESS);
    2079: }

    2081: /*@
    2082:   TSSetStepNumber - Sets the number of steps completed.

    2084:   Logically Collective

    2086:   Input Parameters:
    2087: + ts    - the `TS` context
    2088: - steps - number of steps completed so far

    2090:   Level: developer

    2092:   Note:
    2093:   For most uses of the `TS` solvers the user need not explicitly call
    2094:   `TSSetStepNumber()`, as the step counter is appropriately updated in
    2095:   `TSSolve()`/`TSStep()`/`TSRollBack()`. Power users may call this routine to
    2096:   reinitialize timestepping by setting the step counter to zero (and time
    2097:   to the initial time) to solve a similar problem with different initial
    2098:   conditions or parameters. Other possible use case is to continue
    2099:   timestepping from a previously interrupted run in such a way that `TS`
    2100:   monitors will be called with a initial nonzero step counter.

    2102: .seealso: [](ch_ts), `TS`, `TSGetStepNumber()`, `TSSetTime()`, `TSSetTimeStep()`, `TSSetSolution()`
    2103: @*/
    2104: PetscErrorCode TSSetStepNumber(TS ts, PetscInt steps)
    2105: {
    2106:   PetscFunctionBegin;
    2109:   PetscCheck(steps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Step number must be non-negative");
    2110:   ts->steps = steps;
    2111:   PetscFunctionReturn(PETSC_SUCCESS);
    2112: }

    2114: /*@
    2115:   TSSetTimeStep - Allows one to reset the timestep at any time,
    2116:   useful for simple pseudo-timestepping codes.

    2118:   Logically Collective

    2120:   Input Parameters:
    2121: + ts        - the `TS` context obtained from `TSCreate()`
    2122: - time_step - the size of the timestep

    2124:   Level: intermediate

    2126: .seealso: [](ch_ts), `TS`, `TSPSEUDO`, `TSGetTimeStep()`, `TSSetTime()`
    2127: @*/
    2128: PetscErrorCode TSSetTimeStep(TS ts, PetscReal time_step)
    2129: {
    2130:   PetscFunctionBegin;
    2133:   ts->time_step = time_step;
    2134:   PetscFunctionReturn(PETSC_SUCCESS);
    2135: }

    2137: /*@
    2138:   TSSetExactFinalTime - Determines whether to adapt the final time step to
    2139:   match the exact final time, to interpolate the solution to the exact final time,
    2140:   or to just return at the final time `TS` computed (which may be slightly larger
    2141:   than the requested final time).

    2143:   Logically Collective

    2145:   Input Parameters:
    2146: + ts     - the time-step context
    2147: - eftopt - exact final time option
    2148: .vb
    2149:   TS_EXACTFINALTIME_STEPOVER    - Don't do anything if final time is exceeded, just use it
    2150:   TS_EXACTFINALTIME_INTERPOLATE - Interpolate back to final time if the final time is exceeded
    2151:   TS_EXACTFINALTIME_MATCHSTEP   - Adapt final time step to ensure the computed final time exactly equals the requested final time
    2152: .ve

    2154:   Options Database Key:
    2155: . -ts_exact_final_time <stepover,interpolate,matchstep> - select the final step approach at runtime

    2157:   Level: beginner

    2159:   Note:
    2160:   If you use the option `TS_EXACTFINALTIME_STEPOVER` the solution may be at a very different time
    2161:   then the final time you selected.

    2163: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSGetExactFinalTime()`
    2164: @*/
    2165: PetscErrorCode TSSetExactFinalTime(TS ts, TSExactFinalTimeOption eftopt)
    2166: {
    2167:   PetscFunctionBegin;
    2170:   ts->exact_final_time = eftopt;
    2171:   PetscFunctionReturn(PETSC_SUCCESS);
    2172: }

    2174: /*@
    2175:   TSGetExactFinalTime - Gets the exact final time option set with `TSSetExactFinalTime()`

    2177:   Not Collective

    2179:   Input Parameter:
    2180: . ts - the `TS` context

    2182:   Output Parameter:
    2183: . eftopt - exact final time option

    2185:   Level: beginner

    2187: .seealso: [](ch_ts), `TS`, `TSExactFinalTimeOption`, `TSSetExactFinalTime()`
    2188: @*/
    2189: PetscErrorCode TSGetExactFinalTime(TS ts, TSExactFinalTimeOption *eftopt)
    2190: {
    2191:   PetscFunctionBegin;
    2193:   PetscAssertPointer(eftopt, 2);
    2194:   *eftopt = ts->exact_final_time;
    2195:   PetscFunctionReturn(PETSC_SUCCESS);
    2196: }

    2198: /*@
    2199:   TSGetTimeStep - Gets the current timestep size.

    2201:   Not Collective

    2203:   Input Parameter:
    2204: . ts - the `TS` context obtained from `TSCreate()`

    2206:   Output Parameter:
    2207: . dt - the current timestep size

    2209:   Level: intermediate

    2211: .seealso: [](ch_ts), `TS`, `TSSetTimeStep()`, `TSGetTime()`
    2212: @*/
    2213: PetscErrorCode TSGetTimeStep(TS ts, PetscReal *dt)
    2214: {
    2215:   PetscFunctionBegin;
    2217:   PetscAssertPointer(dt, 2);
    2218:   *dt = ts->time_step;
    2219:   PetscFunctionReturn(PETSC_SUCCESS);
    2220: }

    2222: /*@
    2223:   TSGetSolution - Returns the solution at the present timestep. It
    2224:   is valid to call this routine inside the function that you are evaluating
    2225:   in order to move to the new timestep. This vector not changed until
    2226:   the solution at the next timestep has been calculated.

    2228:   Not Collective, but v returned is parallel if ts is parallel

    2230:   Input Parameter:
    2231: . ts - the `TS` context obtained from `TSCreate()`

    2233:   Output Parameter:
    2234: . v - the vector containing the solution

    2236:   Level: intermediate

    2238:   Note:
    2239:   If you used `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`); this does not return the solution at the requested
    2240:   final time. It returns the solution at the next timestep.

    2242: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetSolutionComponents()`, `TSSetSolutionFunction()`
    2243: @*/
    2244: PetscErrorCode TSGetSolution(TS ts, Vec *v)
    2245: {
    2246:   PetscFunctionBegin;
    2248:   PetscAssertPointer(v, 2);
    2249:   *v = ts->vec_sol;
    2250:   PetscFunctionReturn(PETSC_SUCCESS);
    2251: }

    2253: /*@
    2254:   TSGetSolutionComponents - Returns any solution components at the present
    2255:   timestep, if available for the time integration method being used.
    2256:   Solution components are quantities that share the same size and
    2257:   structure as the solution vector.

    2259:   Not Collective, but v returned is parallel if ts is parallel

    2261:   Input Parameters:
    2262: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
    2263: . n  - If v is `NULL`, then the number of solution components is
    2264:        returned through n, else the n-th solution component is
    2265:        returned in v.
    2266: - v  - the vector containing the n-th solution component
    2267:        (may be `NULL` to use this function to find out
    2268:         the number of solutions components).

    2270:   Level: advanced

    2272: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
    2273: @*/
    2274: PetscErrorCode TSGetSolutionComponents(TS ts, PetscInt *n, Vec *v)
    2275: {
    2276:   PetscFunctionBegin;
    2278:   if (!ts->ops->getsolutioncomponents) *n = 0;
    2279:   else PetscUseTypeMethod(ts, getsolutioncomponents, n, v);
    2280:   PetscFunctionReturn(PETSC_SUCCESS);
    2281: }

    2283: /*@
    2284:   TSGetAuxSolution - Returns an auxiliary solution at the present
    2285:   timestep, if available for the time integration method being used.

    2287:   Not Collective, but v returned is parallel if ts is parallel

    2289:   Input Parameters:
    2290: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
    2291: - v  - the vector containing the auxiliary solution

    2293:   Level: intermediate

    2295: .seealso: [](ch_ts), `TS`, `TSGetSolution()`
    2296: @*/
    2297: PetscErrorCode TSGetAuxSolution(TS ts, Vec *v)
    2298: {
    2299:   PetscFunctionBegin;
    2301:   if (ts->ops->getauxsolution) PetscUseTypeMethod(ts, getauxsolution, v);
    2302:   else PetscCall(VecZeroEntries(*v));
    2303:   PetscFunctionReturn(PETSC_SUCCESS);
    2304: }

    2306: /*@
    2307:   TSGetTimeError - Returns the estimated error vector, if the chosen
    2308:   `TSType` has an error estimation functionality and `TSSetTimeError()` was called

    2310:   Not Collective, but v returned is parallel if ts is parallel

    2312:   Input Parameters:
    2313: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
    2314: . n  - current estimate (n=0) or previous one (n=-1)
    2315: - v  - the vector containing the error (same size as the solution).

    2317:   Level: intermediate

    2319:   Note:
    2320:   MUST call after `TSSetUp()`

    2322: .seealso: [](ch_ts), `TSGetSolution()`, `TSSetTimeError()`
    2323: @*/
    2324: PetscErrorCode TSGetTimeError(TS ts, PetscInt n, Vec *v)
    2325: {
    2326:   PetscFunctionBegin;
    2328:   if (ts->ops->gettimeerror) PetscUseTypeMethod(ts, gettimeerror, n, v);
    2329:   else PetscCall(VecZeroEntries(*v));
    2330:   PetscFunctionReturn(PETSC_SUCCESS);
    2331: }

    2333: /*@
    2334:   TSSetTimeError - Sets the estimated error vector, if the chosen
    2335:   `TSType` has an error estimation functionality. This can be used
    2336:   to restart such a time integrator with a given error vector.

    2338:   Not Collective, but v returned is parallel if ts is parallel

    2340:   Input Parameters:
    2341: + ts - the `TS` context obtained from `TSCreate()` (input parameter).
    2342: - v  - the vector containing the error (same size as the solution).

    2344:   Level: intermediate

    2346: .seealso: [](ch_ts), `TS`, `TSSetSolution()`, `TSGetTimeError()`
    2347: @*/
    2348: PetscErrorCode TSSetTimeError(TS ts, Vec v)
    2349: {
    2350:   PetscFunctionBegin;
    2352:   PetscCheck(ts->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetUp() first");
    2353:   PetscTryTypeMethod(ts, settimeerror, v);
    2354:   PetscFunctionReturn(PETSC_SUCCESS);
    2355: }

    2357: /* ----- Routines to initialize and destroy a timestepper ---- */
    2358: /*@
    2359:   TSSetProblemType - Sets the type of problem to be solved.

    2361:   Not collective

    2363:   Input Parameters:
    2364: + ts   - The `TS`
    2365: - type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
    2366: .vb
    2367:          U_t - A U = 0      (linear)
    2368:          U_t - A(t) U = 0   (linear)
    2369:          F(t,U,U_t) = 0     (nonlinear)
    2370: .ve

    2372:   Level: beginner

    2374: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
    2375: @*/
    2376: PetscErrorCode TSSetProblemType(TS ts, TSProblemType type)
    2377: {
    2378:   PetscFunctionBegin;
    2380:   ts->problem_type = type;
    2381:   if (type == TS_LINEAR) {
    2382:     SNES snes;
    2383:     PetscCall(TSGetSNES(ts, &snes));
    2384:     PetscCall(SNESSetType(snes, SNESKSPONLY));
    2385:   }
    2386:   PetscFunctionReturn(PETSC_SUCCESS);
    2387: }

    2389: /*@
    2390:   TSGetProblemType - Gets the type of problem to be solved.

    2392:   Not collective

    2394:   Input Parameter:
    2395: . ts - The `TS`

    2397:   Output Parameter:
    2398: . type - One of `TS_LINEAR`, `TS_NONLINEAR` where these types refer to problems of the forms
    2399: .vb
    2400:          M U_t = A U
    2401:          M(t) U_t = A(t) U
    2402:          F(t,U,U_t)
    2403: .ve

    2405:   Level: beginner

    2407: .seealso: [](ch_ts), `TSSetUp()`, `TSProblemType`, `TS`
    2408: @*/
    2409: PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type)
    2410: {
    2411:   PetscFunctionBegin;
    2413:   PetscAssertPointer(type, 2);
    2414:   *type = ts->problem_type;
    2415:   PetscFunctionReturn(PETSC_SUCCESS);
    2416: }

    2418: /*
    2419:     Attempt to check/preset a default value for the exact final time option. This is needed at the beginning of TSSolve() and in TSSetUp()
    2420: */
    2421: static PetscErrorCode TSSetExactFinalTimeDefault(TS ts)
    2422: {
    2423:   PetscBool isnone;

    2425:   PetscFunctionBegin;
    2426:   PetscCall(TSGetAdapt(ts, &ts->adapt));
    2427:   PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));

    2429:   PetscCall(PetscObjectTypeCompare((PetscObject)ts->adapt, TSADAPTNONE, &isnone));
    2430:   if (!isnone && ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP;
    2431:   else if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) ts->exact_final_time = TS_EXACTFINALTIME_INTERPOLATE;
    2432:   PetscFunctionReturn(PETSC_SUCCESS);
    2433: }

    2435: /*@
    2436:   TSSetUp - Sets up the internal data structures for the later use of a timestepper.

    2438:   Collective

    2440:   Input Parameter:
    2441: . ts - the `TS` context obtained from `TSCreate()`

    2443:   Level: advanced

    2445:   Note:
    2446:   For basic use of the `TS` solvers the user need not explicitly call
    2447:   `TSSetUp()`, since these actions will automatically occur during
    2448:   the call to `TSStep()` or `TSSolve()`.  However, if one wishes to control this
    2449:   phase separately, `TSSetUp()` should be called after `TSCreate()`
    2450:   and optional routines of the form TSSetXXX(), but before `TSStep()` and `TSSolve()`.

    2452: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSStep()`, `TSDestroy()`, `TSSolve()`
    2453: @*/
    2454: PetscErrorCode TSSetUp(TS ts)
    2455: {
    2456:   DM dm;
    2457:   PetscErrorCode (*func)(SNES, Vec, Vec, void *);
    2458:   PetscErrorCode (*jac)(SNES, Vec, Mat, Mat, void *);
    2459:   TSIFunctionFn   *ifun;
    2460:   TSIJacobianFn   *ijac;
    2461:   TSI2JacobianFn  *i2jac;
    2462:   TSRHSJacobianFn *rhsjac;

    2464:   PetscFunctionBegin;
    2466:   if (ts->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);

    2468:   if (!((PetscObject)ts)->type_name) {
    2469:     PetscCall(TSGetIFunction(ts, NULL, &ifun, NULL));
    2470:     PetscCall(TSSetType(ts, ifun ? TSBEULER : TSEULER));
    2471:   }

    2473:   if (!ts->vec_sol) {
    2474:     PetscCheck(ts->dm, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call TSSetSolution() first");
    2475:     PetscCall(DMCreateGlobalVector(ts->dm, &ts->vec_sol));
    2476:   }

    2478:   if (ts->tspan) {
    2479:     if (!ts->tspan->vecs_sol) PetscCall(VecDuplicateVecs(ts->vec_sol, ts->tspan->num_span_times, &ts->tspan->vecs_sol));
    2480:   }
    2481:   if (!ts->Jacp && ts->Jacprhs) { /* IJacobianP shares the same matrix with RHSJacobianP if only RHSJacobianP is provided */
    2482:     PetscCall(PetscObjectReference((PetscObject)ts->Jacprhs));
    2483:     ts->Jacp = ts->Jacprhs;
    2484:   }

    2486:   if (ts->quadraturets) {
    2487:     PetscCall(TSSetUp(ts->quadraturets));
    2488:     PetscCall(VecDestroy(&ts->vec_costintegrand));
    2489:     PetscCall(VecDuplicate(ts->quadraturets->vec_sol, &ts->vec_costintegrand));
    2490:   }

    2492:   PetscCall(TSGetRHSJacobian(ts, NULL, NULL, &rhsjac, NULL));
    2493:   if (rhsjac == TSComputeRHSJacobianConstant) {
    2494:     Mat  Amat, Pmat;
    2495:     SNES snes;
    2496:     PetscCall(TSGetSNES(ts, &snes));
    2497:     PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL));
    2498:     /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would
    2499:      * have displaced the RHS matrix */
    2500:     if (Amat && Amat == ts->Arhs) {
    2501:       /* 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 */
    2502:       PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat));
    2503:       PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL));
    2504:       PetscCall(MatDestroy(&Amat));
    2505:     }
    2506:     if (Pmat && Pmat == ts->Brhs) {
    2507:       PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat));
    2508:       PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL));
    2509:       PetscCall(MatDestroy(&Pmat));
    2510:     }
    2511:   }

    2513:   PetscCall(TSGetAdapt(ts, &ts->adapt));
    2514:   PetscCall(TSAdaptSetDefaultType(ts->adapt, ts->default_adapt_type));

    2516:   PetscTryTypeMethod(ts, setup);

    2518:   PetscCall(TSSetExactFinalTimeDefault(ts));

    2520:   /* In the case where we've set a DMTSFunction or what have you, we need the default SNESFunction
    2521:      to be set right but can't do it elsewhere due to the overreliance on ctx=ts.
    2522:    */
    2523:   PetscCall(TSGetDM(ts, &dm));
    2524:   PetscCall(DMSNESGetFunction(dm, &func, NULL));
    2525:   if (!func) PetscCall(DMSNESSetFunction(dm, SNESTSFormFunction, ts));

    2527:   /* If the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it.
    2528:      Otherwise, the SNES will use coloring internally to form the Jacobian.
    2529:    */
    2530:   PetscCall(DMSNESGetJacobian(dm, &jac, NULL));
    2531:   PetscCall(DMTSGetIJacobian(dm, &ijac, NULL));
    2532:   PetscCall(DMTSGetI2Jacobian(dm, &i2jac, NULL));
    2533:   PetscCall(DMTSGetRHSJacobian(dm, &rhsjac, NULL));
    2534:   if (!jac && (ijac || i2jac || rhsjac)) PetscCall(DMSNESSetJacobian(dm, SNESTSFormJacobian, ts));

    2536:   /* if time integration scheme has a starting method, call it */
    2537:   PetscTryTypeMethod(ts, startingmethod);

    2539:   ts->setupcalled = PETSC_TRUE;
    2540:   PetscFunctionReturn(PETSC_SUCCESS);
    2541: }

    2543: /*@
    2544:   TSReset - Resets a `TS` context to the state it was in before `TSSetUp()` was called and removes any allocated `Vec` and `Mat` from its data structures

    2546:   Collective

    2548:   Input Parameter:
    2549: . ts - the `TS` context obtained from `TSCreate()`

    2551:   Level: developer

    2553:   Notes:
    2554:   Any options set on the `TS` object, including those set with `TSSetFromOptions()` remain.

    2556:   See also `TSSetResize()` to change the size of the system being integrated (for example by adaptive mesh refinement) during the time integration.

    2558: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetup()`, `TSDestroy()`, `TSSetResize()`
    2559: @*/
    2560: PetscErrorCode TSReset(TS ts)
    2561: {
    2562:   TS_RHSSplitLink ilink = ts->tsrhssplit, next;

    2564:   PetscFunctionBegin;

    2567:   PetscTryTypeMethod(ts, reset);
    2568:   if (ts->snes) PetscCall(SNESReset(ts->snes));
    2569:   if (ts->adapt) PetscCall(TSAdaptReset(ts->adapt));

    2571:   PetscCall(MatDestroy(&ts->Arhs));
    2572:   PetscCall(MatDestroy(&ts->Brhs));
    2573:   PetscCall(VecDestroy(&ts->Frhs));
    2574:   PetscCall(VecDestroy(&ts->vec_sol));
    2575:   PetscCall(VecDestroy(&ts->vec_sol0));
    2576:   PetscCall(VecDestroy(&ts->vec_dot));
    2577:   PetscCall(VecDestroy(&ts->vatol));
    2578:   PetscCall(VecDestroy(&ts->vrtol));
    2579:   PetscCall(VecDestroyVecs(ts->nwork, &ts->work));

    2581:   PetscCall(MatDestroy(&ts->Jacprhs));
    2582:   PetscCall(MatDestroy(&ts->Jacp));
    2583:   if (ts->forward_solve) PetscCall(TSForwardReset(ts));
    2584:   if (ts->quadraturets) {
    2585:     PetscCall(TSReset(ts->quadraturets));
    2586:     PetscCall(VecDestroy(&ts->vec_costintegrand));
    2587:   }
    2588:   while (ilink) {
    2589:     next = ilink->next;
    2590:     PetscCall(TSDestroy(&ilink->ts));
    2591:     PetscCall(PetscFree(ilink->splitname));
    2592:     PetscCall(ISDestroy(&ilink->is));
    2593:     PetscCall(PetscFree(ilink));
    2594:     ilink = next;
    2595:   }
    2596:   ts->tsrhssplit     = NULL;
    2597:   ts->num_rhs_splits = 0;
    2598:   if (ts->tspan) {
    2599:     PetscCall(PetscFree(ts->tspan->span_times));
    2600:     PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
    2601:     PetscCall(PetscFree(ts->tspan));
    2602:   }
    2603:   ts->rhsjacobian.time  = PETSC_MIN_REAL;
    2604:   ts->rhsjacobian.scale = 1.0;
    2605:   ts->ijacobian.shift   = 1.0;
    2606:   ts->setupcalled       = PETSC_FALSE;
    2607:   PetscFunctionReturn(PETSC_SUCCESS);
    2608: }

    2610: static PetscErrorCode TSResizeReset(TS);

    2612: /*@
    2613:   TSDestroy - Destroys the timestepper context that was created
    2614:   with `TSCreate()`.

    2616:   Collective

    2618:   Input Parameter:
    2619: . ts - the `TS` context obtained from `TSCreate()`

    2621:   Level: beginner

    2623: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
    2624: @*/
    2625: PetscErrorCode TSDestroy(TS *ts)
    2626: {
    2627:   PetscFunctionBegin;
    2628:   if (!*ts) PetscFunctionReturn(PETSC_SUCCESS);
    2630:   if (--((PetscObject)*ts)->refct > 0) {
    2631:     *ts = NULL;
    2632:     PetscFunctionReturn(PETSC_SUCCESS);
    2633:   }

    2635:   PetscCall(TSReset(*ts));
    2636:   PetscCall(TSAdjointReset(*ts));
    2637:   if ((*ts)->forward_solve) PetscCall(TSForwardReset(*ts));
    2638:   PetscCall(TSResizeReset(*ts));

    2640:   /* if memory was published with SAWs then destroy it */
    2641:   PetscCall(PetscObjectSAWsViewOff((PetscObject)*ts));
    2642:   PetscTryTypeMethod(*ts, destroy);

    2644:   PetscCall(TSTrajectoryDestroy(&(*ts)->trajectory));

    2646:   PetscCall(TSAdaptDestroy(&(*ts)->adapt));
    2647:   PetscCall(TSEventDestroy(&(*ts)->event));

    2649:   PetscCall(SNESDestroy(&(*ts)->snes));
    2650:   PetscCall(SNESDestroy(&(*ts)->snesrhssplit));
    2651:   PetscCall(DMDestroy(&(*ts)->dm));
    2652:   PetscCall(TSMonitorCancel(*ts));
    2653:   PetscCall(TSAdjointMonitorCancel(*ts));

    2655:   PetscCall(TSDestroy(&(*ts)->quadraturets));
    2656:   PetscCall(PetscHeaderDestroy(ts));
    2657:   PetscFunctionReturn(PETSC_SUCCESS);
    2658: }

    2660: /*@
    2661:   TSGetSNES - Returns the `SNES` (nonlinear solver) associated with
    2662:   a `TS` (timestepper) context. Valid only for nonlinear problems.

    2664:   Not Collective, but snes is parallel if ts is parallel

    2666:   Input Parameter:
    2667: . ts - the `TS` context obtained from `TSCreate()`

    2669:   Output Parameter:
    2670: . snes - the nonlinear solver context

    2672:   Level: beginner

    2674:   Notes:
    2675:   The user can then directly manipulate the `SNES` context to set various
    2676:   options, etc.  Likewise, the user can then extract and manipulate the
    2677:   `KSP`, and `PC` contexts as well.

    2679:   `TSGetSNES()` does not work for integrators that do not use `SNES`; in
    2680:   this case `TSGetSNES()` returns `NULL` in `snes`.

    2682: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`
    2683: @*/
    2684: PetscErrorCode TSGetSNES(TS ts, SNES *snes)
    2685: {
    2686:   PetscFunctionBegin;
    2688:   PetscAssertPointer(snes, 2);
    2689:   if (!ts->snes) {
    2690:     PetscCall(SNESCreate(PetscObjectComm((PetscObject)ts), &ts->snes));
    2691:     PetscCall(PetscObjectSetOptions((PetscObject)ts->snes, ((PetscObject)ts)->options));
    2692:     PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
    2693:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->snes, (PetscObject)ts, 1));
    2694:     if (ts->dm) PetscCall(SNESSetDM(ts->snes, ts->dm));
    2695:     if (ts->problem_type == TS_LINEAR) PetscCall(SNESSetType(ts->snes, SNESKSPONLY));
    2696:   }
    2697:   *snes = ts->snes;
    2698:   PetscFunctionReturn(PETSC_SUCCESS);
    2699: }

    2701: /*@
    2702:   TSSetSNES - Set the `SNES` (nonlinear solver) to be used by the `TS` timestepping context

    2704:   Collective

    2706:   Input Parameters:
    2707: + ts   - the `TS` context obtained from `TSCreate()`
    2708: - snes - the nonlinear solver context

    2710:   Level: developer

    2712:   Note:
    2713:   Most users should have the `TS` created by calling `TSGetSNES()`

    2715: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
    2716: @*/
    2717: PetscErrorCode TSSetSNES(TS ts, SNES snes)
    2718: {
    2719:   PetscErrorCode (*func)(SNES, Vec, Mat, Mat, void *);

    2721:   PetscFunctionBegin;
    2724:   PetscCall(PetscObjectReference((PetscObject)snes));
    2725:   PetscCall(SNESDestroy(&ts->snes));

    2727:   ts->snes = snes;

    2729:   PetscCall(SNESSetFunction(ts->snes, NULL, SNESTSFormFunction, ts));
    2730:   PetscCall(SNESGetJacobian(ts->snes, NULL, NULL, &func, NULL));
    2731:   if (func == SNESTSFormJacobian) PetscCall(SNESSetJacobian(ts->snes, NULL, NULL, SNESTSFormJacobian, ts));
    2732:   PetscFunctionReturn(PETSC_SUCCESS);
    2733: }

    2735: /*@
    2736:   TSGetKSP - Returns the `KSP` (linear solver) associated with
    2737:   a `TS` (timestepper) context.

    2739:   Not Collective, but `ksp` is parallel if `ts` is parallel

    2741:   Input Parameter:
    2742: . ts - the `TS` context obtained from `TSCreate()`

    2744:   Output Parameter:
    2745: . ksp - the nonlinear solver context

    2747:   Level: beginner

    2749:   Notes:
    2750:   The user can then directly manipulate the `KSP` context to set various
    2751:   options, etc.  Likewise, the user can then extract and manipulate the
    2752:   `PC` context as well.

    2754:   `TSGetKSP()` does not work for integrators that do not use `KSP`;
    2755:   in this case `TSGetKSP()` returns `NULL` in `ksp`.

    2757: .seealso: [](ch_ts), `TS`, `SNES`, `KSP`, `TSCreate()`, `TSSetUp()`, `TSSolve()`, `TSGetSNES()`
    2758: @*/
    2759: PetscErrorCode TSGetKSP(TS ts, KSP *ksp)
    2760: {
    2761:   SNES snes;

    2763:   PetscFunctionBegin;
    2765:   PetscAssertPointer(ksp, 2);
    2766:   PetscCheck(((PetscObject)ts)->type_name, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "KSP is not created yet. Call TSSetType() first");
    2767:   PetscCheck(ts->problem_type == TS_LINEAR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Linear only; use TSGetSNES()");
    2768:   PetscCall(TSGetSNES(ts, &snes));
    2769:   PetscCall(SNESGetKSP(snes, ksp));
    2770:   PetscFunctionReturn(PETSC_SUCCESS);
    2771: }

    2773: /* ----------- Routines to set solver parameters ---------- */

    2775: /*@
    2776:   TSSetMaxSteps - Sets the maximum number of steps to use.

    2778:   Logically Collective

    2780:   Input Parameters:
    2781: + ts       - the `TS` context obtained from `TSCreate()`
    2782: - maxsteps - maximum number of steps to use

    2784:   Options Database Key:
    2785: . -ts_max_steps <maxsteps> - Sets maxsteps

    2787:   Level: intermediate

    2789:   Note:
    2790:   Use `PETSC_DETERMINE` to reset the maximum number of steps to the default from when the object's type was set

    2792:   The default maximum number of steps is 5,000

    2794:   Fortran Note:
    2795:   Use `PETSC_DETERMINE_INTEGER`

    2797: .seealso: [](ch_ts), `TS`, `TSGetMaxSteps()`, `TSSetMaxTime()`, `TSSetExactFinalTime()`
    2798: @*/
    2799: PetscErrorCode TSSetMaxSteps(TS ts, PetscInt maxsteps)
    2800: {
    2801:   PetscFunctionBegin;
    2804:   if (maxsteps == PETSC_DETERMINE) {
    2805:     ts->max_steps = ts->default_max_steps;
    2806:   } else {
    2807:     PetscCheck(maxsteps >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of steps must be non-negative");
    2808:     ts->max_steps = maxsteps;
    2809:   }
    2810:   PetscFunctionReturn(PETSC_SUCCESS);
    2811: }

    2813: /*@
    2814:   TSGetMaxSteps - Gets the maximum number of steps to use.

    2816:   Not Collective

    2818:   Input Parameter:
    2819: . ts - the `TS` context obtained from `TSCreate()`

    2821:   Output Parameter:
    2822: . maxsteps - maximum number of steps to use

    2824:   Level: advanced

    2826: .seealso: [](ch_ts), `TS`, `TSSetMaxSteps()`, `TSGetMaxTime()`, `TSSetMaxTime()`
    2827: @*/
    2828: PetscErrorCode TSGetMaxSteps(TS ts, PetscInt *maxsteps)
    2829: {
    2830:   PetscFunctionBegin;
    2832:   PetscAssertPointer(maxsteps, 2);
    2833:   *maxsteps = ts->max_steps;
    2834:   PetscFunctionReturn(PETSC_SUCCESS);
    2835: }

    2837: /*@
    2838:   TSSetMaxTime - Sets the maximum (or final) time for timestepping.

    2840:   Logically Collective

    2842:   Input Parameters:
    2843: + ts      - the `TS` context obtained from `TSCreate()`
    2844: - maxtime - final time to step to

    2846:   Options Database Key:
    2847: . -ts_max_time <maxtime> - Sets maxtime

    2849:   Level: intermediate

    2851:   Notes:
    2852:   Use `PETSC_DETERMINE` to reset the maximum time to the default from when the object's type was set

    2854:   The default maximum time is 5.0

    2856:   Fortran Note:
    2857:   Use `PETSC_DETERMINE_REAL`

    2859: .seealso: [](ch_ts), `TS`, `TSGetMaxTime()`, `TSSetMaxSteps()`, `TSSetExactFinalTime()`
    2860: @*/
    2861: PetscErrorCode TSSetMaxTime(TS ts, PetscReal maxtime)
    2862: {
    2863:   PetscFunctionBegin;
    2866:   if (maxtime == PETSC_DETERMINE) {
    2867:     ts->max_time = ts->default_max_time;
    2868:   } else {
    2869:     ts->max_time = maxtime;
    2870:   }
    2871:   PetscFunctionReturn(PETSC_SUCCESS);
    2872: }

    2874: /*@
    2875:   TSGetMaxTime - Gets the maximum (or final) time for timestepping.

    2877:   Not Collective

    2879:   Input Parameter:
    2880: . ts - the `TS` context obtained from `TSCreate()`

    2882:   Output Parameter:
    2883: . maxtime - final time to step to

    2885:   Level: advanced

    2887: .seealso: [](ch_ts), `TS`, `TSSetMaxTime()`, `TSGetMaxSteps()`, `TSSetMaxSteps()`
    2888: @*/
    2889: PetscErrorCode TSGetMaxTime(TS ts, PetscReal *maxtime)
    2890: {
    2891:   PetscFunctionBegin;
    2893:   PetscAssertPointer(maxtime, 2);
    2894:   *maxtime = ts->max_time;
    2895:   PetscFunctionReturn(PETSC_SUCCESS);
    2896: }

    2898: // PetscClangLinter pragma disable: -fdoc-*
    2899: /*@
    2900:   TSSetInitialTimeStep - Deprecated, use `TSSetTime()` and `TSSetTimeStep()`.

    2902:   Level: deprecated

    2904: @*/
    2905: PetscErrorCode TSSetInitialTimeStep(TS ts, PetscReal initial_time, PetscReal time_step)
    2906: {
    2907:   PetscFunctionBegin;
    2909:   PetscCall(TSSetTime(ts, initial_time));
    2910:   PetscCall(TSSetTimeStep(ts, time_step));
    2911:   PetscFunctionReturn(PETSC_SUCCESS);
    2912: }

    2914: // PetscClangLinter pragma disable: -fdoc-*
    2915: /*@
    2916:   TSGetDuration - Deprecated, use `TSGetMaxSteps()` and `TSGetMaxTime()`.

    2918:   Level: deprecated

    2920: @*/
    2921: PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime)
    2922: {
    2923:   PetscFunctionBegin;
    2925:   if (maxsteps) {
    2926:     PetscAssertPointer(maxsteps, 2);
    2927:     *maxsteps = ts->max_steps;
    2928:   }
    2929:   if (maxtime) {
    2930:     PetscAssertPointer(maxtime, 3);
    2931:     *maxtime = ts->max_time;
    2932:   }
    2933:   PetscFunctionReturn(PETSC_SUCCESS);
    2934: }

    2936: // PetscClangLinter pragma disable: -fdoc-*
    2937: /*@
    2938:   TSSetDuration - Deprecated, use `TSSetMaxSteps()` and `TSSetMaxTime()`.

    2940:   Level: deprecated

    2942: @*/
    2943: PetscErrorCode TSSetDuration(TS ts, PetscInt maxsteps, PetscReal maxtime)
    2944: {
    2945:   PetscFunctionBegin;
    2946:   if (maxsteps != (PetscInt)PETSC_CURRENT) PetscCall(TSSetMaxSteps(ts, maxsteps));
    2947:   if (maxtime != (PetscReal)PETSC_CURRENT) PetscCall(TSSetMaxTime(ts, maxtime));
    2948:   PetscFunctionReturn(PETSC_SUCCESS);
    2949: }

    2951: // PetscClangLinter pragma disable: -fdoc-*
    2952: /*@
    2953:   TSGetTimeStepNumber - Deprecated, use `TSGetStepNumber()`.

    2955:   Level: deprecated

    2957: @*/
    2958: PetscErrorCode TSGetTimeStepNumber(TS ts, PetscInt *steps)
    2959: {
    2960:   return TSGetStepNumber(ts, steps);
    2961: }

    2963: // PetscClangLinter pragma disable: -fdoc-*
    2964: /*@
    2965:   TSGetTotalSteps - Deprecated, use `TSGetStepNumber()`.

    2967:   Level: deprecated

    2969: @*/
    2970: PetscErrorCode TSGetTotalSteps(TS ts, PetscInt *steps)
    2971: {
    2972:   return TSGetStepNumber(ts, steps);
    2973: }

    2975: /*@
    2976:   TSSetSolution - Sets the initial solution vector
    2977:   for use by the `TS` routines.

    2979:   Logically Collective

    2981:   Input Parameters:
    2982: + ts - the `TS` context obtained from `TSCreate()`
    2983: - u  - the solution vector

    2985:   Level: beginner

    2987: .seealso: [](ch_ts), `TS`, `TSSetSolutionFunction()`, `TSGetSolution()`, `TSCreate()`
    2988: @*/
    2989: PetscErrorCode TSSetSolution(TS ts, Vec u)
    2990: {
    2991:   DM dm;

    2993:   PetscFunctionBegin;
    2996:   PetscCall(PetscObjectReference((PetscObject)u));
    2997:   PetscCall(VecDestroy(&ts->vec_sol));
    2998:   ts->vec_sol = u;

    3000:   PetscCall(TSGetDM(ts, &dm));
    3001:   PetscCall(DMShellSetGlobalVector(dm, u));
    3002:   PetscFunctionReturn(PETSC_SUCCESS);
    3003: }

    3005: /*@C
    3006:   TSSetPreStep - Sets the general-purpose function
    3007:   called once at the beginning of each time step.

    3009:   Logically Collective

    3011:   Input Parameters:
    3012: + ts   - The `TS` context obtained from `TSCreate()`
    3013: - func - The function

    3015:   Calling sequence of `func`:
    3016: . ts - the `TS` context

    3018:   Level: intermediate

    3020: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPostStage()`, `TSSetPostStep()`, `TSStep()`, `TSRestartStep()`
    3021: @*/
    3022: PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS ts))
    3023: {
    3024:   PetscFunctionBegin;
    3026:   ts->prestep = func;
    3027:   PetscFunctionReturn(PETSC_SUCCESS);
    3028: }

    3030: /*@
    3031:   TSPreStep - Runs the user-defined pre-step function provided with `TSSetPreStep()`

    3033:   Collective

    3035:   Input Parameter:
    3036: . ts - The `TS` context obtained from `TSCreate()`

    3038:   Level: developer

    3040:   Note:
    3041:   `TSPreStep()` is typically used within time stepping implementations,
    3042:   so most users would not generally call this routine themselves.

    3044: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSPreStage()`, `TSPostStage()`, `TSPostStep()`
    3045: @*/
    3046: PetscErrorCode TSPreStep(TS ts)
    3047: {
    3048:   PetscFunctionBegin;
    3050:   if (ts->prestep) {
    3051:     Vec              U;
    3052:     PetscObjectId    idprev;
    3053:     PetscBool        sameObject;
    3054:     PetscObjectState sprev, spost;

    3056:     PetscCall(TSGetSolution(ts, &U));
    3057:     PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
    3058:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
    3059:     PetscCallBack("TS callback preset", (*ts->prestep)(ts));
    3060:     PetscCall(TSGetSolution(ts, &U));
    3061:     PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
    3062:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
    3063:     if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
    3064:   }
    3065:   PetscFunctionReturn(PETSC_SUCCESS);
    3066: }

    3068: /*@C
    3069:   TSSetPreStage - Sets the general-purpose function
    3070:   called once at the beginning of each stage.

    3072:   Logically Collective

    3074:   Input Parameters:
    3075: + ts   - The `TS` context obtained from `TSCreate()`
    3076: - func - The function

    3078:   Calling sequence of `func`:
    3079: + ts        - the `TS` context
    3080: - stagetime - the stage time

    3082:   Level: intermediate

    3084:   Note:
    3085:   There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
    3086:   The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
    3087:   attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.

    3089: .seealso: [](ch_ts), `TS`, `TSSetPostStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
    3090: @*/
    3091: PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime))
    3092: {
    3093:   PetscFunctionBegin;
    3095:   ts->prestage = func;
    3096:   PetscFunctionReturn(PETSC_SUCCESS);
    3097: }

    3099: /*@C
    3100:   TSSetPostStage - Sets the general-purpose function
    3101:   called once at the end of each stage.

    3103:   Logically Collective

    3105:   Input Parameters:
    3106: + ts   - The `TS` context obtained from `TSCreate()`
    3107: - func - The function

    3109:   Calling sequence of `func`:
    3110: + ts         - the `TS` context
    3111: . stagetime  - the stage time
    3112: . stageindex - the stage index
    3113: - Y          - Array of vectors (of size = total number of stages) with the stage solutions

    3115:   Level: intermediate

    3117:   Note:
    3118:   There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried.
    3119:   The time step number being computed can be queried using `TSGetStepNumber()` and the total size of the step being
    3120:   attempted can be obtained using `TSGetTimeStep()`. The time at the start of the step is available via `TSGetTime()`.

    3122: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
    3123: @*/
    3124: PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y))
    3125: {
    3126:   PetscFunctionBegin;
    3128:   ts->poststage = func;
    3129:   PetscFunctionReturn(PETSC_SUCCESS);
    3130: }

    3132: /*@C
    3133:   TSSetPostEvaluate - Sets the general-purpose function
    3134:   called at the end of each step evaluation.

    3136:   Logically Collective

    3138:   Input Parameters:
    3139: + ts   - The `TS` context obtained from `TSCreate()`
    3140: - func - The function

    3142:   Calling sequence of `func`:
    3143: . ts - the `TS` context

    3145:   Level: intermediate

    3147:   Note:
    3148:   The function set by `TSSetPostEvaluate()` is called after the solution is evaluated, or after the step rollback.
    3149:   Inside the `func` callback, the solution vector can be obtained with `TSGetSolution()`, and modified, if need be.
    3150:   The time step can be obtained with `TSGetTimeStep()`, and the time at the start of the step - via `TSGetTime()`.
    3151:   The potential changes to the solution vector introduced by event handling (`postevent()`) are not relevant for `TSSetPostEvaluate()`,
    3152:   but are relevant for `TSSetPostStep()`, according to the function call scheme in `TSSolve()`, as shown below
    3153: .vb
    3154:   ...
    3155:   Step()
    3156:   PostEvaluate()
    3157:   EventHandling()
    3158:   step_rollback ? PostEvaluate() : PostStep()
    3159:   ...
    3160: .ve
    3161:   where EventHandling() may result in one of the following three outcomes
    3162: .vb
    3163:   (1) | successful step | solution intact
    3164:   (2) | successful step | solution modified by `postevent()`
    3165:   (3) | step_rollback   | solution rolled back
    3166: .ve

    3168: .seealso: [](ch_ts), `TS`, `TSSetPreStage()`, `TSSetPreStep()`, `TSSetPostStep()`, `TSGetApplicationContext()`
    3169: @*/
    3170: PetscErrorCode TSSetPostEvaluate(TS ts, PetscErrorCode (*func)(TS ts))
    3171: {
    3172:   PetscFunctionBegin;
    3174:   ts->postevaluate = func;
    3175:   PetscFunctionReturn(PETSC_SUCCESS);
    3176: }

    3178: /*@
    3179:   TSPreStage - Runs the user-defined pre-stage function set using `TSSetPreStage()`

    3181:   Collective

    3183:   Input Parameters:
    3184: + ts        - The `TS` context obtained from `TSCreate()`
    3185: - stagetime - The absolute time of the current stage

    3187:   Level: developer

    3189:   Note:
    3190:   `TSPreStage()` is typically used within time stepping implementations,
    3191:   most users would not generally call this routine themselves.

    3193: .seealso: [](ch_ts), `TS`, `TSPostStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
    3194: @*/
    3195: PetscErrorCode TSPreStage(TS ts, PetscReal stagetime)
    3196: {
    3197:   PetscFunctionBegin;
    3199:   if (ts->prestage) PetscCallBack("TS callback prestage", (*ts->prestage)(ts, stagetime));
    3200:   PetscFunctionReturn(PETSC_SUCCESS);
    3201: }

    3203: /*@
    3204:   TSPostStage - Runs the user-defined post-stage function set using `TSSetPostStage()`

    3206:   Collective

    3208:   Input Parameters:
    3209: + ts         - The `TS` context obtained from `TSCreate()`
    3210: . stagetime  - The absolute time of the current stage
    3211: . stageindex - Stage number
    3212: - Y          - Array of vectors (of size = total number of stages) with the stage solutions

    3214:   Level: developer

    3216:   Note:
    3217:   `TSPostStage()` is typically used within time stepping implementations,
    3218:   most users would not generally call this routine themselves.

    3220: .seealso: [](ch_ts), `TS`, `TSPreStage()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
    3221: @*/
    3222: PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
    3223: {
    3224:   PetscFunctionBegin;
    3226:   if (ts->poststage) PetscCallBack("TS callback poststage", (*ts->poststage)(ts, stagetime, stageindex, Y));
    3227:   PetscFunctionReturn(PETSC_SUCCESS);
    3228: }

    3230: /*@
    3231:   TSPostEvaluate - Runs the user-defined post-evaluate function set using `TSSetPostEvaluate()`

    3233:   Collective

    3235:   Input Parameter:
    3236: . ts - The `TS` context obtained from `TSCreate()`

    3238:   Level: developer

    3240:   Note:
    3241:   `TSPostEvaluate()` is typically used within time stepping implementations,
    3242:   most users would not generally call this routine themselves.

    3244: .seealso: [](ch_ts), `TS`, `TSSetPostEvaluate()`, `TSSetPreStep()`, `TSPreStep()`, `TSPostStep()`
    3245: @*/
    3246: PetscErrorCode TSPostEvaluate(TS ts)
    3247: {
    3248:   PetscFunctionBegin;
    3250:   if (ts->postevaluate) {
    3251:     Vec              U;
    3252:     PetscObjectState sprev, spost;

    3254:     PetscCall(TSGetSolution(ts, &U));
    3255:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
    3256:     PetscCallBack("TS callback postevaluate", (*ts->postevaluate)(ts));
    3257:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
    3258:     if (sprev != spost) PetscCall(TSRestartStep(ts));
    3259:   }
    3260:   PetscFunctionReturn(PETSC_SUCCESS);
    3261: }

    3263: /*@C
    3264:   TSSetPostStep - Sets the general-purpose function
    3265:   called once at the end of each successful time step.

    3267:   Logically Collective

    3269:   Input Parameters:
    3270: + ts   - The `TS` context obtained from `TSCreate()`
    3271: - func - The function

    3273:   Calling sequence of `func`:
    3274: . ts - the `TS` context

    3276:   Level: intermediate

    3278:   Note:
    3279:   The function set by `TSSetPostStep()` is called after each successful step. If the event handler locates an event at the
    3280:   given step, and `postevent()` modifies the solution vector, the solution vector obtained by `TSGetSolution()` inside `func` will
    3281:   contain the changes. To get the solution without these changes, use `TSSetPostEvaluate()` to set the appropriate callback.
    3282:   The scheme of the relevant function calls in `TSSolve()` is shown below
    3283: .vb
    3284:   ...
    3285:   Step()
    3286:   PostEvaluate()
    3287:   EventHandling()
    3288:   step_rollback ? PostEvaluate() : PostStep()
    3289:   ...
    3290: .ve
    3291:   where EventHandling() may result in one of the following three outcomes
    3292: .vb
    3293:   (1) | successful step | solution intact
    3294:   (2) | successful step | solution modified by `postevent()`
    3295:   (3) | step_rollback   | solution rolled back
    3296: .ve

    3298: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSRestartStep()`
    3299: @*/
    3300: PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS ts))
    3301: {
    3302:   PetscFunctionBegin;
    3304:   ts->poststep = func;
    3305:   PetscFunctionReturn(PETSC_SUCCESS);
    3306: }

    3308: /*@
    3309:   TSPostStep - Runs the user-defined post-step function that was set with `TSSetPostStep()`

    3311:   Collective

    3313:   Input Parameter:
    3314: . ts - The `TS` context obtained from `TSCreate()`

    3316:   Note:
    3317:   `TSPostStep()` is typically used within time stepping implementations,
    3318:   so most users would not generally call this routine themselves.

    3320:   Level: developer

    3322: .seealso: [](ch_ts), `TS`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostEvaluate()`, `TSGetTimeStep()`, `TSGetStepNumber()`, `TSGetTime()`, `TSSetPotsStep()`
    3323: @*/
    3324: PetscErrorCode TSPostStep(TS ts)
    3325: {
    3326:   PetscFunctionBegin;
    3328:   if (ts->poststep) {
    3329:     Vec              U;
    3330:     PetscObjectId    idprev;
    3331:     PetscBool        sameObject;
    3332:     PetscObjectState sprev, spost;

    3334:     PetscCall(TSGetSolution(ts, &U));
    3335:     PetscCall(PetscObjectGetId((PetscObject)U, &idprev));
    3336:     PetscCall(PetscObjectStateGet((PetscObject)U, &sprev));
    3337:     PetscCallBack("TS callback poststep", (*ts->poststep)(ts));
    3338:     PetscCall(TSGetSolution(ts, &U));
    3339:     PetscCall(PetscObjectCompareId((PetscObject)U, idprev, &sameObject));
    3340:     PetscCall(PetscObjectStateGet((PetscObject)U, &spost));
    3341:     if (!sameObject || sprev != spost) PetscCall(TSRestartStep(ts));
    3342:   }
    3343:   PetscFunctionReturn(PETSC_SUCCESS);
    3344: }

    3346: /*@
    3347:   TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval

    3349:   Collective

    3351:   Input Parameters:
    3352: + ts - time stepping context
    3353: - t  - time to interpolate to

    3355:   Output Parameter:
    3356: . U - state at given time

    3358:   Level: intermediate

    3360:   Developer Notes:
    3361:   `TSInterpolate()` and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints.

    3363: .seealso: [](ch_ts), `TS`, `TSSetExactFinalTime()`, `TSSolve()`
    3364: @*/
    3365: PetscErrorCode TSInterpolate(TS ts, PetscReal t, Vec U)
    3366: {
    3367:   PetscFunctionBegin;
    3370:   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);
    3371:   PetscUseTypeMethod(ts, interpolate, t, U);
    3372:   PetscFunctionReturn(PETSC_SUCCESS);
    3373: }

    3375: /*@
    3376:   TSStep - Steps one time step

    3378:   Collective

    3380:   Input Parameter:
    3381: . ts - the `TS` context obtained from `TSCreate()`

    3383:   Level: developer

    3385:   Notes:
    3386:   The public interface for the ODE/DAE solvers is `TSSolve()`, you should almost for sure be using that routine and not this routine.

    3388:   The hook set using `TSSetPreStep()` is called before each attempt to take the step. In general, the time step size may
    3389:   be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages.

    3391:   This may over-step the final time provided in `TSSetMaxTime()` depending on the time-step used. `TSSolve()` interpolates to exactly the
    3392:   time provided in `TSSetMaxTime()`. One can use `TSInterpolate()` to determine an interpolated solution within the final timestep.

    3394: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSSetPostStage()`, `TSInterpolate()`
    3395: @*/
    3396: PetscErrorCode TSStep(TS ts)
    3397: {
    3398:   static PetscBool cite = PETSC_FALSE;
    3399:   PetscReal        ptime;

    3401:   PetscFunctionBegin;
    3403:   PetscCall(PetscCitationsRegister("@article{tspaper,\n"
    3404:                                    "  title         = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n"
    3405:                                    "  author        = {Abhyankar, Shrirang and Brown, Jed and Constantinescu, Emil and Ghosh, Debojyoti and Smith, Barry F. and Zhang, Hong},\n"
    3406:                                    "  journal       = {arXiv e-preprints},\n"
    3407:                                    "  eprint        = {1806.01437},\n"
    3408:                                    "  archivePrefix = {arXiv},\n"
    3409:                                    "  year          = {2018}\n}\n",
    3410:                                    &cite));
    3411:   PetscCall(TSSetUp(ts));
    3412:   PetscCall(TSTrajectorySetUp(ts->trajectory, ts));
    3413:   if (ts->tspan)
    3414:     ts->tspan->worktol = 0; /* In each step of TSSolve() 'tspan->worktol' will be meaningfully defined (later) only once:
    3415:                                                    in TSAdaptChoose() or TSEvent_dt_cap(), and then reused till the end of the step */

    3417:   PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_INT_MAX, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
    3418:   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()");
    3419:   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");

    3421:   if (!ts->vec_sol0) PetscCall(VecDuplicate(ts->vec_sol, &ts->vec_sol0));
    3422:   PetscCall(VecCopy(ts->vec_sol, ts->vec_sol0));
    3423:   ts->time_step0 = ts->time_step;

    3425:   if (!ts->steps) ts->ptime_prev = ts->ptime;
    3426:   ptime = ts->ptime;

    3428:   ts->ptime_prev_rollback = ts->ptime_prev;
    3429:   ts->reason              = TS_CONVERGED_ITERATING;

    3431:   PetscCall(PetscLogEventBegin(TS_Step, ts, 0, 0, 0));
    3432:   PetscUseTypeMethod(ts, step);
    3433:   PetscCall(PetscLogEventEnd(TS_Step, ts, 0, 0, 0));

    3435:   if (ts->reason >= 0) {
    3436:     ts->ptime_prev = ptime;
    3437:     ts->steps++;
    3438:     ts->steprollback = PETSC_FALSE;
    3439:     ts->steprestart  = PETSC_FALSE;
    3440:     ts->stepresize   = PETSC_FALSE;
    3441:   }

    3443:   if (ts->reason < 0 && ts->errorifstepfailed) {
    3444:     PetscCall(TSMonitorCancel(ts));
    3445:     PetscCheck(ts->reason != TS_DIVERGED_NONLINEAR_SOLVE, PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s, increase -ts_max_snes_failures or use unlimited to attempt recovery", TSConvergedReasons[ts->reason]);
    3446:     SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_NOT_CONVERGED, "TSStep has failed due to %s", TSConvergedReasons[ts->reason]);
    3447:   }
    3448:   PetscFunctionReturn(PETSC_SUCCESS);
    3449: }

    3451: /*@
    3452:   TSEvaluateWLTE - Evaluate the weighted local truncation error norm
    3453:   at the end of a time step with a given order of accuracy.

    3455:   Collective

    3457:   Input Parameters:
    3458: + ts        - time stepping context
    3459: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

    3461:   Input/Output Parameter:
    3462: . order - optional, desired order for the error evaluation or `PETSC_DECIDE`;
    3463:            on output, the actual order of the error evaluation

    3465:   Output Parameter:
    3466: . wlte - the weighted local truncation error norm

    3468:   Level: advanced

    3470:   Note:
    3471:   If the timestepper cannot evaluate the error in a particular step
    3472:   (eg. in the first step or restart steps after event handling),
    3473:   this routine returns wlte=-1.0 .

    3475: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`, `TSErrorWeightedNorm()`
    3476: @*/
    3477: PetscErrorCode TSEvaluateWLTE(TS ts, NormType wnormtype, PetscInt *order, PetscReal *wlte)
    3478: {
    3479:   PetscFunctionBegin;
    3483:   if (order) PetscAssertPointer(order, 3);
    3485:   PetscAssertPointer(wlte, 4);
    3486:   PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
    3487:   PetscUseTypeMethod(ts, evaluatewlte, wnormtype, order, wlte);
    3488:   PetscFunctionReturn(PETSC_SUCCESS);
    3489: }

    3491: /*@
    3492:   TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy.

    3494:   Collective

    3496:   Input Parameters:
    3497: + ts    - time stepping context
    3498: . order - desired order of accuracy
    3499: - done  - whether the step was evaluated at this order (pass `NULL` to generate an error if not available)

    3501:   Output Parameter:
    3502: . U - state at the end of the current step

    3504:   Level: advanced

    3506:   Notes:
    3507:   This function cannot be called until all stages have been evaluated.

    3509:   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.

    3511: .seealso: [](ch_ts), `TS`, `TSStep()`, `TSAdapt`
    3512: @*/
    3513: PetscErrorCode TSEvaluateStep(TS ts, PetscInt order, Vec U, PetscBool *done)
    3514: {
    3515:   PetscFunctionBegin;
    3519:   PetscUseTypeMethod(ts, evaluatestep, order, U, done);
    3520:   PetscFunctionReturn(PETSC_SUCCESS);
    3521: }

    3523: /*@C
    3524:   TSGetComputeInitialCondition - Get the function used to automatically compute an initial condition for the timestepping.

    3526:   Not collective

    3528:   Input Parameter:
    3529: . ts - time stepping context

    3531:   Output Parameter:
    3532: . initCondition - The function which computes an initial condition

    3534:   Calling sequence of `initCondition`:
    3535: + ts - The timestepping context
    3536: - u  - The input vector in which the initial condition is stored

    3538:   Level: advanced

    3540: .seealso: [](ch_ts), `TS`, `TSSetComputeInitialCondition()`, `TSComputeInitialCondition()`
    3541: @*/
    3542: PetscErrorCode TSGetComputeInitialCondition(TS ts, PetscErrorCode (**initCondition)(TS ts, Vec u))
    3543: {
    3544:   PetscFunctionBegin;
    3546:   PetscAssertPointer(initCondition, 2);
    3547:   *initCondition = ts->ops->initcondition;
    3548:   PetscFunctionReturn(PETSC_SUCCESS);
    3549: }

    3551: /*@C
    3552:   TSSetComputeInitialCondition - Set the function used to automatically compute an initial condition for the timestepping.

    3554:   Logically collective

    3556:   Input Parameters:
    3557: + ts            - time stepping context
    3558: - initCondition - The function which computes an initial condition

    3560:   Calling sequence of `initCondition`:
    3561: + ts - The timestepping context
    3562: - e  - The input vector in which the initial condition is to be stored

    3564:   Level: advanced

    3566: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSComputeInitialCondition()`
    3567: @*/
    3568: PetscErrorCode TSSetComputeInitialCondition(TS ts, PetscErrorCode (*initCondition)(TS ts, Vec e))
    3569: {
    3570:   PetscFunctionBegin;
    3573:   ts->ops->initcondition = initCondition;
    3574:   PetscFunctionReturn(PETSC_SUCCESS);
    3575: }

    3577: /*@
    3578:   TSComputeInitialCondition - Compute an initial condition for the timestepping using the function previously set with `TSSetComputeInitialCondition()`

    3580:   Collective

    3582:   Input Parameters:
    3583: + ts - time stepping context
    3584: - u  - The `Vec` to store the condition in which will be used in `TSSolve()`

    3586:   Level: advanced

    3588: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
    3589: @*/
    3590: PetscErrorCode TSComputeInitialCondition(TS ts, Vec u)
    3591: {
    3592:   PetscFunctionBegin;
    3595:   PetscTryTypeMethod(ts, initcondition, u);
    3596:   PetscFunctionReturn(PETSC_SUCCESS);
    3597: }

    3599: /*@C
    3600:   TSGetComputeExactError - Get the function used to automatically compute the exact error for the timestepping.

    3602:   Not collective

    3604:   Input Parameter:
    3605: . ts - time stepping context

    3607:   Output Parameter:
    3608: . exactError - The function which computes the solution error

    3610:   Calling sequence of `exactError`:
    3611: + ts - The timestepping context
    3612: . u  - The approximate solution vector
    3613: - e  - The vector in which the error is stored

    3615:   Level: advanced

    3617: .seealso: [](ch_ts), `TS`, `TSComputeExactError()`
    3618: @*/
    3619: PetscErrorCode TSGetComputeExactError(TS ts, PetscErrorCode (**exactError)(TS ts, Vec u, Vec e))
    3620: {
    3621:   PetscFunctionBegin;
    3623:   PetscAssertPointer(exactError, 2);
    3624:   *exactError = ts->ops->exacterror;
    3625:   PetscFunctionReturn(PETSC_SUCCESS);
    3626: }

    3628: /*@C
    3629:   TSSetComputeExactError - Set the function used to automatically compute the exact error for the timestepping.

    3631:   Logically collective

    3633:   Input Parameters:
    3634: + ts         - time stepping context
    3635: - exactError - The function which computes the solution error

    3637:   Calling sequence of `exactError`:
    3638: + ts - The timestepping context
    3639: . u  - The approximate solution vector
    3640: - e  - The  vector in which the error is stored

    3642:   Level: advanced

    3644: .seealso: [](ch_ts), `TS`, `TSGetComputeExactError()`, `TSComputeExactError()`
    3645: @*/
    3646: PetscErrorCode TSSetComputeExactError(TS ts, PetscErrorCode (*exactError)(TS ts, Vec u, Vec e))
    3647: {
    3648:   PetscFunctionBegin;
    3651:   ts->ops->exacterror = exactError;
    3652:   PetscFunctionReturn(PETSC_SUCCESS);
    3653: }

    3655: /*@
    3656:   TSComputeExactError - Compute the solution error for the timestepping using the function previously set with `TSSetComputeExactError()`

    3658:   Collective

    3660:   Input Parameters:
    3661: + ts - time stepping context
    3662: . u  - The approximate solution
    3663: - e  - The `Vec` used to store the error

    3665:   Level: advanced

    3667: .seealso: [](ch_ts), `TS`, `TSGetComputeInitialCondition()`, `TSSetComputeInitialCondition()`, `TSSolve()`
    3668: @*/
    3669: PetscErrorCode TSComputeExactError(TS ts, Vec u, Vec e)
    3670: {
    3671:   PetscFunctionBegin;
    3675:   PetscTryTypeMethod(ts, exacterror, u, e);
    3676:   PetscFunctionReturn(PETSC_SUCCESS);
    3677: }

    3679: /*@C
    3680:   TSSetResize - Sets the resize callbacks.

    3682:   Logically Collective

    3684:   Input Parameters:
    3685: + ts       - The `TS` context obtained from `TSCreate()`
    3686: . rollback - Whether a resize will restart the step
    3687: . setup    - The setup function
    3688: . transfer - The transfer function
    3689: - ctx      - [optional] The user-defined context

    3691:   Calling sequence of `setup`:
    3692: + ts     - the `TS` context
    3693: . step   - the current step
    3694: . time   - the current time
    3695: . state  - the current vector of state
    3696: . resize - (output parameter) `PETSC_TRUE` if need resizing, `PETSC_FALSE` otherwise
    3697: - ctx    - user defined context

    3699:   Calling sequence of `transfer`:
    3700: + ts      - the `TS` context
    3701: . nv      - the number of vectors to be transferred
    3702: . vecsin  - array of vectors to be transferred
    3703: . vecsout - array of transferred vectors
    3704: - ctx     - user defined context

    3706:   Notes:
    3707:   The `setup` function is called inside `TSSolve()` after `TSEventHandler()` or after `TSPostStep()`
    3708:   depending on the `rollback` value: if `rollback` is true, then these callbacks behave as error indicators
    3709:   and will flag the need to remesh and restart the current step. Otherwise, they will simply flag the solver
    3710:   that the size of the discrete problem has changed.
    3711:   In both cases, the solver will collect the needed vectors that will be
    3712:   transferred from the old to the new sizes using the `transfer` callback. These vectors will include the
    3713:   current solution vector, and other vectors needed by the specific solver used.
    3714:   For example, `TSBDF` uses previous solutions vectors to solve for the next time step.
    3715:   Other application specific objects associated with the solver, i.e. Jacobian matrices and `DM`,
    3716:   will be automatically reset if the sizes are changed and they must be specified again by the user
    3717:   inside the `transfer` function.
    3718:   The input and output arrays passed to `transfer` are allocated by PETSc.
    3719:   Vectors in `vecsout` must be created by the user.
    3720:   Ownership of vectors in `vecsout` is transferred to PETSc.

    3722:   Level: advanced

    3724: .seealso: [](ch_ts), `TS`, `TSSetDM()`, `TSSetIJacobian()`, `TSSetRHSJacobian()`
    3725: @*/
    3726: PetscErrorCode TSSetResize(TS ts, PetscBool rollback, PetscErrorCode (*setup)(TS ts, PetscInt step, PetscReal time, Vec state, PetscBool *resize, void *ctx), PetscErrorCode (*transfer)(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx), void *ctx)
    3727: {
    3728:   PetscFunctionBegin;
    3730:   ts->resizerollback = rollback;
    3731:   ts->resizesetup    = setup;
    3732:   ts->resizetransfer = transfer;
    3733:   ts->resizectx      = ctx;
    3734:   PetscFunctionReturn(PETSC_SUCCESS);
    3735: }

    3737: /*
    3738:   TSResizeRegisterOrRetrieve - Register or import vectors transferred with `TSResize()`.

    3740:   Collective

    3742:   Input Parameters:
    3743: + ts   - The `TS` context obtained from `TSCreate()`
    3744: - flg - If `PETSC_TRUE` each TS implementation (e.g. `TSBDF`) will register vectors to be transferred, if `PETSC_FALSE` vectors will be imported from transferred vectors.

    3746:   Level: developer

    3748:   Note:
    3749:   `TSResizeRegisterOrRetrieve()` is declared PETSC_INTERN since it is
    3750:    used within time stepping implementations,
    3751:    so most users would not generally call this routine themselves.

    3753: .seealso: [](ch_ts), `TS`, `TSSetResize()`
    3754: @*/
    3755: static PetscErrorCode TSResizeRegisterOrRetrieve(TS ts, PetscBool flg)
    3756: {
    3757:   PetscFunctionBegin;
    3759:   PetscTryTypeMethod(ts, resizeregister, flg);
    3760:   /* PetscTryTypeMethod(adapt, resizeregister, flg); */
    3761:   PetscFunctionReturn(PETSC_SUCCESS);
    3762: }

    3764: static PetscErrorCode TSResizeReset(TS ts)
    3765: {
    3766:   PetscFunctionBegin;
    3768:   PetscCall(PetscObjectListDestroy(&ts->resizetransferobjs));
    3769:   PetscFunctionReturn(PETSC_SUCCESS);
    3770: }

    3772: static PetscErrorCode TSResizeTransferVecs(TS ts, PetscInt cnt, Vec vecsin[], Vec vecsout[])
    3773: {
    3774:   PetscFunctionBegin;
    3777:   for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPush(vecsin[i]));
    3778:   if (ts->resizetransfer) {
    3779:     PetscCall(PetscInfo(ts, "Transferring %" PetscInt_FMT " vectors\n", cnt));
    3780:     PetscCallBack("TS callback resize transfer", (*ts->resizetransfer)(ts, cnt, vecsin, vecsout, ts->resizectx));
    3781:   }
    3782:   for (PetscInt i = 0; i < cnt; i++) PetscCall(VecLockReadPop(vecsin[i]));
    3783:   PetscFunctionReturn(PETSC_SUCCESS);
    3784: }

    3786: /*@C
    3787:   TSResizeRegisterVec - Register a vector to be transferred with `TSResize()`.

    3789:   Collective

    3791:   Input Parameters:
    3792: + ts   - The `TS` context obtained from `TSCreate()`
    3793: . name - A string identifying the vector
    3794: - vec  - The vector

    3796:   Level: developer

    3798:   Note:
    3799:   `TSResizeRegisterVec()` is typically used within time stepping implementations,
    3800:   so most users would not generally call this routine themselves.

    3802: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRetrieveVec()`
    3803: @*/
    3804: PetscErrorCode TSResizeRegisterVec(TS ts, const char name[], Vec vec)
    3805: {
    3806:   PetscFunctionBegin;
    3808:   PetscAssertPointer(name, 2);
    3810:   PetscCall(PetscObjectListAdd(&ts->resizetransferobjs, name, (PetscObject)vec));
    3811:   PetscFunctionReturn(PETSC_SUCCESS);
    3812: }

    3814: /*@C
    3815:   TSResizeRetrieveVec - Retrieve a vector registered with `TSResizeRegisterVec()`.

    3817:   Collective

    3819:   Input Parameters:
    3820: + ts   - The `TS` context obtained from `TSCreate()`
    3821: . name - A string identifying the vector
    3822: - vec  - The vector

    3824:   Level: developer

    3826:   Note:
    3827:   `TSResizeRetrieveVec()` is typically used within time stepping implementations,
    3828:   so most users would not generally call this routine themselves.

    3830: .seealso: [](ch_ts), `TS`, `TSSetResize()`, `TSResize()`, `TSResizeRegisterVec()`
    3831: @*/
    3832: PetscErrorCode TSResizeRetrieveVec(TS ts, const char name[], Vec *vec)
    3833: {
    3834:   PetscFunctionBegin;
    3836:   PetscAssertPointer(name, 2);
    3837:   PetscAssertPointer(vec, 3);
    3838:   PetscCall(PetscObjectListFind(ts->resizetransferobjs, name, (PetscObject *)vec));
    3839:   PetscFunctionReturn(PETSC_SUCCESS);
    3840: }

    3842: static PetscErrorCode TSResizeGetVecArray(TS ts, PetscInt *nv, const char **names[], Vec *vecs[])
    3843: {
    3844:   PetscInt        cnt;
    3845:   PetscObjectList tmp;
    3846:   Vec            *vecsin  = NULL;
    3847:   const char    **namesin = NULL;

    3849:   PetscFunctionBegin;
    3850:   for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next)
    3851:     if (tmp->obj && tmp->obj->classid == VEC_CLASSID) cnt++;
    3852:   if (names) PetscCall(PetscMalloc1(cnt, &namesin));
    3853:   if (vecs) PetscCall(PetscMalloc1(cnt, &vecsin));
    3854:   for (tmp = ts->resizetransferobjs, cnt = 0; tmp; tmp = tmp->next) {
    3855:     if (tmp->obj && tmp->obj->classid == VEC_CLASSID) {
    3856:       if (vecs) vecsin[cnt] = (Vec)tmp->obj;
    3857:       if (names) namesin[cnt] = tmp->name;
    3858:       cnt++;
    3859:     }
    3860:   }
    3861:   if (nv) *nv = cnt;
    3862:   if (names) *names = namesin;
    3863:   if (vecs) *vecs = vecsin;
    3864:   PetscFunctionReturn(PETSC_SUCCESS);
    3865: }

    3867: /*@
    3868:   TSResize - Runs the user-defined transfer functions provided with `TSSetResize()`

    3870:   Collective

    3872:   Input Parameter:
    3873: . ts - The `TS` context obtained from `TSCreate()`

    3875:   Level: developer

    3877:   Note:
    3878:   `TSResize()` is typically used within time stepping implementations,
    3879:   so most users would not generally call this routine themselves.

    3881: .seealso: [](ch_ts), `TS`, `TSSetResize()`
    3882: @*/
    3883: PetscErrorCode TSResize(TS ts)
    3884: {
    3885:   PetscInt     nv      = 0;
    3886:   const char **names   = NULL;
    3887:   Vec         *vecsin  = NULL;
    3888:   const char  *solname = "ts:vec_sol";

    3890:   PetscFunctionBegin;
    3892:   if (!ts->resizesetup) PetscFunctionReturn(PETSC_SUCCESS);
    3893:   if (ts->resizesetup) {
    3894:     PetscCall(VecLockReadPush(ts->vec_sol));
    3895:     PetscCallBack("TS callback resize setup", (*ts->resizesetup)(ts, ts->steps, ts->ptime, ts->vec_sol, &ts->stepresize, ts->resizectx));
    3896:     PetscCall(VecLockReadPop(ts->vec_sol));
    3897:     if (ts->stepresize) {
    3898:       if (ts->resizerollback) {
    3899:         PetscCall(TSRollBack(ts));
    3900:         ts->time_step = ts->time_step0;
    3901:       }
    3902:       PetscCall(TSResizeRegisterVec(ts, solname, ts->vec_sol));
    3903:       PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_TRUE)); /* specific impls register their own objects */
    3904:     }
    3905:   }

    3907:   PetscCall(TSResizeGetVecArray(ts, &nv, &names, &vecsin));
    3908:   if (nv) {
    3909:     Vec *vecsout, vecsol;

    3911:     /* Reset internal objects */
    3912:     PetscCall(TSReset(ts));

    3914:     /* Transfer needed vectors (users can call SetJacobian, SetDM, etc. here) */
    3915:     PetscCall(PetscCalloc1(nv, &vecsout));
    3916:     PetscCall(TSResizeTransferVecs(ts, nv, vecsin, vecsout));
    3917:     for (PetscInt i = 0; i < nv; i++) {
    3918:       const char *name;
    3919:       char       *oname;

    3921:       PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
    3922:       PetscCall(PetscStrallocpy(name, &oname));
    3923:       PetscCall(TSResizeRegisterVec(ts, names[i], vecsout[i]));
    3924:       if (vecsout[i]) PetscCall(PetscObjectSetName((PetscObject)vecsout[i], oname));
    3925:       PetscCall(PetscFree(oname));
    3926:       PetscCall(VecDestroy(&vecsout[i]));
    3927:     }
    3928:     PetscCall(PetscFree(vecsout));
    3929:     PetscCall(TSResizeRegisterOrRetrieve(ts, PETSC_FALSE)); /* specific impls import the transferred objects */

    3931:     PetscCall(TSResizeRetrieveVec(ts, solname, &vecsol));
    3932:     if (vecsol) PetscCall(TSSetSolution(ts, vecsol));
    3933:     PetscAssert(ts->vec_sol, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_NULL, "Missing TS solution");
    3934:   }

    3936:   PetscCall(PetscFree(names));
    3937:   PetscCall(PetscFree(vecsin));
    3938:   PetscCall(TSResizeReset(ts));
    3939:   PetscFunctionReturn(PETSC_SUCCESS);
    3940: }

    3942: /*@
    3943:   TSSolve - Steps the requested number of timesteps.

    3945:   Collective

    3947:   Input Parameters:
    3948: + ts - the `TS` context obtained from `TSCreate()`
    3949: - u  - the solution vector  (can be null if `TSSetSolution()` was used and `TSSetExactFinalTime`(ts,`TS_EXACTFINALTIME_MATCHSTEP`) was not used,
    3950:        otherwise it must contain the initial conditions and will contain the solution at the final requested time

    3952:   Level: beginner

    3954:   Notes:
    3955:   The final time returned by this function may be different from the time of the internally
    3956:   held state accessible by `TSGetSolution()` and `TSGetTime()` because the method may have
    3957:   stepped over the final time.

    3959: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetSolution()`, `TSStep()`, `TSGetTime()`, `TSGetSolveTime()`
    3960: @*/
    3961: PetscErrorCode TSSolve(TS ts, Vec u)
    3962: {
    3963:   Vec solution;

    3965:   PetscFunctionBegin;

    3969:   PetscCall(TSSetExactFinalTimeDefault(ts));
    3970:   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 */
    3971:     if (!ts->vec_sol || u == ts->vec_sol) {
    3972:       PetscCall(VecDuplicate(u, &solution));
    3973:       PetscCall(TSSetSolution(ts, solution));
    3974:       PetscCall(VecDestroy(&solution)); /* grant ownership */
    3975:     }
    3976:     PetscCall(VecCopy(u, ts->vec_sol));
    3977:     PetscCheck(!ts->forward_solve, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Sensitivity analysis does not support the mode TS_EXACTFINALTIME_INTERPOLATE");
    3978:   } else if (u) PetscCall(TSSetSolution(ts, u));
    3979:   PetscCall(TSSetUp(ts));
    3980:   PetscCall(TSTrajectorySetUp(ts->trajectory, ts));

    3982:   PetscCheck(ts->max_time < PETSC_MAX_REAL || ts->max_steps != PETSC_INT_MAX, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "You must call TSSetMaxTime() or TSSetMaxSteps(), or use -ts_max_time <time> or -ts_max_steps <steps>");
    3983:   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()");
    3984:   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");
    3985:   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");

    3987:   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 */
    3988:     PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[0]));
    3989:     ts->tspan->spanctr = 1;
    3990:   }

    3992:   if (ts->forward_solve) PetscCall(TSForwardSetUp(ts));

    3994:   /* reset number of steps only when the step is not restarted. ARKIMEX
    3995:      restarts the step after an event. Resetting these counters in such case causes
    3996:      TSTrajectory to incorrectly save the output files
    3997:   */
    3998:   /* reset time step and iteration counters */
    3999:   if (!ts->steps) {
    4000:     ts->ksp_its           = 0;
    4001:     ts->snes_its          = 0;
    4002:     ts->num_snes_failures = 0;
    4003:     ts->reject            = 0;
    4004:     ts->steprestart       = PETSC_TRUE;
    4005:     ts->steprollback      = PETSC_FALSE;
    4006:     ts->stepresize        = PETSC_FALSE;
    4007:     ts->rhsjacobian.time  = PETSC_MIN_REAL;
    4008:   }

    4010:   /* make sure initial time step does not overshoot final time or the next point in tspan */
    4011:   if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
    4012:     PetscReal maxdt;
    4013:     PetscReal dt = ts->time_step;

    4015:     if (ts->tspan) maxdt = ts->tspan->span_times[ts->tspan->spanctr] - ts->ptime;
    4016:     else maxdt = ts->max_time - ts->ptime;
    4017:     ts->time_step = dt >= maxdt ? maxdt : (PetscIsCloseAtTol(dt, maxdt, 10 * PETSC_MACHINE_EPSILON, 0) ? maxdt : dt);
    4018:   }
    4019:   ts->reason = TS_CONVERGED_ITERATING;

    4021:   {
    4022:     PetscViewer       viewer;
    4023:     PetscViewerFormat format;
    4024:     PetscBool         flg;
    4025:     static PetscBool  incall = PETSC_FALSE;

    4027:     if (!incall) {
    4028:       /* Estimate the convergence rate of the time discretization */
    4029:       PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ts), ((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_estimate", &viewer, &format, &flg));
    4030:       if (flg) {
    4031:         PetscConvEst conv;
    4032:         DM           dm;
    4033:         PetscReal   *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */
    4034:         PetscInt     Nf;
    4035:         PetscBool    checkTemporal = PETSC_TRUE;

    4037:         incall = PETSC_TRUE;
    4038:         PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_convergence_temporal", &checkTemporal, &flg));
    4039:         PetscCall(TSGetDM(ts, &dm));
    4040:         PetscCall(DMGetNumFields(dm, &Nf));
    4041:         PetscCall(PetscCalloc1(PetscMax(Nf, 1), &alpha));
    4042:         PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)ts), &conv));
    4043:         PetscCall(PetscConvEstUseTS(conv, checkTemporal));
    4044:         PetscCall(PetscConvEstSetSolver(conv, (PetscObject)ts));
    4045:         PetscCall(PetscConvEstSetFromOptions(conv));
    4046:         PetscCall(PetscConvEstSetUp(conv));
    4047:         PetscCall(PetscConvEstGetConvRate(conv, alpha));
    4048:         PetscCall(PetscViewerPushFormat(viewer, format));
    4049:         PetscCall(PetscConvEstRateView(conv, alpha, viewer));
    4050:         PetscCall(PetscViewerPopFormat(viewer));
    4051:         PetscCall(PetscViewerDestroy(&viewer));
    4052:         PetscCall(PetscConvEstDestroy(&conv));
    4053:         PetscCall(PetscFree(alpha));
    4054:         incall = PETSC_FALSE;
    4055:       }
    4056:     }
    4057:   }

    4059:   PetscCall(TSViewFromOptions(ts, NULL, "-ts_view_pre"));

    4061:   if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */
    4062:     PetscUseTypeMethod(ts, solve);
    4063:     if (u) PetscCall(VecCopy(ts->vec_sol, u));
    4064:     ts->solvetime = ts->ptime;
    4065:     solution      = ts->vec_sol;
    4066:   } else { /* Step the requested number of timesteps. */
    4067:     if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
    4068:     else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;

    4070:     if (!ts->steps) {
    4071:       PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
    4072:       PetscCall(TSEventInitialize(ts->event, ts, ts->ptime, ts->vec_sol));
    4073:     }

    4075:     while (!ts->reason) {
    4076:       PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));
    4077:       if (!ts->steprollback || (ts->stepresize && ts->resizerollback)) PetscCall(TSPreStep(ts));
    4078:       PetscCall(TSStep(ts));
    4079:       if (ts->testjacobian) PetscCall(TSRHSJacobianTest(ts, NULL));
    4080:       if (ts->testjacobiantranspose) PetscCall(TSRHSJacobianTestTranspose(ts, NULL));
    4081:       if (ts->quadraturets && ts->costintegralfwd) { /* Must evaluate the cost integral before event is handled. The cost integral value can also be rolled back. */
    4082:         if (ts->reason >= 0) ts->steps--;            /* Revert the step number changed by TSStep() */
    4083:         PetscCall(TSForwardCostIntegral(ts));
    4084:         if (ts->reason >= 0) ts->steps++;
    4085:       }
    4086:       if (ts->forward_solve) {            /* compute forward sensitivities before event handling because postevent() may change RHS and jump conditions may have to be applied */
    4087:         if (ts->reason >= 0) ts->steps--; /* Revert the step number changed by TSStep() */
    4088:         PetscCall(TSForwardStep(ts));
    4089:         if (ts->reason >= 0) ts->steps++;
    4090:       }
    4091:       PetscCall(TSPostEvaluate(ts));
    4092:       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. */
    4093:       if (ts->steprollback) PetscCall(TSPostEvaluate(ts));
    4094:       if (!ts->steprollback && ts->resizerollback) PetscCall(TSResize(ts));
    4095:       /* check convergence */
    4096:       if (!ts->reason) {
    4097:         if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS;
    4098:         else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME;
    4099:       }
    4100:       if (!ts->steprollback) {
    4101:         PetscCall(TSTrajectorySet(ts->trajectory, ts, ts->steps, ts->ptime, ts->vec_sol));
    4102:         PetscCall(TSPostStep(ts));
    4103:         if (!ts->resizerollback) PetscCall(TSResize(ts));

    4105:         if (ts->tspan && ts->tspan->spanctr < ts->tspan->num_span_times) {
    4106:           PetscCheck(ts->tspan->worktol > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Unexpected state !(tspan->worktol > 0) in TSSolve()");
    4107:           if (PetscIsCloseAtTol(ts->ptime, ts->tspan->span_times[ts->tspan->spanctr], ts->tspan->worktol, 0)) PetscCall(VecCopy(ts->vec_sol, ts->tspan->vecs_sol[ts->tspan->spanctr++]));
    4108:         }
    4109:       }
    4110:     }
    4111:     PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));

    4113:     if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) {
    4114:       if (!u) u = ts->vec_sol;
    4115:       PetscCall(TSInterpolate(ts, ts->max_time, u));
    4116:       ts->solvetime = ts->max_time;
    4117:       solution      = u;
    4118:       PetscCall(TSMonitor(ts, -1, ts->solvetime, solution));
    4119:     } else {
    4120:       if (u) PetscCall(VecCopy(ts->vec_sol, u));
    4121:       ts->solvetime = ts->ptime;
    4122:       solution      = ts->vec_sol;
    4123:     }
    4124:   }

    4126:   PetscCall(TSViewFromOptions(ts, NULL, "-ts_view"));
    4127:   PetscCall(VecViewFromOptions(solution, (PetscObject)ts, "-ts_view_solution"));
    4128:   PetscCall(PetscObjectSAWsBlock((PetscObject)ts));
    4129:   if (ts->adjoint_solve) PetscCall(TSAdjointSolve(ts));
    4130:   PetscFunctionReturn(PETSC_SUCCESS);
    4131: }

    4133: /*@
    4134:   TSGetTime - Gets the time of the most recently completed step.

    4136:   Not Collective

    4138:   Input Parameter:
    4139: . ts - the `TS` context obtained from `TSCreate()`

    4141:   Output Parameter:
    4142: . t - the current time. This time may not corresponds to the final time set with `TSSetMaxTime()`, use `TSGetSolveTime()`.

    4144:   Level: beginner

    4146:   Note:
    4147:   When called during time step evaluation (e.g. during residual evaluation or via hooks set using `TSSetPreStep()`,
    4148:   `TSSetPreStage()`, `TSSetPostStage()`, or `TSSetPostStep()`), the time is the time at the start of the step being evaluated.

    4150: .seealso: [](ch_ts), `TS`, `TSGetSolveTime()`, `TSSetTime()`, `TSGetTimeStep()`, `TSGetStepNumber()`
    4151: @*/
    4152: PetscErrorCode TSGetTime(TS ts, PetscReal *t)
    4153: {
    4154:   PetscFunctionBegin;
    4156:   PetscAssertPointer(t, 2);
    4157:   *t = ts->ptime;
    4158:   PetscFunctionReturn(PETSC_SUCCESS);
    4159: }

    4161: /*@
    4162:   TSGetPrevTime - Gets the starting time of the previously completed step.

    4164:   Not Collective

    4166:   Input Parameter:
    4167: . ts - the `TS` context obtained from `TSCreate()`

    4169:   Output Parameter:
    4170: . t - the previous time

    4172:   Level: beginner

    4174: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSGetSolveTime()`, `TSGetTimeStep()`
    4175: @*/
    4176: PetscErrorCode TSGetPrevTime(TS ts, PetscReal *t)
    4177: {
    4178:   PetscFunctionBegin;
    4180:   PetscAssertPointer(t, 2);
    4181:   *t = ts->ptime_prev;
    4182:   PetscFunctionReturn(PETSC_SUCCESS);
    4183: }

    4185: /*@
    4186:   TSSetTime - Allows one to reset the time.

    4188:   Logically Collective

    4190:   Input Parameters:
    4191: + ts - the `TS` context obtained from `TSCreate()`
    4192: - t  - the time

    4194:   Level: intermediate

    4196: .seealso: [](ch_ts), `TS`, `TSGetTime()`, `TSSetMaxSteps()`
    4197: @*/
    4198: PetscErrorCode TSSetTime(TS ts, PetscReal t)
    4199: {
    4200:   PetscFunctionBegin;
    4203:   ts->ptime = t;
    4204:   PetscFunctionReturn(PETSC_SUCCESS);
    4205: }

    4207: /*@
    4208:   TSSetOptionsPrefix - Sets the prefix used for searching for all
    4209:   TS options in the database.

    4211:   Logically Collective

    4213:   Input Parameters:
    4214: + ts     - The `TS` context
    4215: - prefix - The prefix to prepend to all option names

    4217:   Level: advanced

    4219:   Note:
    4220:   A hyphen (-) must NOT be given at the beginning of the prefix name.
    4221:   The first character of all runtime options is AUTOMATICALLY the
    4222:   hyphen.

    4224: .seealso: [](ch_ts), `TS`, `TSSetFromOptions()`, `TSAppendOptionsPrefix()`
    4225: @*/
    4226: PetscErrorCode TSSetOptionsPrefix(TS ts, const char prefix[])
    4227: {
    4228:   SNES snes;

    4230:   PetscFunctionBegin;
    4232:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ts, prefix));
    4233:   PetscCall(TSGetSNES(ts, &snes));
    4234:   PetscCall(SNESSetOptionsPrefix(snes, prefix));
    4235:   PetscFunctionReturn(PETSC_SUCCESS);
    4236: }

    4238: /*@
    4239:   TSAppendOptionsPrefix - Appends to the prefix used for searching for all
    4240:   TS options in the database.

    4242:   Logically Collective

    4244:   Input Parameters:
    4245: + ts     - The `TS` context
    4246: - prefix - The prefix to prepend to all option names

    4248:   Level: advanced

    4250:   Note:
    4251:   A hyphen (-) must NOT be given at the beginning of the prefix name.
    4252:   The first character of all runtime options is AUTOMATICALLY the
    4253:   hyphen.

    4255: .seealso: [](ch_ts), `TS`, `TSGetOptionsPrefix()`, `TSSetOptionsPrefix()`, `TSSetFromOptions()`
    4256: @*/
    4257: PetscErrorCode TSAppendOptionsPrefix(TS ts, const char prefix[])
    4258: {
    4259:   SNES snes;

    4261:   PetscFunctionBegin;
    4263:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)ts, prefix));
    4264:   PetscCall(TSGetSNES(ts, &snes));
    4265:   PetscCall(SNESAppendOptionsPrefix(snes, prefix));
    4266:   PetscFunctionReturn(PETSC_SUCCESS);
    4267: }

    4269: /*@
    4270:   TSGetOptionsPrefix - Sets the prefix used for searching for all
    4271:   `TS` options in the database.

    4273:   Not Collective

    4275:   Input Parameter:
    4276: . ts - The `TS` context

    4278:   Output Parameter:
    4279: . prefix - A pointer to the prefix string used

    4281:   Level: intermediate

    4283:   Fortran Notes:
    4284:   The user should pass in a string 'prefix' of
    4285:   sufficient length to hold the prefix.

    4287: .seealso: [](ch_ts), `TS`, `TSAppendOptionsPrefix()`, `TSSetFromOptions()`
    4288: @*/
    4289: PetscErrorCode TSGetOptionsPrefix(TS ts, const char *prefix[])
    4290: {
    4291:   PetscFunctionBegin;
    4293:   PetscAssertPointer(prefix, 2);
    4294:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ts, prefix));
    4295:   PetscFunctionReturn(PETSC_SUCCESS);
    4296: }

    4298: /*@C
    4299:   TSGetRHSJacobian - Returns the Jacobian J at the present timestep.

    4301:   Not Collective, but parallel objects are returned if ts is parallel

    4303:   Input Parameter:
    4304: . ts - The `TS` context obtained from `TSCreate()`

    4306:   Output Parameters:
    4307: + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t)  (or `NULL`)
    4308: . Pmat - The matrix from which the preconditioner is constructed, usually the same as `Amat`  (or `NULL`)
    4309: . func - Function to compute the Jacobian of the RHS  (or `NULL`)
    4310: - ctx  - User-defined context for Jacobian evaluation routine  (or `NULL`)

    4312:   Level: intermediate

    4314:   Note:
    4315:   You can pass in `NULL` for any return argument you do not need.

    4317: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`

    4319: @*/
    4320: PetscErrorCode TSGetRHSJacobian(TS ts, Mat *Amat, Mat *Pmat, TSRHSJacobianFn **func, void **ctx)
    4321: {
    4322:   DM dm;

    4324:   PetscFunctionBegin;
    4325:   if (Amat || Pmat) {
    4326:     SNES snes;
    4327:     PetscCall(TSGetSNES(ts, &snes));
    4328:     PetscCall(SNESSetUpMatrices(snes));
    4329:     PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
    4330:   }
    4331:   PetscCall(TSGetDM(ts, &dm));
    4332:   PetscCall(DMTSGetRHSJacobian(dm, func, ctx));
    4333:   PetscFunctionReturn(PETSC_SUCCESS);
    4334: }

    4336: /*@C
    4337:   TSGetIJacobian - Returns the implicit Jacobian at the present timestep.

    4339:   Not Collective, but parallel objects are returned if ts is parallel

    4341:   Input Parameter:
    4342: . ts - The `TS` context obtained from `TSCreate()`

    4344:   Output Parameters:
    4345: + Amat - The (approximate) Jacobian of F(t,U,U_t)
    4346: . Pmat - The matrix from which the preconditioner is constructed, often the same as `Amat`
    4347: . f    - The function to compute the matrices
    4348: - ctx  - User-defined context for Jacobian evaluation routine

    4350:   Level: advanced

    4352:   Note:
    4353:   You can pass in `NULL` for any return argument you do not need.

    4355: .seealso: [](ch_ts), `TS`, `TSGetTimeStep()`, `TSGetRHSJacobian()`, `TSGetMatrices()`, `TSGetTime()`, `TSGetStepNumber()`
    4356: @*/
    4357: PetscErrorCode TSGetIJacobian(TS ts, Mat *Amat, Mat *Pmat, TSIJacobianFn **f, void **ctx)
    4358: {
    4359:   DM dm;

    4361:   PetscFunctionBegin;
    4362:   if (Amat || Pmat) {
    4363:     SNES snes;
    4364:     PetscCall(TSGetSNES(ts, &snes));
    4365:     PetscCall(SNESSetUpMatrices(snes));
    4366:     PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL));
    4367:   }
    4368:   PetscCall(TSGetDM(ts, &dm));
    4369:   PetscCall(DMTSGetIJacobian(dm, f, ctx));
    4370:   PetscFunctionReturn(PETSC_SUCCESS);
    4371: }

    4373: #include <petsc/private/dmimpl.h>
    4374: /*@
    4375:   TSSetDM - Sets the `DM` that may be used by some nonlinear solvers or preconditioners under the `TS`

    4377:   Logically Collective

    4379:   Input Parameters:
    4380: + ts - the `TS` integrator object
    4381: - dm - the dm, cannot be `NULL`

    4383:   Level: intermediate

    4385:   Notes:
    4386:   A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`,
    4387:   even when not using interfaces like `DMTSSetIFunction()`.  Use `DMClone()` to get a distinct `DM` when solving
    4388:   different problems using the same function space.

    4390: .seealso: [](ch_ts), `TS`, `DM`, `TSGetDM()`, `SNESSetDM()`, `SNESGetDM()`
    4391: @*/
    4392: PetscErrorCode TSSetDM(TS ts, DM dm)
    4393: {
    4394:   SNES snes;
    4395:   DMTS tsdm;

    4397:   PetscFunctionBegin;
    4400:   PetscCall(PetscObjectReference((PetscObject)dm));
    4401:   if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */
    4402:     if (ts->dm->dmts && !dm->dmts) {
    4403:       PetscCall(DMCopyDMTS(ts->dm, dm));
    4404:       PetscCall(DMGetDMTS(ts->dm, &tsdm));
    4405:       /* Grant write privileges to the replacement DM */
    4406:       if (tsdm->originaldm == ts->dm) tsdm->originaldm = dm;
    4407:     }
    4408:     PetscCall(DMDestroy(&ts->dm));
    4409:   }
    4410:   ts->dm = dm;

    4412:   PetscCall(TSGetSNES(ts, &snes));
    4413:   PetscCall(SNESSetDM(snes, dm));
    4414:   PetscFunctionReturn(PETSC_SUCCESS);
    4415: }

    4417: /*@
    4418:   TSGetDM - Gets the `DM` that may be used by some preconditioners

    4420:   Not Collective

    4422:   Input Parameter:
    4423: . ts - the `TS`

    4425:   Output Parameter:
    4426: . dm - the `DM`

    4428:   Level: intermediate

    4430: .seealso: [](ch_ts), `TS`, `DM`, `TSSetDM()`, `SNESSetDM()`, `SNESGetDM()`
    4431: @*/
    4432: PetscErrorCode TSGetDM(TS ts, DM *dm)
    4433: {
    4434:   PetscFunctionBegin;
    4436:   if (!ts->dm) {
    4437:     PetscCall(DMShellCreate(PetscObjectComm((PetscObject)ts), &ts->dm));
    4438:     if (ts->snes) PetscCall(SNESSetDM(ts->snes, ts->dm));
    4439:   }
    4440:   *dm = ts->dm;
    4441:   PetscFunctionReturn(PETSC_SUCCESS);
    4442: }

    4444: /*@
    4445:   SNESTSFormFunction - Function to evaluate nonlinear residual

    4447:   Logically Collective

    4449:   Input Parameters:
    4450: + snes - nonlinear solver
    4451: . U    - the current state at which to evaluate the residual
    4452: - ctx  - user context, must be a TS

    4454:   Output Parameter:
    4455: . F - the nonlinear residual

    4457:   Level: advanced

    4459:   Note:
    4460:   This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.
    4461:   It is most frequently passed to `MatFDColoringSetFunction()`.

    4463: .seealso: [](ch_ts), `SNESSetFunction()`, `MatFDColoringSetFunction()`
    4464: @*/
    4465: PetscErrorCode SNESTSFormFunction(SNES snes, Vec U, Vec F, void *ctx)
    4466: {
    4467:   TS ts = (TS)ctx;

    4469:   PetscFunctionBegin;
    4474:   PetscCheck(ts->ops->snesfunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesfunction for TS of type %s", ((PetscObject)ts)->type_name);
    4475:   PetscCall((*ts->ops->snesfunction)(snes, U, F, ts));
    4476:   PetscFunctionReturn(PETSC_SUCCESS);
    4477: }

    4479: /*@
    4480:   SNESTSFormJacobian - Function to evaluate the Jacobian

    4482:   Collective

    4484:   Input Parameters:
    4485: + snes - nonlinear solver
    4486: . U    - the current state at which to evaluate the residual
    4487: - ctx  - user context, must be a `TS`

    4489:   Output Parameters:
    4490: + A - the Jacobian
    4491: - B - the preconditioning matrix (may be the same as A)

    4493:   Level: developer

    4495:   Note:
    4496:   This function is not normally called by users and is automatically registered with the `SNES` used by `TS`.

    4498: .seealso: [](ch_ts), `SNESSetJacobian()`
    4499: @*/
    4500: PetscErrorCode SNESTSFormJacobian(SNES snes, Vec U, Mat A, Mat B, void *ctx)
    4501: {
    4502:   TS ts = (TS)ctx;

    4504:   PetscFunctionBegin;
    4510:   PetscCheck(ts->ops->snesjacobian, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No method snesjacobian for TS of type %s", ((PetscObject)ts)->type_name);
    4511:   PetscCall((*ts->ops->snesjacobian)(snes, U, A, B, ts));
    4512:   PetscFunctionReturn(PETSC_SUCCESS);
    4513: }

    4515: /*@C
    4516:   TSComputeRHSFunctionLinear - Evaluate the right-hand side via the user-provided Jacobian, for linear problems Udot = A U only

    4518:   Collective

    4520:   Input Parameters:
    4521: + ts  - time stepping context
    4522: . t   - time at which to evaluate
    4523: . U   - state at which to evaluate
    4524: - ctx - context

    4526:   Output Parameter:
    4527: . F - right-hand side

    4529:   Level: intermediate

    4531:   Note:
    4532:   This function is intended to be passed to `TSSetRHSFunction()` to evaluate the right-hand side for linear problems.
    4533:   The matrix (and optionally the evaluation context) should be passed to `TSSetRHSJacobian()`.

    4535: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSJacobianConstant()`
    4536: @*/
    4537: PetscErrorCode TSComputeRHSFunctionLinear(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
    4538: {
    4539:   Mat Arhs, Brhs;

    4541:   PetscFunctionBegin;
    4542:   PetscCall(TSGetRHSMats_Private(ts, &Arhs, &Brhs));
    4543:   /* undo the damage caused by shifting */
    4544:   PetscCall(TSRecoverRHSJacobian(ts, Arhs, Brhs));
    4545:   PetscCall(TSComputeRHSJacobian(ts, t, U, Arhs, Brhs));
    4546:   PetscCall(MatMult(Arhs, U, F));
    4547:   PetscFunctionReturn(PETSC_SUCCESS);
    4548: }

    4550: /*@C
    4551:   TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent.

    4553:   Collective

    4555:   Input Parameters:
    4556: + ts  - time stepping context
    4557: . t   - time at which to evaluate
    4558: . U   - state at which to evaluate
    4559: - ctx - context

    4561:   Output Parameters:
    4562: + A - pointer to operator
    4563: - B - pointer to preconditioning matrix

    4565:   Level: intermediate

    4567:   Note:
    4568:   This function is intended to be passed to `TSSetRHSJacobian()` to evaluate the Jacobian for linear time-independent problems.

    4570: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSSetRHSJacobian()`, `TSComputeRHSFunctionLinear()`
    4571: @*/
    4572: PetscErrorCode TSComputeRHSJacobianConstant(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
    4573: {
    4574:   PetscFunctionBegin;
    4575:   PetscFunctionReturn(PETSC_SUCCESS);
    4576: }

    4578: /*@C
    4579:   TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only

    4581:   Collective

    4583:   Input Parameters:
    4584: + ts   - time stepping context
    4585: . t    - time at which to evaluate
    4586: . U    - state at which to evaluate
    4587: . Udot - time derivative of state vector
    4588: - ctx  - context

    4590:   Output Parameter:
    4591: . F - left hand side

    4593:   Level: intermediate

    4595:   Notes:
    4596:   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
    4597:   user is required to write their own `TSComputeIFunction()`.
    4598:   This function is intended to be passed to `TSSetIFunction()` to evaluate the left hand side for linear problems.
    4599:   The matrix (and optionally the evaluation context) should be passed to `TSSetIJacobian()`.

    4601:   Note that using this function is NOT equivalent to using `TSComputeRHSFunctionLinear()` since that solves Udot = A U

    4603: .seealso: [](ch_ts), `TS`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIJacobianConstant()`, `TSComputeRHSFunctionLinear()`
    4604: @*/
    4605: PetscErrorCode TSComputeIFunctionLinear(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
    4606: {
    4607:   Mat A, B;

    4609:   PetscFunctionBegin;
    4610:   PetscCall(TSGetIJacobian(ts, &A, &B, NULL, NULL));
    4611:   PetscCall(TSComputeIJacobian(ts, t, U, Udot, 1.0, A, B, PETSC_TRUE));
    4612:   PetscCall(MatMult(A, Udot, F));
    4613:   PetscFunctionReturn(PETSC_SUCCESS);
    4614: }

    4616: /*@C
    4617:   TSComputeIJacobianConstant - Reuses the matrix previously computed with the provided `TSIJacobianFn` for a semi-implicit DAE or ODE

    4619:   Collective

    4621:   Input Parameters:
    4622: + ts    - time stepping context
    4623: . t     - time at which to evaluate
    4624: . U     - state at which to evaluate
    4625: . Udot  - time derivative of state vector
    4626: . shift - shift to apply
    4627: - ctx   - context

    4629:   Output Parameters:
    4630: + A - pointer to operator
    4631: - B - pointer to matrix from which the preconditioner is built (often `A`)

    4633:   Level: advanced

    4635:   Notes:
    4636:   This function is intended to be passed to `TSSetIJacobian()` to evaluate the Jacobian for linear time-independent problems.

    4638:   It is only appropriate for problems of the form

    4640:   $$
    4641:   M \dot{U} = F(U,t)
    4642:   $$

    4644:   where M is constant and F is non-stiff.  The user must pass M to `TSSetIJacobian()`.  The current implementation only
    4645:   works with IMEX time integration methods such as `TSROSW` and `TSARKIMEX`, since there is no support for de-constructing
    4646:   an implicit operator of the form

    4648:   $$
    4649:   shift*M + J
    4650:   $$

    4652:   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
    4653:   a copy of M or reassemble it when requested.

    4655: .seealso: [](ch_ts), `TS`, `TSROSW`, `TSARKIMEX`, `TSSetIFunction()`, `TSSetIJacobian()`, `TSComputeIFunctionLinear()`
    4656: @*/
    4657: PetscErrorCode TSComputeIJacobianConstant(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat A, Mat B, void *ctx)
    4658: {
    4659:   PetscFunctionBegin;
    4660:   PetscCall(MatScale(A, shift / ts->ijacobian.shift));
    4661:   ts->ijacobian.shift = shift;
    4662:   PetscFunctionReturn(PETSC_SUCCESS);
    4663: }

    4665: /*@
    4666:   TSGetEquationType - Gets the type of the equation that `TS` is solving.

    4668:   Not Collective

    4670:   Input Parameter:
    4671: . ts - the `TS` context

    4673:   Output Parameter:
    4674: . equation_type - see `TSEquationType`

    4676:   Level: beginner

    4678: .seealso: [](ch_ts), `TS`, `TSSetEquationType()`, `TSEquationType`
    4679: @*/
    4680: PetscErrorCode TSGetEquationType(TS ts, TSEquationType *equation_type)
    4681: {
    4682:   PetscFunctionBegin;
    4684:   PetscAssertPointer(equation_type, 2);
    4685:   *equation_type = ts->equation_type;
    4686:   PetscFunctionReturn(PETSC_SUCCESS);
    4687: }

    4689: /*@
    4690:   TSSetEquationType - Sets the type of the equation that `TS` is solving.

    4692:   Not Collective

    4694:   Input Parameters:
    4695: + ts            - the `TS` context
    4696: - equation_type - see `TSEquationType`

    4698:   Level: advanced

    4700: .seealso: [](ch_ts), `TS`, `TSGetEquationType()`, `TSEquationType`
    4701: @*/
    4702: PetscErrorCode TSSetEquationType(TS ts, TSEquationType equation_type)
    4703: {
    4704:   PetscFunctionBegin;
    4706:   ts->equation_type = equation_type;
    4707:   PetscFunctionReturn(PETSC_SUCCESS);
    4708: }

    4710: /*@
    4711:   TSGetConvergedReason - Gets the reason the `TS` iteration was stopped.

    4713:   Not Collective

    4715:   Input Parameter:
    4716: . ts - the `TS` context

    4718:   Output Parameter:
    4719: . reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
    4720:             manual pages for the individual convergence tests for complete lists

    4722:   Level: beginner

    4724:   Note:
    4725:   Can only be called after the call to `TSSolve()` is complete.

    4727: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
    4728: @*/
    4729: PetscErrorCode TSGetConvergedReason(TS ts, TSConvergedReason *reason)
    4730: {
    4731:   PetscFunctionBegin;
    4733:   PetscAssertPointer(reason, 2);
    4734:   *reason = ts->reason;
    4735:   PetscFunctionReturn(PETSC_SUCCESS);
    4736: }

    4738: /*@
    4739:   TSSetConvergedReason - Sets the reason for handling the convergence of `TSSolve()`.

    4741:   Logically Collective; reason must contain common value

    4743:   Input Parameters:
    4744: + ts     - the `TS` context
    4745: - reason - negative value indicates diverged, positive value converged, see `TSConvergedReason` or the
    4746:             manual pages for the individual convergence tests for complete lists

    4748:   Level: advanced

    4750:   Note:
    4751:   Can only be called while `TSSolve()` is active.

    4753: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
    4754: @*/
    4755: PetscErrorCode TSSetConvergedReason(TS ts, TSConvergedReason reason)
    4756: {
    4757:   PetscFunctionBegin;
    4759:   ts->reason = reason;
    4760:   PetscFunctionReturn(PETSC_SUCCESS);
    4761: }

    4763: /*@
    4764:   TSGetSolveTime - Gets the time after a call to `TSSolve()`

    4766:   Not Collective

    4768:   Input Parameter:
    4769: . ts - the `TS` context

    4771:   Output Parameter:
    4772: . ftime - the final time. This time corresponds to the final time set with `TSSetMaxTime()`

    4774:   Level: beginner

    4776:   Note:
    4777:   Can only be called after the call to `TSSolve()` is complete.

    4779: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSConvergedReason`
    4780: @*/
    4781: PetscErrorCode TSGetSolveTime(TS ts, PetscReal *ftime)
    4782: {
    4783:   PetscFunctionBegin;
    4785:   PetscAssertPointer(ftime, 2);
    4786:   *ftime = ts->solvetime;
    4787:   PetscFunctionReturn(PETSC_SUCCESS);
    4788: }

    4790: /*@
    4791:   TSGetSNESIterations - Gets the total number of nonlinear iterations
    4792:   used by the time integrator.

    4794:   Not Collective

    4796:   Input Parameter:
    4797: . ts - `TS` context

    4799:   Output Parameter:
    4800: . nits - number of nonlinear iterations

    4802:   Level: intermediate

    4804:   Note:
    4805:   This counter is reset to zero for each successive call to `TSSolve()`.

    4807: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetKSPIterations()`
    4808: @*/
    4809: PetscErrorCode TSGetSNESIterations(TS ts, PetscInt *nits)
    4810: {
    4811:   PetscFunctionBegin;
    4813:   PetscAssertPointer(nits, 2);
    4814:   *nits = ts->snes_its;
    4815:   PetscFunctionReturn(PETSC_SUCCESS);
    4816: }

    4818: /*@
    4819:   TSGetKSPIterations - Gets the total number of linear iterations
    4820:   used by the time integrator.

    4822:   Not Collective

    4824:   Input Parameter:
    4825: . ts - `TS` context

    4827:   Output Parameter:
    4828: . lits - number of linear iterations

    4830:   Level: intermediate

    4832:   Note:
    4833:   This counter is reset to zero for each successive call to `TSSolve()`.

    4835: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `SNESGetKSPIterations()`
    4836: @*/
    4837: PetscErrorCode TSGetKSPIterations(TS ts, PetscInt *lits)
    4838: {
    4839:   PetscFunctionBegin;
    4841:   PetscAssertPointer(lits, 2);
    4842:   *lits = ts->ksp_its;
    4843:   PetscFunctionReturn(PETSC_SUCCESS);
    4844: }

    4846: /*@
    4847:   TSGetStepRejections - Gets the total number of rejected steps.

    4849:   Not Collective

    4851:   Input Parameter:
    4852: . ts - `TS` context

    4854:   Output Parameter:
    4855: . rejects - number of steps rejected

    4857:   Level: intermediate

    4859:   Note:
    4860:   This counter is reset to zero for each successive call to `TSSolve()`.

    4862: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetSNESFailures()`, `TSSetMaxSNESFailures()`, `TSSetErrorIfStepFails()`
    4863: @*/
    4864: PetscErrorCode TSGetStepRejections(TS ts, PetscInt *rejects)
    4865: {
    4866:   PetscFunctionBegin;
    4868:   PetscAssertPointer(rejects, 2);
    4869:   *rejects = ts->reject;
    4870:   PetscFunctionReturn(PETSC_SUCCESS);
    4871: }

    4873: /*@
    4874:   TSGetSNESFailures - Gets the total number of failed `SNES` solves in a `TS`

    4876:   Not Collective

    4878:   Input Parameter:
    4879: . ts - `TS` context

    4881:   Output Parameter:
    4882: . fails - number of failed nonlinear solves

    4884:   Level: intermediate

    4886:   Note:
    4887:   This counter is reset to zero for each successive call to `TSSolve()`.

    4889: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSSetMaxSNESFailures()`
    4890: @*/
    4891: PetscErrorCode TSGetSNESFailures(TS ts, PetscInt *fails)
    4892: {
    4893:   PetscFunctionBegin;
    4895:   PetscAssertPointer(fails, 2);
    4896:   *fails = ts->num_snes_failures;
    4897:   PetscFunctionReturn(PETSC_SUCCESS);
    4898: }

    4900: /*@
    4901:   TSSetMaxStepRejections - Sets the maximum number of step rejections before a time step fails

    4903:   Not Collective

    4905:   Input Parameters:
    4906: + ts      - `TS` context
    4907: - rejects - maximum number of rejected steps, pass `PETSC_UNLIMITED` for unlimited

    4909:   Options Database Key:
    4910: . -ts_max_reject - Maximum number of step rejections before a step fails

    4912:   Level: intermediate

    4914:   Developer Note:
    4915:   The options database name is incorrect.

    4917: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxSNESFailures()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSSetErrorIfStepFails()`, `TSGetConvergedReason()`
    4918: @*/
    4919: PetscErrorCode TSSetMaxStepRejections(TS ts, PetscInt rejects)
    4920: {
    4921:   PetscFunctionBegin;
    4923:   if (rejects == PETSC_UNLIMITED || rejects == -1) {
    4924:     ts->max_reject = PETSC_UNLIMITED;
    4925:   } else {
    4926:     PetscCheck(rejects >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of rejections");
    4927:     ts->max_reject = rejects;
    4928:   }
    4929:   PetscFunctionReturn(PETSC_SUCCESS);
    4930: }

    4932: /*@
    4933:   TSSetMaxSNESFailures - Sets the maximum number of failed `SNES` solves

    4935:   Not Collective

    4937:   Input Parameters:
    4938: + ts    - `TS` context
    4939: - fails - maximum number of failed nonlinear solves, pass `PETSC_UNLIMITED` to allow any number of failures.

    4941:   Options Database Key:
    4942: . -ts_max_snes_failures - Maximum number of nonlinear solve failures

    4944:   Level: intermediate

    4946: .seealso: [](ch_ts), `TS`, `SNES`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `SNESGetConvergedReason()`, `TSGetConvergedReason()`
    4947: @*/
    4948: PetscErrorCode TSSetMaxSNESFailures(TS ts, PetscInt fails)
    4949: {
    4950:   PetscFunctionBegin;
    4952:   if (fails == PETSC_UNLIMITED || fails == -1) {
    4953:     ts->max_snes_failures = PETSC_UNLIMITED;
    4954:   } else {
    4955:     PetscCheck(fails >= 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures");
    4956:     ts->max_snes_failures = fails;
    4957:   }
    4958:   PetscFunctionReturn(PETSC_SUCCESS);
    4959: }

    4961: /*@
    4962:   TSSetErrorIfStepFails - Immediately error if no step succeeds during `TSSolve()`

    4964:   Not Collective

    4966:   Input Parameters:
    4967: + ts  - `TS` context
    4968: - err - `PETSC_TRUE` to error if no step succeeds, `PETSC_FALSE` to return without failure

    4970:   Options Database Key:
    4971: . -ts_error_if_step_fails - Error if no step succeeds

    4973:   Level: intermediate

    4975: .seealso: [](ch_ts), `TS`, `TSGetSNESIterations()`, `TSGetKSPIterations()`, `TSSetMaxStepRejections()`, `TSGetStepRejections()`, `TSGetSNESFailures()`, `TSGetConvergedReason()`
    4976: @*/
    4977: PetscErrorCode TSSetErrorIfStepFails(TS ts, PetscBool err)
    4978: {
    4979:   PetscFunctionBegin;
    4981:   ts->errorifstepfailed = err;
    4982:   PetscFunctionReturn(PETSC_SUCCESS);
    4983: }

    4985: /*@
    4986:   TSGetAdapt - Get the adaptive controller context for the current method

    4988:   Collective if controller has not yet been created

    4990:   Input Parameter:
    4991: . ts - time stepping context

    4993:   Output Parameter:
    4994: . adapt - adaptive controller

    4996:   Level: intermediate

    4998: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSAdaptSetType()`, `TSAdaptChoose()`
    4999: @*/
    5000: PetscErrorCode TSGetAdapt(TS ts, TSAdapt *adapt)
    5001: {
    5002:   PetscFunctionBegin;
    5004:   PetscAssertPointer(adapt, 2);
    5005:   if (!ts->adapt) {
    5006:     PetscCall(TSAdaptCreate(PetscObjectComm((PetscObject)ts), &ts->adapt));
    5007:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ts->adapt, (PetscObject)ts, 1));
    5008:   }
    5009:   *adapt = ts->adapt;
    5010:   PetscFunctionReturn(PETSC_SUCCESS);
    5011: }

    5013: /*@
    5014:   TSSetTolerances - Set tolerances for local truncation error when using an adaptive controller

    5016:   Logically Collective

    5018:   Input Parameters:
    5019: + ts    - time integration context
    5020: . atol  - scalar absolute tolerances
    5021: . vatol - vector of absolute tolerances or `NULL`, used in preference to `atol` if present
    5022: . rtol  - scalar relative tolerances
    5023: - vrtol - vector of relative tolerances or `NULL`, used in preference to `rtol` if present

    5025:   Options Database Keys:
    5026: + -ts_rtol <rtol> - relative tolerance for local truncation error
    5027: - -ts_atol <atol> - Absolute tolerance for local truncation error

    5029:   Level: beginner

    5031:   Notes:
    5032:   `PETSC_CURRENT` or `PETSC_DETERMINE` may be used for `atol` or `rtol` to indicate the current value
    5033:   or the default value from when the object's type was set.

    5035:   With PETSc's implicit schemes for DAE problems, the calculation of the local truncation error
    5036:   (LTE) includes both the differential and the algebraic variables. If one wants the LTE to be
    5037:   computed only for the differential or the algebraic part then this can be done using the vector of
    5038:   tolerances vatol. For example, by setting the tolerance vector with the desired tolerance for the
    5039:   differential part and infinity for the algebraic part, the LTE calculation will include only the
    5040:   differential variables.

    5042:   Fortran Note:
    5043:   Use `PETSC_CURRENT_INTEGER` or `PETSC_DETERMINE_INTEGER`.

    5045: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSGetTolerances()`
    5046: @*/
    5047: PetscErrorCode TSSetTolerances(TS ts, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol)
    5048: {
    5049:   PetscFunctionBegin;
    5050:   if (atol == (PetscReal)PETSC_DETERMINE) {
    5051:     ts->atol = ts->default_atol;
    5052:   } else if (atol != (PetscReal)PETSC_CURRENT) {
    5053:     PetscCheck(atol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)atol);
    5054:     ts->atol = atol;
    5055:   }

    5057:   if (vatol) {
    5058:     PetscCall(PetscObjectReference((PetscObject)vatol));
    5059:     PetscCall(VecDestroy(&ts->vatol));
    5060:     ts->vatol = vatol;
    5061:   }

    5063:   if (rtol == (PetscReal)PETSC_DETERMINE) {
    5064:     ts->rtol = ts->default_rtol;
    5065:   } else if (rtol != (PetscReal)PETSC_CURRENT) {
    5066:     PetscCheck(rtol >= 0.0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative", (double)rtol);
    5067:     ts->rtol = rtol;
    5068:   }

    5070:   if (vrtol) {
    5071:     PetscCall(PetscObjectReference((PetscObject)vrtol));
    5072:     PetscCall(VecDestroy(&ts->vrtol));
    5073:     ts->vrtol = vrtol;
    5074:   }
    5075:   PetscFunctionReturn(PETSC_SUCCESS);
    5076: }

    5078: /*@
    5079:   TSGetTolerances - Get tolerances for local truncation error when using adaptive controller

    5081:   Logically Collective

    5083:   Input Parameter:
    5084: . ts - time integration context

    5086:   Output Parameters:
    5087: + atol  - scalar absolute tolerances, `NULL` to ignore
    5088: . vatol - vector of absolute tolerances, `NULL` to ignore
    5089: . rtol  - scalar relative tolerances, `NULL` to ignore
    5090: - vrtol - vector of relative tolerances, `NULL` to ignore

    5092:   Level: beginner

    5094: .seealso: [](ch_ts), `TS`, `TSAdapt`, `TSErrorWeightedNorm()`, `TSSetTolerances()`
    5095: @*/
    5096: PetscErrorCode TSGetTolerances(TS ts, PetscReal *atol, Vec *vatol, PetscReal *rtol, Vec *vrtol)
    5097: {
    5098:   PetscFunctionBegin;
    5099:   if (atol) *atol = ts->atol;
    5100:   if (vatol) *vatol = ts->vatol;
    5101:   if (rtol) *rtol = ts->rtol;
    5102:   if (vrtol) *vrtol = ts->vrtol;
    5103:   PetscFunctionReturn(PETSC_SUCCESS);
    5104: }

    5106: /*@
    5107:   TSErrorWeightedNorm - compute a weighted norm of the difference between two state vectors based on supplied absolute and relative tolerances

    5109:   Collective

    5111:   Input Parameters:
    5112: + ts        - time stepping context
    5113: . U         - state vector, usually ts->vec_sol
    5114: . Y         - state vector to be compared to U
    5115: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

    5117:   Output Parameters:
    5118: + norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
    5119: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
    5120: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

    5122:   Options Database Key:
    5123: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

    5125:   Level: developer

    5127: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedENorm()`
    5128: @*/
    5129: PetscErrorCode TSErrorWeightedNorm(TS ts, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
    5130: {
    5131:   PetscInt norma_loc, norm_loc, normr_loc;

    5133:   PetscFunctionBegin;
    5135:   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));
    5136:   if (wnormtype == NORM_2) {
    5137:     if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
    5138:     if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
    5139:     if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
    5140:   }
    5141:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
    5142:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
    5143:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
    5144:   PetscFunctionReturn(PETSC_SUCCESS);
    5145: }

    5147: /*@
    5148:   TSErrorWeightedENorm - compute a weighted error norm based on supplied absolute and relative tolerances

    5150:   Collective

    5152:   Input Parameters:
    5153: + ts        - time stepping context
    5154: . E         - error vector
    5155: . U         - state vector, usually ts->vec_sol
    5156: . Y         - state vector, previous time step
    5157: - wnormtype - norm type, either `NORM_2` or `NORM_INFINITY`

    5159:   Output Parameters:
    5160: + norm  - weighted norm, a value of 1.0 achieves a balance between absolute and relative tolerances
    5161: . norma - weighted norm, a value of 1.0 means that the error meets the absolute tolerance set by the user
    5162: - normr - weighted norm, a value of 1.0 means that the error meets the relative tolerance set by the user

    5164:   Options Database Key:
    5165: . -ts_adapt_wnormtype <wnormtype> - 2, INFINITY

    5167:   Level: developer

    5169: .seealso: [](ch_ts), `TS`, `VecErrorWeightedNorms()`, `TSErrorWeightedNorm()`
    5170: @*/
    5171: PetscErrorCode TSErrorWeightedENorm(TS ts, Vec E, Vec U, Vec Y, NormType wnormtype, PetscReal *norm, PetscReal *norma, PetscReal *normr)
    5172: {
    5173:   PetscInt norma_loc, norm_loc, normr_loc;

    5175:   PetscFunctionBegin;
    5177:   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));
    5178:   if (wnormtype == NORM_2) {
    5179:     if (norm_loc) *norm = PetscSqrtReal(PetscSqr(*norm) / norm_loc);
    5180:     if (norma_loc) *norma = PetscSqrtReal(PetscSqr(*norma) / norma_loc);
    5181:     if (normr_loc) *normr = PetscSqrtReal(PetscSqr(*normr) / normr_loc);
    5182:   }
    5183:   PetscCheck(!PetscIsInfOrNanScalar(*norm), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norm");
    5184:   PetscCheck(!PetscIsInfOrNanScalar(*norma), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in norma");
    5185:   PetscCheck(!PetscIsInfOrNanScalar(*normr), PetscObjectComm((PetscObject)ts), PETSC_ERR_FP, "Infinite or not-a-number generated in normr");
    5186:   PetscFunctionReturn(PETSC_SUCCESS);
    5187: }

    5189: /*@
    5190:   TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler

    5192:   Logically Collective

    5194:   Input Parameters:
    5195: + ts      - time stepping context
    5196: - cfltime - maximum stable time step if using forward Euler (value can be different on each process)

    5198:   Note:
    5199:   After calling this function, the global CFL time can be obtained by calling TSGetCFLTime()

    5201:   Level: intermediate

    5203: .seealso: [](ch_ts), `TSGetCFLTime()`, `TSADAPTCFL`
    5204: @*/
    5205: PetscErrorCode TSSetCFLTimeLocal(TS ts, PetscReal cfltime)
    5206: {
    5207:   PetscFunctionBegin;
    5209:   ts->cfltime_local = cfltime;
    5210:   ts->cfltime       = -1.;
    5211:   PetscFunctionReturn(PETSC_SUCCESS);
    5212: }

    5214: /*@
    5215:   TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler

    5217:   Collective

    5219:   Input Parameter:
    5220: . ts - time stepping context

    5222:   Output Parameter:
    5223: . cfltime - maximum stable time step for forward Euler

    5225:   Level: advanced

    5227: .seealso: [](ch_ts), `TSSetCFLTimeLocal()`
    5228: @*/
    5229: PetscErrorCode TSGetCFLTime(TS ts, PetscReal *cfltime)
    5230: {
    5231:   PetscFunctionBegin;
    5232:   if (ts->cfltime < 0) PetscCallMPI(MPIU_Allreduce(&ts->cfltime_local, &ts->cfltime, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
    5233:   *cfltime = ts->cfltime;
    5234:   PetscFunctionReturn(PETSC_SUCCESS);
    5235: }

    5237: /*@
    5238:   TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu

    5240:   Input Parameters:
    5241: + ts - the `TS` context.
    5242: . xl - lower bound.
    5243: - xu - upper bound.

    5245:   Level: advanced

    5247:   Note:
    5248:   If this routine is not called then the lower and upper bounds are set to
    5249:   `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.

    5251: .seealso: [](ch_ts), `TS`
    5252: @*/
    5253: PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu)
    5254: {
    5255:   SNES snes;

    5257:   PetscFunctionBegin;
    5258:   PetscCall(TSGetSNES(ts, &snes));
    5259:   PetscCall(SNESVISetVariableBounds(snes, xl, xu));
    5260:   PetscFunctionReturn(PETSC_SUCCESS);
    5261: }

    5263: /*@
    5264:   TSComputeLinearStability - computes the linear stability function at a point

    5266:   Collective

    5268:   Input Parameters:
    5269: + ts - the `TS` context
    5270: . xr - real part of input argument
    5271: - xi - imaginary part of input argument

    5273:   Output Parameters:
    5274: + yr - real part of function value
    5275: - yi - imaginary part of function value

    5277:   Level: developer

    5279: .seealso: [](ch_ts), `TS`, `TSSetRHSFunction()`, `TSComputeIFunction()`
    5280: @*/
    5281: PetscErrorCode TSComputeLinearStability(TS ts, PetscReal xr, PetscReal xi, PetscReal *yr, PetscReal *yi)
    5282: {
    5283:   PetscFunctionBegin;
    5285:   PetscUseTypeMethod(ts, linearstability, xr, xi, yr, yi);
    5286:   PetscFunctionReturn(PETSC_SUCCESS);
    5287: }

    5289: /*@
    5290:   TSRestartStep - Flags the solver to restart the next step

    5292:   Collective

    5294:   Input Parameter:
    5295: . ts - the `TS` context obtained from `TSCreate()`

    5297:   Level: advanced

    5299:   Notes:
    5300:   Multistep methods like `TSBDF` or Runge-Kutta methods with FSAL property require restarting the solver in the event of
    5301:   discontinuities. These discontinuities may be introduced as a consequence of explicitly modifications to the solution
    5302:   vector (which PETSc attempts to detect and handle) or problem coefficients (which PETSc is not able to detect). For
    5303:   the sake of correctness and maximum safety, users are expected to call `TSRestart()` whenever they introduce
    5304:   discontinuities in callback routines (e.g. prestep and poststep routines, or implicit/rhs function routines with
    5305:   discontinuous source terms).

    5307: .seealso: [](ch_ts), `TS`, `TSBDF`, `TSSolve()`, `TSSetPreStep()`, `TSSetPostStep()`
    5308: @*/
    5309: PetscErrorCode TSRestartStep(TS ts)
    5310: {
    5311:   PetscFunctionBegin;
    5313:   ts->steprestart = PETSC_TRUE;
    5314:   PetscFunctionReturn(PETSC_SUCCESS);
    5315: }

    5317: /*@
    5318:   TSRollBack - Rolls back one time step

    5320:   Collective

    5322:   Input Parameter:
    5323: . ts - the `TS` context obtained from `TSCreate()`

    5325:   Level: advanced

    5327: .seealso: [](ch_ts), `TS`, `TSGetStepRollBack()`, `TSCreate()`, `TSSetUp()`, `TSDestroy()`, `TSSolve()`, `TSSetPreStep()`, `TSSetPreStage()`, `TSInterpolate()`
    5328: @*/
    5329: PetscErrorCode TSRollBack(TS ts)
    5330: {
    5331:   PetscFunctionBegin;
    5333:   PetscCheck(!ts->steprollback, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TSRollBack already called");
    5334:   PetscTryTypeMethod(ts, rollback);
    5335:   PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol));
    5336:   ts->time_step  = ts->ptime - ts->ptime_prev;
    5337:   ts->ptime      = ts->ptime_prev;
    5338:   ts->ptime_prev = ts->ptime_prev_rollback;
    5339:   ts->steps--;
    5340:   ts->steprollback = PETSC_TRUE;
    5341:   PetscFunctionReturn(PETSC_SUCCESS);
    5342: }

    5344: /*@
    5345:   TSGetStepRollBack - Get the internal flag indicating if you are rolling back a step

    5347:   Not collective

    5349:   Input Parameter:
    5350: . ts - the `TS` context obtained from `TSCreate()`

    5352:   Output Parameter:
    5353: . flg - the rollback flag

    5355:   Level: advanced

    5357: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSRollBack()`
    5358: @*/
    5359: PetscErrorCode TSGetStepRollBack(TS ts, PetscBool *flg)
    5360: {
    5361:   PetscFunctionBegin;
    5363:   PetscAssertPointer(flg, 2);
    5364:   *flg = ts->steprollback;
    5365:   PetscFunctionReturn(PETSC_SUCCESS);
    5366: }

    5368: /*@
    5369:   TSGetStepResize - Get the internal flag indicating if the current step is after a resize.

    5371:   Not collective

    5373:   Input Parameter:
    5374: . ts - the `TS` context obtained from `TSCreate()`

    5376:   Output Parameter:
    5377: . flg - the resize flag

    5379:   Level: advanced

    5381: .seealso: [](ch_ts), `TS`, `TSCreate()`, `TSSetResize()`
    5382: @*/
    5383: PetscErrorCode TSGetStepResize(TS ts, PetscBool *flg)
    5384: {
    5385:   PetscFunctionBegin;
    5387:   PetscAssertPointer(flg, 2);
    5388:   *flg = ts->stepresize;
    5389:   PetscFunctionReturn(PETSC_SUCCESS);
    5390: }

    5392: /*@
    5393:   TSGetStages - Get the number of stages and stage values

    5395:   Input Parameter:
    5396: . ts - the `TS` context obtained from `TSCreate()`

    5398:   Output Parameters:
    5399: + ns - the number of stages
    5400: - Y  - the current stage vectors

    5402:   Level: advanced

    5404:   Note:
    5405:   Both `ns` and `Y` can be `NULL`.

    5407: .seealso: [](ch_ts), `TS`, `TSCreate()`
    5408: @*/
    5409: PetscErrorCode TSGetStages(TS ts, PetscInt *ns, Vec **Y)
    5410: {
    5411:   PetscFunctionBegin;
    5413:   if (ns) PetscAssertPointer(ns, 2);
    5414:   if (Y) PetscAssertPointer(Y, 3);
    5415:   if (!ts->ops->getstages) {
    5416:     if (ns) *ns = 0;
    5417:     if (Y) *Y = NULL;
    5418:   } else PetscUseTypeMethod(ts, getstages, ns, Y);
    5419:   PetscFunctionReturn(PETSC_SUCCESS);
    5420: }

    5422: /*@C
    5423:   TSComputeIJacobianDefaultColor - Computes the Jacobian using finite differences and coloring to exploit matrix sparsity.

    5425:   Collective

    5427:   Input Parameters:
    5428: + ts    - the `TS` context
    5429: . t     - current timestep
    5430: . U     - state vector
    5431: . Udot  - time derivative of state vector
    5432: . shift - shift to apply, see note below
    5433: - ctx   - an optional user context

    5435:   Output Parameters:
    5436: + J - Jacobian matrix (not altered in this routine)
    5437: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)

    5439:   Level: intermediate

    5441:   Notes:
    5442:   If F(t,U,Udot)=0 is the DAE, the required Jacobian is

    5444:   dF/dU + shift*dF/dUdot

    5446:   Most users should not need to explicitly call this routine, as it
    5447:   is used internally within the nonlinear solvers.

    5449:   This will first try to get the coloring from the `DM`.  If the `DM` type has no coloring
    5450:   routine, then it will try to get the coloring from the matrix.  This requires that the
    5451:   matrix have nonzero entries precomputed.

    5453: .seealso: [](ch_ts), `TS`, `TSSetIJacobian()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
    5454: @*/
    5455: PetscErrorCode TSComputeIJacobianDefaultColor(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal shift, Mat J, Mat B, void *ctx)
    5456: {
    5457:   SNES          snes;
    5458:   MatFDColoring color;
    5459:   PetscBool     hascolor, matcolor = PETSC_FALSE;

    5461:   PetscFunctionBegin;
    5462:   PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_fd_color_use_mat", &matcolor, NULL));
    5463:   PetscCall(PetscObjectQuery((PetscObject)B, "TSMatFDColoring", (PetscObject *)&color));
    5464:   if (!color) {
    5465:     DM         dm;
    5466:     ISColoring iscoloring;

    5468:     PetscCall(TSGetDM(ts, &dm));
    5469:     PetscCall(DMHasColoring(dm, &hascolor));
    5470:     if (hascolor && !matcolor) {
    5471:       PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
    5472:       PetscCall(MatFDColoringCreate(B, iscoloring, &color));
    5473:       PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode (*)(void))SNESTSFormFunction, (void *)ts));
    5474:       PetscCall(MatFDColoringSetFromOptions(color));
    5475:       PetscCall(MatFDColoringSetUp(B, iscoloring, color));
    5476:       PetscCall(ISColoringDestroy(&iscoloring));
    5477:     } else {
    5478:       MatColoring mc;

    5480:       PetscCall(MatColoringCreate(B, &mc));
    5481:       PetscCall(MatColoringSetDistance(mc, 2));
    5482:       PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
    5483:       PetscCall(MatColoringSetFromOptions(mc));
    5484:       PetscCall(MatColoringApply(mc, &iscoloring));
    5485:       PetscCall(MatColoringDestroy(&mc));
    5486:       PetscCall(MatFDColoringCreate(B, iscoloring, &color));
    5487:       PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode (*)(void))SNESTSFormFunction, (void *)ts));
    5488:       PetscCall(MatFDColoringSetFromOptions(color));
    5489:       PetscCall(MatFDColoringSetUp(B, iscoloring, color));
    5490:       PetscCall(ISColoringDestroy(&iscoloring));
    5491:     }
    5492:     PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)color));
    5493:     PetscCall(PetscObjectDereference((PetscObject)color));
    5494:   }
    5495:   PetscCall(TSGetSNES(ts, &snes));
    5496:   PetscCall(MatFDColoringApply(B, color, U, snes));
    5497:   if (J != B) {
    5498:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    5499:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    5500:   }
    5501:   PetscFunctionReturn(PETSC_SUCCESS);
    5502: }

    5504: /*@C
    5505:   TSSetFunctionDomainError - Set a function that tests if the current state vector is valid

    5507:   Input Parameters:
    5508: + ts   - the `TS` context
    5509: - func - function called within `TSFunctionDomainError()`

    5511:   Calling sequence of `func`:
    5512: + ts     - the `TS` context
    5513: . time   - the current time (of the stage)
    5514: . state  - the state to check if it is valid
    5515: - accept - (output parameter) `PETSC_FALSE` if the state is not acceptable, `PETSC_TRUE` if acceptable

    5517:   Level: intermediate

    5519:   Notes:
    5520:   If an implicit ODE solver is being used then, in addition to providing this routine, the
    5521:   user's code should call `SNESSetFunctionDomainError()` when domain errors occur during
    5522:   function evaluations where the functions are provided by `TSSetIFunction()` or `TSSetRHSFunction()`.
    5523:   Use `TSGetSNES()` to obtain the `SNES` object

    5525:   Developer Notes:
    5526:   The naming of this function is inconsistent with the `SNESSetFunctionDomainError()`
    5527:   since one takes a function pointer and the other does not.

    5529: .seealso: [](ch_ts), `TSAdaptCheckStage()`, `TSFunctionDomainError()`, `SNESSetFunctionDomainError()`, `TSGetSNES()`
    5530: @*/
    5531: PetscErrorCode TSSetFunctionDomainError(TS ts, PetscErrorCode (*func)(TS ts, PetscReal time, Vec state, PetscBool *accept))
    5532: {
    5533:   PetscFunctionBegin;
    5535:   ts->functiondomainerror = func;
    5536:   PetscFunctionReturn(PETSC_SUCCESS);
    5537: }

    5539: /*@
    5540:   TSFunctionDomainError - Checks if the current state is valid

    5542:   Input Parameters:
    5543: + ts        - the `TS` context
    5544: . stagetime - time of the simulation
    5545: - Y         - state vector to check.

    5547:   Output Parameter:
    5548: . accept - Set to `PETSC_FALSE` if the current state vector is valid.

    5550:   Level: developer

    5552:   Note:
    5553:   This function is called by the `TS` integration routines and calls the user provided function (set with `TSSetFunctionDomainError()`)
    5554:   to check if the current state is valid.

    5556: .seealso: [](ch_ts), `TS`, `TSSetFunctionDomainError()`
    5557: @*/
    5558: PetscErrorCode TSFunctionDomainError(TS ts, PetscReal stagetime, Vec Y, PetscBool *accept)
    5559: {
    5560:   PetscFunctionBegin;
    5562:   *accept = PETSC_TRUE;
    5563:   if (ts->functiondomainerror) PetscCall((*ts->functiondomainerror)(ts, stagetime, Y, accept));
    5564:   PetscFunctionReturn(PETSC_SUCCESS);
    5565: }

    5567: /*@
    5568:   TSClone - This function clones a time step `TS` object.

    5570:   Collective

    5572:   Input Parameter:
    5573: . tsin - The input `TS`

    5575:   Output Parameter:
    5576: . tsout - The output `TS` (cloned)

    5578:   Level: developer

    5580:   Notes:
    5581:   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.
    5582:   It will likely be replaced in the future with a mechanism of switching methods on the fly.

    5584:   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
    5585: .vb
    5586:  SNES snes_dup = NULL;
    5587:  TSGetSNES(ts,&snes_dup);
    5588:  TSSetSNES(ts,snes_dup);
    5589: .ve

    5591: .seealso: [](ch_ts), `TS`, `SNES`, `TSCreate()`, `TSSetType()`, `TSSetUp()`, `TSDestroy()`, `TSSetProblemType()`
    5592: @*/
    5593: PetscErrorCode TSClone(TS tsin, TS *tsout)
    5594: {
    5595:   TS     t;
    5596:   SNES   snes_start;
    5597:   DM     dm;
    5598:   TSType type;

    5600:   PetscFunctionBegin;
    5601:   PetscAssertPointer(tsin, 1);
    5602:   *tsout = NULL;

    5604:   PetscCall(PetscHeaderCreate(t, TS_CLASSID, "TS", "Time stepping", "TS", PetscObjectComm((PetscObject)tsin), TSDestroy, TSView));

    5606:   /* General TS description */
    5607:   t->numbermonitors    = 0;
    5608:   t->monitorFrequency  = 1;
    5609:   t->setupcalled       = 0;
    5610:   t->ksp_its           = 0;
    5611:   t->snes_its          = 0;
    5612:   t->nwork             = 0;
    5613:   t->rhsjacobian.time  = PETSC_MIN_REAL;
    5614:   t->rhsjacobian.scale = 1.;
    5615:   t->ijacobian.shift   = 1.;

    5617:   PetscCall(TSGetSNES(tsin, &snes_start));
    5618:   PetscCall(TSSetSNES(t, snes_start));

    5620:   PetscCall(TSGetDM(tsin, &dm));
    5621:   PetscCall(TSSetDM(t, dm));

    5623:   t->adapt = tsin->adapt;
    5624:   PetscCall(PetscObjectReference((PetscObject)t->adapt));

    5626:   t->trajectory = tsin->trajectory;
    5627:   PetscCall(PetscObjectReference((PetscObject)t->trajectory));

    5629:   t->event = tsin->event;
    5630:   if (t->event) t->event->refct++;

    5632:   t->problem_type      = tsin->problem_type;
    5633:   t->ptime             = tsin->ptime;
    5634:   t->ptime_prev        = tsin->ptime_prev;
    5635:   t->time_step         = tsin->time_step;
    5636:   t->max_time          = tsin->max_time;
    5637:   t->steps             = tsin->steps;
    5638:   t->max_steps         = tsin->max_steps;
    5639:   t->equation_type     = tsin->equation_type;
    5640:   t->atol              = tsin->atol;
    5641:   t->rtol              = tsin->rtol;
    5642:   t->max_snes_failures = tsin->max_snes_failures;
    5643:   t->max_reject        = tsin->max_reject;
    5644:   t->errorifstepfailed = tsin->errorifstepfailed;

    5646:   PetscCall(TSGetType(tsin, &type));
    5647:   PetscCall(TSSetType(t, type));

    5649:   t->vec_sol = NULL;

    5651:   t->cfltime          = tsin->cfltime;
    5652:   t->cfltime_local    = tsin->cfltime_local;
    5653:   t->exact_final_time = tsin->exact_final_time;

    5655:   t->ops[0] = tsin->ops[0];

    5657:   if (((PetscObject)tsin)->fortran_func_pointers) {
    5658:     PetscInt i;
    5659:     PetscCall(PetscMalloc((10) * sizeof(void (*)(void)), &((PetscObject)t)->fortran_func_pointers));
    5660:     for (i = 0; i < 10; i++) ((PetscObject)t)->fortran_func_pointers[i] = ((PetscObject)tsin)->fortran_func_pointers[i];
    5661:   }
    5662:   *tsout = t;
    5663:   PetscFunctionReturn(PETSC_SUCCESS);
    5664: }

    5666: static PetscErrorCode RHSWrapperFunction_TSRHSJacobianTest(void *ctx, Vec x, Vec y)
    5667: {
    5668:   TS ts = (TS)ctx;

    5670:   PetscFunctionBegin;
    5671:   PetscCall(TSComputeRHSFunction(ts, 0, x, y));
    5672:   PetscFunctionReturn(PETSC_SUCCESS);
    5673: }

    5675: /*@
    5676:   TSRHSJacobianTest - Compares the multiply routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.

    5678:   Logically Collective

    5680:   Input Parameter:
    5681: . ts - the time stepping routine

    5683:   Output Parameter:
    5684: . flg - `PETSC_TRUE` if the multiply is likely correct

    5686:   Options Database Key:
    5687: . -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view - run the test at each timestep of the integrator

    5689:   Level: advanced

    5691:   Note:
    5692:   This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian

    5694: .seealso: [](ch_ts), `TS`, `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTestTranspose()`
    5695: @*/
    5696: PetscErrorCode TSRHSJacobianTest(TS ts, PetscBool *flg)
    5697: {
    5698:   Mat              J, B;
    5699:   TSRHSJacobianFn *func;
    5700:   void            *ctx;

    5702:   PetscFunctionBegin;
    5703:   PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
    5704:   PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
    5705:   PetscCall(MatShellTestMult(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
    5706:   PetscFunctionReturn(PETSC_SUCCESS);
    5707: }

    5709: /*@
    5710:   TSRHSJacobianTestTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on the `TS` given RHS function.

    5712:   Logically Collective

    5714:   Input Parameter:
    5715: . ts - the time stepping routine

    5717:   Output Parameter:
    5718: . flg - `PETSC_TRUE` if the multiply is likely correct

    5720:   Options Database Key:
    5721: . -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view - run the test at each timestep of the integrator

    5723:   Level: advanced

    5725:   Notes:
    5726:   This only works for problems defined using `TSSetRHSFunction()` and Jacobian NOT `TSSetIFunction()` and Jacobian

    5728: .seealso: [](ch_ts), `TS`, `Mat`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`, `TSRHSJacobianTest()`
    5729: @*/
    5730: PetscErrorCode TSRHSJacobianTestTranspose(TS ts, PetscBool *flg)
    5731: {
    5732:   Mat              J, B;
    5733:   void            *ctx;
    5734:   TSRHSJacobianFn *func;

    5736:   PetscFunctionBegin;
    5737:   PetscCall(TSGetRHSJacobian(ts, &J, &B, &func, &ctx));
    5738:   PetscCall((*func)(ts, 0.0, ts->vec_sol, J, B, ctx));
    5739:   PetscCall(MatShellTestMultTranspose(J, RHSWrapperFunction_TSRHSJacobianTest, ts->vec_sol, ts, flg));
    5740:   PetscFunctionReturn(PETSC_SUCCESS);
    5741: }

    5743: /*@
    5744:   TSSetUseSplitRHSFunction - Use the split RHSFunction when a multirate method is used.

    5746:   Logically Collective

    5748:   Input Parameters:
    5749: + ts                   - timestepping context
    5750: - use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

    5752:   Options Database Key:
    5753: . -ts_use_splitrhsfunction - <true,false>

    5755:   Level: intermediate

    5757:   Note:
    5758:   This is only for multirate methods

    5760: .seealso: [](ch_ts), `TS`, `TSGetUseSplitRHSFunction()`
    5761: @*/
    5762: PetscErrorCode TSSetUseSplitRHSFunction(TS ts, PetscBool use_splitrhsfunction)
    5763: {
    5764:   PetscFunctionBegin;
    5766:   ts->use_splitrhsfunction = use_splitrhsfunction;
    5767:   PetscFunctionReturn(PETSC_SUCCESS);
    5768: }

    5770: /*@
    5771:   TSGetUseSplitRHSFunction - Gets whether to use the split RHSFunction when a multirate method is used.

    5773:   Not Collective

    5775:   Input Parameter:
    5776: . ts - timestepping context

    5778:   Output Parameter:
    5779: . use_splitrhsfunction - `PETSC_TRUE` indicates that the split RHSFunction will be used

    5781:   Level: intermediate

    5783: .seealso: [](ch_ts), `TS`, `TSSetUseSplitRHSFunction()`
    5784: @*/
    5785: PetscErrorCode TSGetUseSplitRHSFunction(TS ts, PetscBool *use_splitrhsfunction)
    5786: {
    5787:   PetscFunctionBegin;
    5789:   *use_splitrhsfunction = ts->use_splitrhsfunction;
    5790:   PetscFunctionReturn(PETSC_SUCCESS);
    5791: }

    5793: /*@
    5794:   TSSetMatStructure - sets the relationship between the nonzero structure of the RHS Jacobian matrix to the IJacobian matrix.

    5796:   Logically  Collective

    5798:   Input Parameters:
    5799: + ts  - the time-stepper
    5800: - str - the structure (the default is `UNKNOWN_NONZERO_PATTERN`)

    5802:   Level: intermediate

    5804:   Note:
    5805:   When the relationship between the nonzero structures is known and supplied the solution process can be much faster

    5807: .seealso: [](ch_ts), `TS`, `MatAXPY()`, `MatStructure`
    5808:  @*/
    5809: PetscErrorCode TSSetMatStructure(TS ts, MatStructure str)
    5810: {
    5811:   PetscFunctionBegin;
    5813:   ts->axpy_pattern = str;
    5814:   PetscFunctionReturn(PETSC_SUCCESS);
    5815: }

    5817: /*@
    5818:   TSSetTimeSpan - sets the time span. The solution will be computed and stored for each time requested in the span

    5820:   Collective

    5822:   Input Parameters:
    5823: + ts         - the time-stepper
    5824: . n          - number of the time points (>=2)
    5825: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.

    5827:   Options Database Key:
    5828: . -ts_time_span <t0,...tf> - Sets the time span

    5830:   Level: intermediate

    5832:   Notes:
    5833:   The elements in tspan must be all increasing. They correspond to the intermediate points for time integration.
    5834:   `TS_EXACTFINALTIME_MATCHSTEP` must be used to make the last time step in each sub-interval match the intermediate points specified.
    5835:   The intermediate solutions are saved in a vector array that can be accessed with `TSGetTimeSpanSolutions()`. Thus using time span may
    5836:   pressure the memory system when using a large number of span points.

    5838: .seealso: [](ch_ts), `TS`, `TSGetTimeSpan()`, `TSGetTimeSpanSolutions()`
    5839:  @*/
    5840: PetscErrorCode TSSetTimeSpan(TS ts, PetscInt n, PetscReal *span_times)
    5841: {
    5842:   PetscFunctionBegin;
    5844:   PetscCheck(n >= 2, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONG, "Minimum time span size is 2 but %" PetscInt_FMT " is provided", n);
    5845:   if (ts->tspan && n != ts->tspan->num_span_times) {
    5846:     PetscCall(PetscFree(ts->tspan->span_times));
    5847:     PetscCall(VecDestroyVecs(ts->tspan->num_span_times, &ts->tspan->vecs_sol));
    5848:     PetscCall(PetscMalloc1(n, &ts->tspan->span_times));
    5849:   }
    5850:   if (!ts->tspan) {
    5851:     TSTimeSpan tspan;
    5852:     PetscCall(PetscNew(&tspan));
    5853:     PetscCall(PetscMalloc1(n, &tspan->span_times));
    5854:     tspan->reltol  = 1e-6;
    5855:     tspan->abstol  = 10 * PETSC_MACHINE_EPSILON;
    5856:     tspan->worktol = 0;
    5857:     ts->tspan      = tspan;
    5858:   }
    5859:   ts->tspan->num_span_times = n;
    5860:   PetscCall(PetscArraycpy(ts->tspan->span_times, span_times, n));
    5861:   PetscCall(TSSetTime(ts, ts->tspan->span_times[0]));
    5862:   PetscCall(TSSetMaxTime(ts, ts->tspan->span_times[n - 1]));
    5863:   PetscFunctionReturn(PETSC_SUCCESS);
    5864: }

    5866: /*@C
    5867:   TSGetTimeSpan - gets the time span set with `TSSetTimeSpan()`

    5869:   Not Collective

    5871:   Input Parameter:
    5872: . ts - the time-stepper

    5874:   Output Parameters:
    5875: + n          - number of the time points (>=2)
    5876: - span_times - array of the time points. The first element and the last element are the initial time and the final time respectively.

    5878:   Level: beginner

    5880:   Note:
    5881:   The values obtained are valid until the `TS` object is destroyed.

    5883:   Both `n` and `span_times` can be `NULL`.

    5885: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`, `TSGetTimeSpanSolutions()`
    5886:  @*/
    5887: PetscErrorCode TSGetTimeSpan(TS ts, PetscInt *n, const PetscReal *span_times[])
    5888: {
    5889:   PetscFunctionBegin;
    5891:   if (n) PetscAssertPointer(n, 2);
    5892:   if (span_times) PetscAssertPointer(span_times, 3);
    5893:   if (!ts->tspan) {
    5894:     if (n) *n = 0;
    5895:     if (span_times) *span_times = NULL;
    5896:   } else {
    5897:     if (n) *n = ts->tspan->num_span_times;
    5898:     if (span_times) *span_times = ts->tspan->span_times;
    5899:   }
    5900:   PetscFunctionReturn(PETSC_SUCCESS);
    5901: }

    5903: /*@
    5904:   TSGetTimeSpanSolutions - Get the number of solutions and the solutions at the time points specified by the time span.

    5906:   Input Parameter:
    5907: . ts - the `TS` context obtained from `TSCreate()`

    5909:   Output Parameters:
    5910: + nsol - the number of solutions
    5911: - Sols - the solution vectors

    5913:   Level: intermediate

    5915:   Notes:
    5916:   Both `nsol` and `Sols` can be `NULL`.

    5918:   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()`.
    5919:   For example, manipulating the step size, especially with a reduced precision, may cause `TS` to step over certain points in the span.

    5921: .seealso: [](ch_ts), `TS`, `TSSetTimeSpan()`
    5922: @*/
    5923: PetscErrorCode TSGetTimeSpanSolutions(TS ts, PetscInt *nsol, Vec **Sols)
    5924: {
    5925:   PetscFunctionBegin;
    5927:   if (nsol) PetscAssertPointer(nsol, 2);
    5928:   if (Sols) PetscAssertPointer(Sols, 3);
    5929:   if (!ts->tspan) {
    5930:     if (nsol) *nsol = 0;
    5931:     if (Sols) *Sols = NULL;
    5932:   } else {
    5933:     if (nsol) *nsol = ts->tspan->spanctr;
    5934:     if (Sols) *Sols = ts->tspan->vecs_sol;
    5935:   }
    5936:   PetscFunctionReturn(PETSC_SUCCESS);
    5937: }

    5939: /*@
    5940:   TSPruneIJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.

    5942:   Collective

    5944:   Input Parameters:
    5945: + ts - the `TS` context
    5946: . J  - Jacobian matrix (not altered in this routine)
    5947: - B  - newly computed Jacobian matrix to use with preconditioner

    5949:   Level: intermediate

    5951:   Notes:
    5952:   This function improves the `MatFDColoring` performance when the Jacobian matrix was over-allocated or contains
    5953:   many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
    5954:   and multiple fields are involved.

    5956:   Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
    5957:   structure. For `MatFDColoring`, the values of nonzero entries are not important. So one can
    5958:   usually call `TSComputeIJacobian()` with randomized input vectors to generate a dummy Jacobian.
    5959:   `TSComputeIJacobian()` should be called before `TSSolve()` but after `TSSetUp()`.

    5961: .seealso: [](ch_ts), `TS`, `MatFDColoring`, `TSComputeIJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
    5962: @*/
    5963: PetscErrorCode TSPruneIJacobianColor(TS ts, Mat J, Mat B)
    5964: {
    5965:   MatColoring   mc            = NULL;
    5966:   ISColoring    iscoloring    = NULL;
    5967:   MatFDColoring matfdcoloring = NULL;

    5969:   PetscFunctionBegin;
    5970:   /* Generate new coloring after eliminating zeros in the matrix */
    5971:   PetscCall(MatEliminateZeros(B, PETSC_TRUE));
    5972:   PetscCall(MatColoringCreate(B, &mc));
    5973:   PetscCall(MatColoringSetDistance(mc, 2));
    5974:   PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
    5975:   PetscCall(MatColoringSetFromOptions(mc));
    5976:   PetscCall(MatColoringApply(mc, &iscoloring));
    5977:   PetscCall(MatColoringDestroy(&mc));
    5978:   /* Replace the old coloring with the new one */
    5979:   PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
    5980:   PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode (*)(void))SNESTSFormFunction, (void *)ts));
    5981:   PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
    5982:   PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
    5983:   PetscCall(PetscObjectCompose((PetscObject)B, "TSMatFDColoring", (PetscObject)matfdcoloring));
    5984:   PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
    5985:   PetscCall(ISColoringDestroy(&iscoloring));
    5986:   PetscFunctionReturn(PETSC_SUCCESS);
    5987: }