Actual source code: ct_vdp_imex.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  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*);


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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

187: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
188: {
189:   User              user = (User)ptr;
190:   PetscScalar       *f;
191:   const PetscScalar *x,*xdot;
192:   PetscErrorCode    ierr;
193:   
195:   VecGetArrayRead(X,&x);
196:   VecGetArrayRead(Xdot,&xdot);
197:   VecGetArray(F,&f);
198:   f[0] = xdot[0];
199:   f[1] = xdot[1]-((1.-x[0]*x[0])*x[1]-x[0])/user->mu;
200:   VecRestoreArrayRead(X,&x);
201:   VecRestoreArrayRead(Xdot,&xdot);
202:   VecRestoreArray(F,&f);
203:   return(0);
204: }

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

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

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

233: /*TEST

235:    test:
236:      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
237:      requires: x

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

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

249:  TEST*/