Internal:IMR 2018 Short Course

Coarse-Grained

Unique cavity-based operator and hierarchical domain partitioning for fast parallel generation of anisotropic meshes,
Loseille A. & Alauzet F & Menier V.

Computer-Aided Design 2017, [1]

Abstract

"Our strategy is based on typical domain splitting algorithm where the initial mesh is split in parts that are then meshed in parallel while the fictitious boundaries between parts are kept unchanged. Then we iterate the procedure to adapt previously unmodified parts of the domain, i.e., the interface mesh. Both the volume and the surface meshes are adapted simultaneously and the efficiency of the method is independent of the complexity of the geometry.

The originality of the method relies on

1. metric-based static load-balancing,
2. hierarchical mesh partitioning techniques to (re)split the (complex) interfaces meshes,
3. a fast, robust and generic sequential cavity-based mesh modification kernel.

In this paper, we follow the coarse-grained parallelization in an adaptive context within the metric-based framework."

Sequential Mesh Generation (essentially is a AF-based method)

Loop the steps :

1. Creation of unit mesh in metric space:
   1. Collapse small edges
   2. Point Creation
   3. Anisotropic Filtering
   4. Point Insertion
2. Mesh Optimization
   1. Edge/face flips using the cavity operator
      
      
3. Surface Approximation base on surrogate geometry ( it is not clear is this is on every loop or at the end)

Hierarchical Partitioning

To enable coarse-grain parallelism a hierarchical partitioning scheme is used.

"At the first level, the initial domain is split and considered for adaptation. From the set of constrained faces (at interface) of the previous level, a new volume mesh (domain) is deduced for the next level. This new domain is then split and adapted in parallel. This process is then applied until convergence."

The new volume is is evaluated using interface meshes (i.e. union of cavities of interface elements )

Domain Decomposition

At each step of the hierarchical partitioning domain decomposition is used to partition the current mesh.

How load is evaluated

"To handle non uniform refinements (in term of sizes and directions), a metric-based static load balancing formula is used to a priori equilibrate the work on each sub-domain."

The unique cavity framework offers constant speed whatever the considered operator.

"This feature allows us to derive an accurate metric-based work that is easily deduced from the input metric-field and input mesh only."

Mesh Partitioning

First Hilbert ordering of vertices is used.

Then BFS with restart and the work load as weight is used to partition the elements.

Load Balancing

Static since previous steps tried to make partition as well balanced as possible.

Performance Data

Start-to-End : speed-up of 20-33 on 120 cores

Parallel anisotropic 3D mesh adaptation by mesh modification,
Frédéric Alauzet & Xiangrong Li & E. Seegyoung Seol & Mark S. Shephard

Engineering with Computers (2006) 21: 247–258, DOI 10.1007/s00366-005-0009-3 [2]

Algorithm

1. Mesh coarsening to eliminate the initial short edges
2. Intelligent mesh refinement that includes operations
   for ensuring the proper geometric approximation of
   the mesh to the geometric domain and proper shape
   to match the metric field, and
3. Mesh optimization to improve the quality of the
   final mesh

" The Flexible distributed Mesh DataBase (FMDB) has been used to support the needed mesh-based parallel
operations on distributed meshes including communications between mesh partitions, migrating mesh entities between the partitions and the dynamic mesh load balancing. "

"The mesh load balancing is also critical for the parallel procedure to achieve desirable performance. The Zoltan library is used to make partition assignment during the dynamic load balancing."

Coarsening

Mesh migration-based parallelization is used:

"More precisely, whenever a mesh modification operation on partition boundaries is desired, we migrate the mesh entities within the polyhedron associated with (i.e., the mesh regions affected by) the operation into one partition and since then we apply the operation as in serial."

The above step is repeated until all coarsening operations are applied.

Process id-based scheme is use to resolve conflicting requests.

Mesh Refinement

"Mesh refinement relies on tagging mesh edges and
applying edge, face and region subdivision templates,
which changes mesh topology but not geometry. "

Steps

1. "Marked mesh edges and their neighboring mesh
   faces on partition boundaries are subdivided using the
   four templates."
2. Boundary Modifications Synchronization
3. Refinement of volume with no communication

