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