Actual source code: ex30.cxx
petsc-3.6.1 2015-08-06
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: // MOAB includes:
22: #if defined (PETSC_HAVE_MOAB)
23: # include <moab/Core.hpp>
24: # include <moab/ReadUtilIface.hpp>
25: # include <MBTagConventions.hpp>
26: #else
27: #error You must have MOAB for this example. Reconfigure using --download-moab
28: #endif
30: typedef struct {
31: PetscScalar u,v;
32: } Field;
34: typedef struct _User *User;
35: struct _User {
36: PetscReal A,B; /* Reaction coefficients */
37: PetscReal alpha; /* Diffusion coefficient */
38: PetscReal uleft,uright; /* Dirichlet boundary conditions */
39: PetscReal vleft,vright; /* Dirichlet boundary conditions */
40: PetscInt npts; /* Number of mesh points */
41: PetscBool io;
43: moab::ParallelComm *pcomm;
44: moab::Interface *mbint;
45: moab::Range *owned_vertexes;
46: moab::Range *owned_edges;
47: moab::Range *all_vertexes;
48: moab::Range *shared_vertexes;
49: moab::Tag unknowns_tag;
50: PetscInt unknowns_tag_size;
51: moab::Tag id_tag;
52: };
54: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
55: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
56: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
58: PetscErrorCode create_app_data(_User& user);
59: PetscErrorCode destroy_app_data(_User& user);
60: PetscErrorCode create_matrix(_User &user, DM &dm, Mat *J);
62: /****************
63: * *
64: * MAIN *
65: * *
66: ****************/
69: int main(int argc,char **argv)
70: {
71: TS ts; /* nonlinear solver */
72: Vec X; /* solution, residual vectors */
73: Mat J; /* Jacobian matrix */
74: PetscInt steps,maxsteps,mx;
75: PetscErrorCode ierr;
76: PetscReal ftime,hx,dt;
77: _User user; /* user-defined work context */
78: TSConvergedReason reason;
79: DM dm;
80: moab::ErrorCode merr;
82: PetscInitialize(&argc,&argv,(char *)0,help);
84: /* Fill in the user defined work context (creates mesh, solution field on mesh) */
85: create_app_data(user);
87: /* Create the DM to manage the mesh */
88: DMMoabCreateMoab(user.pcomm->comm(),user.mbint,user.pcomm,&user.id_tag,user.owned_vertexes,&dm);
89: DMSetUp(dm);
91: /* Create the solution vector: */
92: DMMoabCreateVector(dm,user.unknowns_tag,user.owned_vertexes,PETSC_TRUE,PETSC_FALSE,&X);
93: /* Create the matrix */
94: create_matrix(user, dm, &J);
96: /* Create timestepping solver context */
97: TSCreate(PETSC_COMM_WORLD,&ts);
98: TSSetType(ts,TSARKIMEX);
100: TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);
101: TSSetIFunction(ts,NULL,FormIFunction,&user);
102: TSSetIJacobian(ts,J,J,FormIJacobian,&user);
103: TSSetDM(ts,dm);
105: ftime = 10.0;
106: maxsteps = 10000;
107: TSSetDuration(ts,maxsteps,ftime);
109: /* Set initial conditions */
110: TSSetSolution(ts,X);
111: VecGetSize(X,&mx);
112: hx = 1.0/(PetscReal)(mx/2-1);
113: dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
114: TSSetInitialTimeStep(ts,0.0,dt);
116: /* Set runtime options */
117: TSSetFromOptions(ts);
119: /* Solve nonlinear system */
120: TSSolve(ts,X);
121: TSGetSolveTime(ts,&ftime);
122: TSGetTimeStepNumber(ts,&steps);
123: TSGetConvergedReason(ts,&reason);
124: PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
126: /* Write out the final mesh */
127: if (user.io) {
128: #ifdef MOAB_HDF5_H
129: merr = user.mbint->write_file("ex30.h5m");MBERRNM(merr);
130: #else
131: merr = user.mbint->write_file("ex30.vtk");MBERRNM(merr);
132: #endif
133: }
135: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
137: /* Free work space. */
138: MatDestroy(&J);
139: VecDestroy(&X);
140: TSDestroy(&ts);
141: DMDestroy(&dm);
143: /* Free all MOAB related resources */
144: destroy_app_data(user);
146: PetscFinalize();
147: return 0;
148: }
152: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
153: {
154: User user = (User)ptr;
155: PetscInt i;
156: Field *x,*xdot,*f;
157: PetscReal hx;
158: PetscErrorCode ierr;
159: moab::ErrorCode merr;
160: int *vertex_ids;
163: hx = 1.0/(PetscReal)(user->npts-1);
165: // Get connectivity information:
166: int verts_per_entity;
167: int count,num_edges;
168: moab::EntityHandle *connect;
169: moab::Range::iterator iter = user->owned_edges->begin();
170: merr = user->mbint->connect_iterate(iter,user->owned_edges->end(),connect,verts_per_entity,num_edges);MBERRNM(merr);
172: // Get tag data:
173: moab::Tag x_tag;
174: DMMoabGetVecTag(X,&x_tag);
175: merr = user->mbint->tag_iterate(x_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
176: count,reinterpret_cast<void*&>(x));MBERRNM(merr);
178: moab::Tag xdot_tag;
179: DMMoabGetVecTag(Xdot,&xdot_tag);
180: merr = user->mbint->tag_iterate(xdot_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
181: count,reinterpret_cast<void*&>(xdot));MBERRNM(merr);
183: moab::Tag f_tag;
184: DMMoabGetVecTag(F,&f_tag);
185: merr = user->mbint->tag_iterate(f_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
186: count,reinterpret_cast<void*&>(f));MBERRNM(merr);
188: // Get the global IDs on all vertexes:
189: moab::Tag id_tag;
190: user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);
191: merr = user->mbint->tag_iterate(id_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
192: count,reinterpret_cast<void*&>(vertex_ids));MBERRNM(merr);
194: // Exchange tags that are needed for assembly:
195: merr = user->pcomm->exchange_tags(x_tag,*user->all_vertexes);MBERRNM(merr);
196: merr = user->pcomm->exchange_tags(xdot_tag,*user->all_vertexes);MBERRNM(merr);
198: // Compute f:
199: const Field zero_field = {0.0,0.0};
200: std::fill(f,f+user->all_vertexes->size(),zero_field);
202: const moab::EntityHandle first_vertex = *user->all_vertexes->begin();
204: for (i = 0; i < num_edges; i++) {
205: const moab::EntityHandle idx_left = connect[2*i]-first_vertex;
206: const moab::EntityHandle idx_right = connect[2*i+1]-first_vertex;
208: const int id_left = vertex_ids[idx_left ];
209: const int id_right = vertex_ids[idx_right];
211: if (id_left == 0) {
212: // Apply left BC
213: f[idx_left].u += hx * (x[idx_left].u - user->uleft);
214: f[idx_left].v += hx * (x[idx_left].v - user->vleft);
215: } else {
216: f[idx_left].u += hx * xdot[idx_left].u - user->alpha*(-2*x[idx_left].u + x[idx_right].u)/hx;
217: f[idx_left].v += hx * xdot[idx_left].v - user->alpha*(-2*x[idx_left].v + x[idx_right].v)/hx;
218: }
220: if (id_right == user->npts-1) {
221: // Apply right BC
222: f[idx_right].u += hx * (x[idx_right].u - user->uright);
223: f[idx_right].v += hx * (x[idx_right].v - user->vright);
224: } else {
225: f[idx_right].u -= user->alpha*x[idx_left].u/hx;
226: f[idx_right].v -= user->alpha*x[idx_left].v/hx;
227: }
228: }
230: // Add tags on shared vertexes:
231: merr = user->pcomm->reduce_tags(f_tag,MPI_SUM,*user->shared_vertexes);MBERRNM(merr);
232: return(0);
233: }
237: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
238: {
239: User user = (User)ptr;
240: PetscReal hx;
241: Field *x,*f;
242: PetscErrorCode ierr;
243: moab::Tag id_tag;
244: moab::ErrorCode merr;
247: hx = 1.0/(PetscReal)(user->npts-1);
249: merr = user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr);
251: /* Get pointers to vector data */
252: VecGetArray(X,reinterpret_cast<PetscScalar**>(&x));
253: VecGetArray(F,reinterpret_cast<PetscScalar**>(&f));
255: /* Compute function over the locally owned part of the grid */
256: const moab::EntityHandle first_vertex = *user->owned_vertexes->begin();
257: moab::Range::iterator iter;
258: for(iter = user->owned_vertexes->begin(); iter != user->owned_vertexes->end(); iter++) {
259: moab::EntityHandle i = *iter - first_vertex;
260: PetscScalar u = x[i].u,v = x[i].v;
261: f[i].u = hx*(user->A + u*u*v - (user->B+1)*u);
262: f[i].v = hx*(user->B*u - u*u*v);
263: }
265: /* Restore vectors */
266: VecRestoreArray(X,reinterpret_cast<PetscScalar**>(&x));
267: VecRestoreArray(F,reinterpret_cast<PetscScalar**>(&f));
268: return(0);
269: }
271: /* --------------------------------------------------------------------- */
272: /*
273: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
274: */
277: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr)
278: {
279: User user = (User)ptr;
280: PetscErrorCode ierr;
281: PetscInt i;
282: PetscReal hx;
283: moab::Tag id_tag;
284: moab::ErrorCode merr;
287: hx = 1.0/(PetscReal)(user->npts-1);
289: merr = user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr);
291: /* Compute function over the locally owned part of the grid */
292: moab::Range::iterator iter;
293: for(iter = user->owned_vertexes->begin(); iter != user->owned_vertexes->end(); iter++) {
294: merr = user->mbint->tag_get_data(id_tag,&(*iter),1,&i);MBERRNM(merr);
296: if (i == 0 || i == user->npts-1) {
297: // Boundary conditions...
298: const PetscInt row = i,col = i;
299: const PetscScalar vals[2][2] = {{hx,0},{0,hx}};
300: MatSetValuesBlocked(Jpre,1,&row,1,&col,&vals[0][0],INSERT_VALUES);
301: } else {
302: //
303: const PetscInt row = i,col[] = {i-1,i,i+1};
304: const PetscScalar dxxL = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx;
305: const PetscScalar vals[2][3][2] = {{{dxxL,0},{a*hx+dxx0,0},{dxxR,0}},
306: {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}};
307: MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);
308: }
309: }
311: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
312: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
313: if (J != Jpre) {
314: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
315: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
316: }
317: return(0);
318: }
321: /********************************
322: * *
323: * initialize_moab_mesh *
324: * *
325: ********************************/
329: PetscErrorCode initialize_moab_mesh(moab::ParallelComm* pcomm,int npts,int nghost,moab::Tag &unknowns_tag,PetscInt &unknowns_tag_size,moab::Tag &id_tag)
330: {
331: moab::ErrorCode merr;
332: PetscInt num_procs;
333: PetscInt rank;
336: MPI_Comm_size( pcomm->comm(),&num_procs );
337: MPI_Comm_rank( pcomm->comm(),&rank );
339: // Begin with some error checking:
340: if(npts < 2) {
341: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"npts must be >= 2");
342: }
344: if(num_procs >= npts) {
345: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Num procs must be < npts");
346: }
348: // No errors,proceed with building the mesh:
349: moab::Interface *mbint = pcomm->get_moab();
350: moab::ReadUtilIface* readMeshIface;
351: mbint->query_interface(readMeshIface);
353: // Determine which elements (cells) this process owns:
354: const PetscInt nele = npts-1;
355: PetscInt my_nele = nele / num_procs; // Number of elements owned by this process
356: const PetscInt extra = nele % num_procs;
357: PetscInt vstart = rank * my_nele; // The starting element for this process
358: if(rank < extra) my_nele++;
359: vstart += std::min(extra, rank);
360: const PetscInt my_npts = my_nele + 1;
362: // Create the local portion of the mesh:
363: // Create vertexes and set the coodinate of each vertex:
364: moab::EntityHandle vertex_first;
365: std::vector<double*> vertex_coords;
366: const int sequence_size = (my_nele + 2) + 1;
367: merr = readMeshIface->get_node_coords(3,my_npts,1,vertex_first,vertex_coords,sequence_size);MBERRNM(merr);
369: const double xs = 0.0, xe = 1.0;
370: const double dx = (xe - xs) / nele;
371: for (int i = 0; i < my_npts; i++) {
372: vertex_coords[0][i] = (i+vstart)*dx;
373: vertex_coords[1][i] = vertex_coords[2][i] = 0.0;
374: }
376: moab::Range owned_vertexes;
377: merr = mbint->get_entities_by_type(0,moab::MBVERTEX,owned_vertexes);MBERRNM(merr);
379: // Create elements between mesh points. This is done so that VisIt
380: // will interpret the output as a mesh that can be plotted...
381: moab::EntityHandle edge_first;
382: moab::EntityHandle *connectivity = 0;
383: merr = readMeshIface->get_element_connect(my_nele,2,moab::MBEDGE,1,edge_first,connectivity);MBERRNM(merr);
385: for (int i = 0; i < my_nele; i+=1) {
386: connectivity[2*i] = vertex_first + i;
387: connectivity[2*i+1] = vertex_first + (i+1);
388: }
390: merr = readMeshIface->update_adjacencies(edge_first,my_nele,2,connectivity);MBERRNM(merr);
392: // Set tags on all of the vertexes...
393: // Create tag handle to represent the unknowns, u and v and
394: // initialize them. We will create a single tag whose type is an
395: // array of two doubles and whose name is "unknowns"
396: Field default_value = {0.0,0.0};
397: unknowns_tag_size = sizeof(Field)/sizeof(PetscScalar);
398: merr = mbint->tag_get_handle("unknowns",unknowns_tag_size,moab::MB_TYPE_DOUBLE,unknowns_tag,
399: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT,&default_value);MBERRNM(merr);
401: // Apply the "unknowns" tag to the vertexes with some initial value...
402: std::vector<Field> tag_data(my_npts);
403: for (int i = 0; i < my_npts; i++) {
404: tag_data[i].u = 1+sin(2*PETSC_PI*vertex_coords[0][i]);
405: tag_data[i].v = 3;
406: }
408: merr = mbint->tag_set_data(unknowns_tag,owned_vertexes,tag_data.data());MBERRNM(merr);
410: // Get the global ID tag. The global ID tag is applied to each
411: // vertex. It acts as an global identifier which MOAB uses to
412: // assemble the individual pieces of the mesh:
413: merr = mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr);
415: std::vector<int> global_ids(my_npts);
416: for (int i = 0; i < my_npts; i++) {
417: global_ids[i] = i+vstart;
418: }
420: merr = mbint->tag_set_data(id_tag,owned_vertexes,global_ids.data());MBERRNM(merr);
422: moab::Range owned_edges;
423: merr = mbint->get_entities_by_type(0,moab::MBEDGE,owned_edges);MBERRNM(merr);
424: merr = pcomm->resolve_shared_ents(0,owned_edges,1,0);MBERRNM(merr);
426: // Reassign global IDs on all entities.
427: merr = pcomm->assign_global_ids(0,1,0,false);MBERRNM(merr);
428: merr = pcomm->exchange_ghost_cells( 1,0,nghost,0,true);MBERRNM(merr);
430: // Everything is set up, now just do a tag exchange to update tags
431: // on all of the shared vertexes:
432: merr = pcomm->exchange_tags(id_tag,owned_vertexes);MBERRNM(merr);
433: return(0);
434: }
439: PetscErrorCode create_app_data(_User& user)
440: {
442: moab::ErrorCode merr;
446: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Advection-reaction options","ex30.cxx"); {
447: user.A = 1;
448: user.B = 3;
449: user.alpha = 0.02;
450: user.uleft = 1;
451: user.uright = 1;
452: user.vleft = 3;
453: user.vright = 3;
454: user.npts = 11;
455: user.io = PETSC_FALSE;
456: PetscOptionsReal("-A","Reaction rate","ex30.cxx",user.A,&user.A,NULL);
457: PetscOptionsReal("-B","Reaction rate","ex30.cxx",user.B,&user.B,NULL);
458: PetscOptionsReal("-alpha","Diffusion coefficient","ex30.cxx",user.alpha,&user.alpha,NULL);
459: PetscOptionsReal("-uleft","Dirichlet boundary condition","ex30.cxx",user.uleft,&user.uleft,NULL);
460: PetscOptionsReal("-uright","Dirichlet boundary condition","ex30.cxx",user.uright,&user.uright,NULL);
461: PetscOptionsReal("-vleft","Dirichlet boundary condition","ex30.cxx",user.vleft,&user.vleft,NULL);
462: PetscOptionsReal("-vright","Dirichlet boundary condition","ex30.cxx",user.vright,&user.vright,NULL);
463: PetscOptionsInt("-npts","Number of mesh points","ex30.cxx",user.npts,&user.npts,NULL);
464: PetscOptionsBool("-io", "Write out the solution and mesh data", "ex30.cxx", user.io, &user.io, NULL);
465: } PetscOptionsEnd();
467: user.mbint = new moab::Core;
468: user.pcomm = new moab::ParallelComm(user.mbint,PETSC_COMM_WORLD);
470: initialize_moab_mesh(user.pcomm,user.npts,1,user.unknowns_tag,user.unknowns_tag_size,user.id_tag); // creates moab mesh and field of unknowns on vertices
472: // Set the edge/vertex range once so we don't have to do it
473: // again. To do this we get all of the edges/vertexes then filter so
474: // we have only the owned entities in the ranges:
475: user.owned_edges = new moab::Range;
476: merr = user.mbint->get_entities_by_type(0,moab::MBEDGE,*user.owned_edges);MBERRNM(merr);
478: user.owned_vertexes = new moab::Range;
479: merr = user.mbint->get_entities_by_type(0,moab::MBVERTEX,*user.owned_vertexes);MBERRNM(merr);
480: user.all_vertexes = new moab::Range(*user.owned_vertexes);
482: user.shared_vertexes = new moab::Range;
483: merr = user.mbint->get_entities_by_type(0,moab::MBVERTEX,*user.shared_vertexes);MBERRNM(merr);
485: merr = user.pcomm->filter_pstatus(*user.owned_edges,PSTATUS_NOT_OWNED,PSTATUS_NOT);MBERRNM(merr);
487: merr = user.pcomm->filter_pstatus(*user.owned_vertexes,PSTATUS_NOT_OWNED,PSTATUS_NOT);MBERRNM(merr);
489: merr = user.pcomm->filter_pstatus(*user.shared_vertexes,PSTATUS_SHARED,PSTATUS_OR);MBERRNM(merr);
491: // Do some error checking...make sure that tag_data is in a single
492: // sequence:
493: int count;
494: void *data;
495: moab::Tag tag;
496: merr = user.mbint->tag_get_handle("unknowns",tag);MBERRNM(merr);
497: merr = user.mbint->tag_iterate(tag,user.all_vertexes->begin(),user.all_vertexes->end(),
498: count,data);MBERRNM(merr);
499: if((unsigned int) count != user.all_vertexes->size()) {
500: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Tag data not laid out contiguously %i %i",
501: count,user.all_vertexes->size());
502: }
503: return(0);
504: }
509: PetscErrorCode create_matrix(_User &user, DM &dm, Mat *J)
510: {
512: moab::ErrorCode merr;
513: moab::Tag ltog_tag;
514: const moab::Range *range;
517: DMMoabGetLocalToGlobalTag(dm,<og_tag);
518: DMMoabGetLocalVertices(dm,&range,NULL);
520: // Create the matrix:
521: MatCreateBAIJ(user.pcomm->comm(),2,2*range->size(),2*range->size(),
522: PETSC_DECIDE,PETSC_DECIDE,3, NULL, 3, NULL, J);
524: // Set local to global numbering using the ltog_tag:
525: ISLocalToGlobalMapping ltog;
526: PetscInt *gindices = new PetscInt[range->size()];
527: merr = user.pcomm->get_moab()->tag_get_data(ltog_tag, *range, gindices);MBERRNM(merr);
529: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 2, range->size(), gindices,PETSC_COPY_VALUES, <og);
530: MatSetLocalToGlobalMapping(*J,ltog,ltog);
532: // Clean up:
533: ISLocalToGlobalMappingDestroy(<og);
534: delete [] gindices;
535: return(0);
536: }
540: PetscErrorCode destroy_app_data(_User& user)
541: {
544: delete user.owned_edges;
545: delete user.owned_vertexes;
546: delete user.all_vertexes;
547: delete user.shared_vertexes;
548: delete user.pcomm;
549: delete user.mbint;
551: return(0);
552: }