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