# Presentation: 603061

## PamgenA Parallel Finite-Element Mesh Generation Library

**TUG 2008**
**Monday, October 21, 2008**
**David Hensinger (SNL)**

- Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company,for the United States Department of Energy under contract DE-AC04-94AL85000.

## Parallel Mesh GenerationLibrary Driver & Challenges

**Analysis codes and capability machines are leaving pre-processing tools and shared memory machines behind**
**Serial mesh generation strategies were unable to supply analysts’ demands for:**
**Billions of Elements**
**Quick Turn-around**

**Scale well: multiple thousands of processors and billion + elements**
**Exploit determinism**
**Eschew communication**
**Ease of Use: run different decompositions without modifying input deck**
**Automatic load balancing**
**Consistent topology/geometry**

- Purple, 1532 nodes

- R.S. 12960 Nodes

**Notes:**
This is largely boiler plate.  
The driver was an unstructured ALE code (ALEGRA). 
The original mesh generation capability served ALEGRA’s multi-block structured mesh capability. This proved to be so popular that it was requested on the unstructured side of ALEGRA as well. This adaption to generation of unstructured mesh has subsequently outlived the multi-block structured capability.

Unstructured meshes are frequently constructed interactively as a serial process.
For certain classes of problems the mesh could be tersely described and easily generated. And yet serial machines balked at producing multiple million element meshes.

Goal is to fill capability machines and scale well doing so. 
This will not solve the general meshing problem but will allow problems with nearly trivial geometries to be run at huge sizes.

## Capabilities

**Topologies**
**Cubes**
**Solid Cylinders – Full and partial**
**Hollow Cylinders – Full and partial**
**Geometries **
**Those suggested above plus projection to sphere**
**Arbitrary user-specified geometry transformation**
**Boundary conditions – Call out nodes and element faces on topological faces,edges,corners**
**Decompositions**
**Optimal bisection**
**User controlled**
**Sequential**
**Random**

**Notes:**
This is an overview together with a graphic showing the mesh topology zoo.  Through this discussion it is useful to remember that the topology specifies how elements connect to each other. The topologies are fairly limited. 

The geometries are limited by the difficulty in specifying complex shapes in a text string passed to a parser. However the arbitrary geometry transformation capability takes a user supplied “C” function and applies that function to each node’s coordinates. The original coordinates are available during the calculation. An example of this follows.

Regions for boundary condition application are spelled out in the mesh description passed to the library. An example follows.

The decomposition the mesh generation step that is most limited by the lack of communications. The most useful decomposition is the default. This optimizes the number of slices made through the mesh in each topological direction to produce the required number subdomains. This strategy is most useful wen the number of processors is the product of multiple prime numbers. The next most useful is a user-controlled option that allows the user to specify the number of processors in each topological direction of the mesh.  The product of the processors specified must equal the number of processors used. The sequential distribution deals splits the elements across the processors equal in order by global id. This is a poor decomposition for large numbers of processors. The random decomposition distributres the elements randomly. It is good for testing communication subsystems in unstructured codes.

## Enabling/Limiting Assumptions

**All Processors Have Identical Information**
**No need to communicate**
**All Processors are Identical (except for Ids)**
**Calculations produce identical results**
**No need to communicate**
**Communication is Unavailable**
**Resist the temptation**
**Scale really well**
**All Meshes Will Consist of One or More Structured Blocks**
**Connectivity easily calculated**
**Geometry derives from topologies**

**Notes:**
These assumptions overlap to some degree. Forbidding communication would require all of the others except for the limitation on topology.
The limitation of topology simplifies the calculations and specifications required; but more advanced topologies could be allowed with the addition of substantial complexity.
Forbidding communication significantly helps scaling and development.

## Execution Stages

**Information Distribution – Each processor gets a complete description of the mesh.**
**Decomposition - Run equivalently on each processor such that each processor:**
**Receives a list of its local elements**
**Can calculate the processor of any element**
**“****Serial” Information Generation – Nodes, Elements, Connectivity local to a processor**
**“****Parallel” Information Generation – Inter-processor communication information.**
**Geometric Transformations – User provided subroutines to calculate new nodal coordinates**

**Notes:**
The execution stages are similar are the same as those that might be followed by a serial mesh generation, but the decomposition is moved to the beginning of the mesh generation process instead of being relegated to a post processing step.

Requiring each processor to perform the decomposition for the entire mesh while still maintaining scalability forces reliance on basic strategies. Since the topologies and geometries are also simple the decomposition strategies are sufficient. More about these on a later slide. The default decomposition is an optimized slicing strategy that exploits knowledge of the mesh specification.

Generating “serial” information is equivalent to generating the entire mesh and then throwing away the non-local parts. The actual generation works from the other direction and exploits knowledge of topology to generate local elements and build up node lists then compact them and re-assign local ids.....

The “parallel” information consists of synchronized boundary node lists corresponding to adjacent processors. The lists are generated by visiting elements, calculating their neighboring elements, determining the processor ownership of the neighboring element, and adding the connecting nodes to the appropriate list. The lists are subsequently sorted by a global index.

