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;
173:   PetscCall(VecGetArray(F, &f));
174:   f[0] = x[1];
175:   f[1] = 0.0;
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;
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;
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;
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));

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: