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