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*/