Parallel Vertex Snapping (on curved boundary)

Use direct snapping or collapse operations , if vertex is on partition boundary communicate to ensure validity and/or collect cavity.

Load Balancing

Using Zoltan redistribute the global distributed mesh in (FMDB)

Performance

No time/speedup data are available besides a phrase "In the example, the speedup of the parallel procedure is between 1.2 and 2" pg 255

Parallel Optimization of Tetrahedral Meshes,
D. Benítez, E. Rodríguez, J.M. Escobar, R. Montenegro

6th European Conference on Computational Mechanics (ECCM 6) & 7th European Con-ference on Computational Fluid Dynamics (ECFD 7), 2018 [3]

Summary

A new algorithm on distributed-memory parallel computers for our simultaneous untangling and smoothing of tetrahedral meshes .

Overview

The method is based on

• partitioning a mesh
• optimizing interior vertices
• optimizing boundary vertices of interior partitions, and
• communicating updated coordinates of boundary vertices.

In each mesh sweep, the interior vertices of all partitions are optimized in parallel. For partition boundaries, independent sets of non-fixed vertices are created; each of them is optimized in parallel after all interior vertices have been optimized

Distributed Combinatorial Maps for Parallel Mesh Processing,
Guillaume Damiand, Aldo Gonzalez-Lorenzo, Florence Zara and Florent Dupont

Algorithms 2018, 11(7), 105; doi:10.3390/a11070105 [4]

Summary

"Our main contribution is the definition of distributed combinatorial maps (called n-dmaps), which allow us to
represent the topology of big meshes by splitting them into independent parts. "

"To sum up, this topological representation can be used to describe an nD mesh, while providing
high level operations to iterate through the cells of the mesh, retrieve the incidence and adjacency
relations between the cells, and modify the mesh by using the numerous existing operations. "

Code is part of CGAL

https://doc.cgal.org/latest/Combinatorial_map/index.html#Chapter_Combinatorial_Maps

Application

Adaptive Hexahedral Subdivision (ie octree refinement)

• MPI-based
• Uses a single template (1:8) for splitting the octree
• Shared faces and edges between processes are split independently and asyncronous messages are used to achieve conformity. ( Like PCDM 2D)

Performance Data

Speedup 6.3 on 32 cores

Parallel Anisotropic Tetrahedral Adaptation,
Michael Park & David Darmofal

46th AIAA Aerospace Sciences Meeting and Exhibit, 2008 [5]

Summary

• Coase-grained
• MPI

Mesh Operations

• Node Movement
• Face/Edge Swapping
• Tetrahedra Split and Collapse

Parallelization

Domain Decomposition

• Partition Mesh by assigning  vertices uniquely to each partition.
• Add also ghosted vertices to each node to complete the boundary tetrahedra.

Parallel Execution

• Operations are performed in parallel first, on all purely local nodes.
• Then on boundary by colouring the subdomains and operating on one color each time.
• Changes on boundary are recorded on a transcript that is broadcasted to the neighbors at the end of each pass in order to replicate steps.

Performance Data

No performance Data provided

Parallel anisotropic mesh adaptation with boundary layers for automated viscous flow simulations ,
Onkar Sahni, Aleksandr Ovcharenko, Kedar C. Chitale, Kenneth E. Jansen, Mark S. Shephard

(Engineering with Computers 2017), [6]

Summary ( of parallel code)

• Uses the same template refinement like other papers of this group
• Distributed memory model
• Based on a distributed mesh database

"A partitioned mesh is managed by the parallel mesh database that tracks the mesh  entities residing on inter-part boundaries."

" It supports the ability to update mesh entities on inter-part boundaries if they are modified (e.g., due to an edge split). It also supports the movement or migration of mesh entities from one part to another, which is referred to as a mesh migration step and in which the inter-part boundaries are automatically updated"

"Thus, before applying the mesh modification operation, such cavities[between subdomains] are localized on one part or processor such that the cavity retriangulation can be carried out as in the serial case. In cavity localization, all regions and stacks
of regions involved in the mesh modification operation are migrated onto a single part ."

 

Algorithm

Refinement

