1 static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and uses MOAB.\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 // PETSc includes: 18 #include <petscts.h> 19 #include <petscdmmoab.h> 20 21 typedef struct { 22 PetscScalar u, v; 23 } Field; 24 25 struct pUserCtx { 26 PetscReal A, B; /* Reaction coefficients */ 27 PetscReal alpha; /* Diffusion coefficient */ 28 Field leftbc; /* Dirichlet boundary conditions at left boundary */ 29 Field rightbc; /* Dirichlet boundary conditions at right boundary */ 30 PetscInt n, npts; /* Number of mesh points */ 31 PetscInt ntsteps; /* Number of time steps */ 32 PetscInt nvars; /* Number of variables in the equation system */ 33 PetscBool io; 34 }; 35 typedef pUserCtx *UserCtx; 36 37 PetscErrorCode Initialize_AppContext(UserCtx *puser) { 38 UserCtx user; 39 40 PetscFunctionBegin; 41 PetscCall(PetscNew(&user)); 42 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "ex35.cxx"); 43 { 44 user->nvars = 2; 45 user->A = 1; 46 user->B = 3; 47 user->alpha = 0.02; 48 user->leftbc.u = 1; 49 user->rightbc.u = 1; 50 user->leftbc.v = 3; 51 user->rightbc.v = 3; 52 user->n = 10; 53 user->ntsteps = 10000; 54 user->io = PETSC_FALSE; 55 PetscCall(PetscOptionsReal("-A", "Reaction rate", "ex35.cxx", user->A, &user->A, NULL)); 56 PetscCall(PetscOptionsReal("-B", "Reaction rate", "ex35.cxx", user->B, &user->B, NULL)); 57 PetscCall(PetscOptionsReal("-alpha", "Diffusion coefficient", "ex35.cxx", user->alpha, &user->alpha, NULL)); 58 PetscCall(PetscOptionsScalar("-uleft", "Dirichlet boundary condition", "ex35.cxx", user->leftbc.u, &user->leftbc.u, NULL)); 59 PetscCall(PetscOptionsScalar("-uright", "Dirichlet boundary condition", "ex35.cxx", user->rightbc.u, &user->rightbc.u, NULL)); 60 PetscCall(PetscOptionsScalar("-vleft", "Dirichlet boundary condition", "ex35.cxx", user->leftbc.v, &user->leftbc.v, NULL)); 61 PetscCall(PetscOptionsScalar("-vright", "Dirichlet boundary condition", "ex35.cxx", user->rightbc.v, &user->rightbc.v, NULL)); 62 PetscCall(PetscOptionsInt("-n", "Number of 1-D elements", "ex35.cxx", user->n, &user->n, NULL)); 63 PetscCall(PetscOptionsInt("-ndt", "Number of time steps", "ex35.cxx", user->ntsteps, &user->ntsteps, NULL)); 64 PetscCall(PetscOptionsBool("-io", "Write the mesh and solution output to a file.", "ex35.cxx", user->io, &user->io, NULL)); 65 user->npts = user->n + 1; 66 } 67 PetscOptionsEnd(); 68 69 *puser = user; 70 PetscFunctionReturn(0); 71 } 72 73 PetscErrorCode Destroy_AppContext(UserCtx *user) { 74 PetscFunctionBegin; 75 PetscCall(PetscFree(*user)); 76 PetscFunctionReturn(0); 77 } 78 79 static PetscErrorCode FormInitialSolution(TS, Vec, void *); 80 static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *); 81 static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *); 82 static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 83 84 /**************** 85 * * 86 * MAIN * 87 * * 88 ****************/ 89 int main(int argc, char **argv) { 90 TS ts; /* nonlinear solver */ 91 Vec X; /* solution, residual vectors */ 92 Mat J; /* Jacobian matrix */ 93 PetscInt steps, mx; 94 PetscReal hx, dt, ftime; 95 UserCtx user; /* user-defined work context */ 96 TSConvergedReason reason; 97 DM dm; 98 const char *fields[2] = {"U", "V"}; 99 100 PetscFunctionBeginUser; 101 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 102 103 /* Initialize the user context struct */ 104 PetscCall(Initialize_AppContext(&user)); 105 106 /* Fill in the user defined work context: */ 107 PetscCall(DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, PETSC_FALSE, NULL, user->n, 1, &dm)); 108 PetscCall(DMMoabSetFieldNames(dm, user->nvars, fields)); 109 PetscCall(DMMoabSetBlockSize(dm, user->nvars)); 110 PetscCall(DMSetFromOptions(dm)); 111 112 /* SetUp the data structures for DMMOAB */ 113 PetscCall(DMSetUp(dm)); 114 115 /* Create timestepping solver context */ 116 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 117 PetscCall(TSSetDM(ts, dm)); 118 PetscCall(TSSetType(ts, TSARKIMEX)); 119 PetscCall(TSSetEquationType(ts, TS_EQ_DAE_IMPLICIT_INDEX1)); 120 PetscCall(DMSetMatType(dm, MATBAIJ)); 121 PetscCall(DMCreateMatrix(dm, &J)); 122 123 PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, user)); 124 PetscCall(TSSetIFunction(ts, NULL, FormIFunction, user)); 125 PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, user)); 126 127 ftime = 10.0; 128 PetscCall(TSSetMaxSteps(ts, user->ntsteps)); 129 PetscCall(TSSetMaxTime(ts, ftime)); 130 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 131 132 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 133 Create the solution vector and set the initial conditions 134 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 135 PetscCall(DMCreateGlobalVector(dm, &X)); 136 137 PetscCall(FormInitialSolution(ts, X, user)); 138 PetscCall(TSSetSolution(ts, X)); 139 PetscCall(VecGetSize(X, &mx)); 140 hx = 1.0 / (PetscReal)(mx / 2 - 1); 141 dt = 0.4 * PetscSqr(hx) / user->alpha; /* Diffusive stability limit */ 142 PetscCall(TSSetTimeStep(ts, dt)); 143 144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 145 Set runtime options 146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 147 PetscCall(TSSetFromOptions(ts)); 148 149 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150 Solve nonlinear system 151 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 152 PetscCall(TSSolve(ts, X)); 153 PetscCall(TSGetSolveTime(ts, &ftime)); 154 PetscCall(TSGetStepNumber(ts, &steps)); 155 PetscCall(TSGetConvergedReason(ts, &reason)); 156 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps)); 157 158 if (user->io) { 159 /* Print the numerical solution to screen and then dump to file */ 160 PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD)); 161 162 /* Write out the solution along with the mesh */ 163 PetscCall(DMMoabSetGlobalFieldVector(dm, X)); 164 #ifdef MOAB_HAVE_HDF5 165 PetscCall(DMMoabOutput(dm, "ex35.h5m", "")); 166 #else 167 /* MOAB does not support true parallel writers that aren't HDF5 based 168 And so if you are using VTK as the output format in parallel, 169 the data could be jumbled due to the order in which the processors 170 write out their parts of the mesh and solution tags 171 */ 172 PetscCall(DMMoabOutput(dm, "ex35.vtk", "")); 173 #endif 174 } 175 176 /* Free work space. 177 Free all PETSc related resources: */ 178 PetscCall(MatDestroy(&J)); 179 PetscCall(VecDestroy(&X)); 180 PetscCall(TSDestroy(&ts)); 181 PetscCall(DMDestroy(&dm)); 182 183 /* Free all MOAB related resources: */ 184 PetscCall(Destroy_AppContext(&user)); 185 PetscCall(PetscFinalize()); 186 return 0; 187 } 188 189 /* 190 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot 191 */ 192 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr) { 193 UserCtx user = (UserCtx)ptr; 194 PetscInt dof; 195 PetscReal hx; 196 DM dm; 197 const moab::Range *vlocal; 198 PetscBool vonboundary; 199 200 PetscFunctionBegin; 201 PetscCall(TSGetDM(ts, &dm)); 202 203 /* get the essential MOAB mesh related quantities needed for FEM assembly */ 204 PetscCall(DMMoabGetLocalVertices(dm, &vlocal, NULL)); 205 206 /* compute local element sizes - structured grid */ 207 hx = 1.0 / user->n; 208 209 /* Compute function over the locally owned part of the grid 210 Assemble the operator by looping over edges and computing 211 contribution for each vertex dof */ 212 for (moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) { 213 const moab::EntityHandle vhandle = *iter; 214 215 PetscCall(DMMoabGetDofsBlocked(dm, 1, &vhandle, &dof)); 216 217 /* check if vertex is on the boundary */ 218 PetscCall(DMMoabIsEntityOnBoundary(dm, vhandle, &vonboundary)); 219 220 if (vonboundary) { 221 const PetscScalar bcvals[2][2] = { 222 {hx, 0 }, 223 {0, hx} 224 }; 225 PetscCall(MatSetValuesBlocked(Jpre, 1, &dof, 1, &dof, &bcvals[0][0], INSERT_VALUES)); 226 } else { 227 const PetscInt row = dof, col[] = {dof - 1, dof, dof + 1}; 228 const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx; 229 const PetscScalar vals[2][3][2] = { 230 {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}}, 231 {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}} 232 }; 233 PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES)); 234 } 235 } 236 237 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY)); 238 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY)); 239 if (J != Jpre) { 240 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 241 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 242 } 243 PetscFunctionReturn(0); 244 } 245 246 static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { 247 UserCtx user = (UserCtx)ptr; 248 DM dm; 249 PetscReal hx; 250 const Field *x; 251 Field *f; 252 PetscInt dof; 253 const moab::Range *ownedvtx; 254 255 PetscFunctionBegin; 256 hx = 1.0 / user->n; 257 PetscCall(TSGetDM(ts, &dm)); 258 259 /* Get pointers to vector data */ 260 PetscCall(VecSet(F, 0.0)); 261 262 PetscCall(DMMoabVecGetArrayRead(dm, X, &x)); 263 PetscCall(DMMoabVecGetArray(dm, F, &f)); 264 265 PetscCall(DMMoabGetLocalVertices(dm, &ownedvtx, NULL)); 266 267 /* Compute function over the locally owned part of the grid */ 268 for (moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) { 269 const moab::EntityHandle vhandle = *iter; 270 PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof)); 271 272 PetscScalar u = x[dof].u, v = x[dof].v; 273 f[dof].u = hx * (user->A + u * u * v - (user->B + 1) * u); 274 f[dof].v = hx * (user->B * u - u * u * v); 275 } 276 277 /* Restore vectors */ 278 PetscCall(DMMoabVecRestoreArrayRead(dm, X, &x)); 279 PetscCall(DMMoabVecRestoreArray(dm, F, &f)); 280 PetscFunctionReturn(0); 281 } 282 283 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx) { 284 UserCtx user = (UserCtx)ctx; 285 DM dm; 286 Field *x, *xdot, *f; 287 PetscReal hx; 288 Vec Xloc; 289 PetscInt i, bcindx; 290 PetscBool elem_on_boundary; 291 const moab::Range *vlocal; 292 293 PetscFunctionBegin; 294 hx = 1.0 / user->n; 295 PetscCall(TSGetDM(ts, &dm)); 296 297 /* get the essential MOAB mesh related quantities needed for FEM assembly */ 298 PetscCall(DMMoabGetLocalVertices(dm, &vlocal, NULL)); 299 300 /* reset the residual vector */ 301 PetscCall(VecSet(F, 0.0)); 302 303 PetscCall(DMGetLocalVector(dm, &Xloc)); 304 PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc)); 305 PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc)); 306 307 /* get the local representation of the arrays from Vectors */ 308 PetscCall(DMMoabVecGetArrayRead(dm, Xloc, &x)); 309 PetscCall(DMMoabVecGetArrayRead(dm, Xdot, &xdot)); 310 PetscCall(DMMoabVecGetArray(dm, F, &f)); 311 312 /* loop over local elements */ 313 for (moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) { 314 const moab::EntityHandle vhandle = *iter; 315 316 PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &i)); 317 318 /* check if vertex is on the boundary */ 319 PetscCall(DMMoabIsEntityOnBoundary(dm, vhandle, &elem_on_boundary)); 320 321 if (elem_on_boundary) { 322 PetscCall(DMMoabGetDofsBlocked(dm, 1, &vhandle, &bcindx)); 323 if (bcindx == 0) { /* Apply left BC */ 324 f[i].u = hx * (x[i].u - user->leftbc.u); 325 f[i].v = hx * (x[i].v - user->leftbc.v); 326 } else { /* Apply right BC */ 327 f[i].u = hx * (x[i].u - user->rightbc.u); 328 f[i].v = hx * (x[i].v - user->rightbc.v); 329 } 330 } else { 331 f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx; 332 f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx; 333 } 334 } 335 336 /* Restore data */ 337 PetscCall(DMMoabVecRestoreArrayRead(dm, Xloc, &x)); 338 PetscCall(DMMoabVecRestoreArrayRead(dm, Xdot, &xdot)); 339 PetscCall(DMMoabVecRestoreArray(dm, F, &f)); 340 PetscCall(DMRestoreLocalVector(dm, &Xloc)); 341 PetscFunctionReturn(0); 342 } 343 344 PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx) { 345 UserCtx user = (UserCtx)ctx; 346 PetscReal vpos[3]; 347 DM dm; 348 Field *x; 349 const moab::Range *vowned; 350 PetscInt dof; 351 moab::Range::iterator iter; 352 353 PetscFunctionBegin; 354 PetscCall(TSGetDM(ts, &dm)); 355 356 /* get the essential MOAB mesh related quantities needed for FEM assembly */ 357 PetscCall(DMMoabGetLocalVertices(dm, &vowned, NULL)); 358 359 PetscCall(VecSet(X, 0.0)); 360 361 /* Get pointers to vector data */ 362 PetscCall(DMMoabVecGetArray(dm, X, &x)); 363 364 /* Compute function over the locally owned part of the grid */ 365 for (moab::Range::iterator iter = vowned->begin(); iter != vowned->end(); iter++) { 366 const moab::EntityHandle vhandle = *iter; 367 PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof)); 368 369 /* compute the mid-point of the element and use a 1-point lumped quadrature */ 370 PetscCall(DMMoabGetVertexCoordinates(dm, 1, &vhandle, vpos)); 371 372 PetscReal xi = vpos[0]; 373 x[dof].u = user->leftbc.u * (1. - xi) + user->rightbc.u * xi + PetscSinReal(2. * PETSC_PI * xi); 374 x[dof].v = user->leftbc.v * (1. - xi) + user->rightbc.v * xi; 375 } 376 377 /* Restore vectors */ 378 PetscCall(DMMoabVecRestoreArray(dm, X, &x)); 379 PetscFunctionReturn(0); 380 } 381 382 /*TEST 383 384 build: 385 requires: moab 386 387 test: 388 args: -n 20 -ts_type rosw -ts_rosw_type 2p -ts_dt 5e-2 -ts_adapt_type none 389 390 test: 391 suffix: 2 392 nsize: 2 393 args: -n 50 -ts_type glee -ts_adapt_type none -ts_dt 0.1 -io 394 TODO: 395 396 TEST*/ 397