Actual source code: ex35.cxx
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
8: Initial conditions:
9: u(x,0) = 1 + sin(2 pi x)
10: v(x,0) = 3
12: Boundary conditions:
13: u(0,t) = u(1,t) = 1
14: v(0,t) = v(1,t) = 3
15: */
17: // PETSc includes:
18: #include <petscts.h>
19: #include <petscdmmoab.h>
21: typedef struct {
22: PetscScalar u,v;
23: } Field;
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;
37: PetscErrorCode Initialize_AppContext(UserCtx *puser)
38: {
39: UserCtx user;
40: PetscErrorCode ierr;
43: PetscNew(&user);
45: 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: PetscOptionsReal("-A","Reaction rate","ex35.cxx",user->A,&user->A,NULL);
59: PetscOptionsReal("-B","Reaction rate","ex35.cxx",user->B,&user->B,NULL);
60: PetscOptionsReal("-alpha","Diffusion coefficient","ex35.cxx",user->alpha,&user->alpha,NULL);
61: PetscOptionsScalar("-uleft","Dirichlet boundary condition","ex35.cxx",user->leftbc.u,&user->leftbc.u,NULL);
62: PetscOptionsScalar("-uright","Dirichlet boundary condition","ex35.cxx",user->rightbc.u,&user->rightbc.u,NULL);
63: PetscOptionsScalar("-vleft","Dirichlet boundary condition","ex35.cxx",user->leftbc.v,&user->leftbc.v,NULL);
64: PetscOptionsScalar("-vright","Dirichlet boundary condition","ex35.cxx",user->rightbc.v,&user->rightbc.v,NULL);
65: PetscOptionsInt("-n","Number of 1-D elements","ex35.cxx",user->n,&user->n,NULL);
66: PetscOptionsInt("-ndt","Number of time steps","ex35.cxx",user->ntsteps,&user->ntsteps,NULL);
67: 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: PetscOptionsEnd();
72: *puser = user;
73: return(0);
74: }
76: PetscErrorCode Destroy_AppContext(UserCtx *user)
77: {
81: PetscFree(*user);
82: return(0);
83: }
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*);
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;
106: DM dm;
107: const char *fields[2] = {"U","V"};
109: PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;
111: /* Initialize the user context struct */
112: Initialize_AppContext(&user);
114: /* Fill in the user defined work context: */
115: DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, PETSC_FALSE, NULL, user->n, 1, &dm);
116: DMMoabSetFieldNames(dm, user->nvars, fields);
117: DMMoabSetBlockSize(dm, user->nvars);
118: DMSetFromOptions(dm);
120: /* SetUp the data structures for DMMOAB */
121: DMSetUp(dm);
123: /* Create timestepping solver context */
124: TSCreate(PETSC_COMM_WORLD,&ts);
125: TSSetDM(ts, dm);
126: TSSetType(ts,TSARKIMEX);
127: TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);
128: DMSetMatType(dm,MATBAIJ);
129: DMCreateMatrix(dm,&J);
131: TSSetRHSFunction(ts,NULL,FormRHSFunction,user);
132: TSSetIFunction(ts,NULL,FormIFunction,user);
133: TSSetIJacobian(ts,J,J,FormIJacobian,user);
135: ftime = 10.0;
136: TSSetMaxSteps(ts,user->ntsteps);
137: TSSetMaxTime(ts,ftime);
138: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Create the solution vector and set the initial conditions
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: DMCreateGlobalVector(dm, &X);
145: FormInitialSolution(ts,X,user);
146: TSSetSolution(ts,X);
147: VecGetSize(X,&mx);
148: hx = 1.0/(PetscReal)(mx/2-1);
149: dt = 0.4 * PetscSqr(hx) / user->alpha; /* Diffusive stability limit */
150: TSSetTimeStep(ts,dt);
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Set runtime options
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155: TSSetFromOptions(ts);
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Solve nonlinear system
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160: TSSolve(ts,X);
161: TSGetSolveTime(ts,&ftime);
162: TSGetStepNumber(ts,&steps);
163: TSGetConvergedReason(ts,&reason);
164: PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],ftime,steps);
166: if (user->io) {
167: /* Print the numerical solution to screen and then dump to file */
168: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
170: /* Write out the solution along with the mesh */
171: DMMoabSetGlobalFieldVector(dm, X);
172: #ifdef MOAB_HAVE_HDF5
173: DMMoabOutput(dm, "ex35.h5m", "");
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: DMMoabOutput(dm, "ex35.vtk", "");
181: #endif
182: }
184: /* Free work space.
185: Free all PETSc related resources: */
186: MatDestroy(&J);
187: VecDestroy(&X);
188: TSDestroy(&ts);
189: DMDestroy(&dm);
191: /* Free all MOAB related resources: */
192: Destroy_AppContext(&user);
193: PetscFinalize();
194: return ierr;
195: }
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;
212: TSGetDM(ts, &dm);
214: /* get the essential MOAB mesh related quantities needed for FEM assembly */
215: DMMoabGetLocalVertices(dm, &vlocal, NULL);
217: /* compute local element sizes - structured grid */
218: hx = 1.0/user->n;
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;
226: DMMoabGetDofsBlocked(dm, 1, &vhandle, &dof);
228: /* check if vertex is on the boundary */
229: DMMoabIsEntityOnBoundary(dm,vhandle,&vonboundary);
231: if (vonboundary) {
232: const PetscScalar bcvals[2][2] = {{hx,0},{0,hx}};
233: MatSetValuesBlocked(Jpre,1,&dof,1,&dof,&bcvals[0][0],INSERT_VALUES);
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: MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);
241: }
242: }
244: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
245: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
246: if (J != Jpre) {
247: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
248: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
249: }
250: return(0);
251: }
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;
266: hx = 1.0/user->n;
267: TSGetDM(ts,&dm);
269: /* Get pointers to vector data */
270: VecSet(F,0.0);
272: DMMoabVecGetArrayRead(dm, X, &x);
273: DMMoabVecGetArray(dm, F, &f);
275: DMMoabGetLocalVertices(dm, &ownedvtx, NULL);
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: DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof);
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: }
287: /* Restore vectors */
288: DMMoabVecRestoreArrayRead(dm, X, &x);
289: DMMoabVecRestoreArray(dm, F, &f);
290: return(0);
291: }
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;
307: hx = 1.0/user->n;
308: TSGetDM(ts, &dm);
310: /* get the essential MOAB mesh related quantities needed for FEM assembly */
311: DMMoabGetLocalVertices(dm, &vlocal, NULL);
313: /* reset the residual vector */
314: VecSet(F,0.0);
316: DMGetLocalVector(dm,&Xloc);
317: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
318: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
320: /* get the local representation of the arrays from Vectors */
321: DMMoabVecGetArrayRead(dm, Xloc, &x);
322: DMMoabVecGetArrayRead(dm, Xdot, &xdot);
323: DMMoabVecGetArray(dm, F, &f);
325: /* loop over local elements */
326: for (moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) {
327: const moab::EntityHandle vhandle = *iter;
329: DMMoabGetDofsBlockedLocal(dm,1,&vhandle,&i);
331: /* check if vertex is on the boundary */
332: DMMoabIsEntityOnBoundary(dm,vhandle,&elem_on_boundary);
334: if (elem_on_boundary) {
335: DMMoabGetDofsBlocked(dm, 1, &vhandle, &bcindx);
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: }
350: /* Restore data */
351: DMMoabVecRestoreArrayRead(dm, Xloc, &x);
352: DMMoabVecRestoreArrayRead(dm, Xdot, &xdot);
353: DMMoabVecRestoreArray(dm, F, &f);
354: DMRestoreLocalVector(dm, &Xloc);
355: return(0);
356: }
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;
371: TSGetDM(ts, &dm);
373: /* get the essential MOAB mesh related quantities needed for FEM assembly */
374: DMMoabGetLocalVertices(dm, &vowned, NULL);
376: VecSet(X, 0.0);
378: /* Get pointers to vector data */
379: DMMoabVecGetArray(dm, X, &x);
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: DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof);
386: /* compute the mid-point of the element and use a 1-point lumped quadrature */
387: DMMoabGetVertexCoordinates(dm,1,&vhandle,vpos);
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: }
394: /* Restore vectors */
395: DMMoabVecRestoreArray(dm, X, &x);
396: return(0);
397: }
399: /*TEST
401: build:
402: requires: moab
404: test:
405: args: -n 20 -ts_type rosw -ts_rosw_type 2p -ts_dt 5e-2 -ts_adapt_type none
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:
413: TEST*/