1. Tag Long edges/faces for refinement
2. For each face /edge
   1. Split  face/edge
   2. If on shared boundary copy data to buffer
3. end for
4. send.receive modifications

After communication new elements point to remote neighbors correctly. i.e. the distributed mesh is properly updated.

Coarsening and Swapping

1. Tag vertices/elements that need modification
2. While there is more work
   1. For each vertex/element
   2. If the whole cavity is present perform operation
   3. Else add missing parts to a buffer
   4. end for
   5. Send/receive buffers
3. end while

Special rules are used for the boundary layers. Main idea is that all the elements of the layer need to be local if one is to be split, since splits propagate from upper to lower layers. Similar with collapse and/or swap.

 

Performance Data

Strong scaling

Using 512 cores as base case  there is 75% efficiency on 1024 cores and 30% efficiency on 8192 cores

Weak scaling

Using 256 cores as base case  there is 97% efficiency on 1536 cores and 86% efficiency on 9216 cores

Note : in weak efficiency evaluations uses actual work increase instead of core increase 

Parallel tetrahedral mesh adaptation with dynamic load balancing,
Leonid Oliker, Rupak Biswas, Harold N Gabow

( Parallel Computing 2000), [7]

Summary

• Main problem tackled : How to load balance mesh adaptation which by the physics on the problem requires adaptation only in small regions which may be in their entirety in a subdomain
• Edge-based adaptation
• MPI
• Coarse-grained
• Utilizes PLUM see relevant paper

Adaptation Procedure

In every iteration:

1. A solution is acquired
2. Mesh adaptation based on error of the solver
3. The mesh is (re)partitioned using a collection of metrics in PLUM
4. back to 1 until convergence

Tetrahedral Adaptation

Uses three templates (1:8,1:4,1:2)

1. Mark (in parallel) all edges until to edge can be split further
2. Communicate changes with subdomain neighbors
3. Split edges in parallel with no communication

Performance Data

Efficiency (strong scaling ?)

76% on 64 processors with load balancing off

82% on 64 processors with load balancing on

Fine-Grained

Mesh adaptation for moving objects on shared memory hardware,
Dan Ibanez & Mark Shephard

Research Note on 25th IMR , 2016 (Suggested by Dan As reference paper for Omega_h) [8]

Summary

Parallel shared memory adaptation on GPUs and threads.

Single code-base for threads and GPUs using Kokkos.

Mesh operators

1. Edge split (at midpoint)
2. Edge Collapse
3. Edge Swap (using n -> (2n-4) flips up to n = 7  )

Independent set selection

Use a modified version of Luby's algorithm to select a maximal (vs maximum) independent set of cavities . The cavities can be then performed in parallel with no communication/synchronization. This procedure is used throughout the mesh adaptation process.

Adaptation Steps

1. Size Optimization
   1. Refinement :
      1. Edges longer than desired threshold are marked.
      2. From the previous set if a split will produce low quality tetrahedra then the edge is unmarked.
      3. An independent set is selected and they are split in parallel.
   2. Coarsening :
      1. Vertices attached to edges shorter than a desired threshold are marked.
      2. For each vertex the best quality collapse is selected.
      3. An independent set is selected and splits are performed in parallel.
2. Quality Optimization
   1. Edge Swapping :
      1. Low quality elements are marked together with four "layers" (neighbors of neighbors 4 level deep)  of elements around them.
      2. Mark all edges of all elements
      3. For each edge calculate best flip and keep only those than improve quality
      4. An independent set is selected and flips are performed in parallel.
   2. Coarsening for Shape Correction
      1. Vertices of low quality elements are marked
      2. Continue similar of steps 2-3 of coarsening above

It is interesting to note that mesh is readonly and thus is rebuild (in parallel) after every each pass.

Performance Data

Speedup of 49 on 64 threads (Intel Knight Landing)

Speedup of 3.6 on 8 threads (Intel Xeon 2620 v4)

Speedup of 7  using a GPU ( NVidia GTX 980 Ti)

Thread-Parallel Anisotropic Mesh Adaptation,
Gerard J. Gorman, Georgios Rokos, James Southern, Paul H. J. Kelly

New Challenges in Grid Generation and Adaptivity for Scientific Computing, (2015) [9]

