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 PetscErrorCode Brusselator(int, char **, PetscInt); 45 46 int main(int argc, char **argv) 47 { 48 int ierr = MPI_Init(&argc, &argv); 49 if (ierr) return ierr; 50 for (PetscInt cycle = 0; cycle < 4; cycle++) { 51 if (Brusselator(argc, argv, cycle)) return 1; 52 } 53 return MPI_Finalize(); 54 } 55 56 PetscErrorCode Brusselator(int argc, char **argv, PetscInt cycle) 57 { 58 TS ts; /* nonlinear solver */ 59 Vec X; /* solution, residual vectors */ 60 Mat J; /* Jacobian matrix */ 61 PetscInt steps, mx; 62 DM da; 63 PetscReal ftime, hx, dt, xmax, xmin; 64 struct _User user; /* user-defined work context */ 65 TSConvergedReason reason; 66 67 PetscFunctionBeginUser; 68 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 69 70 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 71 Create distributed array (DMDA) to manage parallel grid and vectors 72 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 73 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da)); 74 PetscCall(DMSetFromOptions(da)); 75 PetscCall(DMSetUp(da)); 76 77 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 78 Extract global vectors from DMDA; 79 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 80 PetscCall(DMCreateGlobalVector(da, &X)); 81 82 /* Initialize user application context */ 83 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", ""); 84 { 85 user.A = 1; 86 user.B = 3; 87 user.alpha = 0.1; 88 user.uleft = 1; 89 user.uright = 1; 90 user.vleft = 3; 91 user.vright = 3; 92 PetscCall(PetscOptionsReal("-A", "Reaction rate", "", user.A, &user.A, NULL)); 93 PetscCall(PetscOptionsReal("-B", "Reaction rate", "", user.B, &user.B, NULL)); 94 PetscCall(PetscOptionsReal("-alpha", "Diffusion coefficient", "", user.alpha, &user.alpha, NULL)); 95 PetscCall(PetscOptionsReal("-uleft", "Dirichlet boundary condition", "", user.uleft, &user.uleft, NULL)); 96 PetscCall(PetscOptionsReal("-uright", "Dirichlet boundary condition", "", user.uright, &user.uright, NULL)); 97 PetscCall(PetscOptionsReal("-vleft", "Dirichlet boundary condition", "", user.vleft, &user.vleft, NULL)); 98 PetscCall(PetscOptionsReal("-vright", "Dirichlet boundary condition", "", user.vright, &user.vright, NULL)); 99 } 100 PetscOptionsEnd(); 101 102 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 103 Create timestepping solver context 104 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 105 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 106 PetscCall(TSSetDM(ts, da)); 107 PetscCall(TSSetType(ts, TSARKIMEX)); 108 PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user)); 109 PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user)); 110 PetscCall(DMSetMatType(da, MATAIJ)); 111 PetscCall(DMCreateMatrix(da, &J)); 112 PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user)); 113 114 ftime = 1.0; 115 PetscCall(TSSetMaxTime(ts, ftime)); 116 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 117 118 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119 Set initial conditions 120 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 121 PetscCall(FormInitialSolution(ts, X, &user)); 122 PetscCall(TSSetSolution(ts, X)); 123 PetscCall(VecGetSize(X, &mx)); 124 hx = 1.0 / (PetscReal)(mx / 2 - 1); 125 dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */ 126 dt *= PetscPowRealInt(0.2, cycle); /* Shrink the time step in convergence study. */ 127 PetscCall(TSSetTimeStep(ts, dt)); 128 PetscCall(TSSetTolerances(ts, 1e-3 * PetscPowRealInt(0.5, cycle), NULL, 1e-3 * PetscPowRealInt(0.5, cycle), NULL)); 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(TSGetConvergedReason(ts, &reason)); 142 PetscCall(VecMin(X, NULL, &xmin)); 143 PetscCall(VecMax(X, NULL, &xmax)); 144 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)); 145 146 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147 Free work space. 148 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 149 PetscCall(MatDestroy(&J)); 150 PetscCall(VecDestroy(&X)); 151 PetscCall(TSDestroy(&ts)); 152 PetscCall(DMDestroy(&da)); 153 PetscCall(PetscFinalize()); 154 return PETSC_SUCCESS; 155 } 156 157 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr) 158 { 159 User user = (User)ptr; 160 DM da; 161 DMDALocalInfo info; 162 PetscInt i; 163 Field *x, *xdot, *f; 164 PetscReal hx; 165 Vec Xloc; 166 167 PetscFunctionBeginUser; 168 PetscCall(TSGetDM(ts, &da)); 169 PetscCall(DMDAGetLocalInfo(da, &info)); 170 hx = 1.0 / (PetscReal)(info.mx - 1); 171 172 /* 173 Scatter ghost points to local vector,using the 2-step process 174 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 175 By placing code between these two statements, computations can be 176 done while messages are in transition. 177 */ 178 PetscCall(DMGetLocalVector(da, &Xloc)); 179 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); 180 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc)); 181 182 /* Get pointers to vector data */ 183 PetscCall(DMDAVecGetArrayRead(da, Xloc, &x)); 184 PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot)); 185 PetscCall(DMDAVecGetArray(da, F, &f)); 186 187 /* Compute function over the locally owned part of the grid */ 188 for (i = info.xs; i < info.xs + info.xm; i++) { 189 if (i == 0) { 190 f[i].u = hx * (x[i].u - user->uleft); 191 f[i].v = hx * (x[i].v - user->vleft); 192 } else if (i == info.mx - 1) { 193 f[i].u = hx * (x[i].u - user->uright); 194 f[i].v = hx * (x[i].v - user->vright); 195 } else { 196 f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx; 197 f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx; 198 } 199 } 200 201 /* Restore vectors */ 202 PetscCall(DMDAVecRestoreArrayRead(da, Xloc, &x)); 203 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot)); 204 PetscCall(DMDAVecRestoreArray(da, F, &f)); 205 PetscCall(DMRestoreLocalVector(da, &Xloc)); 206 PetscFunctionReturn(PETSC_SUCCESS); 207 } 208 209 static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) 210 { 211 User user = (User)ptr; 212 DM da; 213 DMDALocalInfo info; 214 PetscInt i; 215 PetscReal hx; 216 Field *x, *f; 217 218 PetscFunctionBeginUser; 219 PetscCall(TSGetDM(ts, &da)); 220 PetscCall(DMDAGetLocalInfo(da, &info)); 221 hx = 1.0 / (PetscReal)(info.mx - 1); 222 223 /* Get pointers to vector data */ 224 PetscCall(DMDAVecGetArrayRead(da, X, &x)); 225 PetscCall(DMDAVecGetArray(da, F, &f)); 226 227 /* Compute function over the locally owned part of the grid */ 228 for (i = info.xs; i < info.xs + info.xm; i++) { 229 PetscScalar u = x[i].u, v = x[i].v; 230 f[i].u = hx * (user->A + u * u * v - (user->B + 1) * u); 231 f[i].v = hx * (user->B * u - u * u * v); 232 } 233 234 /* Restore vectors */ 235 PetscCall(DMDAVecRestoreArrayRead(da, X, &x)); 236 PetscCall(DMDAVecRestoreArray(da, F, &f)); 237 PetscFunctionReturn(PETSC_SUCCESS); 238 } 239 240 /* --------------------------------------------------------------------- */ 241 /* 242 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot 243 */ 244 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr) 245 { 246 User user = (User)ptr; 247 DMDALocalInfo info; 248 PetscInt i; 249 PetscReal hx; 250 DM da; 251 Field *x, *xdot; 252 253 PetscFunctionBeginUser; 254 PetscCall(TSGetDM(ts, &da)); 255 PetscCall(DMDAGetLocalInfo(da, &info)); 256 hx = 1.0 / (PetscReal)(info.mx - 1); 257 258 /* Get pointers to vector data */ 259 PetscCall(DMDAVecGetArrayRead(da, X, &x)); 260 PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot)); 261 262 /* Compute function over the locally owned part of the grid */ 263 for (i = info.xs; i < info.xs + info.xm; i++) { 264 if (i == 0 || i == info.mx - 1) { 265 const PetscInt row = i, col = i; 266 const PetscScalar vals[2][2] = { 267 {hx, 0 }, 268 {0, hx} 269 }; 270 PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 1, &col, &vals[0][0], INSERT_VALUES)); 271 } else { 272 const PetscInt row = i, col[] = {i - 1, i, i + 1}; 273 const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx; 274 const PetscScalar vals[2][3][2] = { 275 {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}}, 276 {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}} 277 }; 278 PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES)); 279 } 280 } 281 282 /* Restore vectors */ 283 PetscCall(DMDAVecRestoreArrayRead(da, X, &x)); 284 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot)); 285 286 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY)); 287 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY)); 288 if (J != Jpre) { 289 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 290 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 291 } 292 PetscFunctionReturn(PETSC_SUCCESS); 293 } 294 295 PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx) 296 { 297 User user = (User)ctx; 298 DM da; 299 PetscInt i; 300 DMDALocalInfo info; 301 Field *x; 302 PetscReal hx; 303 304 PetscFunctionBeginUser; 305 PetscCall(TSGetDM(ts, &da)); 306 PetscCall(DMDAGetLocalInfo(da, &info)); 307 hx = 1.0 / (PetscReal)(info.mx - 1); 308 309 /* Get pointers to vector data */ 310 PetscCall(DMDAVecGetArray(da, X, &x)); 311 312 /* Compute function over the locally owned part of the grid */ 313 for (i = info.xs; i < info.xs + info.xm; i++) { 314 PetscReal xi = i * hx; 315 x[i].u = user->uleft * (1. - xi) + user->uright * xi + PetscSinReal(2. * PETSC_PI * xi); 316 x[i].v = user->vleft * (1. - xi) + user->vright * xi; 317 } 318 PetscCall(DMDAVecRestoreArray(da, X, &x)); 319 PetscFunctionReturn(PETSC_SUCCESS); 320 } 321 322 /*TEST 323 324 test: 325 args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3 326 327 test: 328 suffix: 2 329 args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3 330 331 TEST*/ 332