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:  * Processors:1
 11:  */

 13: /*
 14:  * This program solves the van der Pol equation
 15:  * y' = z                               (1)
 16:  * z' = (((1-y^2)*z-y)/eps              (2)
 17:  * on the domain 0<=x<=0.5, with the initial conditions
 18:  * y(0) = 2,
 19:  * z(0) = -2/3 + 10/81*eps - 292/2187*eps^2-1814/19683*eps^3
 20:  * IMEX schemes are applied to the splitted equation
 21:  * [y'] = [z]  + [0                 ]
 22:  * [z']   [0]    [(((1-y^2)*z-y)/eps]
 23:  *
 24:  * F(x)= [z]
 25:  *       [0]
 26:  *
 27:  * G(x)= [y'] -   [0                 ]
 28:  *       [z']     [(((1-y^2)*z-y)/eps]
 29:  *
 30:  * JG(x) =  G_x + a G_xdot
 31:  */

 33: #include <petscdmda.h>
 34: #include <petscts.h>

 36: typedef struct _User *User;
 37: struct _User {
 38:   PetscReal mu;  /*stiffness control coefficient: epsilon*/
 39: };

 41: static PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
 42: static PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*);
 43: static PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

 45: int main(int argc, char **argv)
 46: {
 47:   TS                ts;
 48:   Vec               x; /*solution vector*/
 49:   Mat               A; /*Jacobian*/
 50:   PetscInt          steps,mx,eimex_rowcol[2],two;
 51:   PetscErrorCode    ierr;
 52:   PetscScalar       *x_ptr;
 53:   PetscReal         ftime,dt,norm;
 54:   Vec               ref;
 55:   struct _User      user;       /* user-defined work context */
 56:   PetscViewer       viewer;

 58:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
 59:   /* Initialize user application context */
 60:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"van der Pol options","");
 61:   user.mu      = 1e0;
 62:   PetscOptionsReal("-eps","Stiffness controller","",user.mu,&user.mu,NULL);
 63:   PetscOptionsEnd();

 65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 66:    Set runtime options
 67:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 68:   /*
 69:    PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
 70:    */

 72:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 73:    Create necessary matrix and vectors, solve same ODE on every process
 74:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 75:   MatCreate(PETSC_COMM_WORLD,&A);
 76:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
 77:   MatSetFromOptions(A);
 78:   MatSetUp(A);
 79:   MatCreateVecs(A,&x,NULL);

 81:   MatCreateVecs(A,&ref,NULL);
 82:   VecGetArray(ref,&x_ptr);
 83:   /*
 84:    * [0,1], mu=10^-3
 85:    */
 86:   /*
 87:    x_ptr[0] = -1.8881254106283;
 88:    x_ptr[1] =  0.7359074233370;*/

 90:   /*
 91:    * [0,0.5],mu=10^-3
 92:    */
 93:   /*
 94:    x_ptr[0] = 1.596980778659137;
 95:    x_ptr[1] = -1.029103015879544;
 96:    */
 97:   /*
 98:    * [0,0.5],mu=1
 99:    */
100:   x_ptr[0] = 1.619084329683235;
101:   x_ptr[1] = -0.803530465176385;

103:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104:    Create timestepping solver context
105:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106:   TSCreate(PETSC_COMM_WORLD,&ts);
107:   TSSetType(ts,TSEIMEX);
108:   TSSetRHSFunction(ts,NULL,RHSFunction,&user);
109:   TSSetIFunction(ts,NULL,IFunction,&user);
110:   TSSetIJacobian(ts,A,A,IJacobian,&user);

112:   dt    = 0.00001;
113:   ftime = 1.1;
114:   TSSetTimeStep(ts,dt);
115:   TSSetMaxTime(ts,ftime);
116:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:    Set initial conditions
119:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120:   VecGetArray(x,&x_ptr);
121:   x_ptr[0] = 2.;
122:   x_ptr[1] = -2./3. + 10./81.*(user.mu) - 292./2187.* (user.mu) * (user.mu)
123:     -1814./19683.*(user.mu)*(user.mu)*(user.mu);
124:   TSSetSolution(ts,x);
125:   VecGetSize(x,&mx);

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:    Set runtime options
129:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   TSSetFromOptions(ts);

132:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133:    Solve nonlinear system
134:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135:   TSSolve(ts,x);
136:   TSGetTime(ts,&ftime);
137:   TSGetStepNumber(ts,&steps);

139:   VecAXPY(x,-1.0,ref);
140:   VecNorm(x,NORM_2,&norm);
141:   TSGetTimeStep(ts,&dt);

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

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

158:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159:    Free work space.
160:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161:   MatDestroy(&A);
162:   VecDestroy(&x);
163:   VecDestroy(&ref);
164:   TSDestroy(&ts);
165:   PetscFinalize();
166:   return ierr;
167: }

169: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
170: {
171:   PetscErrorCode    ierr;
172:   PetscScalar       *f;
173:   const PetscScalar *x;

176:   VecGetArrayRead(X,&x);
177:   VecGetArray(F,&f);
178:   f[0] = x[1];
179:   f[1] = 0.0;
180:   VecRestoreArrayRead(X,&x);
181:   VecRestoreArray(F,&f);
182:   return(0);
183: }

185: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
186: {
187:   User              user = (User)ptr;
188:   PetscScalar       *f;
189:   const PetscScalar *x,*xdot;
190:   PetscErrorCode    ierr;

193:   VecGetArrayRead(X,&x);
194:   VecGetArrayRead(Xdot,&xdot);
195:   VecGetArray(F,&f);
196:   f[0] = xdot[0];
197:   f[1] = xdot[1]-((1.-x[0]*x[0])*x[1]-x[0])/user->mu;
198:   VecRestoreArrayRead(X,&x);
199:   VecRestoreArrayRead(Xdot,&xdot);
200:   VecRestoreArray(F,&f);
201:   return(0);
202: }

204: static PetscErrorCode IJacobian(TS  ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ptr)
205: {
206:   PetscErrorCode    ierr;
207:   User              user = (User)ptr;
208:   PetscReal         mu = user->mu;
209:   PetscInt          rowcol[] = {0,1};
210:   PetscScalar       J[2][2];
211:   const PetscScalar *x;

214:   VecGetArrayRead(X,&x);
215:   J[0][0] = a;
216:   J[0][1] = 0;
217:   J[1][0] = (2.*x[0]*x[1]+1.)/mu;
218:   J[1][1] = a - (1. - x[0]*x[0])/mu;
219:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
220:   VecRestoreArrayRead(X,&x);

222:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
223:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
224:   if (A != B) {
225:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
226:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
227:   }
228:   return(0);
229: }

231: /*TEST

233:    test:
234:      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
235:      requires: x

237:    test:
238:      suffix: adapt
239:      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
240:      requires: x

242:    test:
243:      suffix: loop
244:      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}
245:      requires: x

247:  TEST*/