(this is the first reference paper for pragmatic)

Summary

• Shared Memory , OpenMP
• Features independent tasks through graph coloring and deferred updates
• Template refinement

It is interesting to note that point insertion is done in two steps first insert all midpoints on edges and then in a separate parallel loop "refine" the element ie connect new point with old vertices of the element.

Colouring

Colouring is used to find independent cavities.

The parallel coloring algorithm for maximal sets by Gebremedhin is used.

Parallel Updates of data structures

Data structures are arrays.

There are two types of updates

1. Deferred updates :
   Modifications on the arrays are recorded on the temporary buffer of a thread ( tid = i % nthreads, where i the entry to update) and applied in parallel at the end of each pass. see Algorithm 3 pg 15
2. Worklist updates
   Used for adding new elements in the arrays. Implemented using a global index of last used position that is increased atomically

 

Parallel Refinement

1. Split edges (but not element) in parallel.
2. Copy new vertices to global buffer
3. Barrier
4. Update point-adjacency in parallel
5. Barrier
6. Refine elements in parallel
7. Barrier
8. Copy new element to global list in parallel

Parallel Coarsening

1. Colour vertices
2. Identify those that can be collapsed in parallel exit if none
3. Evaluate a maximal independent set of the previous set
4. Collapse edges using deferred operations in parallel
5. Commit the deferred operations
6. Repair colouring
7. Go to 2

Parallel Swapping

Uses a modified version of the parallel Coarsening

Parallel Smoothing

1. Color vertices
2. Smooth one color at a time in parallel

Performance Data

Between 5-9 speedup (strong scaling) on 16 cores for the different phase of adaptation

A Two-Level Multithreaded Delaunay Kernel,
Jean-François Remacle, Vincent Bertrand,Christophe Geuzaine

24th International Meshing Roundtable (2015), [10]

Summary

• It is only Delaunay Triangulation
• OpenMP

Algorithm

1. Sort points in space using a Hilbert curve
2. Divide the curve into nthreads equal parts
3. Loop until there are no points to insert
   1. Calculate M2 cavities per thread
   2. Barrier
   3. For each of the M2 cavities insert the point if there is no overlap (if it is the lower thread-id advances and the rest of the points will be reinserted later)
   4. Barrier

Performance Data

Strong Scaling Speedup

~5 on 8 threads

~12 on 32 threads

Journal version (CAD 2017) of the paper reports speedup of 18 on 32 threads

One machine, one minute, three billion tetrahedra,
Célestin Marot & Jeanne Pellerin & Jean-François Remacle

preprint [11]

Summary (only of the parallel part)

• Shared memory (OpenMP)
• Delaunay-based
• Create parititions using space-filling curve to reduce chance of overlap (and thus rollbacks)

Triangulation Algorithm

• Points to be inserted are partitioned into nthreads partitions using a space-filling curve
• Mesh is partitioned into nthread partitions  using space-filling curves
• If cavity or mesh walk crosses partition it is aborted
• Re-partition by shifting and rotating the curve and try again to insert rejected points
• Decrease the number of threads if success rate is low

Parallel Mesh Refinement Algorithm

1. Boundary points insertion (in parallel)
2. Boundary Recovery
3. While tets need to be refined
   1. Evaluate circumcenters of all bad tets (in parallel)
   2. Filter out very close points (in parallel)
   3. Insert points using the triangulation algorithm

Performance Data

Mesh refinement

Small case 123 million

(Strong scaling) speedup 1.9 on 8 threads

Big case 600 million tets

(Strong scaling) speedup 4.5 on 127 threads

A Parallel Local Reconnection Approach for Tetrahedral Mesh Improvement,
Mengmeng Shang, Chaoyan Zhu, Jianjun Chen, Zhoufang Xiao, Yao Zheng

(IMR 2016), [12]

Summary

• OpenMP code
• Uses Hilbert curve to decrease chance of rollback in flips
• Uses coloring for smoothing

Parallel Flip Algorithm

1. Sort using a Hilbert curve all the bad elements
2. For each thread evaluate M cavities (elements/points are marked by vertex id)
3. BARRIER
4. For each thread
   1. For each cavity check if there is an overlap higher thread id gets priority else proceed
