1 /*
2 * ex_vdp.c
3 *
4 * Created on: Jun 1, 2012
5 * Author: Hong Zhang
6 */
7 static char help[] = "Solves the van der Pol equation. \n Input parameters include:\n";
8
9 /*
10 * This program solves the van der Pol equation
11 * y' = z (1)
12 * z' = (((1-y^2)*z-y)/eps (2)
13 * on the domain 0<=x<=0.5, with the initial conditions
14 * y(0) = 2,
15 * z(0) = -2/3 + 10/81*eps - 292/2187*eps^2-1814/19683*eps^3
16 * IMEX schemes are applied to the split equation
17 * [y'] = [z] + [0 ]
18 * [z'] [0] [(((1-y^2)*z-y)/eps]
19 *
20 * F(x)= [z]
21 * [0]
22 *
23 * G(x)= [y'] - [0 ]
24 * [z'] [(((1-y^2)*z-y)/eps]
25 *
26 * JG(x) = G_x + a G_xdot
27 */
28
29 #include <petscdmda.h>
30 #include <petscts.h>
31
32 typedef struct _User *User;
33 struct _User {
34 PetscReal mu; /*stiffness control coefficient: epsilon*/
35 };
36
37 static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
38 static PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *);
39 static PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
40
main(int argc,char ** argv)41 int main(int argc, char **argv)
42 {
43 TS ts;
44 Vec x; /* solution vector */
45 Mat A; /* Jacobian */
46 PetscInt steps, mx, eimex_rowcol[2], two;
47 PetscScalar *x_ptr;
48 PetscReal ftime, dt, norm;
49 Vec ref;
50 struct _User user; /* user-defined work context */
51 PetscViewer viewer;
52
53 PetscFunctionBeginUser;
54 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
55 /* Initialize user application context */
56 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "van der Pol options", "");
57 user.mu = 1e0;
58 PetscCall(PetscOptionsReal("-eps", "Stiffness controller", "", user.mu, &user.mu, NULL));
59 PetscOptionsEnd();
60
61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62 Set runtime options
63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64 /*
65 PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL));
66 */
67
68 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69 Create necessary matrix and vectors, solve same ODE on every process
70 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
72 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
73 PetscCall(MatSetFromOptions(A));
74 PetscCall(MatSetUp(A));
75 PetscCall(MatCreateVecs(A, &x, NULL));
76
77 PetscCall(MatCreateVecs(A, &ref, NULL));
78 PetscCall(VecGetArray(ref, &x_ptr));
79 /*
80 * [0,1], mu=10^-3
81 */
82 /*
83 x_ptr[0] = -1.8881254106283;
84 x_ptr[1] = 0.7359074233370;*/
85
86 /*
87 * [0,0.5],mu=10^-3
88 */
89 /*
90 x_ptr[0] = 1.596980778659137;
91 x_ptr[1] = -1.029103015879544;
92 */
93 /*
94 * [0,0.5],mu=1
95 */
96 x_ptr[0] = 1.619084329683235;
97 x_ptr[1] = -0.803530465176385;
98
99 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100 Create timestepping solver context
101 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
103 PetscCall(TSSetType(ts, TSEIMEX));
104 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user));
105 PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
106 PetscCall(TSSetIJacobian(ts, A, A, IJacobian, &user));
107
108 dt = 0.00001;
109 ftime = 1.1;
110 PetscCall(TSSetTimeStep(ts, dt));
111 PetscCall(TSSetMaxTime(ts, ftime));
112 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114 Set initial conditions
115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116 PetscCall(VecGetArray(x, &x_ptr));
117 x_ptr[0] = 2.;
118 x_ptr[1] = -2. / 3. + 10. / 81. * (user.mu) - 292. / 2187. * (user.mu) * (user.mu) - 1814. / 19683. * (user.mu) * (user.mu) * (user.mu);
119 PetscCall(TSSetSolution(ts, x));
120 PetscCall(VecGetSize(x, &mx));
121
122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123 Set runtime options
124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125 PetscCall(TSSetFromOptions(ts));
126
127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128 Solve nonlinear system
129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130 PetscCall(TSSolve(ts, x));
131 PetscCall(TSGetTime(ts, &ftime));
132 PetscCall(TSGetStepNumber(ts, &steps));
133
134 PetscCall(VecAXPY(x, -1.0, ref));
135 PetscCall(VecNorm(x, NORM_2, &norm));
136 PetscCall(TSGetTimeStep(ts, &dt));
137
138 eimex_rowcol[0] = 0;
139 eimex_rowcol[1] = 0;
140 two = 2;
141 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-ts_eimex_row_col", eimex_rowcol, &two, NULL));
142 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "order %11s %18s %37s\n", "dt", "norm", "final solution components 0 and 1"));
143 PetscCall(VecGetArray(x, &x_ptr));
144 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(%" PetscInt_FMT ",%" PetscInt_FMT ") %10.8f %18.15f %18.15f %18.15f\n", eimex_rowcol[0], eimex_rowcol[1], (double)dt, (double)norm, (double)PetscRealPart(x_ptr[0]), (double)PetscRealPart(x_ptr[1])));
145 PetscCall(VecRestoreArray(x, &x_ptr));
146
147 /* Write line in convergence log */
148 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
149 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII));
150 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_APPEND));
151 PetscCall(PetscViewerFileSetName(viewer, "eimex_nonstiff_vdp.txt"));
152 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %10.8f %18.15f\n", eimex_rowcol[0], eimex_rowcol[1], (double)dt, (double)norm));
153 PetscCall(PetscViewerDestroy(&viewer));
154
155 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156 Free work space.
157 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158 PetscCall(MatDestroy(&A));
159 PetscCall(VecDestroy(&x));
160 PetscCall(VecDestroy(&ref));
161 PetscCall(TSDestroy(&ts));
162 PetscCall(PetscFinalize());
163 return 0;
164 }
165
RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)166 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
167 {
168 PetscScalar *f;
169 const PetscScalar *x;
170
171 PetscFunctionBegin;
172 PetscCall(VecGetArrayRead(X, &x));
173 PetscCall(VecGetArray(F, &f));
174 f[0] = x[1];
175 f[1] = 0.0;
176 PetscCall(VecRestoreArrayRead(X, &x));
177 PetscCall(VecRestoreArray(F, &f));
178 PetscFunctionReturn(PETSC_SUCCESS);
179 }
180
IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void * ptr)181 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
182 {
183 User user = (User)ptr;
184 PetscScalar *f;
185 const PetscScalar *x, *xdot;
186
187 PetscFunctionBegin;
188 PetscCall(VecGetArrayRead(X, &x));
189 PetscCall(VecGetArrayRead(Xdot, &xdot));
190 PetscCall(VecGetArray(F, &f));
191 f[0] = xdot[0];
192 f[1] = xdot[1] - ((1. - x[0] * x[0]) * x[1] - x[0]) / user->mu;
193 PetscCall(VecRestoreArrayRead(X, &x));
194 PetscCall(VecRestoreArrayRead(Xdot, &xdot));
195 PetscCall(VecRestoreArray(F, &f));
196 PetscFunctionReturn(PETSC_SUCCESS);
197 }
198
IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void * ptr)199 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ptr)
200 {
201 User user = (User)ptr;
202 PetscReal mu = user->mu;
203 PetscInt rowcol[] = {0, 1};
204 PetscScalar J[2][2];
205 const PetscScalar *x;
206
207 PetscFunctionBegin;
208 PetscCall(VecGetArrayRead(X, &x));
209 J[0][0] = a;
210 J[0][1] = 0;
211 J[1][0] = (2. * x[0] * x[1] + 1.) / mu;
212 J[1][1] = a - (1. - x[0] * x[0]) / mu;
213 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
214 PetscCall(VecRestoreArrayRead(X, &x));
215
216 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
217 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
218 if (A != B) {
219 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
220 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
221 }
222 PetscFunctionReturn(PETSC_SUCCESS);
223 }
224
225 /*TEST
226
227 test:
228 args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_time_step 0.01 -ts_max_time 10 -ts_eimex_row_col 3,3 -ts_monitor_lg_solution
229 requires: x
230
231 test:
232 suffix: adapt
233 args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_time_step 0.01 -ts_max_time 10 -ts_eimex_order_adapt -ts_eimex_max_rows 7 -ts_monitor_lg_solution
234 requires: x
235
236 test:
237 suffix: loop
238 args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_time_step {{0.005 0.001 0.0005}separate output} -ts_max_steps {{100 500 1000}separate output} -ts_eimex_row_col {{1,1 2,1 3,1 2,2 3,2 3,3}separate output}
239 requires: x
240
241 TEST*/
242