1 static char help[] = "Small ODE to test TS accuracy.\n";
2
3 /*
4 The ODE
5 u1_t = cos(t),
6 u2_t = sin(u2)
7 with analytical solution
8 u1(t) = sin(t),
9 u2(t) = 2 * atan(exp(t) * tan(0.5))
10 is used to test the accuracy of TS schemes.
11 */
12
13 #include <petscts.h>
14
15 /*
16 Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
17 */
RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void * s)18 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *s)
19 {
20 PetscScalar *f;
21 const PetscScalar *u;
22
23 PetscFunctionBeginUser;
24 PetscCall(VecGetArrayRead(U, &u));
25 PetscCall(VecGetArray(F, &f));
26
27 f[0] = PetscCosReal(t);
28 f[1] = PetscSinScalar(u[1]);
29
30 PetscCall(VecRestoreArrayRead(U, &u));
31 PetscCall(VecRestoreArray(F, &f));
32 PetscFunctionReturn(PETSC_SUCCESS);
33 }
34
35 /*
36 Defines the exact solution.
37 */
ExactSolution(PetscReal t,Vec U)38 static PetscErrorCode ExactSolution(PetscReal t, Vec U)
39 {
40 PetscScalar *u;
41
42 PetscFunctionBeginUser;
43 PetscCall(VecGetArray(U, &u));
44
45 u[0] = PetscSinReal(t);
46 u[1] = 2 * PetscAtanReal(PetscExpReal(t) * PetscTanReal(0.5));
47
48 PetscCall(VecRestoreArray(U, &u));
49 PetscFunctionReturn(PETSC_SUCCESS);
50 }
51
main(int argc,char ** argv)52 int main(int argc, char **argv)
53 {
54 TS ts; /* ODE integrator */
55 Vec U; /* numerical solution will be stored here */
56 Vec Uex; /* analytical (exact) solution will be stored here */
57 PetscMPIInt size;
58 PetscInt n = 2;
59 PetscScalar *u;
60 PetscReal t, final_time = 1.0, dt = 0.25;
61 PetscReal error;
62 TSAdapt adapt;
63
64 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65 Initialize program
66 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67 PetscFunctionBeginUser;
68 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
69 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
70 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
71
72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73 Create timestepping solver context
74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
76 PetscCall(TSSetType(ts, TSROSW));
77
78 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79 Set ODE routines
80 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
82 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
83
84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85 Set initial conditions
86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87 PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
88 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
89 PetscCall(VecSetUp(U));
90 PetscCall(VecGetArray(U, &u));
91 u[0] = 0.0;
92 u[1] = 1.0;
93 PetscCall(VecRestoreArray(U, &u));
94 PetscCall(TSSetSolution(ts, U));
95
96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 Set solver options
98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99 PetscCall(TSSetSaveTrajectory(ts));
100 PetscCall(TSSetMaxTime(ts, final_time));
101 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
102 PetscCall(TSSetTimeStep(ts, dt));
103 /* The adaptive time step controller is forced to take constant time steps. */
104 PetscCall(TSGetAdapt(ts, &adapt));
105 PetscCall(TSAdaptSetType(adapt, TSADAPTNONE));
106
107 PetscCall(TSSetFromOptions(ts));
108
109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110 Run timestepping solver and compute error
111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112 PetscCall(TSSolve(ts, U));
113 PetscCall(TSGetTime(ts, &t));
114
115 if (PetscAbsReal(t - final_time) > 100 * PETSC_MACHINE_EPSILON) {
116 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Note: There is a difference of %g between the prescribed final time %g and the actual final time.\n", (double)(final_time - t), (double)final_time));
117 }
118 PetscCall(VecDuplicate(U, &Uex));
119 PetscCall(ExactSolution(t, Uex));
120
121 PetscCall(VecAYPX(Uex, -1.0, U));
122 PetscCall(VecNorm(Uex, NORM_2, &error));
123 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error at final time: %.2E\n", (double)error));
124
125 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126 Free work space. All PETSc objects should be destroyed when they are no longer needed.
127 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128 PetscCall(VecDestroy(&U));
129 PetscCall(VecDestroy(&Uex));
130 PetscCall(TSDestroy(&ts));
131
132 PetscCall(PetscFinalize());
133 return 0;
134 }
135
136 /*TEST
137
138 test:
139 suffix: 3bs
140 args: -ts_type rk -ts_rk_type 3bs
141 requires: !single
142
143 test:
144 suffix: 5bs
145 args: -ts_type rk -ts_rk_type 5bs
146 requires: !single
147
148 test:
149 suffix: 5dp
150 args: -ts_type rk -ts_rk_type 5dp
151 requires: !single
152
153 test:
154 suffix: 6vr
155 args: -ts_type rk -ts_rk_type 6vr
156 requires: !single
157
158 test:
159 suffix: 7vr
160 args: -ts_type rk -ts_rk_type 7vr
161 requires: !single
162
163 test:
164 suffix: 8vr
165 args: -ts_type rk -ts_rk_type 8vr
166 requires: !single
167
168 TEST*/
169