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
Initialize_AppContext(UserCtx * puser)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(PETSC_SUCCESS);
72 }
73
Destroy_AppContext(UserCtx * user)74 PetscErrorCode Destroy_AppContext(UserCtx *user)
75 {
76 PetscFunctionBegin;
77 PetscCall(PetscFree(*user));
78 PetscFunctionReturn(PETSC_SUCCESS);
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 ****************/
main(int argc,char ** argv)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, nullptr, 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 */
FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void * ptr)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] = {
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 }
240
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 }
249
FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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;
259
260 PetscFunctionBegin;
261 hx = 1.0 / user->n;
262 PetscCall(TSGetDM(ts, &dm));
263
264 /* Get pointers to vector data */
265 PetscCall(VecSet(F, 0.0));
266
267 PetscCall(DMMoabVecGetArrayRead(dm, X, &x));
268 PetscCall(DMMoabVecGetArray(dm, F, &f));
269
270 PetscCall(DMMoabGetLocalVertices(dm, &ownedvtx, NULL));
271
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));
276
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 }
281
282 /* Restore vectors */
283 PetscCall(DMMoabVecRestoreArrayRead(dm, X, &x));
284 PetscCall(DMMoabVecRestoreArray(dm, F, &f));
285 PetscFunctionReturn(PETSC_SUCCESS);
286 }
287
FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,PetscCtx ctx)288 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx 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;
298
299 PetscFunctionBegin;
300 hx = 1.0 / user->n;
301 PetscCall(TSGetDM(ts, &dm));
302
303 /* get the essential MOAB mesh related quantities needed for FEM assembly */
304 PetscCall(DMMoabGetLocalVertices(dm, &vlocal, NULL));
305
306 /* reset the residual vector */
307 PetscCall(VecSet(F, 0.0));
308
309 PetscCall(DMGetLocalVector(dm, &Xloc));
310 PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
311 PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));
312
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));
317
318 /* loop over local elements */
319 for (moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) {
320 const moab::EntityHandle vhandle = *iter;
321
322 PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &i));
323
324 /* check if vertex is on the boundary */
325 PetscCall(DMMoabIsEntityOnBoundary(dm, vhandle, &elem_on_boundary));
326
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 }
341
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 }
349
FormInitialSolution(TS ts,Vec X,PetscCtx ctx)350 PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx 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;
359
360 PetscFunctionBegin;
361 PetscCall(TSGetDM(ts, &dm));
362
363 /* get the essential MOAB mesh related quantities needed for FEM assembly */
364 PetscCall(DMMoabGetLocalVertices(dm, &vowned, NULL));
365
366 PetscCall(VecSet(X, 0.0));
367
368 /* Get pointers to vector data */
369 PetscCall(DMMoabVecGetArray(dm, X, &x));
370
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));
375
376 /* compute the mid-point of the element and use a 1-point lumped quadrature */
377 PetscCall(DMMoabGetVertexCoordinates(dm, 1, &vhandle, vpos));
378
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 }
383
384 /* Restore vectors */
385 PetscCall(DMMoabVecRestoreArray(dm, X, &x));
386 PetscFunctionReturn(PETSC_SUCCESS);
387 }
388
389 /*TEST
390
391 build:
392 requires: moab !complex
393
394 test:
395 args: -n 20 -ts_type rosw -ts_rosw_type 2p -ts_time_step 5e-2 -ts_adapt_type none
396
397 test:
398 suffix: 2
399 nsize: 2
400 args: -n 50 -ts_type glee -ts_adapt_type none -ts_time_step 0.1 -io
401 TODO:
402
403 TEST*/
404