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: 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();
70: *puser = user;
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: PetscErrorCode Destroy_AppContext(UserCtx *user)
75: {
76: PetscFunctionBegin;
77: PetscCall(PetscFree(*user));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
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 *);
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"};
103: PetscFunctionBeginUser;
104: PetscCall(PetscInitialize(&argc, &argv, nullptr, help));
106: /* Initialize the user context struct */
107: PetscCall(Initialize_AppContext(&user));
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));
115: /* SetUp the data structures for DMMOAB */
116: PetscCall(DMSetUp(dm));
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));
126: PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, user));
127: PetscCall(TSSetIFunction(ts, NULL, FormIFunction, user));
128: PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, user));
130: ftime = 10.0;
131: PetscCall(TSSetMaxSteps(ts, user->ntsteps));
132: PetscCall(TSSetMaxTime(ts, ftime));
133: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Create the solution vector and set the initial conditions
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: PetscCall(DMCreateGlobalVector(dm, &X));
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));
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Set runtime options
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: PetscCall(TSSetFromOptions(ts));
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));
161: if (user->io) {
162: /* Print the numerical solution to screen and then dump to file */
163: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
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: }
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));
186: /* Free all MOAB related resources: */
187: PetscCall(Destroy_AppContext(&user));
188: PetscCall(PetscFinalize());
189: return 0;
190: }
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;
204: PetscFunctionBegin;
205: PetscCall(TSGetDM(ts, &dm));
207: /* get the essential MOAB mesh related quantities needed for FEM assembly */
208: PetscCall(DMMoabGetLocalVertices(dm, &vlocal, NULL));
210: /* compute local element sizes - structured grid */
211: hx = 1.0 / user->n;
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;
219: PetscCall(DMMoabGetDofsBlocked(dm, 1, &vhandle, &dof));
221: /* check if vertex is on the boundary */
222: PetscCall(DMMoabIsEntityOnBoundary(dm, vhandle, &vonboundary));
224: if (vonboundary) {
225: const PetscScalar bcvals[2][2] = {
226: {hx, 0 },
227: {0, hx}
228: };
229: PetscCall(MatSetValuesBlocked(Jpre, 1, &dof, 1, &dof, &bcvals[0][0], INSERT_VALUES));
230: } else {
231: const PetscInt row = dof, col[] = {dof - 1, dof, dof + 1};
232: const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx;
233: const PetscScalar vals[2][3][2] = {
234: {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}},
235: {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}}
236: };
237: PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES));
238: }
239: }
241: PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
242: PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
243: if (J != Jpre) {
244: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
245: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
246: }
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
251: {
252: UserCtx user = (UserCtx)ptr;
253: DM dm;
254: PetscReal hx;
255: const Field *x;
256: Field *f;
257: PetscInt dof;
258: const moab::Range *ownedvtx;
260: PetscFunctionBegin;
261: hx = 1.0 / user->n;
262: PetscCall(TSGetDM(ts, &dm));
264: /* Get pointers to vector data */
265: PetscCall(VecSet(F, 0.0));
267: PetscCall(DMMoabVecGetArrayRead(dm, X, &x));
268: PetscCall(DMMoabVecGetArray(dm, F, &f));
270: PetscCall(DMMoabGetLocalVertices(dm, &ownedvtx, NULL));
272: /* Compute function over the locally owned part of the grid */
273: for (moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) {
274: const moab::EntityHandle vhandle = *iter;
275: PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof));
277: PetscScalar u = x[dof].u, v = x[dof].v;
278: f[dof].u = hx * (user->A + u * u * v - (user->B + 1) * u);
279: f[dof].v = hx * (user->B * u - u * u * v);
280: }
282: /* Restore vectors */
283: PetscCall(DMMoabVecRestoreArrayRead(dm, X, &x));
284: PetscCall(DMMoabVecRestoreArray(dm, F, &f));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
289: {
290: UserCtx user = (UserCtx)ctx;
291: DM dm;
292: Field *x, *xdot, *f;
293: PetscReal hx;
294: Vec Xloc;
295: PetscInt i, bcindx;
296: PetscBool elem_on_boundary;
297: const moab::Range *vlocal;
299: PetscFunctionBegin;
300: hx = 1.0 / user->n;
301: PetscCall(TSGetDM(ts, &dm));
303: /* get the essential MOAB mesh related quantities needed for FEM assembly */
304: PetscCall(DMMoabGetLocalVertices(dm, &vlocal, NULL));
306: /* reset the residual vector */
307: PetscCall(VecSet(F, 0.0));
309: PetscCall(DMGetLocalVector(dm, &Xloc));
310: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
311: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
313: /* get the local representation of the arrays from Vectors */
314: PetscCall(DMMoabVecGetArrayRead(dm, Xloc, &x));
315: PetscCall(DMMoabVecGetArrayRead(dm, Xdot, &xdot));
316: PetscCall(DMMoabVecGetArray(dm, F, &f));
318: /* loop over local elements */
319: for (moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) {
320: const moab::EntityHandle vhandle = *iter;
322: PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &i));
324: /* check if vertex is on the boundary */
325: PetscCall(DMMoabIsEntityOnBoundary(dm, vhandle, &elem_on_boundary));
327: if (elem_on_boundary) {
328: PetscCall(DMMoabGetDofsBlocked(dm, 1, &vhandle, &bcindx));
329: if (bcindx == 0) { /* Apply left BC */
330: f[i].u = hx * (x[i].u - user->leftbc.u);
331: f[i].v = hx * (x[i].v - user->leftbc.v);
332: } else { /* Apply right BC */
333: f[i].u = hx * (x[i].u - user->rightbc.u);
334: f[i].v = hx * (x[i].v - user->rightbc.v);
335: }
336: } else {
337: f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx;
338: f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx;
339: }
340: }
342: /* Restore data */
343: PetscCall(DMMoabVecRestoreArrayRead(dm, Xloc, &x));
344: PetscCall(DMMoabVecRestoreArrayRead(dm, Xdot, &xdot));
345: PetscCall(DMMoabVecRestoreArray(dm, F, &f));
346: PetscCall(DMRestoreLocalVector(dm, &Xloc));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
351: {
352: UserCtx user = (UserCtx)ctx;
353: PetscReal vpos[3];
354: DM dm;
355: Field *x;
356: const moab::Range *vowned;
357: PetscInt dof;
358: moab::Range::iterator iter;
360: PetscFunctionBegin;
361: PetscCall(TSGetDM(ts, &dm));
363: /* get the essential MOAB mesh related quantities needed for FEM assembly */
364: PetscCall(DMMoabGetLocalVertices(dm, &vowned, NULL));
366: PetscCall(VecSet(X, 0.0));
368: /* Get pointers to vector data */
369: PetscCall(DMMoabVecGetArray(dm, X, &x));
371: /* Compute function over the locally owned part of the grid */
372: for (moab::Range::iterator iter = vowned->begin(); iter != vowned->end(); iter++) {
373: const moab::EntityHandle vhandle = *iter;
374: PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof));
376: /* compute the mid-point of the element and use a 1-point lumped quadrature */
377: PetscCall(DMMoabGetVertexCoordinates(dm, 1, &vhandle, vpos));
379: PetscReal xi = vpos[0];
380: x[dof].u = user->leftbc.u * (1. - xi) + user->rightbc.u * xi + PetscSinReal(2. * PETSC_PI * xi);
381: x[dof].v = user->leftbc.v * (1. - xi) + user->rightbc.v * xi;
382: }
384: /* Restore vectors */
385: PetscCall(DMMoabVecRestoreArray(dm, X, &x));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: /*TEST
391: build:
392: requires: moab !complex
394: test:
395: args: -n 20 -ts_type rosw -ts_rosw_type 2p -ts_dt 5e-2 -ts_adapt_type none
397: test:
398: suffix: 2
399: nsize: 2
400: args: -n 50 -ts_type glee -ts_adapt_type none -ts_dt 0.1 -io
401: TODO:
403: TEST*/