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 CHKERRQ(PetscNew(&user)); 44 45 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Advection-reaction options","ex35.cxx");CHKERRQ(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 CHKERRQ(PetscOptionsReal("-A","Reaction rate","ex35.cxx",user->A,&user->A,NULL)); 59 CHKERRQ(PetscOptionsReal("-B","Reaction rate","ex35.cxx",user->B,&user->B,NULL)); 60 CHKERRQ(PetscOptionsReal("-alpha","Diffusion coefficient","ex35.cxx",user->alpha,&user->alpha,NULL)); 61 CHKERRQ(PetscOptionsScalar("-uleft","Dirichlet boundary condition","ex35.cxx",user->leftbc.u,&user->leftbc.u,NULL)); 62 CHKERRQ(PetscOptionsScalar("-uright","Dirichlet boundary condition","ex35.cxx",user->rightbc.u,&user->rightbc.u,NULL)); 63 CHKERRQ(PetscOptionsScalar("-vleft","Dirichlet boundary condition","ex35.cxx",user->leftbc.v,&user->leftbc.v,NULL)); 64 CHKERRQ(PetscOptionsScalar("-vright","Dirichlet boundary condition","ex35.cxx",user->rightbc.v,&user->rightbc.v,NULL)); 65 CHKERRQ(PetscOptionsInt("-n","Number of 1-D elements","ex35.cxx",user->n,&user->n,NULL)); 66 CHKERRQ(PetscOptionsInt("-ndt","Number of time steps","ex35.cxx",user->ntsteps,&user->ntsteps,NULL)); 67 CHKERRQ(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();CHKERRQ(ierr); 71 72 *puser = user; 73 PetscFunctionReturn(0); 74 } 75 76 PetscErrorCode Destroy_AppContext(UserCtx *user) 77 { 78 PetscFunctionBegin; 79 CHKERRQ(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 CHKERRQ(PetscInitialize(&argc,&argv,(char *)0,help)); 106 107 /* Initialize the user context struct */ 108 CHKERRQ(Initialize_AppContext(&user)); 109 110 /* Fill in the user defined work context: */ 111 CHKERRQ(DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, PETSC_FALSE, NULL, user->n, 1, &dm)); 112 CHKERRQ(DMMoabSetFieldNames(dm, user->nvars, fields)); 113 CHKERRQ(DMMoabSetBlockSize(dm, user->nvars)); 114 CHKERRQ(DMSetFromOptions(dm)); 115 116 /* SetUp the data structures for DMMOAB */ 117 CHKERRQ(DMSetUp(dm)); 118 119 /* Create timestepping solver context */ 120 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 121 CHKERRQ(TSSetDM(ts, dm)); 122 CHKERRQ(TSSetType(ts,TSARKIMEX)); 123 CHKERRQ(TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1)); 124 CHKERRQ(DMSetMatType(dm,MATBAIJ)); 125 CHKERRQ(DMCreateMatrix(dm,&J)); 126 127 CHKERRQ(TSSetRHSFunction(ts,NULL,FormRHSFunction,user)); 128 CHKERRQ(TSSetIFunction(ts,NULL,FormIFunction,user)); 129 CHKERRQ(TSSetIJacobian(ts,J,J,FormIJacobian,user)); 130 131 ftime = 10.0; 132 CHKERRQ(TSSetMaxSteps(ts,user->ntsteps)); 133 CHKERRQ(TSSetMaxTime(ts,ftime)); 134 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 135 136 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 137 Create the solution vector and set the initial conditions 138 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 139 CHKERRQ(DMCreateGlobalVector(dm, &X)); 140 141 CHKERRQ(FormInitialSolution(ts,X,user)); 142 CHKERRQ(TSSetSolution(ts,X)); 143 CHKERRQ(VecGetSize(X,&mx)); 144 hx = 1.0/(PetscReal)(mx/2-1); 145 dt = 0.4 * PetscSqr(hx) / user->alpha; /* Diffusive stability limit */ 146 CHKERRQ(TSSetTimeStep(ts,dt)); 147 148 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 149 Set runtime options 150 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 151 CHKERRQ(TSSetFromOptions(ts)); 152 153 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 154 Solve nonlinear system 155 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 156 CHKERRQ(TSSolve(ts,X)); 157 CHKERRQ(TSGetSolveTime(ts,&ftime)); 158 CHKERRQ(TSGetStepNumber(ts,&steps)); 159 CHKERRQ(TSGetConvergedReason(ts,&reason)); 160 CHKERRQ(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 CHKERRQ(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); 165 166 /* Write out the solution along with the mesh */ 167 CHKERRQ(DMMoabSetGlobalFieldVector(dm, X)); 168 #ifdef MOAB_HAVE_HDF5 169 CHKERRQ(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 CHKERRQ(DMMoabOutput(dm, "ex35.vtk", "")); 177 #endif 178 } 179 180 /* Free work space. 181 Free all PETSc related resources: */ 182 CHKERRQ(MatDestroy(&J)); 183 CHKERRQ(VecDestroy(&X)); 184 CHKERRQ(TSDestroy(&ts)); 185 CHKERRQ(DMDestroy(&dm)); 186 187 /* Free all MOAB related resources: */ 188 CHKERRQ(Destroy_AppContext(&user)); 189 CHKERRQ(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 CHKERRQ(TSGetDM(ts, &dm)); 207 208 /* get the essential MOAB mesh related quantities needed for FEM assembly */ 209 CHKERRQ(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 CHKERRQ(DMMoabGetDofsBlocked(dm, 1, &vhandle, &dof)); 221 222 /* check if vertex is on the boundary */ 223 CHKERRQ(DMMoabIsEntityOnBoundary(dm,vhandle,&vonboundary)); 224 225 if (vonboundary) { 226 const PetscScalar bcvals[2][2] = {{hx,0},{0,hx}}; 227 CHKERRQ(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 CHKERRQ(MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES)); 235 } 236 } 237 238 CHKERRQ(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY)); 239 CHKERRQ(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY)); 240 if (J != Jpre) { 241 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 242 CHKERRQ(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 CHKERRQ(TSGetDM(ts,&dm)); 260 261 /* Get pointers to vector data */ 262 CHKERRQ(VecSet(F,0.0)); 263 264 CHKERRQ(DMMoabVecGetArrayRead(dm, X, &x)); 265 CHKERRQ(DMMoabVecGetArray(dm, F, &f)); 266 267 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMMoabVecRestoreArrayRead(dm, X, &x)); 281 CHKERRQ(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 CHKERRQ(TSGetDM(ts, &dm)); 299 300 /* get the essential MOAB mesh related quantities needed for FEM assembly */ 301 CHKERRQ(DMMoabGetLocalVertices(dm, &vlocal, NULL)); 302 303 /* reset the residual vector */ 304 CHKERRQ(VecSet(F,0.0)); 305 306 CHKERRQ(DMGetLocalVector(dm,&Xloc)); 307 CHKERRQ(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc)); 308 CHKERRQ(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc)); 309 310 /* get the local representation of the arrays from Vectors */ 311 CHKERRQ(DMMoabVecGetArrayRead(dm, Xloc, &x)); 312 CHKERRQ(DMMoabVecGetArrayRead(dm, Xdot, &xdot)); 313 CHKERRQ(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 CHKERRQ(DMMoabGetDofsBlockedLocal(dm,1,&vhandle,&i)); 320 321 /* check if vertex is on the boundary */ 322 CHKERRQ(DMMoabIsEntityOnBoundary(dm,vhandle,&elem_on_boundary)); 323 324 if (elem_on_boundary) { 325 CHKERRQ(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 CHKERRQ(DMMoabVecRestoreArrayRead(dm, Xloc, &x)); 342 CHKERRQ(DMMoabVecRestoreArrayRead(dm, Xdot, &xdot)); 343 CHKERRQ(DMMoabVecRestoreArray(dm, F, &f)); 344 CHKERRQ(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 CHKERRQ(TSGetDM(ts, &dm)); 360 361 /* get the essential MOAB mesh related quantities needed for FEM assembly */ 362 CHKERRQ(DMMoabGetLocalVertices(dm, &vowned, NULL)); 363 364 CHKERRQ(VecSet(X, 0.0)); 365 366 /* Get pointers to vector data */ 367 CHKERRQ(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 CHKERRQ(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof)); 373 374 /* compute the mid-point of the element and use a 1-point lumped quadrature */ 375 CHKERRQ(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 CHKERRQ(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