Actual source code: ex26.c

  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 *);

  9: int main(int argc, char **argv)
 10: {
 11:   TS  ts;
 12:   Vec x;
 13:   Vec f;
 14:   Mat A;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 19:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
 20:   PetscCall(TSSetEquationType(ts, TS_EQ_ODE_IMPLICIT));
 21:   PetscCall(VecCreate(PETSC_COMM_WORLD, &f));
 22:   PetscCall(VecSetSizes(f, 1, PETSC_DECIDE));
 23:   PetscCall(VecSetFromOptions(f));
 24:   PetscCall(VecSetUp(f));
 25:   PetscCall(TSSetIFunction(ts, f, IFunction, NULL));
 26:   PetscCall(VecDestroy(&f));

 28:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 29:   PetscCall(MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE));
 30:   PetscCall(MatSetFromOptions(A));
 31:   PetscCall(MatSetUp(A));
 32:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 33:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 34:   /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
 35:   PetscCall(MatShift(A, (PetscReal)1));
 36:   PetscCall(MatShift(A, (PetscReal)-1));
 37:   PetscCall(TSSetIJacobian(ts, A, A, IJacobian, NULL));
 38:   PetscCall(MatDestroy(&A));

 40:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 41:   PetscCall(VecSetSizes(x, 1, PETSC_DECIDE));
 42:   PetscCall(VecSetFromOptions(x));
 43:   PetscCall(VecSetUp(x));
 44:   PetscCall(TSSetSolution(ts, x));
 45:   PetscCall(VecDestroy(&x));
 46:   PetscCall(TSSetFromOptions(ts));

 48:   PetscCall(TSSetStepNumber(ts, 0));
 49:   PetscCall(TSSetTimeStep(ts, 1));
 50:   PetscCall(TSSetTime(ts, 0));
 51:   PetscCall(TSSetMaxTime(ts, PETSC_MAX_REAL));
 52:   PetscCall(TSSetMaxSteps(ts, 3));

 54:   /*
 55:       When an ARKIMEX scheme with an explicit stage is used this will error with a message informing the user it is not possible to use
 56:       a non-trivial mass matrix with ARKIMEX schemes with explicit stages.
 57:   */
 58:   PetscCall(TSSolve(ts, NULL));

 60:   PetscCall(TSDestroy(&ts));
 61:   PetscCall(PetscFinalize());
 62:   return 0;
 63: }

 65: PetscErrorCode IFunction(TS ts, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx)
 66: {
 67:   PetscFunctionBeginUser;
 68:   PetscCall(VecCopy(xdot, f));
 69:   PetscCall(VecScale(f, 2.0));
 70:   PetscCall(VecShift(f, -1.0));
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

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

 78:   PetscFunctionBeginUser;
 79:   j = shift * 2.0;
 80:   PetscCall(MatSetValue(B, 0, 0, j, INSERT_VALUES));
 81:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 82:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
 83:   PetscFunctionReturn(PETSC_SUCCESS);
 84: }

 86: /*TEST

 88:     test:
 89:       suffix: arkimex_explicit_stage
 90:       requires: !defined(PETSCTEST_VALGRIND) defined(PETSC_USE_DEBUG)
 91:       args: -ts_type arkimex -petsc_ci_portable_error_output -error_output_stdout
 92:       filter: grep -E -v "(options_left|memory block|leaked context|not freed before MPI_Finalize|Could be the program crashed)"

 94:     test:
 95:       suffix: arkimex_implicit_stage
 96:       args: -ts_type arkimex -ts_arkimex_type l2 -ts_monitor_solution -ts_monitor

 98: TEST*/