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 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: */
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;
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();
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Set runtime options
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: /*
65: PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL));
66: */
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));
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;*/
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: 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));
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));
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Set runtime options
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: PetscCall(TSSetFromOptions(ts));
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Solve nonlinear system
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: PetscCall(TSSolve(ts, x));
131: PetscCall(TSGetTime(ts, &ftime));
132: PetscCall(TSGetStepNumber(ts, &steps));
134: PetscCall(VecAXPY(x, -1.0, ref));
135: PetscCall(VecNorm(x, NORM_2, &norm));
136: PetscCall(TSGetTimeStep(ts, &dt));
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));
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));
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: }
166: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
167: {
168: PetscScalar *f;
169: const PetscScalar *x;
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: }
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;
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: }
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;
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));
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: }
225: /*TEST
227: test:
228: 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
229: requires: x
231: test:
232: suffix: adapt
233: 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
234: requires: x
236: test:
237: suffix: loop
238: 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}
239: requires: x
241: TEST*/