Towards High-Performance
Unstructured-MeshComputations

Xing Cai
Simula Research Laboratory&University of Oslo, Norway

November 29, 2023

Background

Modern processors are fast with number crunching, but the memory speeds
lag behind: The “memory wall” problem
Most scientific computing applications are memory-traffic bound
Unstructured meshes are needed for modeling realistic problems, butirregular memory accesses are inevitable
How to alleviate the challenges of achieving performance for unstructuredmesh computations?

This talk presents some recent research activities at Simula

Re-ordering of mesh entities
Automated code generation for GPU computing
Re-purposing an ML-specific processor for unstructured mesh computation
Physics-guided mesh partitioning
Detailed modeling of heterogeneous point-to-point MPI communication
Joint work with J. Trotter, A. Thune, L. Burchard, K. Hustad, J. Langguth, S. Funke, A. Rustad,
S.-A. Reinemo, T. Skeie

Re-ordering of mesh entities may improve memory
performance

Numerical discretization (FEM, FVM) over unstructured meshes will
inevitably lead to irregular memory accesses
Proper re-ordering of the mesh entities may improve data reuse in thecaches, thus reducing memory trafficRe-ordering example 1: sparse matrix-vector multiplication

y=Ax#pragma omp parallel forfor (int i=0; i<num_rows; i++) {double tmp = 0.0;
for (int j=irows[i]; j<irows[i+1]; j++)tmp += A_values[j] * x[jcols[j]];y[i] = tmp;}

Matrix A is sparse, stored in the CSR format (irows, jcols, A_values)
Several re-ordering strategies may improve memory performance of SpMV
•Reverse Cuthill-McKee re-ordering•

Graph partitioning-based re-ordering•

Nested dissection•

OthersRe-ordering example 1: sparse matrix-vector multiplication

Before re-ordering

After re-ordering

A new publication at SC23

Re-ordering example 2: finite element assembly procedure
1. Gather cell coordinatesand coefficients

2. Compute element
vector and matrix

3. Update global
vector and matrix

1. Reverse Cuthill-McKee re-ordering of themesh nodes
2. Re-ordering of the mesh cells in ascending
lexicographic order according to their nodeindicesRe-ordering example 2: finite element assembly procedure

Automated code generation can simplify GPU programming

GPU programming for numerical computations can be challenging
Automated code generation can alleviate the programming challenge
Combination of high-level domain-specific language and special compiler
The FEniCS framework has used automated code generation to deploy
finite element computation to CPU clusters in a user friendly way
We recently enabled automated GPU computing for FEniCS•

GPU offloading of the finite element assembly procedure•

Seamless coupling with GPU-capable linear algebra backends

Automated generation of GPU-accelerated assembly code
1. Mathematical equations! "∇u ⋅∇v dx = ! fv dx
2. High-level user code (Python)Form language (UFL)

ParserCode generator

3. Form compiler

GPU finite element code

123456789101112elem = FiniteElement("Lagrange",tetrahedron,1)
coords = VectorElement("Lagrange",tetrahedron,1)
mesh = Mesh(coords)
V = FunctionSpace(mesh, elem)u = TrialFunction(V)v = TestFunction(V)f = Coefficient(V)kappa = Constant(mesh)a = kappa * inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx

Automated GPU acceleration of finite element assemblyBest CPU performance

Best GPU performance

Automated GPU acceleration of entire FEM computation

Re-purposing Graphcore intelligence processing units (IPUs)• Massively parallel architecture1472 tiny cores per chip
Each core has its private SRAM
No chip-level shared memory• Special programming style

The computation needs to be formulated as a“dataflow” graphNodes: small computational tasks for the cores
Edges: flow of data between the nodesInter-core communication is “implied”• IPUs: originally designed for ML
Simula ported several scientific computing
applications to IPUs:
heart simulation, graph analytics,
sequence alignment, etc.

A simple model of cardiac electrophysiology• The monodomain model of cardiac electrophysiology• Operator splitting results in a “PDE” part and an “ODE part
• PDE part: a 3D diffusion equation• ODE part: a system of nonlinear ODEs at every mesh entity• Also subject of the JHPCN project between Simula & U. Tokyo

Porting a simple cardiac simulator to Graph IPUs
In comparison with using Nvidia’s A100 GPUs
• The PDE part of the cardiac simulator runs
faster on IPUs
• The ODE part is slower on IPUs

Physics-guided mesh partitioning
Mesh partitioning is the first step of parallel computing• Partitioning an unstructured
computational mesh is non-trivial• Mesh partitioning affects theparallelization overhead, may also
impact the numerical performance• When mesh entities haveheterogeneous connectivity
strength, the partitioning problem
needs special care

Example: Parallel solution of the “black-oil” reservoir model• Large-scale reservoir simulations require parallelization• Parallel preconditioners are essential for the iterative linear solvers, but
their effectiveness is sensitive to mesh partitioning• Mesh partitioning must therefore balance between parallel efficiency
and numerical efficiency

Transmissibility as a measure of numerical connectivity• Cell-centered finite volumediscretization• Transmissibility across theboundary of two computational
cells is a good measure of numericalconnectivity• Tightly connected cells should
ideally not be divided between two
subdomains

Three graph partitioners used for mesh partitioning!!

