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