## Library Interface

**Create a “mesh” within the library**
**	****	****int Create_Pamgen_Mesh(char * ****	****	****mesh_description, int dimension, int ****	****rank, int num_procs);**
**Query the library to build up representation in client code – the same as may be done with a file interface API**
**	****	****Im_ex_get_info(...)**
**	****	****Im_ex_get_connectivit(...)**
**After queries are completed the library memory can be cleared with a delete function.**

**Notes:**
The construction occurs with a single function call that takes a “C” string that is then parsed within the library. The string is a terse description of the mesh. The construction of the strings is covered in the SAND report that documents the library.

The API used to query the library is a direct descendent of the NEMESIS and EXODUSII APIs. These  API’s are used to query data files and load mesh information into application memory. The meshing library is meant as a partial replacement for an external mesh file and strives to emulate the existing interfaces.  The API functions used by Pamgen are generally the same as those used by EXODUS II and NEMESIS but are prefixed by im_.

Within the ALEGRA code the replacement of file access for mesh loading with library memory access was straightforward. One significant difference is that no disk-staging is required for the mesh library access. There is a downside that the mesh memory is duplicated in both the lilbrary and client code, but he library memory can be subsequently freed. 

## ExampleA Brick mesh

**mesh**
**  ****brick**
**  ****numz 2**
**    ****zblock 1 2. interval 5**
**    ****zblock 2 8. interval 4**
**  ****numx 2**
**    ****xblock 1 5.0 interval 5**
**    ****xblock 2 5.0 interval 5**
**  ****numy 2**
**    ****yblock 1 10. first size 1. last size .1**
**    ****yblock 2 10. first size .1 last size 1.**
**  ****end**
**end**

## ExampleA Partial Cylinder with Node Sets and Side Sets

**mesh**
**  ****radial trisection**
**    ****trisection blocks, 2**
**    ****zmin -0.00075**
**    ****numz 1**
**      ****zblock 1 1. interval 4**
**    ****numr 3**
**      ****rblock 1 2.0 interval 4**
**      ****rblock 2 3.0 interval 4**
**      ****rblock 3 4.0 interval 4**
**    ****numa 1**
**      ****ablock 1 90. interval 12**
**    ****end**
**  ****set assign**
**    ****nodeset, ilo, 100 **
**    ****block sideset, ilo, 35, 2**
**    ****block sideset, ihi, 45, 2**
**  ****end**
**end**

## Example A Brick With Random Decomposition for 2 Processors

**mesh**
**  ****brick**
**  ****numz 1**
**    ****zblock 1 2.0 interval 15**
**  ****numx 1**
**    ****xblock 1 2.0 interval 15**
**numy 1**
**    ****yblock 1 2.0 interval 15**
**  ****end**
**  ****decomposition stratgy**
**    ****random**
**  ****end**
**end**

## ExampleA Portion of a Solid Cylinder

**mesh**
**  ****radial trisection**
**  ****trisection blocks, 2**
**  ****zmin -0.00075**
**  ****numz 1**
**    ****zblock 1 1. interval 4**
**  ****numr 3**
**    ****rblock 1 2.0 interval 4**
**    ****rblock 2 3.0 interval 4**
**    ****rblock 3 4.0 interval 4**
**  ****numa 1**
**    ****ablock 1 90. interval 12**
**  ****end**
**end**

## ExampleA 2D Block of Mesh with Geometry Transformation

**mesh**
**  ****rectilinear**
**    ****nx = 10**
**    ****ny = 10**
**    ****bx = 3**
**    ****by = 3**
**    ****gmin = -1.0 -1.0**
**   ****gmax = 1.0 1.0**
**  ****end**
**  ****user defined geometry transformation**
**  ****"**
**    ****double r = sqrt(inxcoord*inxcoord +inycoord*inycoord);**
**    ****double theta = atan2(inycoord,inxcoord);**
**    ****if(r > 0.5)**
**    ****{**
**    ****theta = theta + (3.14159 / 4.0)*((r-    0.5)/0.5);**
**    ****outxcoord = r*cos(theta);**
**    ****outycoord = r*sin(theta);**
**   ****}**
**  ****"**
**  ****end**
**end**

## ExampleA 3D Solid Cylinder Decomposed for Eight Processors

**mesh**
**  ****radial trisection**
**    ****trisection blocks, 4**
**    ****numz 1**
**      ****zblock 1 4.0 interval 1**
**    ****numr 3**
**      ****rblock 1 2. interval 4**
**      ****rblock 2 3. interval 4**
**      ****rblock 3 5. interval 4**
**    ****numa 1**
**      ****ablock 1 360. interval 32**
**  ****end**
**  ****decomposition strategy**
**    ****numprocs j, 8**
**  ****end**
**end**

## Availability and Distribution

**Available under GNU Lesser General Public License (LGPL).**
**Distributed as a component of Trilinos**
**Documented in report: SAND 2008-1933**
**http://trilinos.sandia.gov/packages/pamgen**