!#!#! !##• The computational mesh is first!
translated to a graph, each cell!
becomes a vertex• If two cells share a face, the twocorresponding vertices areconnected by an edge in the graph
1. Each edge has uniform weight
2. Each edge is weighted by thetransmissibility itself
3. Each edge is weighted by thelogarithmic of transmissibility

#"

$!&!)
!' !(!#%!#*!"
!%

!!

"⊂$!"

Comparing the edge-weight schemes using 64 MPI processesComparing the edge-weight schemes using 64 MPI processesComparing the edge-weight schemes using 64 MPI processes

Detailed modeling of communication overhead• The inter-connect on a parallel system
is often heterogeneous• The actual process-to-process
communication is also heterogeneous• Detailed understanding of thecommunication overhead is important• State-of-the-art models haveweaknesses• We have developed new models

The state-of-the-art models have several weaknesses
Postal model:! "&= $ + '(The max-rate model [Gropp et al. 2016]:! ", *+,&=$+-./(+,'(12 ,'(345 )Our new model is based on a “staircase” principle

Staircase model&'(!!"#$% = )* !(&) ,"#"#$%!!"#$% = !!&'+(&! )(+! &+!"#),-.$%((&!)!! = # + max(*!"#$% , *!&#'( )N!"#$1!"#$ (1)2!"#$ (2)3!"#$ (3)4!"#$ (4)5!"#$ (5)Our new model is based on a “staircase” principle

Staircase model&'(!!"#$% = )* !(&) ,"#"#$%!!"#$% = !!&'+(&! )(+! &+!"#),-.$%((&!)!! = # + max(*!"#$% , *!&#'( )N!"#$1!"#$ (1)2!"#$ (2)3!"#$ (3)4!"#$ (4)5!"#$ (5)Our new model is based on a “staircase” principle

Staircase model&'(!!"#$% = )* !(&) ,"#"#$%!!"#$% = !!&'+(&! )(+! &+!"#),-.$%((&!)!! = # + max(*!"#$% , *!&#'( )N!"#$

1!"#$ (1)2!"#$ (2)3!"#$ (3)4!"#$ (4)5!"#$ (5)Our new model is based on a “staircase” principle

Staircase model&'(!!"#$% = )* !(&) ,"#"#$%!!"#$% = !!&'+(&! )(+! &+!"#),-.$%((&!)!! = # + max(*!"#$% , *!&#'( )0123,67 (0))*+,-. = 45

N!"#$1!"#$ (1)2!"#$ (2)3!"#$ (3)4!"#$ (4)5%&'( (5)Our new model is based on a “staircase” principle

Staircase model&'(!!"#$% = )* !(&) ,"#"#$%!!"#$% = !!&'+(&! )(+! &+!"#),-.$%((&!)!! = # + max(*!"#$% , *!&#'( )=247,;< (=))0+,-. = 9:)*+,-. = )0+,-. +12(45647)
9:;<(1),N!"#$1!"#$ (1)2!"#$ (2)3!"#$ (3)4%&'( (4)5!"#$ (5)Our new model is based on a “staircase” principle

Staircase model&'(!!"#$% = )* !(&) ,"#"#$%!!"#$% = !!&'+(&! )(+! &+!"#),-.$%((&!)!! = # + max(*!"#$% , *!&#'( )@24?
,;< (@))=+,-. = 9:)0+,-. = )=+,-. +>2(4764?)
9:;<(>)
12(4 64 ),7)*+,-. = )0+,-. + 9: 5 (1),;<N!"#$1!"#$ (1)2!"#$ (2)3%&'( (3)4!"#$ (4)5!"#$ (5)Our new model is based on a “staircase” principleStaircase model&'(!!"#$% = )* !(&) ,"#"#$%!!"#$% = !!&'+(&! )(+! &+!"#),-.$%((&!)!! = # + max(*!"#$% , *!&#'( )A13@,:; (A))>+,-. = 89)<+,-. = )>+,-. +?1(3=53@)
89:;(?),*1(3 53 )=)0+,-. = )<+,-. + 89 6 (*),:;)*+,-. = )0+,-. +01(34536)
89:;(0),N!"#$1!"#$ (1)2%&'( (2)3!"#$ (3)4!"#$ (4)5!"#$ (5)Our new model is based on a “staircase” principle

Staircase model&'(!!"#$% = )* !(&) ,"#"#$%!!"#$% = !!&'+(&! )(+! &+!"#),-.$%((&!)!! = # + max(*!"#$% , *!&#'( )B=2A,9: (B))@+,-. = 78);+,-. = )@+,-. +*=(2?42A)
789:(*),0=(2 42 )?
)<+,-. = );+,-. + 78 > (0),9:)0+,-. = )<+,-. +<=(2542>)
789:(<),(2 425),9: (;))*+,-. = )0+,-. + 783

N!"#$

1%&'( (1)2!"#$ (2)3!"#$ (3)4!"#$ (4)5!"#$ (5)The Staircase model works in more general cases when messages are
mixed intra-socket, inter-socket, and inter-node

Mixed intra-node bandwidth estimate:&*+&())#$%(,!"!"($, &$ ) = '! !"!" ($) + ' ! !"!"($),

General case estimate:!"#$#%& = !""(#)%*($+, + !""(#,)*($+,&$ =-!"#"$$-!"# .-!

Example of detailed modeling of communication overhead

Concluding remarks
New research needed for mesh reordering & partitioning
Further developments of automated code
generation for accelerated computing
Important international collaborations• JHPCN project with U. Tokyo• SparCity project (EuroHPC)