5. BARRIER
6. Go to 2

Parallel Smoothing Algorithm

1. Color vertices using a greedy heuristic in parallel
2. For each color
   1. smooth in parallel ( no synchronization needed)

Performance Data

Strong scaling

8.62 speedup on 16 cores for smoothing

6.9 speedup on 16 cores for Flips

Other

Thread-parallel mesh improvement using face and edge swapping and vertex insertion,
Reza Zangeneh & Carl F. Ollivier-Gooch

Computational Geometry Journal, 2018 [13]

Summary

• OpenMP
• Delaunay-based
• Parallel swap algorithm using optimistic approach
• Parallel point insertion using PDR theory with an nested adaptive octree

Parallel face/edge swap algorithm

1. Divide all faces/edges into threads queues
2. Face/edge migration for load balancing
3. While there are faces/edges in thread queues
   1. Check if flips is needed and evaluate cavity (in parallel- no communication needed)
   2. Face/edge migration for load balancing
   3. use an optimistic approach to lock cavities in each thread queue (in parallel with atomics)
   4. Face/edge migration for load balancing
   5. Perform flips in parallel with no communication

Parallel Point insertion algorithm

"We also exploit nested parallelization features in OpenMP to create separate quad/oct tree structures within each leaf
assigned to a thread. In this way, each leaf of the initial tree can be the root node of finer trees which will be assigned
to new spawned threads if more refinement is required"

There is one (static) octree A .

For each leaf of A there is a second octree B which is refined progressively.

1. Initialize octree A and assign tetrahedra to leaves
2. Assign leaves of A to threads
3. While there are unrefined tetrahedra
   1. If any thread has no leaves to refine split the leaves of the thread local octree B to create new leaves.
   2. Find independent set of leaves I among all thread local B octrees
   3. for all leaves in I
      1. spawn a thread and a task and refine it up to the current refinement level
   4. increase the refinement level

Performance Data

(Strong scaling)

Parallel Swap Speedup ~35 on 64 cores

Parallel Insertion Speedup ~18 on 64 cores

(Weak scaling)

Parallel Swap Speedup ~21 on 64 cores

Parallel Insertion Speedup ~17 on 64 cores

A parallel advancing front grid generation scheme,
Rainald Löhner

Int. J. Numer. Meth. Engng. 2001 [14]

Summary

Decompose Front using octree, advance from front "far enough". Shift octree to remesh interface boundaries.

"A parallel advancing front scheme has been developed. The domain to be gridded is first subdivided spatially using a relatively coarse octree. Boxes are then identified and gridded in parallel. A scheme that resembles closely the advancing front technique on scalar machines is recovered by only considering the boxes of the active front that generate small elements."

Grid generation

The algorithm uses many passes.

In each pass the faces of the front are divided using an octree. From faces in "far enough " octree leaves points are advanced simutaneusly

only the boxes covering the active front in regions where the smallest
new elements are being introduced are considered.

Shared Faces between leaves are eliminated shifting the octree leaves and remeshing the front.

Work balancing is achieved using overdecomposition of the octree.

Mesh Improvement

Swapping

Checks are done in parallel , actual flipping sequentially

Removal of bad elements (Edge Contraction)

Indentify edges around bad elements in parallel , contract them either by collapsing an endpoint onto the other or both onto a new point in the middle.

Laplacian Smoothing

Color edges/elements into independent sets and run in parallel. When negative volume elements are created, remove them together with adjacent elements and remesh the hole.

For surface meshing (preprocessing) each surface patch is meshed in parallel.

Results

Shared memory.

Performance more than 50% efficiency on 16 cores

Recent Advances in Parallel Advancing Front Grid Generation,
Rainald Löhner

Archives of Computational Methods in Engineering 2014,(there is also ashorter version in IMR 2013) [15]

Summary

Two Level Parallelism OpenMP + MPI

Coarse grain in both levels ? (for openmp kernel see 2001 paper)

Parallel Algorithm For mesh generation

