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 splitted 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 TS ts; 43 Vec x; /* solution vector */ 44 Mat A; /* Jacobian */ 45 PetscInt steps, mx, eimex_rowcol[2], two; 46 PetscScalar *x_ptr; 47 PetscReal ftime, dt, norm; 48 Vec ref; 49 struct _User user; /* user-defined work context */ 50 PetscViewer viewer; 51 52 PetscFunctionBeginUser; 53 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 54 /* Initialize user application context */ 55 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "van der Pol options", ""); 56 user.mu = 1e0; 57 PetscCall(PetscOptionsReal("-eps", "Stiffness controller", "", user.mu, &user.mu, NULL)); 58 PetscOptionsEnd(); 59 60 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 61 Set runtime options 62 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 63 /* 64 PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 65 */ 66 67 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 68 Create necessary matrix and vectors, solve same ODE on every process 69 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 70 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 71 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2)); 72 PetscCall(MatSetFromOptions(A)); 73 PetscCall(MatSetUp(A)); 74 PetscCall(MatCreateVecs(A, &x, NULL)); 75 76 PetscCall(MatCreateVecs(A, &ref, NULL)); 77 PetscCall(VecGetArray(ref, &x_ptr)); 78 /* 79 * [0,1], mu=10^-3 80 */ 81 /* 82 x_ptr[0] = -1.8881254106283; 83 x_ptr[1] = 0.7359074233370;*/ 84 85 /* 86 * [0,0.5],mu=10^-3 87 */ 88 /* 89 x_ptr[0] = 1.596980778659137; 90 x_ptr[1] = -1.029103015879544; 91 */ 92 /* 93 * [0,0.5],mu=1 94 */ 95 x_ptr[0] = 1.619084329683235; 96 x_ptr[1] = -0.803530465176385; 97 98 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 99 Create timestepping solver context 100 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 101 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 102 PetscCall(TSSetType(ts, TSEIMEX)); 103 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user)); 104 PetscCall(TSSetIFunction(ts, NULL, IFunction, &user)); 105 PetscCall(TSSetIJacobian(ts, A, A, IJacobian, &user)); 106 107 dt = 0.00001; 108 ftime = 1.1; 109 PetscCall(TSSetTimeStep(ts, dt)); 110 PetscCall(TSSetMaxTime(ts, ftime)); 111 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 113 Set initial conditions 114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 115 PetscCall(VecGetArray(x, &x_ptr)); 116 x_ptr[0] = 2.; 117 x_ptr[1] = -2. / 3. + 10. / 81. * (user.mu) - 292. / 2187. * (user.mu) * (user.mu) - 1814. / 19683. * (user.mu) * (user.mu) * (user.mu); 118 PetscCall(TSSetSolution(ts, x)); 119 PetscCall(VecGetSize(x, &mx)); 120 121 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 122 Set runtime options 123 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 124 PetscCall(TSSetFromOptions(ts)); 125 126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 Solve nonlinear system 128 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 129 PetscCall(TSSolve(ts, x)); 130 PetscCall(TSGetTime(ts, &ftime)); 131 PetscCall(TSGetStepNumber(ts, &steps)); 132 133 PetscCall(VecAXPY(x, -1.0, ref)); 134 PetscCall(VecNorm(x, NORM_2, &norm)); 135 PetscCall(TSGetTimeStep(ts, &dt)); 136 137 eimex_rowcol[0] = 0; 138 eimex_rowcol[1] = 0; 139 two = 2; 140 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-ts_eimex_row_col", eimex_rowcol, &two, NULL)); 141 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "order %11s %18s %37s\n", "dt", "norm", "final solution components 0 and 1")); 142 PetscCall(VecGetArray(x, &x_ptr)); 143 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]))); 144 PetscCall(VecRestoreArray(x, &x_ptr)); 145 146 /* Write line in convergence log */ 147 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer)); 148 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII)); 149 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_APPEND)); 150 PetscCall(PetscViewerFileSetName(viewer, "eimex_nonstiff_vdp.txt")); 151 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %10.8f %18.15f\n", eimex_rowcol[0], eimex_rowcol[1], (double)dt, (double)norm)); 152 PetscCall(PetscViewerDestroy(&viewer)); 153 154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 155 Free work space. 156 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 157 PetscCall(MatDestroy(&A)); 158 PetscCall(VecDestroy(&x)); 159 PetscCall(VecDestroy(&ref)); 160 PetscCall(TSDestroy(&ts)); 161 PetscCall(PetscFinalize()); 162 return 0; 163 } 164 165 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { 166 PetscScalar *f; 167 const PetscScalar *x; 168 169 PetscFunctionBegin; 170 PetscCall(VecGetArrayRead(X, &x)); 171 PetscCall(VecGetArray(F, &f)); 172 f[0] = x[1]; 173 f[1] = 0.0; 174 PetscCall(VecRestoreArrayRead(X, &x)); 175 PetscCall(VecRestoreArray(F, &f)); 176 PetscFunctionReturn(0); 177 } 178 179 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr) { 180 User user = (User)ptr; 181 PetscScalar *f; 182 const PetscScalar *x, *xdot; 183 184 PetscFunctionBegin; 185 PetscCall(VecGetArrayRead(X, &x)); 186 PetscCall(VecGetArrayRead(Xdot, &xdot)); 187 PetscCall(VecGetArray(F, &f)); 188 f[0] = xdot[0]; 189 f[1] = xdot[1] - ((1. - x[0] * x[0]) * x[1] - x[0]) / user->mu; 190 PetscCall(VecRestoreArrayRead(X, &x)); 191 PetscCall(VecRestoreArrayRead(Xdot, &xdot)); 192 PetscCall(VecRestoreArray(F, &f)); 193 PetscFunctionReturn(0); 194 } 195 196 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ptr) { 197 User user = (User)ptr; 198 PetscReal mu = user->mu; 199 PetscInt rowcol[] = {0, 1}; 200 PetscScalar J[2][2]; 201 const PetscScalar *x; 202 203 PetscFunctionBegin; 204 PetscCall(VecGetArrayRead(X, &x)); 205 J[0][0] = a; 206 J[0][1] = 0; 207 J[1][0] = (2. * x[0] * x[1] + 1.) / mu; 208 J[1][1] = a - (1. - x[0] * x[0]) / mu; 209 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES)); 210 PetscCall(VecRestoreArrayRead(X, &x)); 211 212 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 213 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 214 if (A != B) { 215 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 216 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 217 } 218 PetscFunctionReturn(0); 219 } 220 221 /*TEST 222 223 test: 224 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 225 requires: x 226 227 test: 228 suffix: adapt 229 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 230 requires: x 231 232 test: 233 suffix: loop 234 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} 235 requires: x 236 237 TEST*/ 238