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