1 static char help[] = "Solves the van der Pol equation.\n\ 2 Input parameters include:\n"; 3 4 /* ------------------------------------------------------------------------ 5 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) = - 2/3 +10/(81*\mu) - 292/(2187*\mu^2), 11 and 12 \mu = 10^6 ( y'(0) ~ -0.6666665432100101). 13 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. 14 15 Notes: 16 This code demonstrates the TS solver interface to an ODE -- RHSFunction for explicit form and IFunction for implicit form. 17 18 ------------------------------------------------------------------------- */ 19 20 #include <petscts.h> 21 22 typedef struct _n_User *User; 23 struct _n_User { 24 PetscReal mu; 25 PetscReal next_output; 26 }; 27 28 /* 29 User-defined routines 30 */ 31 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, PetscCtx ctx) 32 { 33 User user = (User)ctx; 34 PetscScalar *f; 35 const PetscScalar *x; 36 37 PetscFunctionBeginUser; 38 PetscCall(VecGetArrayRead(X, &x)); 39 PetscCall(VecGetArray(F, &f)); 40 f[0] = x[1]; 41 f[1] = user->mu * (1. - x[0] * x[0]) * x[1] - x[0]; 42 PetscCall(VecRestoreArrayRead(X, &x)); 43 PetscCall(VecRestoreArray(F, &f)); 44 PetscFunctionReturn(PETSC_SUCCESS); 45 } 46 47 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx) 48 { 49 User user = (User)ctx; 50 const PetscScalar *x, *xdot; 51 PetscScalar *f; 52 53 PetscFunctionBeginUser; 54 PetscCall(VecGetArrayRead(X, &x)); 55 PetscCall(VecGetArrayRead(Xdot, &xdot)); 56 PetscCall(VecGetArray(F, &f)); 57 f[0] = xdot[0] - x[1]; 58 f[1] = xdot[1] - user->mu * ((1.0 - x[0] * x[0]) * x[1] - x[0]); 59 PetscCall(VecRestoreArrayRead(X, &x)); 60 PetscCall(VecRestoreArrayRead(Xdot, &xdot)); 61 PetscCall(VecRestoreArray(F, &f)); 62 PetscFunctionReturn(PETSC_SUCCESS); 63 } 64 65 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, PetscCtx ctx) 66 { 67 User user = (User)ctx; 68 PetscInt rowcol[] = {0, 1}; 69 const PetscScalar *x; 70 PetscScalar J[2][2]; 71 72 PetscFunctionBeginUser; 73 PetscCall(VecGetArrayRead(X, &x)); 74 J[0][0] = a; 75 J[0][1] = -1.0; 76 J[1][0] = user->mu * (2.0 * x[0] * x[1] + 1.0); 77 J[1][1] = a - user->mu * (1.0 - x[0] * x[0]); 78 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES)); 79 PetscCall(VecRestoreArrayRead(X, &x)); 80 81 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 82 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 83 if (A != B) { 84 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 85 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 86 } 87 PetscFunctionReturn(PETSC_SUCCESS); 88 } 89 90 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */ 91 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, PetscCtx ctx) 92 { 93 const PetscScalar *x; 94 PetscReal tfinal, dt; 95 User user = (User)ctx; 96 Vec interpolatedX; 97 98 PetscFunctionBeginUser; 99 PetscCall(TSGetTimeStep(ts, &dt)); 100 PetscCall(TSGetMaxTime(ts, &tfinal)); 101 102 while (user->next_output <= t && user->next_output <= tfinal) { 103 PetscCall(VecDuplicate(X, &interpolatedX)); 104 PetscCall(TSInterpolate(ts, user->next_output, interpolatedX)); 105 PetscCall(VecGetArrayRead(interpolatedX, &x)); 106 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]))); 107 PetscCall(VecRestoreArrayRead(interpolatedX, &x)); 108 PetscCall(VecDestroy(&interpolatedX)); 109 user->next_output += 0.1; 110 } 111 PetscFunctionReturn(PETSC_SUCCESS); 112 } 113 114 int main(int argc, char **argv) 115 { 116 TS ts; /* nonlinear solver */ 117 Vec x; /* solution, residual vectors */ 118 Mat A; /* Jacobian matrix */ 119 PetscInt steps; 120 PetscReal ftime = 0.5; 121 PetscBool monitor = PETSC_FALSE, implicitform = PETSC_TRUE; 122 PetscScalar *x_ptr; 123 PetscMPIInt size; 124 struct _n_User user; 125 126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 Initialize program 128 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 129 PetscFunctionBeginUser; 130 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 131 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 132 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 133 134 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 135 Set runtime options 136 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 137 user.next_output = 0.0; 138 user.mu = 1.0e3; 139 PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL)); 140 PetscCall(PetscOptionsGetBool(NULL, NULL, "-implicitform", &implicitform, NULL)); 141 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Physical parameters", NULL); 142 PetscCall(PetscOptionsReal("-mu", "Stiffness parameter", "<1.0e6>", user.mu, &user.mu, NULL)); 143 PetscOptionsEnd(); 144 145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 Create necessary matrix and vectors, solve same ODE on every process 147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 148 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 149 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2)); 150 PetscCall(MatSetFromOptions(A)); 151 PetscCall(MatSetUp(A)); 152 153 PetscCall(MatCreateVecs(A, &x, NULL)); 154 155 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 156 Create timestepping solver context 157 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 158 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 159 if (implicitform) { 160 PetscCall(TSSetIFunction(ts, NULL, IFunction, &user)); 161 PetscCall(TSSetIJacobian(ts, A, A, IJacobian, &user)); 162 PetscCall(TSSetType(ts, TSBEULER)); 163 } else { 164 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user)); 165 PetscCall(TSSetType(ts, TSRK)); 166 } 167 PetscCall(TSSetMaxTime(ts, ftime)); 168 PetscCall(TSSetTimeStep(ts, 0.001)); 169 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 170 if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL)); 171 172 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 173 Set initial conditions 174 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 175 PetscCall(VecGetArray(x, &x_ptr)); 176 x_ptr[0] = 2.0; 177 x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user.mu) - 292.0 / (2187.0 * user.mu * user.mu); 178 PetscCall(VecRestoreArray(x, &x_ptr)); 179 180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 181 Set runtime options 182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 183 PetscCall(TSSetFromOptions(ts)); 184 185 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 186 Solve nonlinear system 187 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 188 PetscCall(TSSolve(ts, x)); 189 PetscCall(TSGetSolveTime(ts, &ftime)); 190 PetscCall(TSGetStepNumber(ts, &steps)); 191 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "steps %" PetscInt_FMT ", ftime %g\n", steps, (double)ftime)); 192 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Free work space. All PETSc objects should be destroyed when they 196 are no longer needed. 197 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 198 PetscCall(MatDestroy(&A)); 199 PetscCall(VecDestroy(&x)); 200 PetscCall(TSDestroy(&ts)); 201 202 PetscCall(PetscFinalize()); 203 return 0; 204 } 205 206 /*TEST 207 208 test: 209 requires: !single 210 args: -mu 1e6 211 212 test: 213 requires: !single 214 suffix: 2 215 args: -implicitform false -ts_type rk -ts_rk_type 5dp -ts_adapt_type dsp 216 217 test: 218 requires: !single 219 suffix: 3 220 args: -implicitform false -ts_type rk -ts_rk_type 5dp -ts_adapt_type dsp -ts_adapt_dsp_filter H0312 221 222 TEST*/ 223