Actual source code: ex49.c
1: static char help[] = "Solves the van der Pol equation.\n\
2: Input parameters include:\n";
4: /* ------------------------------------------------------------------------
6: This program solves the van der Pol DAE ODE equivalent
7: y' = z (1)
8: z' = mu[(1-y^2)z-y]
9: on the domain 0 <= x <= 1, with the boundary conditions
10: y(0) = 2, y'(0) = -6.6e-01,
11: and
12: mu = 10^6.
13: This is a nonlinear equation.
15: This is a copy and modification of ex20.c to exactly match a test
16: problem that comes with the Radau5 integrator package.
18: ------------------------------------------------------------------------- */
20: #include <petscts.h>
22: typedef struct _n_User *User;
23: struct _n_User {
24: PetscReal mu;
25: PetscReal next_output;
26: };
28: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
29: {
30: User user = (User)ctx;
31: const PetscScalar *x, *xdot;
32: PetscScalar *f;
34: PetscFunctionBeginUser;
35: PetscCall(VecGetArrayRead(X, &x));
36: PetscCall(VecGetArrayRead(Xdot, &xdot));
37: PetscCall(VecGetArray(F, &f));
38: f[0] = xdot[0] - x[1];
39: f[1] = xdot[1] - user->mu * ((1.0 - x[0] * x[0]) * x[1] - x[0]);
40: PetscCall(VecRestoreArrayRead(X, &x));
41: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
42: PetscCall(VecRestoreArray(F, &f));
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
47: {
48: User user = (User)ctx;
49: PetscInt rowcol[] = {0, 1};
50: const PetscScalar *x;
51: PetscScalar J[2][2];
53: PetscFunctionBeginUser;
54: PetscCall(VecGetArrayRead(X, &x));
55: J[0][0] = a;
56: J[0][1] = -1.0;
57: J[1][0] = user->mu * (1.0 + 2.0 * x[0] * x[1]);
58: J[1][1] = a - user->mu * (1.0 - x[0] * x[0]);
59: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
60: PetscCall(VecRestoreArrayRead(X, &x));
62: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
63: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
64: if (A != B) {
65: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
66: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
67: }
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: int main(int argc, char **argv)
72: {
73: TS ts; /* nonlinear solver */
74: Vec x; /* solution, residual vectors */
75: Mat A; /* Jacobian matrix */
76: PetscInt steps;
77: PetscReal ftime = 2;
78: PetscScalar *x_ptr;
79: PetscMPIInt size;
80: struct _n_User user;
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Initialize program
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: PetscFunctionBeginUser;
86: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
87: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
88: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Set runtime options
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: user.next_output = 0.0;
94: user.mu = 1.0e6;
95: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Physical parameters", NULL);
96: PetscCall(PetscOptionsReal("-mu", "Stiffness parameter", "<1.0e6>", user.mu, &user.mu, NULL));
97: PetscOptionsEnd();
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create necessary matrix and vectors, solve same ODE on every process
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
103: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
104: PetscCall(MatSetFromOptions(A));
105: PetscCall(MatSetUp(A));
107: PetscCall(MatCreateVecs(A, &x, NULL));
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Create timestepping solver context
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
113: PetscCall(TSSetType(ts, TSBEULER));
114: PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
115: PetscCall(TSSetIJacobian(ts, A, A, IJacobian, &user));
117: PetscCall(TSSetMaxTime(ts, ftime));
118: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
119: PetscCall(TSSetTolerances(ts, 1.e-4, NULL, 1.e-4, NULL));
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Set initial conditions
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: PetscCall(VecGetArray(x, &x_ptr));
124: x_ptr[0] = 2.0;
125: x_ptr[1] = -6.6e-01;
126: PetscCall(VecRestoreArray(x, &x_ptr));
127: PetscCall(TSSetTimeStep(ts, .000001));
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Set runtime options
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: PetscCall(TSSetFromOptions(ts));
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Solve nonlinear system
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: PetscCall(TSSolve(ts, x));
138: PetscCall(TSGetSolveTime(ts, &ftime));
139: PetscCall(TSGetStepNumber(ts, &steps));
140: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "steps %" PetscInt_FMT ", ftime %g\n", steps, (double)ftime));
141: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Free work space. All PETSc objects should be destroyed when they
145: are no longer needed.
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: PetscCall(MatDestroy(&A));
148: PetscCall(VecDestroy(&x));
149: PetscCall(TSDestroy(&ts));
151: PetscCall(PetscFinalize());
152: return 0;
153: }
155: /*TEST
157: build:
158: requires: double !complex !defined(PETSC_USE_64BIT_INDICES) radau5
160: test:
161: args: -ts_monitor_solution -ts_type radau5
163: TEST*/