Actual source code: ex35.cxx
petsc-3.7.5 2017-01-01
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;
39: PetscErrorCode Initialize_AppContext(UserCtx *puser)
40: {
41: UserCtx user;
42: PetscErrorCode ierr;
45: PetscNew(&user);
47: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Advection-reaction options","ex35.cxx");
48: {
49: user->nvars = 2;
50: user->A = 1;
51: user->B = 3;
52: user->alpha = 0.02;
53: user->leftbc.u = 1;
54: user->rightbc.u = 1;
55: user->leftbc.v = 3;
56: user->rightbc.v = 3;
57: user->n = 10;
58: user->ntsteps = 10000;
59: user->io = PETSC_FALSE;
60: PetscOptionsReal("-A","Reaction rate","ex35.cxx",user->A,&user->A,NULL);
61: PetscOptionsReal("-B","Reaction rate","ex35.cxx",user->B,&user->B,NULL);
62: PetscOptionsReal("-alpha","Diffusion coefficient","ex35.cxx",user->alpha,&user->alpha,NULL);
63: PetscOptionsScalar("-uleft","Dirichlet boundary condition","ex35.cxx",user->leftbc.u,&user->leftbc.u,NULL);
64: PetscOptionsScalar("-uright","Dirichlet boundary condition","ex35.cxx",user->rightbc.u,&user->rightbc.u,NULL);
65: PetscOptionsScalar("-vleft","Dirichlet boundary condition","ex35.cxx",user->leftbc.v,&user->leftbc.v,NULL);
66: PetscOptionsScalar("-vright","Dirichlet boundary condition","ex35.cxx",user->rightbc.v,&user->rightbc.v,NULL);
67: PetscOptionsInt("-n","Number of 1-D elements","ex35.cxx",user->n,&user->n,NULL);
68: PetscOptionsInt("-ndt","Number of time steps","ex35.cxx",user->ntsteps,&user->ntsteps,NULL);
69: PetscOptionsBool("-io","Write the mesh and solution output to a file.","ex35.cxx",user->io,&user->io,NULL);
70: user->npts = user->n+1;
71: }
72: PetscOptionsEnd();
74: *puser = user;
75: return(0);
76: }
80: PetscErrorCode Destroy_AppContext(UserCtx *user)
81: {
85: PetscFree(*user);
86: return(0);
87: }
89: static PetscErrorCode FormInitialSolution(TS,Vec,void*);
90: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
92: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
93: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
95: /****************
96: * *
97: * MAIN *
98: * *
99: ****************/
102: int main(int argc,char **argv)
103: {
104: TS ts; /* nonlinear solver */
105: Vec X; /* solution, residual vectors */
106: Mat J; /* Jacobian matrix */
107: PetscInt steps,mx;
108: PetscErrorCode ierr;
109: PetscReal hx,dt,ftime;
110: UserCtx user; /* user-defined work context */
111: TSConvergedReason reason;
113: DM dm;
114: moab::ErrorCode merr;
115: moab::Interface* mbImpl;
116: moab::Tag solndofs;
117: const moab::Range *ownedvtx;
118: const PetscReal bounds[2] = {0.0,1.0};
119: const char *fields[2] = {"U","V"};
120: PetscScalar deflt[2]={0.0,0.0};
122: PetscInitialize(&argc,&argv,(char *)0,help);
124: /* Initialize the user context struct */
125: Initialize_AppContext(&user);
127: /* Fill in the user defined work context: */
128: DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, PETSC_FALSE, bounds, user->n, 1, &dm);
129: DMMoabSetBlockSize(dm, user->nvars);
130: DMMoabSetFieldNames(dm, user->nvars, fields);
131: DMSetMatType(dm,MATBAIJ);
132: DMSetFromOptions(dm);
134: /* SetUp the data structures for DMMOAB */
135: DMSetUp(dm);
137: DMMoabGetInterface(dm, &mbImpl);
139: /* Create timestepping solver context */
140: TSCreate(PETSC_COMM_WORLD,&ts);
141: TSSetDM(ts, dm);
142: TSSetType(ts,TSARKIMEX);
143: TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);
144: TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE);
145: TSSetRHSFunction(ts,NULL,FormRHSFunction,user);
146: TSSetIFunction(ts,NULL,FormIFunction,user);
147: DMCreateMatrix(dm,&J);
148: TSSetIJacobian(ts,J,J,FormIJacobian,user);
150: ftime = 10.0;
151: TSSetDuration(ts,user->ntsteps,ftime);
152: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
153:
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Create the solution vector and set the initial conditions
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: /* Use the call to DMMoabCreateVector for creating a named global MOAB Vec object.
158: Alternately, use the following call to DM for creating an unnamed (anonymous) global
159: MOAB Vec object.
161: DMCreateGlobalVector(dm, &X);
162: */
163: DMMoabGetLocalVertices(dm, &ownedvtx, NULL);
164: /* create a tag to store the unknown fields in the problem */
165: merr = mbImpl->tag_get_handle("UNKNOWNS",2,moab::MB_TYPE_DOUBLE,solndofs,
166: moab::MB_TAG_DENSE|moab::MB_TAG_CREAT,deflt);MBERRNM(merr);
168: DMMoabCreateVector(dm, solndofs, ownedvtx, PETSC_TRUE, PETSC_FALSE, &X);
170: FormInitialSolution(ts,X,user);
171: TSSetSolution(ts,X);
172: VecGetSize(X,&mx);
173: hx = 1.0/(PetscReal)(mx/2-1);
174: dt = 0.4 * PetscSqr(hx) / user->alpha; /* Diffusive stability limit */
175: TSSetInitialTimeStep(ts,0.0,dt);
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Set runtime options
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: TSSetFromOptions(ts);
182: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: Solve nonlinear system
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185: TSSolve(ts,X);
186: TSGetSolveTime(ts,&ftime);
187: TSGetTimeStepNumber(ts,&steps);
188: TSGetConvergedReason(ts,&reason);
189: PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],ftime,steps);
191: if (user->io) {
192: /* Print the numerical solution to screen and then dump to file */
193: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
195: /* Write out the solution along with the mesh */
196: DMMoabSetGlobalFieldVector(dm, X);
197: #ifdef MOAB_HDF5_H
198: DMMoabOutput(dm, "ex35.h5m", "");
199: #else
200: /* MOAB does not support true parallel writers that aren't HDF5 based
201: And so if you are using VTK as the output format in parallel,
202: the data could be jumbled due to the order in which the processors
203: write out their parts of the mesh and solution tags
204: */
205: DMMoabOutput(dm, "ex35.vtk", "");
206: #endif
207: }
209: /* Free work space.
210: Free all PETSc related resources: */
211: MatDestroy(&J);
212: VecDestroy(&X);
213: TSDestroy(&ts);
214: DMDestroy(&dm);
216: /* Free all MOAB related resources: */
217: Destroy_AppContext(&user);
219: PetscFinalize();
220: return 0;
221: }
225: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
226: {
227: UserCtx user = (UserCtx)ptr;
228: DM dm;
229: PetscReal hx;
230: const Field *x;
231: Field *f;
232: PetscInt dof_index;
233: const moab::Range *ownedvtx;
234: PetscErrorCode ierr;
237: hx = 1.0/user->n;
238: TSGetDM(ts,&dm);
240: /* Get pointers to vector data */
241: VecSet(F,0.0);
243: DMMoabVecGetArrayRead(dm, X, &x);
244: DMMoabVecGetArray(dm, F, &f);
246: DMMoabGetLocalVertices(dm, &ownedvtx, NULL);
248: /* Compute function over the locally owned part of the grid */
249: for(moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) {
250: const moab::EntityHandle vhandle = *iter;
251: DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof_index);
253: const Field& xx = x[dof_index];
254: f[dof_index].u = hx*(user->A + xx.u*xx.u*xx.v - (user->B+1)*xx.u);
255: f[dof_index].v = hx*(user->B*xx.u - xx.u*xx.u*xx.v);
256: }
258: /* Restore vectors */
259: DMMoabVecRestoreArrayRead(dm, X, &x);
260: DMMoabVecRestoreArray(dm, F, &f);
261: return(0);
262: }
265: /*
266: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
267: */
270: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr)
271: {
272: UserCtx user = (UserCtx)ptr;
273: PetscErrorCode ierr;
274: const moab::EntityHandle *connect;
275: PetscInt vpere=2;
276: PetscReal hx;
277: DM dm;
278: moab::Interface* mbImpl;
279: const moab::Range *elocal;
280: PetscInt dof_indices[2];
281: PetscBool elem_on_boundary;
284: TSGetDM(ts, &dm);
286: /* get the essential MOAB mesh related quantities needed for FEM assembly */
287: DMMoabGetInterface(dm, &mbImpl);
288: DMMoabGetLocalElements(dm, &elocal);
290: /* zero out the discrete operator */
291: MatZeroEntries(Jpre);
293: /* compute local element sizes - structured grid */
294: hx = 1.0/user->n;
296: const int& idl = dof_indices[0];
297: const int& idr = dof_indices[1];
298: const PetscScalar dxxL = user->alpha/hx,dxxR = -user->alpha/hx;
299: const PetscScalar bcvals[2][2] = {{hx,0},{0,hx}};
300: const PetscScalar e_vals[2][2][2] = {{{a *hx/2+dxxL,0},{dxxR,0}},
301: {{0,a*hx/2+dxxL},{0,dxxR}}};
303: /* Compute function over the locally owned part of the grid
304: Assemble the operator by looping over edges and computing
305: contribution for each vertex dof */
306: for(moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
307: const moab::EntityHandle ehandle = *iter;
309: // Get connectivity information in canonical order
310: DMMoabGetElementConnectivity(dm, ehandle, &vpere, &connect);
311: if (vpere != 2) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only EDGE2 element bases are supported in the current example. n(Connectivity)=%D.\n", vpere);
313: DMMoabGetDofsBlocked(dm, vpere, connect, dof_indices);
315: const PetscInt lcols[] = {idl,idr}, rcols[] = {idr, idl};
317: /* check if element is on the boundary */
318: DMMoabIsEntityOnBoundary(dm,ehandle,&elem_on_boundary);
320: if (elem_on_boundary) {
321: if (idl == 0) {
322: // Left Boundary conditions...
323: MatSetValuesBlocked(Jpre,1,&idl,1,&idl,&bcvals[0][0],ADD_VALUES);
324: MatSetValuesBlocked(Jpre,1,&idr,2,rcols,&e_vals[0][0][0],ADD_VALUES);
325: }
326: else {
327: // Right Boundary conditions...
328: MatSetValuesBlocked(Jpre,1,&idr,1,&idr,&bcvals[0][0],ADD_VALUES);
329: MatSetValuesBlocked(Jpre,1,&idl,2,lcols,&e_vals[0][0][0],ADD_VALUES);
330: }
331: }
332: else {
333: MatSetValuesBlocked(Jpre,1,&idr,2,rcols,&e_vals[0][0][0],ADD_VALUES);
334: MatSetValuesBlocked(Jpre,1,&idl,2,lcols,&e_vals[0][0][0],ADD_VALUES);
335: }
336: }
338: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
339: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
340: if (J != Jpre) {
341: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
342: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
343: }
344: return(0);
345: }
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: PetscErrorCode ierr;
357: moab::Interface* mbImpl;
358: const moab::Range *vowned;
359: PetscInt dof_index;
360: moab::Range::iterator iter;
363: TSGetDM(ts, &dm);
364:
365: /* get the essential MOAB mesh related quantities needed for FEM assembly */
366: DMMoabGetInterface(dm, &mbImpl);
368: DMMoabGetLocalVertices(dm, &vowned, NULL);
370: VecSet(X, 0.0);
372: /* Get pointers to vector data */
373: DMMoabVecGetArray(dm, X, &x);
375: /* Compute function over the locally owned part of the grid */
376: for(moab::Range::iterator iter = vowned->begin(); iter != vowned->end(); iter++) {
377: const moab::EntityHandle vhandle = *iter;
378: DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof_index);
380: /* compute the mid-point of the element and use a 1-point lumped quadrature */
381: DMMoabGetVertexCoordinates(dm,1,&vhandle,vpos);
383: PetscReal xi = vpos[0];
384: x[dof_index].u = user->leftbc.u*(1.-xi) + user->rightbc.u*xi + sin(2.*PETSC_PI*xi);
385: x[dof_index].v = user->leftbc.v*(1.-xi) + user->rightbc.v*xi;
386: }
388: /* Restore vectors */
389: DMMoabVecRestoreArray(dm, X, &x);
390: return(0);
391: }
396: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
397: {
398: UserCtx user = (UserCtx)ctx;
399: DM dm;
400: Field *x,*xdot,*f;
401: PetscReal hx,vpos[2*3];
402: PetscErrorCode ierr;
403: PetscInt dof_indices[2],bc_indices[2];
404: const moab::EntityHandle *connect;
405: PetscInt vpere=2,nloc,ngh;
406: PetscBool elem_on_boundary;
407: const int& idx_left = dof_indices[0];
408: const int& idx_right = dof_indices[1];
409: moab::Interface* mbImpl;
410: const moab::Range *elocal;
413: hx = 1.0/user->n;
414: TSGetDM(ts, &dm);
416: /* get the essential MOAB mesh related quantities needed for FEM assembly */
417: DMMoabGetInterface(dm, &mbImpl);
418: DMMoabGetLocalElements(dm, &elocal);
419: DMMoabGetLocalSize(dm, NULL, NULL, &nloc, &ngh);
421: /* reset the residual vector */
422: VecSet(F,0.0);
424: /* get the local representation of the arrays from Vectors */
425: DMMoabVecGetArrayRead(dm, X, &x);
426: DMMoabVecGetArrayRead(dm, Xdot, &xdot);
427: DMMoabVecGetArray(dm, F, &f);
429: /* loop over local elements */
430: for(moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
431: const moab::EntityHandle ehandle = *iter;
433: // Get connectivity information in canonical order
434: DMMoabGetElementConnectivity(dm, ehandle, &vpere, &connect);
435: if (vpere != 2) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only EDGE2 element bases are supported in the current example. n(Connectivity)=%D.\n", vpere);
437: /* compute the mid-point of the element and use a 1-point lumped quadrature */
438: DMMoabGetVertexCoordinates(dm,vpere,connect,vpos);
440: DMMoabGetDofsBlockedLocal(dm, vpere, connect, dof_indices);
442: /* check if element is on the boundary */
443: DMMoabIsEntityOnBoundary(dm,ehandle,&elem_on_boundary);
445: if (elem_on_boundary) {
446: DMMoabGetDofsBlocked(dm, vpere, connect, bc_indices);
447: if (bc_indices[0] == 0) { /* Apply left BC */
448: f[idx_left].u = hx * (x[idx_left].u - user->leftbc.u);
449: f[idx_left].v = hx * (x[idx_left].v - user->leftbc.v);
450: f[idx_right].u += user->alpha*(x[idx_right].u-x[idx_left].u)/hx;
451: f[idx_right].v += user->alpha*(x[idx_right].v-x[idx_left].v)/hx;
452: }
453: else { /* Apply right BC */
454: f[idx_left].u += hx * xdot[idx_left].u + user->alpha*(x[idx_left].u - x[idx_right].u)/hx;
455: f[idx_left].v += hx * xdot[idx_left].v + user->alpha*(x[idx_left].v - x[idx_right].v)/hx;
456: f[idx_right].u = hx * (x[idx_right].u - user->rightbc.u);
457: f[idx_right].v = hx * (x[idx_right].v - user->rightbc.v);
458: }
459: }
460: else {
461: f[idx_left].u += hx * xdot[idx_left].u + user->alpha*(x[idx_left].u - x[idx_right].u)/hx;
462: f[idx_left].v += hx * xdot[idx_left].v + user->alpha*(x[idx_left].v - x[idx_right].v)/hx;
463: f[idx_right].u += user->alpha*(x[idx_right].u-x[idx_left].u)/hx;
464: f[idx_right].v += user->alpha*(x[idx_right].v-x[idx_left].v)/hx;
465: }
466: }
468: /* Restore data */
469: DMMoabVecRestoreArrayRead(dm, X, &x);
470: DMMoabVecRestoreArrayRead(dm, Xdot, &xdot);
471: DMMoabVecRestoreArray(dm, F, &f);
472: return(0);
473: }