1. Scale up  distribution function and create an initial coarse mesh called and initial domain-defining grid (DDG)
2. Decompose/load balance using a combination of greedy algorithms, recursive bisection and spacefilling curves
3. Mesh, in parallel, the zones inside the subdivision domains
4. Modify partitions by adding 1-2 extra layers of elements and redistribute active front (source subdomain id < dest subdomain id) and mesh new interior zones in parallel
5. Modify partitions by adding 1-2 extra layers of elements and redistribute active front (source subdomain id > dest subdomain id) and mesh new interior zones in parallel
6. If required, mesh the remaining regions on processor 1.
    

Mesh Improvement

Mesh improvement uses

• Diagonal Swapping
• Removal of Bad Elements
• Smoothing and
• Selective mesh movement in regions of very low quality

This operations are perform within each subdomain with no communication since modifications on the subdomain boundary are not allowed in this stage

Perfomance Data

Plenty of timing results but no comparison with sequential code.

Top speed is 403,354 elms/sec

With base case 8 cores ( 8 MPI procs) the strong scaling speedup is ~5 on 256 cores generating 121M elements

Surveys

"A decade of progress on anisotropic mesh adaptation for computational fluid dynamics" Frédéric Alauzet, Adrien Loseille, (CAD 2016), [16]

TO DO

1. "A parallel unstructured dynamic mesh adaptation algorithm for 3‐D unsteady flows"

   Young Min Park,Oh Joon Kwon, (International Journal for Numerical Methods in Fluids 2005) [17]

2. Parallel adaptation of general three-dimensional hybrid meshes Christos Kavouklis, Yannis Kallinderis [18]
3. "Scalable generation of large-scale unstructured meshes by a novel domain decomposition approach"

   Jianjun Chen, Zhoufang Xiao, Yao Zheng, Jianfeng Zou, Dawei Zhao, Yufeng Yao, (Advances in Engineering Software, 2018), [19]

4. '"A distributed-memory parallel approach for the generation of multibillion element tetrahedral meshes"

   Xiaoqing Wang, Xianlong Jin,(SNPD 2016), [20]

5. "Parallel Unstructured Mesh Adaptation Method for Moving Body Applications"

   Peter A. Cavallo; Neeraj Sinha; Gregory M. Feldman,(AIAAJ 2005), [21]

6. "A parallel log-barrier method for mesh quality improvement and untangling"

   Shankar P. Sastry, Suzanne M. Shontz, (Engineering with Computers 2014), [22]

7. "Generating Multibillion Element Unstructured Meshes on Distributed Memory Parallel Machines"

   Seren Soner and Can Ozturan, (Scientific Programming 2015), [23]

8. "A Practical Demonstration of Scalable, Parallel Mesh Generation"

   William Dawes, Simon Harvey,Simon Fellows, Neil Eccles, D Jaeggi, Will Kellar, Cambridge Flow Solutions, ( AIAA-2009-0981 ), [24]

9. "Automatic parallel generation of tetrahedral grids by using a domain decomposition approach"

   H. Andrä, O. N. Gluchshenko, E. G. Ivanov, A. N. Kudryavtsev, (Computational Mathematics and Mathematical Physics 2008), [25]

10. "Parallel remeshing of unstructured volume grids for CFD applications"

    U. Tremel, K. A. Sørensen, S. Hitzel, H. Rieger, Oubay Hassan, Nigel P. Weatherill, ( International Journal for Numerical Methods in Fluids 2007), [26]

11. "Parallel refinement and coarsening of tetrahedral meshes"

    de Cougny HL, Shephard MS ( Comput Methods Appl Mech Eng, 1999),[27]

Not used

1. "p4est: Scalable Algorithms for Parallel Adaptive Mesh Refinement on Forests of Octrees"

   Carsten Burstedde, Lucas C. Wilcox, and Omar Ghattas, (SIAM J. Sci. Comput. 2011), [28], (no use of tetrahedral refinement)

2. Löhner R, Ramamurti R (1995) A load balancing algorithm for unstructured grids. Comput Fluid Dyn 5:39–58 (no mesh generation results)
3. "Parallel adaptive mesh refinement method based on bubble-type local mesh generation"

   Wei Guo, Yufeng Nie, Weiwei Zhang, (JPDC 2018), [29] ( it is only 2D and on simple inputs)