Actual source code: ex16.c


  2: static char help[] = "Solves the van der Pol equation and demonstrate IMEX.\n\
  3: Input parameters include:\n\
  4:       -mu : stiffness parameter\n\n";

  6: /*
  7:    Concepts: TS^time-dependent nonlinear problems
  8:    Concepts: TS^van der Pol equation
  9:    Processors: 1
 10: */
 11: /* ------------------------------------------------------------------------

 13:    This program solves the van der Pol equation
 14:        y'' - \mu ((1-y^2)*y' - y) = 0        (1)
 15:    on the domain 0 <= x <= 1, with the boundary conditions
 16:        y(0) = 2, y'(0) = - 2/3 +10/(81*\mu) - 292/(2187*\mu^2),
 17:    This is a nonlinear equation. The well prepared initial condition gives errors that are not dominated by the first few steps of the method when \mu is large.

 19:    Notes:
 20:    This code demonstrates the TS solver interface to two variants of
 21:    linear problems, u_t = f(u,t), namely turning (1) into a system of
 22:    first order differential equations,

 24:    [ y' ] = [          z            ]
 25:    [ z' ]   [ \mu ((1 - y^2) z - y) ]

 27:    which then we can write as a vector equation

 29:    [ u_1' ] = [             u_2           ]  (2)
 30:    [ u_2' ]   [ \mu (1 - u_1^2) u_2 - u_1 ]

 32:    which is now in the desired form of u_t = f(u,t). One way that we
 33:    can split f(u,t) in (2) is to split by component,

 35:    [ u_1' ] = [ u_2 ] + [            0                ]
 36:    [ u_2' ]   [  0  ]   [ \mu ((1 - u_1^2) u_2 - u_1) ]

 38:    where

 40:    [ G(u,t) ] = [ u_2 ]
 41:                 [  0  ]

 43:    and

 45:    [ F(u',u,t) ] = [ u_1' ] - [            0                ]
 46:                    [ u_2' ]   [ \mu ((1 - u_1^2) u_2 - u_1) ]

 48:    Using the definition of the Jacobian of F (from the PETSc user manual),
 49:    in the equation F(u',u,t) = G(u,t),

 51:               dF   dF
 52:    J(F) = a * -- - --
 53:               du'  du

 55:    where d is the partial derivative. In this example,

 57:    dF   [ 1 ; 0 ]
 58:    -- = [       ]
 59:    du'  [ 0 ; 1 ]

 61:    dF   [       0             ;         0        ]
 62:    -- = [                                        ]
 63:    du   [ -\mu (2*u_1*u_2 + 1);  \mu (1 - u_1^2) ]

 65:    Hence,

 67:           [      a             ;          0          ]
 68:    J(F) = [                                          ]
 69:           [ \mu (2*u_1*u_2 + 1); a - \mu (1 - u_1^2) ]

 71:   ------------------------------------------------------------------------- */

 73: #include <petscts.h>

 75: typedef struct _n_User *User;
 76: struct _n_User {
 77:   PetscReal mu;
 78:   PetscBool imex;
 79:   PetscReal next_output;
 80: };

 82: /*
 83: *  User-defined routines
 84: */
 85: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx)
 86: {
 87:   PetscErrorCode    ierr;
 88:   User              user = (User)ctx;
 89:   PetscScalar       *f;
 90:   const PetscScalar *x;

 93:   VecGetArrayRead(X,&x);
 94:   VecGetArray(F,&f);
 95:   f[0] = (user->imex ? x[1] : 0);
 96:   f[1] = 0.0;
 97:   VecRestoreArrayRead(X,&x);
 98:   VecRestoreArray(F,&f);
 99:   return(0);
100: }

102: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
103: {
104:   PetscErrorCode    ierr;
105:   User              user = (User)ctx;
106:   const PetscScalar *x,*xdot;
107:   PetscScalar       *f;

110:   VecGetArrayRead(X,&x);
111:   VecGetArrayRead(Xdot,&xdot);
112:   VecGetArray(F,&f);
113:   f[0] = xdot[0] + (user->imex ? 0 : x[1]);
114:   f[1] = xdot[1] - user->mu*((1. - x[0]*x[0])*x[1] - x[0]);
115:   VecRestoreArrayRead(X,&x);
116:   VecRestoreArrayRead(Xdot,&xdot);
117:   VecRestoreArray(F,&f);
118:   return(0);
119: }

