1 static char help[] = "Solves the trivial ODE 2 du/dt = 1, u(0) = 0. \n\n";
2
3 #include <petscts.h>
4 #include <petscpc.h>
5
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 *);
9
main(int argc,char ** argv)10 int main(int argc, char **argv)
11 {
12 TS ts;
13 Vec x;
14 Mat A;
15 PetscBool flg = PETSC_FALSE, usingimex;
16
17 PetscFunctionBeginUser;
18 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
19
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));
25
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));
34
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));
42
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));
55
56 PetscCall(TSSolve(ts, NULL));
57
58 PetscCall(TSDestroy(&ts));
59 PetscCall(PetscFinalize());
60 return 0;
61 }
62
RHSFunction(TS ts,PetscReal t,Vec x,Vec f,PetscCtx ctx)63 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, PetscCtx ctx)
64 {
65 PetscBool usingimex = *(PetscBool *)ctx;
66
67 PetscFunctionBeginUser;
68 PetscCall(VecSet(f, usingimex ? 0.5 : 1));
69 PetscFunctionReturn(PETSC_SUCCESS);
70 }
71
IFunction(TS ts,PetscReal t,Vec x,Vec xdot,Vec f,PetscCtx ctx)72 PetscErrorCode IFunction(TS ts, PetscReal t, Vec x, Vec xdot, Vec f, PetscCtx ctx)
73 {
74 PetscFunctionBeginUser;
75 PetscCall(VecCopy(xdot, f));
76 PetscCall(VecScale(f, 2.0));
77 PetscFunctionReturn(PETSC_SUCCESS);
78 }
79
IJacobian(TS ts,PetscReal t,Vec x,Vec xdot,PetscReal shift,Mat A,Mat B,PetscCtx ctx)80 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec x, Vec xdot, PetscReal shift, Mat A, Mat B, PetscCtx ctx)
81 {
82 PetscScalar j;
83
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 }
91
92 /*TEST
93
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)"
99
100 test:
101 suffix: arkimex_implicit_stage
102 args: -ts_type arkimex -ts_arkimex_type {{3 l2}} -ts_monitor_solution -ts_monitor
103
104 TEST*/
105