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