1 /* 2 * ex_vdp.c 3 * 4 * Created on: Jun 1, 2012 5 * Author: Hong Zhang 6 */ 7 static char help[] = "Solves the van der Pol equation. \n Input parameters include:\n"; 8 9 /* 10 * This program solves the van der Pol equation 11 * y' = z (1) 12 * z' = (((1-y^2)*z-y)/eps (2) 13 * on the domain 0<=x<=0.5, with the initial conditions 14 * y(0) = 2, 15 * z(0) = -2/3 + 10/81*eps - 292/2187*eps^2-1814/19683*eps^3 16 * IMEX schemes are applied to the split equation 17 * [y'] = [z] + [0 ] 18 * [z'] [0] [(((1-y^2)*z-y)/eps] 19 * 20 * F(x)= [z] 21 * [0] 22 * 23 * G(x)= [y'] - [0 ] 24 * [z'] [(((1-y^2)*z-y)/eps] 25 * 26 * JG(x) = G_x + a G_xdot 27 */ 28 29 #include <petscdmda.h> 30 #include <petscts.h> 31 32 typedef struct _User *User; 33 struct _User { 34 PetscReal mu; /*stiffness control coefficient: epsilon*/ 35 }; 36 37 static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *); 38 static PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *); 39 static PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 40 41 int main(int argc, char **argv) 42 { 43 TS ts; 44 Vec x; /* solution vector */ 45 Mat A; /* Jacobian */ 46 PetscInt steps, mx, eimex_rowcol[2], two; 47 PetscScalar *x_ptr; 48 PetscReal ftime, dt, norm; 49 Vec ref; 50 struct _User user; /* user-defined work context */ 51 PetscViewer viewer; 52 53 PetscFunctionBeginUser; 54 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 55 /* Initialize user application context */ 56 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "van der Pol options", ""); 57 user.mu = 1e0; 58 PetscCall(PetscOptionsReal("-eps", "Stiffness controller", "", user.mu, &user.mu, NULL)); 59 PetscOptionsEnd(); 60 61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 62 Set runtime options 63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 64 /* 65 PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 66 */ 67 68 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 69 Create necessary matrix and vectors, solve same ODE on every process 70 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 71 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 72 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2)); 73 PetscCall(MatSetFromOptions(A)); 74 PetscCall(MatSetUp(A)); 75 PetscCall(MatCreateVecs(A, &x, NULL)); 76 77 PetscCall(MatCreateVecs(A, &ref, NULL)); 78 PetscCall(VecGetArray(ref, &x_ptr)); 79 /* 80 * [0,1], mu=10^-3 81 */ 82 /* 83 x_ptr[0] = -1.8881254106283; 84 x_ptr[1] = 0.7359074233370;*/ 85 86 /* 87 * [0,0.5],mu=10^-3 88 */ 89 /* 90 x_ptr[0] = 1.596980778659137; 91 x_ptr[1] = -1.029103015879544; 92 */ 93 /* 94 * [0,0.5],mu=1 95 */ 96 x_ptr[0] = 1.619084329683235; 97 x_ptr[1] = -0.803530465176385; 98 99 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 100 Create timestepping solver context 101 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 102 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 103 PetscCall(TSSetType(ts, TSEIMEX)); 104 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user)); 105 PetscCall(TSSetIFunction(ts, NULL, IFunction, &user)); 106 PetscCall(TSSetIJacobian(ts, A, A, IJacobian, &user)); 107 108 dt = 0.00001; 109 ftime = 1.1; 110 PetscCall(TSSetTimeStep(ts, dt)); 111 PetscCall(TSSetMaxTime(ts, ftime)); 112 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114 Set initial conditions 115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 116 PetscCall(VecGetArray(x, &x_ptr)); 117 x_ptr[0] = 2.; 118 x_ptr[1] = -2. / 3. + 10. / 81. * (user.mu) - 292. / 2187. * (user.mu) * (user.mu) - 1814. / 19683. * (user.mu) * (user.mu) * (user.mu); 119 PetscCall(TSSetSolution(ts, x)); 120 PetscCall(VecGetSize(x, &mx)); 121 122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 123 Set runtime options 124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 125 PetscCall(TSSetFromOptions(ts)); 126 127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128 Solve nonlinear system 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 PetscCall(TSSolve(ts, x)); 131 PetscCall(TSGetTime(ts, &ftime)); 132 PetscCall(TSGetStepNumber(ts, &steps)); 133 134 PetscCall(VecAXPY(x, -1.0, ref)); 135 PetscCall(VecNorm(x, NORM_2, &norm)); 136 PetscCall(TSGetTimeStep(ts, &dt)); 137 138 eimex_rowcol[0] = 0; 139 eimex_rowcol[1] = 0; 140 two = 2; 141 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-ts_eimex_row_col", eimex_rowcol, &two, NULL)); 142 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "order %11s %18s %37s\n", "dt", "norm", "final solution components 0 and 1")); 143 PetscCall(VecGetArray(x, &x_ptr)); 144 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(%" PetscInt_FMT ",%" PetscInt_FMT ") %10.8f %18.15f %18.15f %18.15f\n", eimex_rowcol[0], eimex_rowcol[1], (double)dt, (double)norm, (double)PetscRealPart(x_ptr[0]), (double)PetscRealPart(x_ptr[1]))); 145 PetscCall(VecRestoreArray(x, &x_ptr)); 146 147 /* Write line in convergence log */ 148 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer)); 149 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII)); 150 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_APPEND)); 151 PetscCall(PetscViewerFileSetName(viewer, "eimex_nonstiff_vdp.txt")); 152 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %10.8f %18.15f\n", eimex_rowcol[0], eimex_rowcol[1], (double)dt, (double)norm)); 153 PetscCall(PetscViewerDestroy(&viewer)); 154 155 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 156 Free work space. 157 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 158 PetscCall(MatDestroy(&A)); 159 PetscCall(VecDestroy(&x)); 160 PetscCall(VecDestroy(&ref)); 161 PetscCall(TSDestroy(&ts)); 162 PetscCall(PetscFinalize()); 163 return 0; 164 } 165 166 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) 167 { 168 PetscScalar *f; 169 const PetscScalar *x; 170 171 PetscFunctionBegin; 172 PetscCall(VecGetArrayRead(X, &x)); 173 PetscCall(VecGetArray(F, &f)); 174 f[0] = x[1]; 175 f[1] = 0.0; 176 PetscCall(VecRestoreArrayRead(X, &x)); 177 PetscCall(VecRestoreArray(F, &f)); 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180 181 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr) 182 { 183 User user = (User)ptr; 184 PetscScalar *f; 185 const PetscScalar *x, *xdot; 186 187 PetscFunctionBegin; 188 PetscCall(VecGetArrayRead(X, &x)); 189 PetscCall(VecGetArrayRead(Xdot, &xdot)); 190 PetscCall(VecGetArray(F, &f)); 191 f[0] = xdot[0]; 192 f[1] = xdot[1] - ((1. - x[0] * x[0]) * x[1] - x[0]) / user->mu; 193 PetscCall(VecRestoreArrayRead(X, &x)); 194 PetscCall(VecRestoreArrayRead(Xdot, &xdot)); 195 PetscCall(VecRestoreArray(F, &f)); 196 PetscFunctionReturn(PETSC_SUCCESS); 197 } 198 199 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ptr) 200 { 201 User user = (User)ptr; 202 PetscReal mu = user->mu; 203 PetscInt rowcol[] = {0, 1}; 204 PetscScalar J[2][2]; 205 const PetscScalar *x; 206 207 PetscFunctionBegin; 208 PetscCall(VecGetArrayRead(X, &x)); 209 J[0][0] = a; 210 J[0][1] = 0; 211 J[1][0] = (2. * x[0] * x[1] + 1.) / mu; 212 J[1][1] = a - (1. - x[0] * x[0]) / mu; 213 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES)); 214 PetscCall(VecRestoreArrayRead(X, &x)); 215 216 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 217 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 218 if (A != B) { 219 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 220 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 221 } 222 PetscFunctionReturn(PETSC_SUCCESS); 223 } 224 225 /*TEST 226 227 test: 228 args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_row_col 3,3 -ts_monitor_lg_solution 229 requires: x 230 231 test: 232 suffix: adapt 233 args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt 0.01 -ts_max_time 10 -ts_eimex_order_adapt -ts_eimex_max_rows 7 -ts_monitor_lg_solution 234 requires: x 235 236 test: 237 suffix: loop 238 args: -ts_type eimex -ts_adapt_type none -pc_type lu -ts_dt {{0.005 0.001 0.0005}separate output} -ts_max_steps {{100 500 1000}separate output} -ts_eimex_row_col {{1,1 2,1 3,1 2,2 3,2 3,3}separate output} 239 requires: x 240 241 TEST*/ 242