Actual source code: ex16.c

  2: /* Program usage:  ./ex16 [-help] [all PETSc options] */

  4: static char help[] = "Solves the van der Pol equation.\n\
  5: Input parameters include:\n\
  6:       -mu : stiffness parameter\n\n";

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

 15:    This program solves the van der Pol equation
 16:        y'' - \mu (1-y^2)*y' + y = 0        (1)
 17:    on the domain 0 <= x <= 1, with the boundary conditions
 18:        y(0) = 2, y'(0) = 0,
 19:    This is a nonlinear equation.

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

 26:    [ y' ] = [          z          ]
 27:    [ z' ]   [ \mu (1 - y^2) z - y ]

 29:    which then we can write as a vector equation

 31:    [ u_1' ] = [             u_2           ]  (2)
 32:    [ u_2' ]   [ \mu (1 - u_1^2) u_2 - u_1 ]

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

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

 40:    where

 42:    [ F(u,t) ] = [ u_2 ]
 43:                 [  0  ]

 45:    and

 47:    [ G(u',u,t) ] = [ u_1' ] - [            0              ]
 48:                    [ u_2' ]   [ \mu (1 - u_1^2) u_2 - u_1 ]

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

 53:               dG   dG
 54:    J(G) = a * -- - --
 55:               du'  du

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

 59:    dG   [ 1 ; 0 ]
 60:    -- = [       ]
 61:    du'  [ 0 ; 1 ]

 63:    dG   [       0       ;         0        ]
 64:    -- = [                                  ]
 65:    du   [ -2 \mu u_1 - 1;  \mu (1 - u_1^2) ]

 67:    Hence,

 69:           [      a       ;          0          ]
 70:    J(G) = [                                    ]
 71:           [ 2 \mu u_1 + 1; a - \mu (1 - u_1^2) ]

 73:   ------------------------------------------------------------------------- */

 75: #include <petscts.h>

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

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

 96:   VecGetArray(X,&x);
 97:   VecGetArray(F,&f);
 98:   f[0] = (user->imex ? x[1] : 0);
 99:   f[1] = 0.0;
100:   VecRestoreArray(X,&x);
101:   VecRestoreArray(F,&f);
102:   return(0);
103: }

107: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
108: {
110:   User user = (User)ctx;
111:   PetscScalar *x,*xdot,*f;

114:   VecGetArray(X,&x);
115:   VecGetArray(Xdot,&xdot);
116:   VecGetArray(F,&f);
117:   f[0] = xdot[0] + (user->imex ? 0 : x[1]);
118:   f[1] = xdot[1] - user->mu*(1. - x[0]*x[0])*x[1] + x[0];
119:   VecRestoreArray(X,&x);
120:   VecRestoreArray(Xdot,&xdot);
121:   VecRestoreArray(F,&f);
122:   return(0);
123: }

127: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *flag,void *ctx)
128: {
130:   User user = (User)ctx;
131:   PetscReal mu = user->mu;
132:   PetscInt rowcol[] = {0,1};
133:   PetscScalar *x,J[2][2];

136:   VecGetArray(X,&x);
137:   J[0][0] = a;                    J[0][1] = (user->imex ? 0 : 1.);
138:   J[1][0] = 2.*mu*x[0]*x[1]+1.;   J[1][1] = a - mu*(1. - x[0]*x[0]);
139:   MatSetValues(*B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
140:   VecRestoreArray(X,&x);

142:   MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
143:   MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
144:   if (*A != *B) {
145:     MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
146:     MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
147:   }
148:   *flag = SAME_NONZERO_PATTERN;
149:   return(0);
150: }

154: static PetscErrorCode RegisterMyARK2(void)
155: {

159:   {
160:     const PetscReal
161:       A[3][3] = {{0,0,0},
162:                  {0.41421356237309504880,0,0},
163:                  {0.75,0.25,0}},
164:       At[3][3] = {{0,0,0},
165:                   {0.12132034355964257320,0.29289321881345247560,0},
166:                   {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}};
167:       TSARKIMEXRegister("myark2",2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,0,PETSC_NULL,PETSC_NULL);
168:   }
169:   return(0);
170: }

174: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
175: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
176: {
178:   const PetscScalar *x;
179:   PetscReal tfinal, dt;
180:   User user = (User)ctx;
181:   Vec interpolatedX;

184:   TSGetTimeStep(ts,&dt);
185:   TSGetDuration(ts,PETSC_NULL,&tfinal);

187:   while (user->next_output <= t && user->next_output <= tfinal) {
188:     VecDuplicate(X,&interpolatedX);
189:     TSInterpolate(ts,user->next_output,interpolatedX);
190:     VecGetArrayRead(interpolatedX,&x);
191:     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]));
192:     VecRestoreArrayRead(interpolatedX,&x);
193:     VecDestroy(&interpolatedX);
194:     user->next_output += 0.1;
195:   }
196:   return(0);
197: }

201: int main(int argc,char **argv)
202: {
203:   TS              ts;           /* nonlinear solver */
204:   Vec             x;            /* solution, residual vectors */
205:   Mat             A;            /* Jacobian matrix */
206:   PetscInt        steps;
207:   PetscReal       ftime=0.5;
208:   PetscBool       monitor = PETSC_FALSE;
209:   PetscScalar     *x_ptr;
210:   PetscMPIInt     size;
211:   struct _n_User  user;
212:   PetscErrorCode  ierr;

214:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215:      Initialize program
216:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217:   PetscInitialize(&argc,&argv,PETSC_NULL,help);

219:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
220:   if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");

222:   RegisterMyARK2();

224:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225:     Set runtime options
226:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
227:   user.mu = 1000;
228:   user.imex = PETSC_TRUE;
229:   user.next_output = 0.0;
230:   PetscOptionsGetReal(PETSC_NULL,"-mu",&user.mu,PETSC_NULL);
231:   PetscOptionsGetBool(PETSC_NULL,"-imex",&user.imex,PETSC_NULL);
232:   PetscOptionsGetBool(PETSC_NULL,"-monitor",&monitor,PETSC_NULL);

234:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235:     Create necessary matrix and vectors, solve same ODE on every process
236:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
237:   MatCreate(PETSC_COMM_WORLD,&A);
238:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
239:   MatSetFromOptions(A);

241:   MatGetVecs(A,&x,PETSC_NULL);

243:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244:      Create timestepping solver context
245:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
246:   TSCreate(PETSC_COMM_WORLD,&ts);
247:   TSSetType(ts,TSBEULER);
248:   TSSetRHSFunction(ts,PETSC_NULL,RHSFunction,&user);
249:   TSSetIFunction(ts,PETSC_NULL,IFunction,&user);
250:   TSSetIJacobian(ts,A,A,IJacobian,&user);
251:   TSSetDuration(ts,PETSC_DEFAULT,ftime);
252:   if (monitor) {
253:     TSMonitorSet(ts,Monitor,&user,PETSC_NULL);
254:   }

256:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257:      Set initial conditions
258:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
259:   VecGetArray(x,&x_ptr);
260:   x_ptr[0] = 2;   x_ptr[1] = 0.66666654321;
261:   VecRestoreArray(x,&x_ptr);
262:   TSSetInitialTimeStep(ts,0.0,.001);

264:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
265:      Set runtime options
266:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
267:   TSSetFromOptions(ts);

269:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270:      Solve nonlinear system
271:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
272:   TSSolve(ts,x,&ftime);
273:   TSGetTimeStepNumber(ts,&steps);
274:   PetscPrintf(PETSC_COMM_WORLD,"mu %G, steps %D, ftime %G\n",user.mu,steps,ftime);
275:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

277:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
278:      Free work space.  All PETSc objects should be destroyed when they
279:      are no longer needed.
280:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
281:   MatDestroy(&A);
282:   VecDestroy(&x);
283:   TSDestroy(&ts);

285:   PetscFinalize();
286:   return(0);
287: }