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) = -6.6e-01, 11 and 12 mu = 10^6. 13 This is a nonlinear equation. 14 15 This is a copy and modification of ex20.c to exactly match a test 16 problem that comes with the Radau5 integrator package. 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 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx) 29 { 30 User user = (User)ctx; 31 const PetscScalar *x, *xdot; 32 PetscScalar *f; 33 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 } 45 46 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, PetscCtx ctx) 47 { 48 User user = (User)ctx; 49 PetscInt rowcol[] = {0, 1}; 50 const PetscScalar *x; 51 PetscScalar J[2][2]; 52 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)); 61 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 } 70 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; 81 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!"); 89 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(); 98 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)); 106 107 PetscCall(MatCreateVecs(A, &x, NULL)); 108 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)); 116 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)); 128 129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130 Set runtime options 131 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 132 PetscCall(TSSetFromOptions(ts)); 133 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)); 142 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)); 150 151 PetscCall(PetscFinalize()); 152 return 0; 153 } 154 155 /*TEST 156 157 build: 158 requires: double !complex !defined(PETSC_USE_64BIT_INDICES) radau5 159 160 test: 161 args: -ts_monitor_solution -ts_type radau5 162 163 TEST*/ 164