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;

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

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

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

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

 48:   TSSetStepNumber(ts, 0);
 49:   TSSetTimeStep(ts, 1);
 50:   TSSetTime(ts, 0);
 51:   TSSetMaxTime(ts, PETSC_MAX_REAL);
 52:   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:   TSSolve(ts, NULL);

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

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

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

 79:   j = shift * 2.0;
 80:   MatSetValue(B, 0, 0, j, INSERT_VALUES);
 81:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 82:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 83:   return 0;
 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*/