Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019jh190022-NAJ

High performance simulations using FreeFem++ on mixed distributedplus shared-memory architecture
Atsushi Suzuki（Cybermedia Ceter, Osaka Univeristy）
Supercomputers used to perform numerical simulations nowadays typically consist of clusters of multi-coreprocessing units nodes. In such a cluster, there are two different kinds of workload distribution: through sharedmemory, e.g., using POSIX threads or OpenMP, and through distributed memory, e.g., using MPI (MessagePassing Interface). For numerical simulations of physical three-dimensional phenomena described by partial
differential equations, it is essential to solve large linear and nonlinear systems usually stored as sparse matrices.
A domain specific language, “FreeFem++” has capability to describe various multi-physics problems withvariational formulation and to generate sparse linear systems. The aim of this project is to combine the highefficiency of Dissection direct sparse solver inside one node through shared-memory parallelism and thenumerical robustness of HPDDM iterative solver over distributed nodes to tackle large-scale FEM problems.

1.

Basic Inform ation

clusters of multi-core processing units nodes. In(1) Collaborating JHPCN Centerssuch a cluster, there are two different kinds ofCybermedia Center, Osaka University

workload distribution: through shared-memory,(2) Research Arease.g., using POSIX threads or OpenMP, and☒Very large-scale numerical computation

through distributed memory, e.g., using MPI.

oVery large-scale data processing

Foro

Very large capacity network technology

three-dimensional (3D) phenomena describedo

Very large-scale information systems

by partial differential equations (PDE), it is(3) Roles of Project M em bers

numericalofphysicalessential to solve large linear and nonlinear- Atsushi Suzuki (Cybermedia Center, Osakasystems usually stored as sparse matrices,University): Implementation of Dissection

whose dimensions are within the 10 to 100

direct solver in HPDDM with GenEO

million range. In this project, we are focused on- Pierre Jolivet (CNRS, Toulouse Institutethe finite element method (FEM), which is,of Compter Science Research, France):

among other methods such as finite differences

Improvement of interface of HPDDM toor finite volumes, a widely used method forFreeFem++

discretizing PDEs. One of its advantages is that- Daisuke Furihata (Cybermedia Center,the FEM may handle irregular 3D geometries

Osaka University): Principal contact of theused as computing domains. Using a union ofresearch

arbitrarily refined tetrahedral, one may pick a
discretization

2.

simulations

PurposeandSignificanceofcomputational

Supercomputers used to perform numerical
nowadaystypically

consistaccordingtothesmoothness of the solution and the availableResearch
simulations

accuracy

decadesofF.

Jacques-Louis
1resources.

Hecht

etLions

Duringal.

thepastatLaboratoire(LJLL),

Sorbonne

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019

Universite, have been developing a domain

obtained from discretization of PDEs. However,specific language, “FreeFem++”. It provides thethis project targets solution of linear sparsemost important and powerful tools to perform

matrices generated in FreeFem++. Thanks tofinite element analyses. This is done through aFreeFem++ and its ability to describe complex

simple description of the weak formulation of amathematical models, real-world three dimen-

PDE using keywords very close to the actual

sional computations using a mixed shared- andmathematical formulation. It also provides

distributed-memory parallelism will be feasible

automatic mesh decomposition for discreti-

in different fields such as shape optimization or

zation, and operators to access or modify sparse

fluid–structureand dense matrices or vectors. Through thisdirectly applicable to various problems in thelanguage,manufacturing industry.

linearsystems

fromweakinteraction.

Thesewill

beformulations with selected finite element space

In Japan, there is the "ADVENTURE" project,

are automatically generated. The user only

started in late 1990 to provide a finite element

needs to describe nonlinear iterations and/orsimulation software with domain decomposition

time evolution algorithms.

methods. The software consists of several

Various kind offinite elements are available, e.g., standard

separatedpiecewise polynomials in a tetrahedron forsimulation of specific physical phenomena: heat

elasticity and fluid problems, vector-valued

flow, incompressible flow, elasticity, electro-finite elements for oil reservoir simulation andmagnetic,

semi-conductor problems, where satisfying apredefined physical parameters to fit physicalconservation principle is important, and edge

models. This employs a completely different

elements designed specially for electromagnetic

approach than the one of FreeFem++ withproblems.

HPDDM, where the finite element matrices are

The aim of this project is to combine the highgenerated from a given weak formulation.

efficiency of direct solver whose name isAn LU factorization of one million unknowns“Dissection”

byinsideofonenodethroughpackages

wheremulti-frontal

totheperformuser

methods,setse.g.,

numerical

valueofPARDISO,

shared-memory parallelism and the numerical

MUMPS, and Dissection, requires about 60GB

robustness of iterative solver named “HPDDM”

memory and one million unknowns stays in aover distributed parallelism to tackle largekind of sweat spot of the computational-scale problems.

complexity with efficient use of cache memory
to compensate low factor of byte/FLOP in3.

Significanceas

JH PCN

Joint

modern multi-core architectures.

Research Project

Cybermedia center of Osaka University has

There are several software engineering projects

been providing vector- and scalar-architecture

toalgebra

clusters. Both hardware resources are ideal

algorithm to solve linear systems or to solve

platforms with multi- and many-core CPU toeigenvalue problems as a general purpose

demonstrate effectiveness and robustness of thelibrary, but almost of them result in solvers for

combined strategy of direct and iterative solvers

problems consisting of dense matrices, and areenhanced by eigenvalue analysis in complex

implement

numerical

linearnot fit to solve sparse matrices that are

actual engineering and multi-physics problems
2Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019that are currently only solvable by a monolithic

https://github.com/hpddm/hpddm

direct solver.

