Actual source code: ex20.c


  2: static char help[] = "Solves the van der Pol equation.\n\
  3: Input parameters include:\n";

  5: /* ------------------------------------------------------------------------

  7:    This program solves the van der Pol DAE ODE equivalent
  8:        y' = z                 (1)
  9:        z' = \mu ((1-y^2)z-y)
 10:    on the domain 0 <= x <= 1, with the boundary conditions
 11:        y(0) = 2, y'(0) = - 2/3 +10/(81*\mu) - 292/(2187*\mu^2),
 12:    and
 13:        \mu = 10^6 ( y'(0) ~ -0.6666665432100101).
 14:    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.

 16:    Notes:
 17:    This code demonstrates the TS solver interface to an ODE -- RHSFunction for explicit form and IFunction for implicit form.

 19:   ------------------------------------------------------------------------- */

 21: #include <petscts.h>

 23: typedef struct _n_User *User;
 24: struct _n_User {
 25:   PetscReal mu;
 26:   PetscReal next_output;
 27: };

 29: /*
 30:    User-defined routines
 31: */
 32: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ctx)
 33: {
 34:   User               user = (User)ctx;
 35:   PetscScalar       *f;
 36:   const PetscScalar *x;

 38:   PetscFunctionBeginUser;
 39:   PetscCall(VecGetArrayRead(X, &x));
 40:   PetscCall(VecGetArray(F, &f));
 41:   f[0] = x[1];
 42:   f[1] = user->mu * (1. - x[0] * x[0]) * x[1] - x[0];
 43:   PetscCall(VecRestoreArrayRead(X, &x));
 44:   PetscCall(VecRestoreArray(F, &f));
 45:   PetscFunctionReturn(PETSC_SUCCESS);
 46: }

 48: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
 49: {
 50:   User               user = (User)ctx;
 51:   const PetscScalar *x, *xdot;
 52:   PetscScalar       *f;

 54:   PetscFunctionBeginUser;
 55:   PetscCall(VecGetArrayRead(X, &x));
 56:   PetscCall(VecGetArrayRead(Xdot, &xdot));
 57:   PetscCall(VecGetArray(F, &f));
 58:   f[0] = xdot[0] - x[1];
 59:   f[1] = xdot[1] - user->mu * ((1.0 - x[0] * x[0]) * x[1] - x[0]);
 60:   PetscCall(VecRestoreArrayRead(X, &x));
 61:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
 62:   PetscCall(VecRestoreArray(F, &f));
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
 67: {
 68:   User               user     = (User)ctx;
 69:   PetscInt           rowcol[] = {0, 1};
 70:   const PetscScalar *x;
 71:   PetscScalar        J[2][2];

 73:   PetscFunctionBeginUser;
 74:   PetscCall(VecGetArrayRead(X, &x));
 75:   J[0][0] = a;
 76:   J[0][1] = -1.0;
 77:   J[1][0] = user->mu * (2.0 * x[0] * x[1] + 1.0);
 78:   J[1][1] = a - user->mu * (1.0 - x[0] * x[0]);
 79:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
 80:   PetscCall(VecRestoreArrayRead(X, &x));

 82:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 83:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 84:   if (A != B) {
 85:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 86:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
 87:   }
 88:   PetscFunctionReturn(PETSC_SUCCESS);
 89: }

 91: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
 92: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, void *ctx)
 93: {
 94:   const PetscScalar *x;
 95:   PetscReal          tfinal, dt;
 96:   User               user = (User)ctx;
 97:   Vec                interpolatedX;

 99:   PetscFunctionBeginUser;
100:   PetscCall(TSGetTimeStep(ts, &dt));
101:   PetscCall(TSGetMaxTime(ts, &tfinal));

103:   while (user->next_output <= t && user->next_output <= tfinal) {
104:     PetscCall(VecDuplicate(X, &interpolatedX));
105:     PetscCall(TSInterpolate(ts, user->next_output, interpolatedX));
106:     PetscCall(VecGetArrayRead(interpolatedX, &x));
107:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%.1f] %" PetscInt_FMT " TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(x[0]), (double)PetscRealPart(x[1])));
108:     PetscCall(VecRestoreArrayRead(interpolatedX, &x));
109:     PetscCall(VecDestroy(&interpolatedX));
110:     user->next_output += 0.1;
111:   }
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }

115: int main(int argc, char **argv)
116: {
117:   TS             ts; /* nonlinear solver */
118:   Vec            x;  /* solution, residual vectors */
119:   Mat            A;  /* Jacobian matrix */
120:   PetscInt       steps;
121:   PetscReal      ftime   = 0.5;
122:   PetscBool      monitor = PETSC_FALSE, implicitform = PETSC_TRUE;
123:   PetscScalar   *x_ptr;
124:   PetscMPIInt    size;
125:   struct _n_User user;

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Initialize program
129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   PetscFunctionBeginUser;
131:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
132:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
133:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:     Set runtime options
137:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138:   user.next_output = 0.0;
139:   user.mu          = 1.0e3;
140:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
141:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-implicitform", &implicitform, NULL));
142:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Physical parameters", NULL);
143:   PetscCall(PetscOptionsReal("-mu", "Stiffness parameter", "<1.0e6>", user.mu, &user.mu, NULL));
144:   PetscOptionsEnd();

146:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147:     Create necessary matrix and vectors, solve same ODE on every process
148:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
150:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
151:   PetscCall(MatSetFromOptions(A));
152:   PetscCall(MatSetUp(A));

154:   PetscCall(MatCreateVecs(A, &x, NULL));

156:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157:      Create timestepping solver context
158:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
160:   if (implicitform) {
161:     PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
162:     PetscCall(TSSetIJacobian(ts, A, A, IJacobian, &user));
163:     PetscCall(TSSetType(ts, TSBEULER));
164:   } else {
165:     PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user));
166:     PetscCall(TSSetType(ts, TSRK));
167:   }
168:   PetscCall(TSSetMaxTime(ts, ftime));
169:   PetscCall(TSSetTimeStep(ts, 0.001));
170:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
171:   if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL));

173:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174:      Set initial conditions
175:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176:   PetscCall(VecGetArray(x, &x_ptr));
177:   x_ptr[0] = 2.0;
178:   x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user.mu) - 292.0 / (2187.0 * user.mu * user.mu);
179:   PetscCall(VecRestoreArray(x, &x_ptr));

181:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182:      Set runtime options
183:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184:   PetscCall(TSSetFromOptions(ts));

186:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187:      Solve nonlinear system
188:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189:   PetscCall(TSSolve(ts, x));
190:   PetscCall(TSGetSolveTime(ts, &ftime));
191:   PetscCall(TSGetStepNumber(ts, &steps));
192:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "steps %" PetscInt_FMT ", ftime %g\n", steps, (double)ftime));
193:   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));

195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196:      Free work space.  All PETSc objects should be destroyed when they
197:      are no longer needed.
198:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199:   PetscCall(MatDestroy(&A));
200:   PetscCall(VecDestroy(&x));
201:   PetscCall(TSDestroy(&ts));

203:   PetscCall(PetscFinalize());
204:   return (0);
205: }

207: /*TEST

209:     test:
210:       requires: !single
211:       args: -mu 1e6

213:     test:
214:       requires: !single
215:       suffix: 2
216:       args: -implicitform false -ts_type rk -ts_rk_type 5dp -ts_adapt_type dsp

218:     test:
219:       requires: !single
220:       suffix: 3
221:       args: -implicitform false -ts_type rk -ts_rk_type 5dp -ts_adapt_type dsp -ts_adapt_dsp_filter H0312

223: TEST*/