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*/