Actual source code: ct_vdp_imex.c
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";
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 splitted 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: */
29: #include <petscdmda.h>
30: #include <petscts.h>
32: typedef struct _User *User;
33: struct _User {
34: PetscReal mu; /*stiffness control coefficient: epsilon*/
35: };
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 *);
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;
54: 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: PetscOptionsReal("-eps", "Stiffness controller", "", user.mu, &user.mu, NULL);
59: PetscOptionsEnd();
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Set runtime options
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: /*
65: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
66: */
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create necessary matrix and vectors, solve same ODE on every process
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: MatCreate(PETSC_COMM_WORLD, &A);
72: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2);
73: MatSetFromOptions(A);
74: MatSetUp(A);
75: MatCreateVecs(A, &x, NULL);
77: MatCreateVecs(A, &ref, NULL);
78: 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;*/
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;
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create timestepping solver context
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: TSCreate(PETSC_COMM_WORLD, &ts);
103: TSSetType(ts, TSEIMEX);
104: TSSetRHSFunction(ts, NULL, RHSFunction, &user);
105: TSSetIFunction(ts, NULL, IFunction, &user);
106: TSSetIJacobian(ts, A, A, IJacobian, &user);
108: dt = 0.00001;
109: ftime = 1.1;
110: TSSetTimeStep(ts, dt);
111: TSSetMaxTime(ts, ftime);
112: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Set initial conditions
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: 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: TSSetSolution(ts, x);
120: VecGetSize(x, &mx);
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Set runtime options
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: TSSetFromOptions(ts);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Solve nonlinear system
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: TSSolve(ts, x);
131: TSGetTime(ts, &ftime);
132: TSGetStepNumber(ts, &steps);
134: VecAXPY(x, -1.0, ref);
135: VecNorm(x, NORM_2, &norm);
136: TSGetTimeStep(ts, &dt);
138: eimex_rowcol[0] = 0;
139: eimex_rowcol[1] = 0;
140: two = 2;
141: PetscOptionsGetIntArray(NULL, NULL, "-ts_eimex_row_col", eimex_rowcol, &two, NULL);
142: PetscPrintf(PETSC_COMM_WORLD, "order %11s %18s %37s\n", "dt", "norm", "final solution components 0 and 1");
143: VecGetArray(x, &x_ptr);
144: 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: VecRestoreArray(x, &x_ptr);
147: /* Write line in convergence log */
148: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
149: PetscViewerSetType(viewer, PETSCVIEWERASCII);
150: PetscViewerFileSetMode(viewer, FILE_MODE_APPEND);
151: PetscViewerFileSetName(viewer, "eimex_nonstiff_vdp.txt");
152: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %10.8f %18.15f\n", eimex_rowcol[0], eimex_rowcol[1], (double)dt, (double)norm);
153: PetscViewerDestroy(&viewer);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Free work space.
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: MatDestroy(&A);
159: VecDestroy(&x);
160: VecDestroy(&ref);
161: TSDestroy(&ts);
162: PetscFinalize();
163: return 0;
164: }
166: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
167: {
168: PetscScalar *f;
169: const PetscScalar *x;
171: VecGetArrayRead(X, &x);
172: VecGetArray(F, &f);
173: f[0] = x[1];
174: f[1] = 0.0;
175: VecRestoreArrayRead(X, &x);
176: VecRestoreArray(F, &f);
177: return 0;
178: }
180: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
181: {
182: User user = (User)ptr;
183: PetscScalar *f;
184: const PetscScalar *x, *xdot;
186: VecGetArrayRead(X, &x);
187: VecGetArrayRead(Xdot, &xdot);
188: VecGetArray(F, &f);
189: f[0] = xdot[0];
190: f[1] = xdot[1] - ((1. - x[0] * x[0]) * x[1] - x[0]) / user->mu;
191: VecRestoreArrayRead(X, &x);
192: VecRestoreArrayRead(Xdot, &xdot);
193: VecRestoreArray(F, &f);
194: return 0;
195: }
197: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ptr)
198: {
199: User user = (User)ptr;
200: PetscReal mu = user->mu;
201: PetscInt rowcol[] = {0, 1};
202: PetscScalar J[2][2];
203: const PetscScalar *x;
205: VecGetArrayRead(X, &x);
206: J[0][0] = a;
207: J[0][1] = 0;
208: J[1][0] = (2. * x[0] * x[1] + 1.) / mu;
209: J[1][1] = a - (1. - x[0] * x[0]) / mu;
210: MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
211: VecRestoreArrayRead(X, &x);
213: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
214: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
215: if (A != B) {
216: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
217: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
218: }
219: return 0;
220: }
222: /*TEST
224: test:
225: args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_row_col 3,3 -ts_monitor_lg_solution
226: requires: x
228: test:
229: suffix: adapt
230: args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_order_adapt -ts_eimex_max_rows 7 -ts_monitor_lg_solution
231: requires: x
233: test:
234: suffix: loop
235: args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt {{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}
236: requires: x
238: TEST*/