Actual source code: ex51.c
1: static char help[] = "Small ODE to test TS accuracy.\n";
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: */
13: #include <petscts.h>
15: /*
16: Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
17: */
18: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *s)
19: {
20: PetscScalar *f;
21: const PetscScalar *u;
24: VecGetArrayRead(U, &u);
25: VecGetArray(F, &f);
27: f[0] = PetscCosReal(t);
28: f[1] = PetscSinReal(u[1]);
30: VecRestoreArrayRead(U, &u);
31: VecRestoreArray(F, &f);
32: return 0;
33: }
35: /*
36: Defines the exact solution.
37: */
38: static PetscErrorCode ExactSolution(PetscReal t, Vec U)
39: {
40: PetscScalar *u;
43: VecGetArray(U, &u);
45: u[0] = PetscSinReal(t);
46: u[1] = 2 * PetscAtanReal(PetscExpReal(t) * PetscTanReal(0.5));
48: VecRestoreArray(U, &u);
49: return 0;
50: }
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;
64: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: Initialize program
66: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68: PetscInitialize(&argc, &argv, (char *)0, help);
69: MPI_Comm_size(PETSC_COMM_WORLD, &size);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Create timestepping solver context
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: TSCreate(PETSC_COMM_WORLD, &ts);
76: TSSetType(ts, TSROSW);
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Set ODE routines
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: TSSetProblemType(ts, TS_NONLINEAR);
82: TSSetRHSFunction(ts, NULL, RHSFunction, NULL);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Set initial conditions
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: VecCreate(PETSC_COMM_WORLD, &U);
88: VecSetSizes(U, n, PETSC_DETERMINE);
89: VecSetUp(U);
90: VecGetArray(U, &u);
91: u[0] = 0.0;
92: u[1] = 1.0;
93: VecRestoreArray(U, &u);
94: TSSetSolution(ts, U);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Set solver options
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: TSSetSaveTrajectory(ts);
100: TSSetMaxTime(ts, final_time);
101: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
102: TSSetTimeStep(ts, dt);
103: /* The adaptive time step controller is forced to take constant time steps. */
104: TSGetAdapt(ts, &adapt);
105: TSAdaptSetType(adapt, TSADAPTNONE);
107: TSSetFromOptions(ts);
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Run timestepping solver and compute error
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: TSSolve(ts, U);
113: TSGetTime(ts, &t);
115: if (PetscAbsReal(t - final_time) > 100 * PETSC_MACHINE_EPSILON) {
116: 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: VecDuplicate(U, &Uex);
119: ExactSolution(t, Uex);
121: VecAYPX(Uex, -1.0, U);
122: VecNorm(Uex, NORM_2, &error);
123: PetscPrintf(PETSC_COMM_WORLD, "Error at final time: %.2E\n", (double)error);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Free work space. All PETSc objects should be destroyed when they are no longer needed.
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: VecDestroy(&U);
129: VecDestroy(&Uex);
130: TSDestroy(&ts);
132: PetscFinalize();
133: return 0;
134: }
136: /*TEST
138: test:
139: suffix: 3bs
140: args: -ts_type rk -ts_rk_type 3bs
141: requires: !single
143: test:
144: suffix: 5bs
145: args: -ts_type rk -ts_rk_type 5bs
146: requires: !single
148: test:
149: suffix: 5dp
150: args: -ts_type rk -ts_rk_type 5dp
151: requires: !single
153: test:
154: suffix: 6vr
155: args: -ts_type rk -ts_rk_type 6vr
156: requires: !single
158: test:
159: suffix: 7vr
160: args: -ts_type rk -ts_rk_type 7vr
161: requires: !single
163: test:
164: suffix: 8vr
165: args: -ts_type rk -ts_rk_type 8vr
166: requires: !single
168: TEST*/