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