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:   PetscErrorCode    ierr;
 48:   PetscScalar       *x_ptr;
 49:   PetscReal         ftime,dt,norm;
 50:   Vec               ref;
 51:   struct _User      user;       /* user-defined work context */
 52:   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)
119:     -1814./19683.*(user.mu)*(user.mu)*(user.mu);
120:   TSSetSolution(ts,x);
121:   VecGetSize(x,&mx);

123:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124:    Set runtime options
125:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126:   TSSetFromOptions(ts);

128:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:    Solve nonlinear system
130:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131:   TSSolve(ts,x);
132:   TSGetTime(ts,&ftime);
133:   TSGetStepNumber(ts,&steps);

135:   VecAXPY(x,-1.0,ref);
136:   VecNorm(x,NORM_2,&norm);
137:   TSGetTimeStep(ts,&dt);

139:   eimex_rowcol[0] = 0; eimex_rowcol[1] = 0; two = 2;
140:   PetscOptionsGetIntArray(NULL,NULL,"-ts_eimex_row_col",eimex_rowcol,&two,NULL);
141:   PetscPrintf(PETSC_COMM_WORLD,"order %11s %18s %37s\n","dt","norm","final solution components 0 and 1");
142:   VecGetArray(x,&x_ptr);
143:   PetscPrintf(PETSC_COMM_WORLD,"(%D,%D) %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]));
144:   VecRestoreArray(x,&x_ptr);

146:   /* Write line in convergence log */
147:   PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
148:   PetscViewerSetType(viewer,PETSCVIEWERASCII);
149:   PetscViewerFileSetMode(viewer,FILE_MODE_APPEND);
150:   PetscViewerFileSetName(viewer,"eimex_nonstiff_vdp.txt");
151:   PetscViewerASCIIPrintf(viewer,"%D %D %10.8f %18.15f\n",eimex_rowcol[0],eimex_rowcol[1],(double)dt,(double)norm);
152:   PetscViewerDestroy(&viewer);

154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:    Free work space.
156:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157:   MatDestroy(&A);
158:   VecDestroy(&x);
159:   VecDestroy(&ref);
160:   TSDestroy(&ts);
161:   PetscFinalize();
162:   return 0;
163: }

165: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
166: {
167:   PetscScalar       *f;
168:   const PetscScalar *x;

170:   VecGetArrayRead(X,&x);
171:   VecGetArray(F,&f);
172:   f[0] = x[1];
173:   f[1] = 0.0;
174:   VecRestoreArrayRead(X,&x);
175:   VecRestoreArray(F,&f);
176:   return 0;
177: }

179: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
180: {
181:   User              user = (User)ptr;
182:   PetscScalar       *f;
183:   const PetscScalar *x,*xdot;

185:   VecGetArrayRead(X,&x);
186:   VecGetArrayRead(Xdot,&xdot);
187:   VecGetArray(F,&f);
188:   f[0] = xdot[0];
189:   f[1] = xdot[1]-((1.-x[0]*x[0])*x[1]-x[0])/user->mu;
190:   VecRestoreArrayRead(X,&x);
191:   VecRestoreArrayRead(Xdot,&xdot);
192:   VecRestoreArray(F,&f);
193:   return 0;
194: }

196: static PetscErrorCode IJacobian(TS  ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ptr)
197: {
198:   User              user = (User)ptr;
199:   PetscReal         mu = user->mu;
200:   PetscInt          rowcol[] = {0,1};
201:   PetscScalar       J[2][2];
202:   const PetscScalar *x;

204:   VecGetArrayRead(X,&x);
205:   J[0][0] = a;
206:   J[0][1] = 0;
207:   J[1][0] = (2.*x[0]*x[1]+1.)/mu;
208:   J[1][1] = a - (1. - x[0]*x[0])/mu;
209:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
210:   VecRestoreArrayRead(X,&x);

212:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
213:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
214:   if (A != B) {
215:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
216:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
217:   }
218:   return 0;
219: }

221: /*TEST

223:    test:
224:      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
225:      requires: x

227:    test:
228:      suffix: adapt
229:      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
230:      requires: x

232:    test:
233:      suffix: loop
234:      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}
235:      requires: x

237:  TEST*/