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