121: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
122: {
123:   PetscErrorCode    ierr;
124:   User              user     = (User)ctx;
125:   PetscReal         mu       = user->mu;
126:   PetscInt          rowcol[] = {0,1};
127:   const PetscScalar *x;
128:   PetscScalar       J[2][2];

131:   VecGetArrayRead(X,&x);
132:   J[0][0] = a;                    J[0][1] = (user->imex ? 0 : 1.);
133:   J[1][0] = mu*(2.*x[0]*x[1]+1.);   J[1][1] = a - mu*(1. - x[0]*x[0]);
134:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
135:   VecRestoreArrayRead(X,&x);

137:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
138:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
139:   if (A != B) {
140:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
141:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
142:   }
143:   return(0);
144: }

146: static PetscErrorCode RegisterMyARK2(void)
147: {

151:   {
152:     const PetscReal
153:       A[3][3] = {{0,0,0},
154:                  {0.41421356237309504880,0,0},
155:                  {0.75,0.25,0}},
156:       At[3][3] = {{0,0,0},
157:                   {0.12132034355964257320,0.29289321881345247560,0},
158:                   {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}},
159:       *bembedt = NULL,*bembed = NULL;
160:     TSARKIMEXRegister("myark2",2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembed,0,NULL,NULL);
161:   }
162:   return(0);
163: }

165: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
166: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
167: {
168:   PetscErrorCode    ierr;
169:   const PetscScalar *x;
170:   PetscReal         tfinal, dt;
171:   User              user = (User)ctx;
172:   Vec               interpolatedX;

175:   TSGetTimeStep(ts,&dt);
176:   TSGetMaxTime(ts,&tfinal);

178:   while (user->next_output <= t && user->next_output <= tfinal) {
179:     VecDuplicate(X,&interpolatedX);
180:     TSInterpolate(ts,user->next_output,interpolatedX);
181:     VecGetArrayRead(interpolatedX,&x);
182:     PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n",user->next_output,step,t,dt,(double)PetscRealPart(x[0]),(double)PetscRealPart(x[1]));
183:     VecRestoreArrayRead(interpolatedX,&x);
184:     VecDestroy(&interpolatedX);

186:     user->next_output += 0.1;
187:   }
188:   return(0);
189: }

191: int main(int argc,char **argv)
192: {
193:   TS             ts;            /* nonlinear solver */
194:   Vec            x;             /* solution, residual vectors */
195:   Mat            A;             /* Jacobian matrix */
196:   PetscInt       steps;
197:   PetscReal      ftime = 0.5;
198:   PetscBool      monitor = PETSC_FALSE;
199:   PetscScalar    *x_ptr;
200:   PetscMPIInt    size;
201:   struct _n_User user;

204:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205:      Initialize program
206:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
207:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
208:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
209:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");

211:   RegisterMyARK2();

213:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214:     Set runtime options
215:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216:   user.mu          = 1000.0;
217:   user.imex        = PETSC_TRUE;
218:   user.next_output = 0.0;

220:   PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
221:   PetscOptionsGetBool(NULL,NULL,"-imex",&user.imex,NULL);
222:   PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);

224:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225:     Create necessary matrix and vectors, solve same ODE on every process
226:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
227:   MatCreate(PETSC_COMM_WORLD,&A);
228:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
229:   MatSetFromOptions(A);
230:   MatSetUp(A);
231:   MatCreateVecs(A,&x,NULL);

233:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234:      Create timestepping solver context
235:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
236:   TSCreate(PETSC_COMM_WORLD,&ts);
237:   TSSetType(ts,TSBEULER);
238:   TSSetRHSFunction(ts,NULL,RHSFunction,&user);
239:   TSSetIFunction(ts,NULL,IFunction,&user);
240:   TSSetIJacobian(ts,A,A,IJacobian,&user);
241:   TSSetMaxTime(ts,ftime);
242:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
243:   if (monitor) {
244:     TSMonitorSet(ts,Monitor,&user,NULL);
245:   }

247:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
248:      Set initial conditions
249:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
250:   VecGetArray(x,&x_ptr);
251:   x_ptr[0] = 2.0;
252:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
253:   VecRestoreArray(x,&x_ptr);
254:   TSSetTimeStep(ts,0.01);

256:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257:      Set runtime options
258:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
259:   TSSetFromOptions(ts);

261:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
262:      Solve nonlinear system
263:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
264:   TSSolve(ts,x);
265:   TSGetSolveTime(ts,&ftime);
266:   TSGetStepNumber(ts,&steps);
267:   PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %D, ftime %g\n",(double)user.mu,steps,(double)ftime);
268:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

270:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
271:      Free work space.  All PETSc objects should be destroyed when they
272:      are no longer needed.
273:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274:   MatDestroy(&A);
275:   VecDestroy(&x);
276:   TSDestroy(&ts);

278:   PetscFinalize();
279:   return ierr;
280: }

282: /*TEST

284:     test:
285:       args: -ts_type arkimex -ts_arkimex_type myark2 -ts_adapt_type none
286:       requires: !single

288: TEST*/