1: static char help[] = "Solves the trivial ODE 2 du/dt = 1, u(0) = 0. \n\n";

  3: #include <petscts.h>
  4: #include <petscpc.h>

  6: PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *);
  7: PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
  8: PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);

 10: int main(int argc, char **argv)
 11: {
 12:   TS        ts;
 13:   Vec       x;
 14:   Mat       A;
 15:   PetscBool flg = PETSC_FALSE, usingimex;

 17:   PetscFunctionBeginUser;
 18:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 20:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
 21:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-set_implicit", &flg, NULL));
 22:   if (flg) PetscCall(TSSetEquationType(ts, TS_EQ_ODE_IMPLICIT));
 23:   PetscCall(TSSetIFunction(ts, NULL, IFunction, NULL));
 24:   PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &usingimex));

 26:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 27:   PetscCall(MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE));
 28:   PetscCall(MatSetFromOptions(A));
 29:   PetscCall(MatSetUp(A));
 30:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 31:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 32:   PetscCall(TSSetIJacobian(ts, A, A, IJacobian, NULL));
 33:   PetscCall(MatDestroy(&A));

 35:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 36:   PetscCall(VecSetSizes(x, 1, PETSC_DECIDE));
 37:   PetscCall(VecSetFromOptions(x));
 38:   PetscCall(VecSetUp(x));
 39:   PetscCall(TSSetSolution(ts, x));
 40:   PetscCall(VecDestroy(&x));
 41:   PetscCall(TSSetFromOptions(ts));

 43:   /* Need to know if we are using an IMEX scheme to decide on the form
 44:      of the RHS function */
 45:   PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSARKIMEX, &usingimex));
 46:   if (usingimex) {
 47:     PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg));
 48:     if (flg) usingimex = PETSC_FALSE;
 49:   }
 50:   PetscCall(TSSetStepNumber(ts, 0));
 51:   PetscCall(TSSetTimeStep(ts, 1));
 52:   PetscCall(TSSetTime(ts, 0));
 53:   PetscCall(TSSetMaxTime(ts, PETSC_MAX_REAL));
 54:   PetscCall(TSSetMaxSteps(ts, 3));

 56:   PetscCall(TSSolve(ts, NULL));

 58:   PetscCall(TSDestroy(&ts));
 59:   PetscCall(PetscFinalize());
 60:   return 0;
 61: }

 63: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
 64: {
 65:   PetscBool usingimex = *(PetscBool *)ctx;

 67:   PetscFunctionBeginUser;
 68:   PetscCall(VecSet(f, usingimex ? 0.5 : 1));
 69:   PetscFunctionReturn(PETSC_SUCCESS);
 70: }

 72: PetscErrorCode IFunction(TS ts, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx)
 73: {
 74:   PetscFunctionBeginUser;
 75:   PetscCall(VecCopy(xdot, f));
 76:   PetscCall(VecScale(f, 2.0));
 77:   PetscFunctionReturn(PETSC_SUCCESS);
 78: }

 80: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec x, Vec xdot, PetscReal shift, Mat A, Mat B, void *ctx)
 81: {
 82:   PetscScalar j;

 84:   PetscFunctionBeginUser;
 85:   j = shift * 2.0;
 86:   PetscCall(MatSetValue(B, 0, 0, j, INSERT_VALUES));
 87:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 88:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: /*TEST

 94:     test:
 95:       suffix: arkimex_explicit_stage
 96:       requires: !defined(PETSCTEST_VALGRIND) defined(PETSC_USE_DEBUG) !defined(PETSC_HAVE_SANITIZER)
 97:       args: -ts_type arkimex -petsc_ci_portable_error_output -error_output_stdout -set_implicit
 98:       filter: grep -E -v "(memory block|leaked context|not freed before MPI_Finalize|Could be the program crashed)"

100:     test:
101:       suffix: arkimex_implicit_stage
102:       args: -ts_type arkimex -ts_arkimex_type {{3 l2}} -ts_monitor_solution -ts_monitor

104: TEST*/