The second step is for a preconditioner in theiterative substructuring method with enhanced

4.

Outline of Research Achievem ents

coarse space, which is known as the balancing

up to FY2018

domaindecomposition(BDD)

method.

Construction of coarse space is realized byThis project started in FY 2019.

collecting zero energy mode of the local stiffness5.

DetailsofFY2019

matrix. Especially for the linear elasticity

Researchproblem, the zero energy mode is explicitly

Achievem ents
HPDDM provides a unified framework for

written as set of rigid body vectors. Therefore,

domaindecomposition

including

HPDDM has an interface to receive such

additive

Schwarz

iterative

vectors through FreeFem++, where user can

substructuring methods. Those methods will bedefine rigid body vectors based on geometricalenhanced

has

information with finite element meshes. Incomparable performance to multigrid methods

contrast to HPDDM, ADVENTURE also uses

such as GAMG or BoomerAMG.

BDDFor linear system with one million unknowns,

employed to avoid singular matrices, which may

direct solvers based on the LU factorization

perturb the convergence of a preconditioned

may reach up to 40% of the peak performance of

iterative method when condition numbers areshared-memory architecture. The multi-frontal

high. HPDDM follows mathematical setting ofmethod and blocked operations lead to the task

BDD without modification and archives fastof dense matrix-matrix products, which isand robust convergence. However difficulty

efficient on modern CPU.

comes for nonlinear elasticity problems. Newton

A careful combination of Dissection anditeration is a common tool to solve nonlinear

HPDDM is expected to deliver a high parallelsystem, but for finite element approximation toefficiency, which could be applicable to a wide

therange of problem, easily prototyped throughperturbation to the zero energy mode by

FreeFem++.

interpolation. Therefore, it is impossible toThere are three steps, depending on the levelpredict exact eigenvectors corresponding to zeroof sophistication of preconditioner for domain

eigenvalue of the stiffness matrix explicitly.

methods

methods

through

GenEO.

andHPDDM

butanonlinear

Dissection

decomposition iteration.

regularization

elasticity

hasparameter

introduces

implemented

a

is

some

powerful

The first step is for the domain decomposition

algorithm for rank-deficient problem, which cansolver with both overlapping/non-overlapping

numerically detect dimension of the kernel

subdomains

space of the large sparse matrix.

without

preconditioner.

This

implementation is rather straightforward toIn addition to existing routines to pass data

follow the existing routines to send and receiveon sparse matrix and solution vector between

sparse matrices data and solutions of dense

Dissection local solver inside of the node andvector, for other director solvers, PARDISO andHPDDMMUMPS. The improved interface is accessible

allocation of kernel vectors which depends onover

distributed

nodes,

dynamic

characteristics of local matrix can handle

through open source repository of HPDDM,
3

Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures
Final Report for JHPCN Joint Research of FY 2019

automatic generation of the coarse space only

creation of Schur complement matrix in theby algebraic procedure without physical and/or

GenEO preconditioner. Dissection has better

mathematical knowledge of the finite element

implementation to compute solution of linearmatrix.

system with multiple RHS of dense vectors.

ThethirdsophisticatedThis procedure has been extended to sparse

preconditioner named GenEO, which is also

vector by analyzing nonzero pattern of thedeveloped in LJLL. The main idea of GenEO issparsity. Thanks to the binary tree structure ofan extension of the coarse space consisting of

nested-dissection

zero eigenvalue to of smaller ones. GenEO

forward substitution could be at least half. Moreutilizes information retrieved by performing

precise measurement on performance update

concurrent eigenvalue analysis of each local

will be a future issue.

problem.

stepThis

isacomplexityofdifferent

The other progress is on implementation oflinear

algebra

FreeFem++ on vector architecture provided aslibraries. For additive Schwarz methods withNEC SX-ACE with proprietary C++ compiler.

non-verlapping subdomains, the local problems

FreeFem++ is written with the standard C++11

consist of a sparse matrix and the local

features and template library. It was necessary

generalized eigenvalue problems are resolved

to introduce slight modification of C++ template

by ARPACK, where a series of matrix-vector

instantiation by adding explicit declarations.

products are performed as part of the Arnoldi

However, it is not possible to implement aprocess.

furthermore

useful feature of FreeFem++, dynamic loading

preconditioned, and thus it is essential to solve

facility due to lack of shared library system oflinear systems with a sparse matrix, which

SX-ACE operating system.

requires usage of Dissection within ARPACK.

7.

For

theGenEO

This

numerical

process

iterative
requiresrealizedordering,

bystate-of-the-art

isfor

is

substructuringadenseofPublications

andconference

PapersPresentations

method,eigensolver,List(1) International

as(Non-refereed)

implemented in LAPACK using a tridiagonalization process. The generalized eigenvalue[1] "Dissection sparse direct solver andproblems consist of a local Schur complement

parallel task management", A. Suzuki,with Dirichlet boundary conditions, and a

HPC Asia 2020, Fukuoka, Japan

restriction of the global Schur complement.

[2] "An industrial application of a freeHere, it is necessary to compute all elements ofboundary problem with contact angle

the Schur complement using solutions of theand volume constraint", A. Suzuki,

direct solver with sparse right-hand sides

11th FreeFEM days, Paris, France(RHS).

(2) Presentations at domestic conference

6.

Progress during FY2019 and Future(Non-refereed)

Prospects[3] “Finite element solution of a free

The key ingredient to enhance performance ofboundary problem of the liquid phasecombined packages “Dissection” and “HPDDM”

with contact angle and volume constraint”

is developed, which is realized by improvement

A. Suzuki, 65th NCTAM, Sapporo, Japanof sparse multiple RHS solver of Dissection for4