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