1 static char help[] = "Solves the van der Pol DAE.\n\
2 Input parameters include:\n";
3
4 /* ------------------------------------------------------------------------
5
6 This program solves the van der Pol DAE
7 y' = -z = f(y,z) (1)
8 0 = y-(z^3/3 - z) = g(y,z)
9 on the domain 0 <= x <= 1, with the boundary conditions
10 y(0) = -2, y'(0) = -2.355301397608119909925287735864250951918
11 This is a nonlinear equation.
12
13 Notes:
14 This code demonstrates the TS solver interface with the Van der Pol DAE,
15 namely it is the case when there is no RHS (meaning the RHS == 0), and the
16 equations are converted to two variants of linear problems, u_t = f(u,t),
17 namely turning (1) into a vector equation in terms of u,
18
19 [ y' + z ] = [ 0 ]
20 [ (z^3/3 - z) - y ] [ 0 ]
21
22 which then we can write as a vector equation
23
24 [ u_1' + u_2 ] = [ 0 ] (2)
25 [ (u_2^3/3 - u_2) - u_1 ] [ 0 ]
26
27 which is now in the desired form of u_t = f(u,t). As this is a DAE, and
28 there is no u_2', there is no need for a split,
29
30 so
31
32 [ F(u',u,t) ] = [ u_1' ] + [ u_2 ]
33 [ 0 ] [ (u_2^3/3 - u_2) - u_1 ]
34
35 Using the definition of the Jacobian of F (from the PETSc user manual),
36 in the equation F(u',u,t) = G(u,t),
37
38 dF dF
39 J(F) = a * -- - --
40 du' du
41
42 where d is the partial derivative. In this example,
43
44 dF [ 1 ; 0 ]
45 -- = [ ]
46 du' [ 0 ; 0 ]
47
48 dF [ 0 ; 1 ]
49 -- = [ ]
50 du [ -1 ; 1 - u_2^2 ]
51
52 Hence,
53
54 [ a ; -1 ]
55 J(F) = [ ]
56 [ 1 ; u_2^2 - 1 ]
57
58 ------------------------------------------------------------------------- */
59
60 #include <petscts.h>
61
62 typedef struct _n_User *User;
63 struct _n_User {
64 PetscReal next_output;
65 };
66
67 /*
68 User-defined routines
69 */
70
IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,PetscCtx ctx)71 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx)
72 {
73 PetscScalar *f;
74 const PetscScalar *x, *xdot;
75
76 PetscFunctionBeginUser;
77 PetscCall(VecGetArrayRead(X, &x));
78 PetscCall(VecGetArrayRead(Xdot, &xdot));
79 PetscCall(VecGetArray(F, &f));
80 f[0] = xdot[0] + x[1];
81 f[1] = (x[1] * x[1] * x[1] / 3.0 - x[1]) - x[0];
82 PetscCall(VecRestoreArrayRead(X, &x));
83 PetscCall(VecRestoreArrayRead(Xdot, &xdot));
84 PetscCall(VecRestoreArray(F, &f));
85 PetscFunctionReturn(PETSC_SUCCESS);
86 }
87
IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,PetscCtx ctx)88 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
89 {
90 PetscInt rowcol[] = {0, 1};
91 PetscScalar J[2][2];
92 const PetscScalar *x;
93
94 PetscFunctionBeginUser;
95 PetscCall(VecGetArrayRead(X, &x));
96 J[0][0] = a;
97 J[0][1] = -1.;
98 J[1][0] = 1.;
99 J[1][1] = -1. + x[1] * x[1];
100 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
101 PetscCall(VecRestoreArrayRead(X, &x));
102
103 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
104 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
105 if (A != B) {
106 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
107 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
108 }
109 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111
RegisterMyARK2(void)112 static PetscErrorCode RegisterMyARK2(void)
113 {
114 PetscFunctionBeginUser;
115 {
116 const PetscReal A[3][3] =
117 {
118 {0, 0, 0},
119 {0.41421356237309504880, 0, 0},
120 {0.75, 0.25, 0}
121 },
122 At[3][3] = {{0, 0, 0}, {0.12132034355964257320, 0.29289321881345247560, 0}, {0.20710678118654752440, 0.50000000000000000000, 0.29289321881345247560}}, *bembedt = NULL, *bembed = NULL;
123 PetscCall(TSARKIMEXRegister("myark2", 2, 3, &At[0][0], NULL, NULL, &A[0][0], NULL, NULL, bembedt, bembed, 0, NULL, NULL));
124 }
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127
128 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
Monitor(TS ts,PetscInt step,PetscReal t,Vec X,PetscCtx ctx)129 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, PetscCtx ctx)
130 {
131 const PetscScalar *x;
132 PetscReal tfinal, dt;
133 User user = (User)ctx;
134 Vec interpolatedX;
135
136 PetscFunctionBeginUser;
137 PetscCall(TSGetTimeStep(ts, &dt));
138 PetscCall(TSGetMaxTime(ts, &tfinal));
139
140 while (user->next_output <= t && user->next_output <= tfinal) {
141 PetscCall(VecDuplicate(X, &interpolatedX));
142 PetscCall(TSInterpolate(ts, user->next_output, interpolatedX));
143 PetscCall(VecGetArrayRead(interpolatedX, &x));
144 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%.1f] %3" PetscInt_FMT " TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(x[0]), (double)PetscRealPart(x[1])));
145 PetscCall(VecRestoreArrayRead(interpolatedX, &x));
146 PetscCall(VecDestroy(&interpolatedX));
147 user->next_output += PetscRealConstant(0.1);
148 }
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151
main(int argc,char ** argv)152 int main(int argc, char **argv)
153 {
154 TS ts; /* nonlinear solver */
155 Vec x; /* solution, residual vectors */
156 Mat A; /* Jacobian matrix */
157 PetscInt steps;
158 PetscReal ftime = 0.5;
159 PetscBool monitor = PETSC_FALSE;
160 PetscScalar *x_ptr;
161 PetscMPIInt size;
162 struct _n_User user;
163
164 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165 Initialize program
166 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167 PetscFunctionBeginUser;
168 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
169 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
170 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
171
172 PetscCall(RegisterMyARK2());
173
174 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175 Set runtime options
176 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177
178 user.next_output = 0.0;
179 PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
180
181 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182 Create necessary matrix and vectors, solve same ODE on every process
183 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
185 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
186 PetscCall(MatSetFromOptions(A));
187 PetscCall(MatSetUp(A));
188 PetscCall(MatCreateVecs(A, &x, NULL));
189
190 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191 Create timestepping solver context
192 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
193 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
194 PetscCall(TSSetType(ts, TSBEULER));
195 PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
196 PetscCall(TSSetIJacobian(ts, A, A, IJacobian, &user));
197 PetscCall(TSSetMaxTime(ts, ftime));
198 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
199 if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL));
200
201 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202 Set initial conditions
203 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204 PetscCall(VecGetArray(x, &x_ptr));
205 x_ptr[0] = -2;
206 x_ptr[1] = -2.355301397608119909925287735864250951918;
207 PetscCall(VecRestoreArray(x, &x_ptr));
208 PetscCall(TSSetTimeStep(ts, .001));
209
210 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211 Set runtime options
212 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
213 PetscCall(TSSetFromOptions(ts));
214
215 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216 Solve nonlinear system
217 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
218 PetscCall(TSSolve(ts, x));
219 PetscCall(TSGetSolveTime(ts, &ftime));
220 PetscCall(TSGetStepNumber(ts, &steps));
221 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "steps %3" PetscInt_FMT ", ftime %g\n", steps, (double)ftime));
222 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
223
224 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225 Free work space. All PETSc objects should be destroyed when they
226 are no longer needed.
227 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
228 PetscCall(MatDestroy(&A));
229 PetscCall(VecDestroy(&x));
230 PetscCall(TSDestroy(&ts));
231
232 PetscCall(PetscFinalize());
233 return 0;
234 }
235
236 /*TEST
237
238 test:
239 requires: !single
240 suffix: a
241 args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp
242 output_file: output/ex19_pi42.out
243
244 test:
245 requires: !single
246 suffix: b
247 args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp -ts_adapt_dsp_filter PI42
248 output_file: output/ex19_pi42.out
249
250 test:
251 requires: !single
252 suffix: c
253 args: -monitor -ts_type bdf -ts_rtol 1e-6 -ts_atol 1e-6 -ts_view -ts_adapt_type dsp -ts_adapt_dsp_pid 0.4,0.2
254 output_file: output/ex19_pi42.out
255
256 test:
257 requires: !single
258 suffix: bdf_reject
259 args: -ts_type bdf -ts_time_step 0.5 -ts_max_steps 1 -ts_max_step_rejections {{0 1 2}separate_output} -ts_error_if_step_fails false -ts_adapt_monitor
260
261 TEST*/
262