1 static const char help[] = "Call PetscInitialize multiple times.\n"; 2 /* 3 This example is based on the Brusselator tutorial of the same name, but tests multiple calls to PetscInitialize(). 4 This is a bad "convergence study" because it only compares min and max values of the solution rather than comparing 5 norms of the errors. For convergence studies, we recommend invoking PetscInitialize() only once and comparing norms 6 of errors (perhaps estimated using an accurate reference solution). 7 8 Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and multiple solves. 9 10 u_t - alpha u_xx = A + u^2 v - (B+1) u 11 v_t - alpha v_xx = B u - u^2 v 12 0 < x < 1; 13 A = 1, B = 3, alpha = 1/10 14 15 Initial conditions: 16 u(x,0) = 1 + sin(2 pi x) 17 v(x,0) = 3 18 19 Boundary conditions: 20 u(0,t) = u(1,t) = 1 21 v(0,t) = v(1,t) = 3 22 */ 23 24 #include <petscdm.h> 25 #include <petscdmda.h> 26 #include <petscts.h> 27 28 typedef struct { 29 PetscScalar u, v; 30 } Field; 31 32 typedef struct _User *User; 33 struct _User { 34 PetscReal A, B; /* Reaction coefficients */ 35 PetscReal alpha; /* Diffusion coefficient */ 36 PetscReal uleft, uright; /* Dirichlet boundary conditions */ 37 PetscReal vleft, vright; /* Dirichlet boundary conditions */ 38 }; 39 40 static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *); 41 static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *); 42 static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 43 static PetscErrorCode FormInitialSolution(TS, Vec, void *); 44 static int Brusselator(int, char **, PetscInt); 45 46 int main(int argc, char **argv) 47 { 48 PetscInt cycle; 49 PetscErrorCode ierr; 50 51 ierr = MPI_Init(&argc, &argv); 52 if (ierr) return ierr; 53 for (cycle = 0; cycle < 4; cycle++) { 54 ierr = Brusselator(argc, argv, cycle); 55 if (ierr) return 1; 56 } 57 ierr = MPI_Finalize(); 58 return ierr; 59 } 60 61 PetscErrorCode Brusselator(int argc, char **argv, PetscInt cycle) 62 { 63 TS ts; /* nonlinear solver */ 64 Vec X; /* solution, residual vectors */ 65 Mat J; /* Jacobian matrix */ 66 PetscInt steps, mx; 67 DM da; 68 PetscReal ftime, hx, dt, xmax, xmin; 69 struct _User user; /* user-defined work context */ 70 TSConvergedReason reason; 71 72 PetscFunctionBeginUser; 73 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 74 75 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 76 Create distributed array (DMDA) to manage parallel grid and vectors 77 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 78 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da)); 79 PetscCall(DMSetFromOptions(da)); 80 PetscCall(DMSetUp(da)); 81 82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 83 Extract global vectors from DMDA; 84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 85 PetscCall(DMCreateGlobalVector(da, &X)); 86 87 /* Initialize user application context */ 88 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", ""); 89 { 90 user.A = 1; 91 user.B = 3; 92 user.alpha = 0.1; 93 user.uleft = 1; 94 user.uright = 1; 95 user.vleft = 3; 96 user.vright = 3; 97 PetscCall(PetscOptionsReal("-A", "Reaction rate", "", user.A, &user.A, NULL)); 98 PetscCall(PetscOptionsReal("-B", "Reaction rate", "", user.B, &user.B, NULL)); 99 PetscCall(PetscOptionsReal("-alpha", "Diffusion coefficient", "", user.alpha, &user.alpha, NULL)); 100 PetscCall(PetscOptionsReal("-uleft", "Dirichlet boundary condition", "", user.uleft, &user.uleft, NULL)); 101 PetscCall(PetscOptionsReal("-uright", "Dirichlet boundary condition", "", user.uright, &user.uright, NULL)); 102 PetscCall(PetscOptionsReal("-vleft", "Dirichlet boundary condition", "", user.vleft, &user.vleft, NULL)); 103 PetscCall(PetscOptionsReal("-vright", "Dirichlet boundary condition", "", user.vright, &user.vright, NULL)); 104 } 105 PetscOptionsEnd(); 106 107 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 108 Create timestepping solver context 109 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 110 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 111 PetscCall(TSSetDM(ts, da)); 112 PetscCall(TSSetType(ts, TSARKIMEX)); 113 PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user)); 114 PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user)); 115 PetscCall(DMSetMatType(da, MATAIJ)); 116 PetscCall(DMCreateMatrix(da, &J)); 117 PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user)); 118 119 ftime = 1.0; 120 PetscCall(TSSetMaxTime(ts, ftime)); 121 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 122 123 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 124 Set initial conditions 125 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 126 PetscCall(FormInitialSolution(ts, X, &user)); 127 PetscCall(TSSetSolution(ts, X)); 128 PetscCall(VecGetSize(X, &mx)); 129 hx = 1.0 / (PetscReal)(mx / 2 - 1); 130 dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */ 131 dt *= PetscPowRealInt(0.2, cycle); /* Shrink the time step in convergence study. */ 132 PetscCall(TSSetTimeStep(ts, dt)); 133 PetscCall(TSSetTolerances(ts, 1e-3 * PetscPowRealInt(0.5, cycle), NULL, 1e-3 * PetscPowRealInt(0.5, cycle), NULL)); 134 135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 136 Set runtime options 137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 138 PetscCall(TSSetFromOptions(ts)); 139 140 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141 Solve nonlinear system 142 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 143 PetscCall(TSSolve(ts, X)); 144 PetscCall(TSGetSolveTime(ts, &ftime)); 145 PetscCall(TSGetStepNumber(ts, &steps)); 146 PetscCall(TSGetConvergedReason(ts, &reason)); 147 PetscCall(VecMin(X, NULL, &xmin)); 148 PetscCall(VecMax(X, NULL, &xmax)); 149 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after % 3" PetscInt_FMT " steps. Range [%6.4f,%6.4f]\n", TSConvergedReasons[reason], (double)ftime, steps, (double)xmin, (double)xmax)); 150 151 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 152 Free work space. 153 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 154 PetscCall(MatDestroy(&J)); 155 PetscCall(VecDestroy(&X)); 156 PetscCall(TSDestroy(&ts)); 157 PetscCall(DMDestroy(&da)); 158 PetscCall(PetscFinalize()); 159 return 0; 160 } 161 162 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr) 163 { 164 User user = (User)ptr; 165 DM da; 166 DMDALocalInfo info; 167 PetscInt i; 168 Field *x, *xdot, *f; 169 PetscReal hx; 170 Vec Xloc; 171 172 PetscFunctionBeginUser; 173 PetscCall(TSGetDM(ts, &da)); 174 PetscCall(DMDAGetLocalInfo(da, &info)); 175 hx = 1.0 / (PetscReal)(info.mx - 1); 176 177 /* 178 Scatter ghost points to local vector,using the 2-step process 179 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 180 By placing code between these two statements, computations can be 181 done while messages are in transition. 182 */ 183 PetscCall(DMGetLocalVector(da, &Xloc)); 184 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 185 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 186 187 /* Get pointers to vector data */ 188 PetscCall(DMDAVecGetArrayRead(da, Xloc, &x)); 189 PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot)); 190 PetscCall(DMDAVecGetArray(da, F, &f)); 191 192 /* Compute function over the locally owned part of the grid */ 193 for (i = info.xs; i < info.xs + info.xm; i++) { 194 if (i == 0) { 195 f[i].u = hx * (x[i].u - user->uleft); 196 f[i].v = hx * (x[i].v - user->vleft); 197 } else if (i == info.mx - 1) { 198 f[i].u = hx * (x[i].u - user->uright); 199 f[i].v = hx * (x[i].v - user->vright); 200 } else { 201 f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx; 202 f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx; 203 } 204 } 205 206 /* Restore vectors */ 207 PetscCall(DMDAVecRestoreArrayRead(da, Xloc, &x)); 208 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot)); 209 PetscCall(DMDAVecRestoreArray(da, F, &f)); 210 PetscCall(DMRestoreLocalVector(da, &Xloc)); 211 PetscFunctionReturn(0); 212 } 213 214 static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) 215 { 216 User user = (User)ptr; 217 DM da; 218 DMDALocalInfo info; 219 PetscInt i; 220 PetscReal hx; 221 Field *x, *f; 222 223 PetscFunctionBeginUser; 224 PetscCall(TSGetDM(ts, &da)); 225 PetscCall(DMDAGetLocalInfo(da, &info)); 226 hx = 1.0 / (PetscReal)(info.mx - 1); 227 228 /* Get pointers to vector data */ 229 PetscCall(DMDAVecGetArrayRead(da, X, &x)); 230 PetscCall(DMDAVecGetArray(da, F, &f)); 231 232 /* Compute function over the locally owned part of the grid */ 233 for (i = info.xs; i < info.xs + info.xm; i++) { 234 PetscScalar u = x[i].u, v = x[i].v; 235 f[i].u = hx * (user->A + u * u * v - (user->B + 1) * u); 236 f[i].v = hx * (user->B * u - u * u * v); 237 } 238 239 /* Restore vectors */ 240 PetscCall(DMDAVecRestoreArrayRead(da, X, &x)); 241 PetscCall(DMDAVecRestoreArray(da, F, &f)); 242 PetscFunctionReturn(0); 243 } 244 245 /* --------------------------------------------------------------------- */ 246 /* 247 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot 248 */ 249 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr) 250 { 251 User user = (User)ptr; 252 DMDALocalInfo info; 253 PetscInt i; 254 PetscReal hx; 255 DM da; 256 Field *x, *xdot; 257 258 PetscFunctionBeginUser; 259 PetscCall(TSGetDM(ts, &da)); 260 PetscCall(DMDAGetLocalInfo(da, &info)); 261 hx = 1.0 / (PetscReal)(info.mx - 1); 262 263 /* Get pointers to vector data */ 264 PetscCall(DMDAVecGetArrayRead(da, X, &x)); 265 PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot)); 266 267 /* Compute function over the locally owned part of the grid */ 268 for (i = info.xs; i < info.xs + info.xm; i++) { 269 if (i == 0 || i == info.mx - 1) { 270 const PetscInt row = i, col = i; 271 const PetscScalar vals[2][2] = { 272 {hx, 0 }, 273 {0, hx} 274 }; 275 PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 1, &col, &vals[0][0], INSERT_VALUES)); 276 } else { 277 const PetscInt row = i, col[] = {i - 1, i, i + 1}; 278 const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx; 279 const PetscScalar vals[2][3][2] = { 280 {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}}, 281 {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}} 282 }; 283 PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES)); 284 } 285 } 286 287 /* Restore vectors */ 288 PetscCall(DMDAVecRestoreArrayRead(da, X, &x)); 289 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot)); 290 291 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY)); 292 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY)); 293 if (J != Jpre) { 294 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 295 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 296 } 297 PetscFunctionReturn(0); 298 } 299 300 PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx) 301 { 302 User user = (User)ctx; 303 DM da; 304 PetscInt i; 305 DMDALocalInfo info; 306 Field *x; 307 PetscReal hx; 308 309 PetscFunctionBeginUser; 310 PetscCall(TSGetDM(ts, &da)); 311 PetscCall(DMDAGetLocalInfo(da, &info)); 312 hx = 1.0 / (PetscReal)(info.mx - 1); 313 314 /* Get pointers to vector data */ 315 PetscCall(DMDAVecGetArray(da, X, &x)); 316 317 /* Compute function over the locally owned part of the grid */ 318 for (i = info.xs; i < info.xs + info.xm; i++) { 319 PetscReal xi = i * hx; 320 x[i].u = user->uleft * (1. - xi) + user->uright * xi + PetscSinReal(2. * PETSC_PI * xi); 321 x[i].v = user->vleft * (1. - xi) + user->vright * xi; 322 } 323 PetscCall(DMDAVecRestoreArray(da, X, &x)); 324 PetscFunctionReturn(0); 325 } 326 327 /*TEST 328 329 test: 330 args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3 331 332 test: 333 suffix: 2 334 args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3 335 336 TEST*/ 337