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:   User              user = (User)ctx;
 88:   PetscScalar       *f;
 89:   const PetscScalar *x;

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

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

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

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

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

134:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
135:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
136:   if (A != B) {
137:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
138:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
139:   }
140:   return 0;
141: }

143: static PetscErrorCode RegisterMyARK2(void)
144: {
146:   {
147:     const PetscReal
148:       A[3][3] = {{0,0,0},
149:                  {0.41421356237309504880,0,0},
150:                  {0.75,0.25,0}},
151:       At[3][3] = {{0,0,0},
152:                   {0.12132034355964257320,0.29289321881345247560,0},
153:                   {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}},
154:       *bembedt = NULL,*bembed = NULL;
155:     TSARKIMEXRegister("myark2",2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembed,0,NULL,NULL);
156:   }
157:   return 0;
158: }

160: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
161: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
162: {
163:   const PetscScalar *x;
164:   PetscReal         tfinal, dt;
165:   User              user = (User)ctx;
166:   Vec               interpolatedX;

169:   TSGetTimeStep(ts,&dt);
170:   TSGetMaxTime(ts,&tfinal);

172:   while (user->next_output <= t && user->next_output <= tfinal) {
173:     VecDuplicate(X,&interpolatedX);
174:     TSInterpolate(ts,user->next_output,interpolatedX);
175:     VecGetArrayRead(interpolatedX,&x);
176:     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]));
177:     VecRestoreArrayRead(interpolatedX,&x);
178:     VecDestroy(&interpolatedX);

180:     user->next_output += 0.1;
181:   }
182:   return 0;
183: }

185: int main(int argc,char **argv)
186: {
187:   TS             ts;            /* nonlinear solver */
188:   Vec            x;             /* solution, residual vectors */
189:   Mat            A;             /* Jacobian matrix */
190:   PetscInt       steps;
191:   PetscReal      ftime = 0.5;
192:   PetscBool      monitor = PETSC_FALSE;
193:   PetscScalar    *x_ptr;
194:   PetscMPIInt    size;
195:   struct _n_User user;

197:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198:      Initialize program
199:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200:   PetscInitialize(&argc,&argv,NULL,help);
201:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

204:   RegisterMyARK2();

206:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207:     Set runtime options
208:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
209:   user.mu          = 1000.0;
210:   user.imex        = PETSC_TRUE;
211:   user.next_output = 0.0;

213:   PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
214:   PetscOptionsGetBool(NULL,NULL,"-imex",&user.imex,NULL);
215:   PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);

217:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218:     Create necessary matrix and vectors, solve same ODE on every process
219:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
220:   MatCreate(PETSC_COMM_WORLD,&A);
221:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
222:   MatSetFromOptions(A);
223:   MatSetUp(A);
224:   MatCreateVecs(A,&x,NULL);

226:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227:      Create timestepping solver context
228:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
229:   TSCreate(PETSC_COMM_WORLD,&ts);
230:   TSSetType(ts,TSBEULER);
231:   TSSetRHSFunction(ts,NULL,RHSFunction,&user);
232:   TSSetIFunction(ts,NULL,IFunction,&user);
233:   TSSetIJacobian(ts,A,A,IJacobian,&user);
234:   TSSetMaxTime(ts,ftime);
235:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
236:   if (monitor) {
237:     TSMonitorSet(ts,Monitor,&user,NULL);
238:   }

240:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241:      Set initial conditions
242:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
243:   VecGetArray(x,&x_ptr);
244:   x_ptr[0] = 2.0;
245:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
246:   VecRestoreArray(x,&x_ptr);
247:   TSSetTimeStep(ts,0.01);

249:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
250:      Set runtime options
251:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
252:   TSSetFromOptions(ts);

254:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255:      Solve nonlinear system
256:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
257:   TSSolve(ts,x);
258:   TSGetSolveTime(ts,&ftime);
259:   TSGetStepNumber(ts,&steps);
260:   PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %D, ftime %g\n",(double)user.mu,steps,(double)ftime);
261:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

263:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264:      Free work space.  All PETSc objects should be destroyed when they
265:      are no longer needed.
266:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
267:   MatDestroy(&A);
268:   VecDestroy(&x);
269:   TSDestroy(&ts);

271:   PetscFinalize();
272:   return 0;
273: }

275: /*TEST

277:     test:
278:       args: -ts_type arkimex -ts_arkimex_type myark2 -ts_adapt_type none
279:       requires: !single

281: TEST*/