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 /* 199 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot 200 */ 201 PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr) 202 { 203 UserCtx user = (UserCtx)ptr; 204 PetscErrorCode ierr; 205 PetscInt dof; 206 PetscReal hx; 207 DM dm; 208 const moab::Range *vlocal; 209 PetscBool vonboundary; 210 211 PetscFunctionBegin; 212 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 213 214 /* get the essential MOAB mesh related quantities needed for FEM assembly */ 215 ierr = DMMoabGetLocalVertices(dm, &vlocal, NULL);CHKERRQ(ierr); 216 217 /* compute local element sizes - structured grid */ 218 hx = 1.0/user->n; 219 220 /* Compute function over the locally owned part of the grid 221 Assemble the operator by looping over edges and computing 222 contribution for each vertex dof */ 223 for(moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) { 224 const moab::EntityHandle vhandle = *iter; 225 226 ierr = DMMoabGetDofsBlocked(dm, 1, &vhandle, &dof);CHKERRQ(ierr); 227 228 /* check if vertex is on the boundary */ 229 ierr = DMMoabIsEntityOnBoundary(dm,vhandle,&vonboundary);CHKERRQ(ierr); 230 231 if (vonboundary) { 232 const PetscScalar bcvals[2][2] = {{hx,0},{0,hx}}; 233 ierr = MatSetValuesBlocked(Jpre,1,&dof,1,&dof,&bcvals[0][0],INSERT_VALUES);CHKERRQ(ierr); 234 } 235 else { 236 const PetscInt row = dof,col[] = {dof-1,dof,dof+1}; 237 const PetscScalar dxxL = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx; 238 const PetscScalar vals[2][3][2] = {{{dxxL,0},{a *hx+dxx0,0},{dxxR,0}}, 239 {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}}; 240 ierr = MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);CHKERRQ(ierr); 241 } 242 } 243 244 ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 245 ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 246 if (J != Jpre) { 247 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 248 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 249 } 250 PetscFunctionReturn(0); 251 } 252 253 254 static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 255 { 256 UserCtx user = (UserCtx)ptr; 257 DM dm; 258 PetscReal hx; 259 const Field *x; 260 Field *f; 261 PetscInt dof; 262 const moab::Range *ownedvtx; 263 PetscErrorCode ierr; 264 265 PetscFunctionBegin; 266 hx = 1.0/user->n; 267 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 268 269 /* Get pointers to vector data */ 270 ierr = VecSet(F,0.0);CHKERRQ(ierr); 271 272 ierr = DMMoabVecGetArrayRead(dm, X, &x);CHKERRQ(ierr); 273 ierr = DMMoabVecGetArray(dm, F, &f);CHKERRQ(ierr); 274 275 ierr = DMMoabGetLocalVertices(dm, &ownedvtx, NULL);CHKERRQ(ierr); 276 277 /* Compute function over the locally owned part of the grid */ 278 for(moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) { 279 const moab::EntityHandle vhandle = *iter; 280 ierr = DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof);CHKERRQ(ierr); 281 282 PetscScalar u = x[dof].u,v = x[dof].v; 283 f[dof].u = hx*(user->A + u*u*v - (user->B+1)*u); 284 f[dof].v = hx*(user->B*u - u*u*v); 285 } 286 287 /* Restore vectors */ 288 ierr = DMMoabVecRestoreArrayRead(dm, X, &x);CHKERRQ(ierr); 289 ierr = DMMoabVecRestoreArray(dm, F, &f);CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 294 static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 295 { 296 UserCtx user = (UserCtx)ctx; 297 DM dm; 298 Field *x,*xdot,*f; 299 PetscReal hx; 300 Vec Xloc; 301 PetscErrorCode ierr; 302 PetscInt i,bcindx; 303 PetscBool elem_on_boundary; 304 const moab::Range *vlocal; 305 306 PetscFunctionBegin; 307 hx = 1.0/user->n; 308 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 309 310 /* get the essential MOAB mesh related quantities needed for FEM assembly */ 311 ierr = DMMoabGetLocalVertices(dm, &vlocal, NULL);CHKERRQ(ierr); 312 313 /* reset the residual vector */ 314 ierr = VecSet(F,0.0);CHKERRQ(ierr); 315 316 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 317 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 318 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 319 320 /* get the local representation of the arrays from Vectors */ 321 ierr = DMMoabVecGetArrayRead(dm, Xloc, &x);CHKERRQ(ierr); 322 ierr = DMMoabVecGetArrayRead(dm, Xdot, &xdot);CHKERRQ(ierr); 323 ierr = DMMoabVecGetArray(dm, F, &f);CHKERRQ(ierr); 324 325 /* loop over local elements */ 326 for(moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) { 327 const moab::EntityHandle vhandle = *iter; 328 329 ierr = DMMoabGetDofsBlockedLocal(dm,1,&vhandle,&i);CHKERRQ(ierr); 330 331 /* check if vertex is on the boundary */ 332 ierr = DMMoabIsEntityOnBoundary(dm,vhandle,&elem_on_boundary);CHKERRQ(ierr); 333 334 if (elem_on_boundary) { 335 ierr = DMMoabGetDofsBlocked(dm, 1, &vhandle, &bcindx);CHKERRQ(ierr); 336 if (bcindx == 0) { /* Apply left BC */ 337 f[i].u = hx * (x[i].u - user->leftbc.u); 338 f[i].v = hx * (x[i].v - user->leftbc.v); 339 } else { /* Apply right BC */ 340 f[i].u = hx * (x[i].u - user->rightbc.u); 341 f[i].v = hx * (x[i].v - user->rightbc.v); 342 } 343 } 344 else { 345 f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx; 346 f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx; 347 } 348 } 349 350 /* Restore data */ 351 ierr = DMMoabVecRestoreArrayRead(dm, Xloc, &x);CHKERRQ(ierr); 352 ierr = DMMoabVecRestoreArrayRead(dm, Xdot, &xdot);CHKERRQ(ierr); 353 ierr = DMMoabVecRestoreArray(dm, F, &f);CHKERRQ(ierr); 354 ierr = DMRestoreLocalVector(dm, &Xloc);CHKERRQ(ierr); 355 PetscFunctionReturn(0); 356 } 357 358 359 PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx) 360 { 361 UserCtx user = (UserCtx)ctx; 362 PetscReal vpos[3]; 363 DM dm; 364 Field *x; 365 PetscErrorCode ierr; 366 const moab::Range *vowned; 367 PetscInt dof; 368 moab::Range::iterator iter; 369 370 PetscFunctionBegin; 371 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 372 373 /* get the essential MOAB mesh related quantities needed for FEM assembly */ 374 ierr = DMMoabGetLocalVertices(dm, &vowned, NULL);CHKERRQ(ierr); 375 376 ierr = VecSet(X, 0.0);CHKERRQ(ierr); 377 378 /* Get pointers to vector data */ 379 ierr = DMMoabVecGetArray(dm, X, &x);CHKERRQ(ierr); 380 381 /* Compute function over the locally owned part of the grid */ 382 for(moab::Range::iterator iter = vowned->begin(); iter != vowned->end(); iter++) { 383 const moab::EntityHandle vhandle = *iter; 384 ierr = DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof);CHKERRQ(ierr); 385 386 /* compute the mid-point of the element and use a 1-point lumped quadrature */ 387 ierr = DMMoabGetVertexCoordinates(dm,1,&vhandle,vpos);CHKERRQ(ierr); 388 389 PetscReal xi = vpos[0]; 390 x[dof].u = user->leftbc.u*(1.-xi) + user->rightbc.u*xi + PetscSinReal(2.*PETSC_PI*xi); 391 x[dof].v = user->leftbc.v*(1.-xi) + user->rightbc.v*xi; 392 } 393 394 /* Restore vectors */ 395 ierr = DMMoabVecRestoreArray(dm, X, &x);CHKERRQ(ierr); 396 PetscFunctionReturn(0); 397 } 398 399 /*TEST 400 401 build: 402 requires: moab 403 404 test: 405 args: -n 20 -ts_type rosw -ts_rosw_type 2p -ts_dt 5e-2 -ts_adapt_type none 406 407 test: 408 suffix: 2 409 nsize: 2 410 args: -n 50 -ts_type glee -ts_adapt_type none -ts_dt 0.1 -io 411 TODO: 412 413 TEST*/ 414