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