EN FR
EN FR

2024Activity reportProject-TeamALPINES

RNSR: 201321071B
  • Research center Inria Paris Centre at Sorbonne University
  • In partnership with:CNRS, Sorbonne Université
  • Team name: Algorithms and parallel tools for integrated numerical simulations
  • In collaboration with:Laboratoire Jacques-Louis Lions (LJLL)
  • Domain:Networks, Systems and Services, Distributed Computing
  • Theme:Distributed and High Performance Computing

Keywords

Computer Science and Digital Science

  • A6.1.1. Continuous Modeling (PDE, ODE)
  • A6.1.4. Multiscale modeling
  • A6.1.5. Multiphysics modeling
  • A6.2.1. Numerical analysis of PDE and ODE
  • A6.2.5. Numerical Linear Algebra
  • A6.2.7. High performance computing
  • A6.3. Computation-data interaction
  • A6.3.1. Inverse problems
  • A7.1. Algorithms

Other Research Topics and Application Domains

  • B3.3.1. Earth and subsoil
  • B9.5.2. Mathematics
  • B9.5.3. Physics

1 Team members, visitors, external collaborators

Research Scientists

  • Frédéric Nataf [Team leader, CNRS, Senior Researcher]
  • Sever Hirstoaga [INRIA, Researcher]
  • Emile Parolin [INRIA, Researcher]
  • Pierre-Henri Tournier [CNRS, Researcher]

Faculty Member

  • Frédéric Hecht [Sorbonne University, Professor]

Post-Doctoral Fellow

  • Lukas Spies [INRIA, Post-Doctoral Fellow, from Oct 2024]

PhD Students

  • Brieuc Antoine Dit Urban [INRIA, from Oct 2024]
  • Siwar Badreddine [INRIA, until Feb 2024]
  • Filippo Brunelli [SORBONNE UNIVERSITE, from Nov 2024]
  • Tom Caruso [INRIA, from Nov 2024]
  • Jean-Guillaume De Damas [INRIA]
  • Nicola Galante [INRIA]
  • Edouard Timsit [INRIA]

Technical Staff

  • Daniel Torres Gonzales [SORBONNE UNIVERSITE, Engineer]

Interns and Apprentices

  • Filippo Brunelli [INRIA, Intern, from Apr 2024 until Aug 2024]
  • Elarif Djabir [INRIA, Intern, from May 2024 until Aug 2024]
  • Vincent Robert [INRIA, Intern, from Apr 2024 until Sep 2024]
  • Chris Teku [INRIA, Intern, from May 2024 until Jul 2024]

Administrative Assistant

  • Laurence Bourcier [INRIA]

2 Overall objectives

2.1 Introduction

The focus of our research is on the development of novel parallel numerical algorithms and tools appropriate for state-of-the-art mathematical models used in complex scientific applications, and in particular numerical simulations. The proposed research program is by nature multi-disciplinary, interweaving aspects of applied mathematics, computer science, as well as those of several specific applications, as porous media flows, elasticity, wave propagation in multi-scale media, molecular simulations, and inverse problems.

Our first objective is to develop numerical methods and tools for complex scientific and industrial applications that will enhance their scalable execution on the emergent heterogeneous hierarchical models of massively parallel machines. Our second objective is to integrate the novel numerical algorithms into a middle-layer that will hide as much as possible the complexity of massively parallel machines from the users of these machines.

3 Research program

3.1 Overview

The research described here is directly relevant to several steps of the numerical simulation chain. Given a numerical simulation that was expressed as a set of differential equations, our research focuses on mesh generation methods for parallel computation, novel numerical algorithms for linear algebra, as well as algorithms and tools for their efficient and scalable implementation on high performance computers. The validation and the exploitation of the results is performed with collaborators from applications and is based on the usage of existing tools. In summary, the topics studied in our group are the following:

  • Numerical methods and algorithms
    • Mesh generation for parallel computation
    • Solvers for numerical linear algebra:

      domain decomposition methods, preconditioning for iterative methods

    • Computational kernels for numerical linear algebra
    • Tensor computations for high dimensional problems
  • Validation on numerical simulations and other numerical applications

3.2 Domain specific language - parallel FreeFem

In the engineering, researchers, and teachers communities, there is a strong demand for simulation frameworks that are simple to install and use, efficient, sustainable, and that solve efficiently and accurately complex problems for which there are no dedicated tools or codes available. In our group we develop FreeFem++, a user-dedicated language for solving PDEs. The goal of FreeFem++ is not to be a substitute for complex numerical codes, but rather to provide an efficient and relatively generic tool for:

  • getting a quick answer to a specific problem,
  • prototyping the resolution of a new complex problem.

The current users of FreeFem++ are mathematicians, engineers, university professors, and students. In general for these users the installation of public libraries as MPI, MUMPS, Ipopt, Blas, lapack, OpenGL, fftw, scotch, PETSc, SLEPc is a very difficult problem. For this reason, the authors of FreeFem++ have created a user friendly language, and over years have enriched its capabilities and provided tools for compiling FreeFem++ such that the users do not need to have special knowledge of computer science. This leads to an important work on porting the software on different emerging architectures.

Today, the main components of parallel FreeFem++ are:

  • definition of a coarse grid,
  • splitting of the coarse grid,
  • mesh generation of all subdomains of the coarse grid, and construction of parallel data structures for vectors and sparse matrices from the mesh of the subdomain,
  • discretization by a chosen numerical method,
  • call to a linear solver,
  • analysis of the result.

All these components are parallel, except for the last point which is not in the focus of our research. However for the moment, the parallel mesh generation algorithm is very simple and not sufficient, for example it addresses only polygonal geometries. Having a better parallel mesh generation algorithm is one of the goals of our project. In addition, in the current version of FreeFem++, the parallelism is not hidden from the user, it is done through direct calls to MPI. Our goal is also to hide all the MPI calls in the specific language part of FreeFem++. In addition to these in-house domain decomposition methods, FreeFem++ is also linked to PETSc solvers which enables an easy use of third parties parallel multigrid methods.

3.3 Solvers for numerical linear algebra

Iterative methods are widely used in industrial applications, and preconditioning is the most important research subject here. Our research considers domain decomposition methods 6 and iterative methods and its goal is to develop solvers that are suitable for parallelism and that exploit the fact that the matrices are arising from the discretization of a system of PDEs on unstructured grids.

One of the main challenges that we address is the lack of robustness and scalability of existing methods as incomplete LU factorizations or Schwarz-based approaches, for which the number of iterations increases significantly with the problem size or with the number of processors. This is often due to the presence of several low frequency modes that hinder the convergence of the iterative method. To address this problem, we study different approaches for dealing with the low frequency modes as coarse space correction in domain decomposition or deflation techniques. We minimize the memory footprint and communication avoiding algorithms by leveraging mixed precision arithmetic.

We also focus on developing boundary integral equation methods that would be adapted to the simulation of wave propagation in complex physical situations, and that would lend themselves to the use of parallel architectures. The final objective is to bring the state of the art on boundary integral equations closer to contemporary industrial needs. From this perspective, we investigate domain decomposition strategies in conjunction with boundary element method as well as acceleration techniques (H-matrices, FMM and the like) that would appear relevant in multi-material and/or multi-domain configurations. Our work on this topic also includes numerical implementation on large scale problems, which appears as a challenge due to the peculiarities of boundary integral equations.

3.4 Computational kernels for numerical linear and multilinear algebra

The design of new numerical methods that are robust and that have well proven convergence properties is one of the challenges addressed in Alpines. Another important challenge is the design of parallel algorithms for the novel numerical methods and the underlying building blocks from numerical linear algebra. The goal is to enable their efficient execution on a diverse set of node architectures and their scaling to emerging high-performance clusters with an increasing number of nodes.

Increased communication cost is one of the main challenges in high performance computing that we address in our research by investigating algorithms that minimize communication, as communication avoiding algorithms. The communication avoiding algorithmic design is an approach originally developed in our group since more than ten years (initially in collaboration with researchers from UC Berkeley and CU Denver). While our first results concerned direct methods of factorization as LU or QR factorizations, our recent work focuses on designing robust algorithms for computing the low rank approximation of a matrix or a tensor. We focus on both deterministic and randomized approaches.

Our research also focuses on solving problems of large size that feature high dimensions as in molecular simulations. The data in this case is represented by objects called tensors, or multilinear arrays. The goal is to design novel tensor techniques to allow their effective compression, i.e. their representation by simpler objects in small dimensions, while controlling the loss of information. The algorithms are aiming to being highly parallel to allow to deal with the large number of dimensions and large data sets, while preserving the required information for obtaining the solution of the problem.

4 Application domains

Our work is by design as general purpose as possible in the field of scientific computing based on partial differential equations modeling including multiphysics and multiscale modeling. Moreover, since our methods are available in standalone open sources libraries and in the free finite element specific language FreeFEM which is developed in our team, they have many potential users that we are not necessarily aware of. We give here some recent works performed in collaboration with team members.

4.1 Radiative transfer in the atmosphere

The Hierarchical matrix acceleration techniques implemented in our parallel Htool 7.1.3 library, which is interfaced with FreeFEM, allows to efficiently solve the radiative transfer equations with complexity n log(n), where n is the number of vertices of the physical domain. The method can be used to study the change of temperature in the atmosphere, and in particular the effect of clouds and greenhouse gases 7, 13.

4.2 Inverse problems

4.2.1 Time Reversal techniques

We focus on methods related to the blend of time reversal techniques and absorbing boundary conditions (ABC) used in a non standard way. Since the seminal paper by [M. Fink et al., Imaging through inhomogeneous media using time reversal mirrors. Ultrasonic Imaging, 13(2):199, 1991.], time reversal is a subject of very active research. The principle is to back-propagate signals to the sources that emitted them. The initial experiment was to refocus, very precisely, a recorded signal after passing through a barrier consisting of randomly distributed metal rods. In [de Rosny and Fink. Overcoming the diffraction limit in wave physics using a time-reversal mirror and a novel acoustic sink. Phys. Rev. Lett., 89 (12), 2002], the source that created the signal is time reversed in order to have a perfect time reversal experiment. In 29, we improve this result from a numerical point of view by showing that it can be done numerically without knowing the source. This is done at the expense of not being able to recover the signal in the vicinity of the source. In 30, we have extended these methods to non uniform media.

4.2.2 Seismic imaging

In the context of seismic imaging, frequency-domain full-waveform inversion (FWI) is suitable for long-offset stationary-recording acquisition, since reliable subsurface models can be reconstructed with a few frequencies and attenuation is easily implemented without computational overhead. In the frequency domain, wave modeling is a Helmholtz-type boundary-value problem which requires to solve a large and sparse system of linear equations per frequency with multiple right-hand sides (sources).

For sparse node acquisitions covering large domains (more than 50 millions of unknowns), an iterative solver should be the method of choice. However, fast convergence remains challenging in the high frequency regime due to the non-definiteness of the operator and to the discretization constraints needed to minimize dispersion, hence requiring efficient preconditioners. For such wave propagation problems, our Schwarz domain decomposition preconditioners are good candidates to accelerate the time to solution and to provide robustness with respect to heterogeneities and frequency, for both finite difference and finite element discretizations 34.

4.2.3 Electroencephalography

The inverse source problem in electroencephalography (EEG) is usually solved using a numerical head model composed of different layers of tissues (scalp, skull, brain, ...) having constant electric conductivity. However, the conductivity of the skull can hardly be assumed to be constant: (i) hard and spongy bones, especially for neonates (fontanels), (ii) conductivity models vary between individuals, and through time for a single person.

In 5, we consider a new numerical method for the EEG inverse source problem which takes into account the heterogeneity of the skull. In particular, in the first step of source reconstruction, the so-called "cortical mapping" procedure which deals with (i) data completion on the scalp knowing partial measurements of the electric potential from electrode measurements and (ii) data transmission from the scalp to the brain surface, we use the quasi-reversibility method which provides a regularized solution of Cauchy problems in a bounded domain. The problem is solved using the finite element method in the scalp and skull layers.

4.2.4 Medical imaging with microwave tomography

Microwave tomography is a novel imaging modality with a large number of potential attractive medical applications, and is based on the difference between the dielectric properties of normal and diseased tissues. Microwave tomography features rapid data acquisition time, and together with rapid tomographic reconstructions allows detecting, identifying and monitoring injuries continuously.

From a computational point of view, microwave imaging requires the solution of an inverse problem based on a minimisation algorithm. Reconstruction algorithms are computationally intensive with successive solutions of the forward problem needing efficient numerical modelling and high-performance parallel computing. This methodology involves distinct research fields: optimisation, inverse problems, approximation and solution methods for the simulation of the forward problem modelled by Maxwell’s equations. The latter is challenging in itself as the modelling must accurately take account of the high heterogeneity and complexity of the different tissues.

Our numerical framework for microwave imaging is based on the FreeFEM 7.1.1 finite element software, and the forward problem is solved efficiently in parallel using Domain Decomposition methods implemented in the HPDDM 7.1.2 library.

Brain stroke detection and monitoring

Stroke, or cerebrovascular accident (CVA), is a major cause of disability and death worldwide. About 85% of CVAs are ischemic due to cerebral infarction, caused by an interruption of the blood supply to some part of the brain, and 15% are haemorrhagic. Differentiating between ischemic and haemorrhagic CVAs is an essential part of the initial workup of the patient, and rapid and accurate diagnosis is crucial for patient survival; here, neuroimaging plays a vital role. Computed Tomography (CT) and Magnetic Resonance Imaging (MRI) are the ‘gold’ standards, but their use is not well suited to efficient medical care of CVAS, as they are bulky diagnostic instruments and cannot be used in continuous brain monitoring. A non-invasive and transportable/portable device for the characterisation and monitoring of CVAs would have clear clinical applications, beginning with the very first instance of patient care in an ambulance and extending to continuous patient monitoring at the hospital.

Microwave tomography holds great promise for rapid brain stroke diagnosis. The raw data acquired by the microwave imaging system can be wirelessly transferred to a remote computing center, where the tomographic images will be computed. The images can then be quickly transferred to the hospital.

In 33, 32, we demonstrate on synthetic data the feasibility of a microwave imaging technique for the characterisation and monitoring of strokes. Using high performance computing, we are able to obtain a tomographic reconstruction of the brain in less than two minutes.

Detection and imaging of rotator cuff tears

One of the most challenging shoulder injuries is rotator cuff tear, which increases with aging and particularly happens among athletes. These tears cause pain and highly affect the functionality of the shoulder.

In  3, 2, we propose to study the feasibility of microwave tomographic imaging to detect these tears. We introduce a wearable imaging system design, adapted to the real shoulder structure and designed to surround it partially. We investigate the efficiency of the imaging system numerically using realistic CAD models for shoulder profile and bones, including humerus and scapula, and taking into account the dielectric properties of the different tissues (bone, tendon, muscle, skin, synovial fluid). The reconstructed images of the shoulder joint obtained using noisy synthetic data show that the proposed imaging system is capable of accurately detecting and localizing large (5mL) and partial (1mL) rotator cuff tears.

4.3 Numerical methods for time harmonic wave propagation problems

Acoustic, electromagnetic and elastic linear waves are ubiquitous phenomena in science and engineering. The numerical simulation of their propagation and interaction is a core task in areas like medical imaging, seismic remote sensing, design of electronic devices, atmospheric particle scattering, radar and sonar modelling, etc. The design of appropriate numerical schemes is challenging because of approximation issues arising from the fact that the solutions are highly oscillatory, which usually require a large computational effort to produce a meaningful numerical approximation. In addition, common formulations suffer from stability issues, such as the numerical dispersion (also known as pollution effect), and sign-indefiniteness, which typically limit the achievable accuracy. The variety of different applications and the challenges of their numerical simulations produces a great research effort in trying to design new efficient schemes: new approximation spaces using degrees of freedom sparingly (wave based and Trefftz methods), preconditionners for solving large linear systems (e.g. by domain decomposition), new formulations, acceleration techniques, etc.

4.4 Molecular simulations

Molecular simulation is one of the most dynamic areas of scientific computing. Its field of application is very broad, ranging from theoretical chemistry and drug design to materials science and nanotechnology. It provides many challenging problems to mathematicians, and computer scientists.

In the context of the ERC Synergy Grant EMC2 we address several important limitations of state of the art molecular simulation. In particular, the simulation of very large molecular systems, or smaller systems in which electrons interact strongly with each other, remains out of reach today. In an interdisciplinary collaboration between chemists, mathematicians and computer scientists, we focus on developing a new generation of reliable molecular simulation algorithms and software.

5 Social and environmental responsibility

5.1 Impact of research results

Our activities in scientific computing help in the design and optimization of non fossil energy resources. Here are three examples we are aware of and there are many others:

  • the startup Airthium specialized in decarbonized energy uses our free parallel finite element software FreeFEM.
  • The French electricity company EDF has integrated our domain decomposition methods GenEO via the open source library HPDDM in its simulation software Salomé in order to be able to perform large scale computations related to nuclear reactors.
  • In the context of climate change, our work on radiative transfer studies the change of temperature in the atmosphere, in particular due to the effect of clouds and greenhouse gases.

6 Highlights of the year

6.1 Awards

ERC Synergy 2024

PSINumScat (Phase Space Inspired Numerical Methods for High-Frequency Scattering)

This Synergy project has been granted to Pierre-Henri Tournier, (the other leaders are Euan A. Spence, professor at the University of Bath, and Jeffrey Galkowski, professor at University College London), aims to design, analyze and implement in free software new methods for the numerical resolution of high-frequency acoustic and electromagnetic wave propagation problems. These methods will be faster, more reliable, and more widely applicable than current state-of-the-art methods, leveraging techniques from fundamental mathematics specifically designed to study high-frequency problems. These techniques come from the field of semi-classical analysis and are based on the study of functions in phase space (the combination of physical space and frequency space). Since numerical analysis and semi-classical analysis of high-frequency wave propagation problems are mature areas of mathematics, it is surprising that these two areas operate largely in isolation from each other. None of the extensive knowledge of high-frequency waves gained from semi-classical analysis has been used in the design of numerical methods for these problems. The project will take advantage of this untapped potential and make the resulting methods available to the scientific community in the free software for multi-physics numerical simulation FreeFEM developed at the Laboratory J.L. Lions-INRIA team Alpines.

7 New software, platforms, open data

7.1 New software

7.1.1 FreeFem++

  • Name:
    FreeFrem
  • Keywords:
    Scientific computing, High performance computing, Boundary element method
  • Scientific Description:

    FreeFem++ is a partial differential equation solver. It has its own language. freefem scripts can solve multiphysics non linear systems in 2D and 3D.

    Problems involving PDE (2d, 3d) from several branches of physics such as fluid-structure interactions require interpolations of data on several meshes and their manipulation within one program. FreeFem++ includes a fast 2d̂-tree-based interpolation algorithm and a language for the manipulation of data on multiple meshes (as a follow up of bamg (now a part of FreeFem++ ).

    FreeFem++ is written in C++ and the FreeFem++ language is a C++ idiom. It runs on Macs, Windows, Unix machines. FreeFem++ replaces the older freefem and freefem+.

  • Functional Description:
    FreeFem is a PDE (partial differential equation) solver based on a flexible language that allows a large number of problems to be expressed (elasticity, fluids, etc) with different finite element approximations on different meshes.
  • URL:
  • Contact:
    Frédéric Hecht

7.1.2 HPDDM

  • Keywords:
    High performance computing, Parallel computing
  • Scientific Description:
    HPDDM is an efficient implementation of various domain decomposition methods (DDM) such as one- and two-level Restricted Additive Schwarz methods, the Finite Element Tearing and Interconnecting (FETI) method, and the Balancing Domain Decomposition (BDD) method. This code has been proven to be efficient for solving various elliptic problems such as scalar diffusion equations, the system of linear elasticity, but also frequency domain problems like the Helmholtz equation. A comparison with modern multigrid methods can be found in the thesis of Pierre Jolivet.
  • Functional Description:
    HPDDM is an efficient implementation of various domain decomposition methods (DDM) such as one- and two-level Restricted Additive Schwarz methods, the Finite Element Tearing and Interconnecting (FETI) method, and the Balancing Domain Decomposition (BDD) method. It is interfaced, among others, via PETSc, FEM, Feel++.
  • URL:
  • Contact:
    Pierre Jolivet
  • Participants:
    Frédéric Nataf, Pierre Jolivet

7.1.3 HTool

  • Keyword:
    Hierarchical matrices
  • Functional Description:
    HTool is a C++ header-only library implementing compression techniques (e.g. Adaptive Cross Approximation) using hierarchical matrices (H-matrices). The library uses MPI and OpenMP for parallelism, and is interfaced with HPDDM for the solution of linear systems.
  • URL:
  • Contact:
    Pierre Marchand

7.1.4 scatensor

  • Name:
    Scalable and Communication Avoiding Tensor Library
  • Keywords:
    Tensor decomposition, Linear algebra, Matrix calculation, Distributed computing, MPI, Communication avoiding
  • Scientific Description:
    Matrix and tensor decomposition in parallel, to obtain a compressed version.
  • Functional Description:
    - Reading distributed matrices from a file. - Generating random matrices. - Distributed QR decomposition with tournament pivoting (QRTP) of a matrix. - Reading/Generating a tensor. - Tucker decomposition of a tensor, using QRTP.
  • Publications:
  • Contact:
    Laura Grigori
  • Participants:
    David Frenkiel, Laura Grigori, Matthias Beaupere

8 New results

Participants: Alpines members.

8.1 A robust and adaptive GenEO-type domain decomposition preconditioner for H(curl) problems in general non-convex three-dimensional geometries

In 1, we develop and analyse domain decomposition methods for linear systems of equations arising from conforming finite element discretisations of positive Maxwell-type equations, namely for H(curl) problems. It is well known that convergence of domain decomposition methods rely heavily on the efficiency of the coarse space used in the second level. We design adaptive coarse spaces that complement a near-kernel space made from the gradient of scalar functions. The new class of preconditioner is inspired by the idea of subspace decomposition, but based on spectral coarse spaces, and is specially designed for curl-conforming discretisations of Maxwell's equations in heterogeneous media on general domains which may have holes. Our approach has wider applicability and theoretical justification than the well-known Hiptmair-Xu auxiliary space preconditioner, with results extending to the variable coefficient case and non-convex domains at the expense of a larger coarse space.

8.2 Coarse spaces for non-symmetric two-level preconditioners based on local generalized eigenproblems

In 27, we present an analysis of adaptive coarse spaces which is quite general since it applies to symmetric and non symmetric problems, to symmetric preconditioners such the additive Schwarz method (ASM) and to the non-symmetric preconditioner restricted additive Schwarz (RAS), as well as to exact or inexact subdomain solves. The coarse space is built by solving generalized eigenvalues in the subdomains and applying a well-chosen operator to the selected eigenvectors.

8.3 A Robust Two-Level Schwarz Preconditioner For Sparse Matrices

In 24 daas:hal-04381509, we introduce a fully algebraic two-level additive Schwarz preconditioner for general sparse large-scale matrices. The preconditioner is analyzed for symmetric positive definite (SPD) matrices. For those matrices, the coarse space is constructed based on approximating two local subspaces in each subdomain. These subspaces are obtained by approximating a number of eigenvectors corresponding to dominant eigenvalues of two judiciously posed generalized eigenvalue problems. The number of eigenvectors can be chosen to control the condition number. For general sparse matrices, the coarse space is constructed by approximating the image of a local operator that can be defined from information in the coefficient matrix. The connection between the coarse spaces for SPD and general matrices is also discussed. Numerical experiments show the great effectiveness of the proposed preconditioners on matrices arising from a wide range of applications. The set of matrices includes SPD, symmetric indefinite, nonsymmetric, and saddle-point matrices. In addition, we compare the proposed preconditioners to the state-of-the-art domain decomposition preconditioners.

8.4 Microwave tomographic imaging of shoulder injury

One of the most challenging shoulder injuries is rotator cuff tear, which increases with aging and particularly happens among athletes. These tears cause pain and highly affect the functionality of the shoulder. The motivation of this work 17 is to detect these tears by microwave tomographic imaging. This imaging method requires the solution of an inverse problem based on a minimization algorithm, with successive solutions of a direct problem. The direct problem consists in solving the time-harmonic Maxwell’s equations discretized with a Nédélec edge finite element method, and we make use of an ORAS Schwarz domain decomposition preconditioner to solve it efficiently in parallel. We propose a wearable imaging system design with 96 antennas distributed over two fully-circular and two half-circular layers, adapted to the real shoulder structure and designed to surround it partially. We test the imaging system numerically using realistic CAD models for shoulder profile and bones, including humerus and scapula, and taking into account the dielectric properties of the different tissues (bone, tendon, muscle, skin, synovial fluid). The reconstructed images of the shoulder joint obtained using noisy synthetic data show that the proposed imaging system is capable of accurately detecting and localizing large (5mL) and partial (1mL) rotator cuff tears. The reconstruction takes 20 minutes on 480 computing cores and shows great promise for rapid diagnosis or medical monitoring. In a second step, we take advantage of numerical modeling to try and optimize the number of antennas in the imaging system, achieving a drastic reduction from 96 to 32 antennas while still being able to detect the injury in the difficult partial tear case, reducing the computing time from 20 to 12 minutes.

In 16, 21, we use the previous forward numerical modeling techniques to generate a generalizable dataset of scattering parameters in order to bypass real-world data collection challenges. We then explore the use of machine learning (ML) algorithms for the fast detection of shoulder tendon injuries. The corresponding data of various healthy and injured models are categorized into two classes. We use a support vector machine (SVM) to differentiate between injured and healthy shoulder models. This approach is more efficient in terms of required memory resources and computing time compared with traditional imaging methods, and we manage to achieve an accuracy of 100%.

8.5 New developments in FreeFEM

Work has been started to ease the parallelization of a user's sequential code. In this setting, the parallelization of the assembly and solution steps is entirely hidden to the user, and a sequential code can be parallelized with only very few changes in the script. For example, PETSc can easily be used to solve the linear system stemming from a coupled problem with a fieldsplit preconditioner.

The FreeFEM software now includes the GenEO for saddle-point 11 examples with PCHPDDM in PETSc.

FreeFEM can now run Markdown (.md) files as well as .edp files. The Markdown syntax can be used to document FreeFEM scripts. When running a .md file, FreeFEM will execute the code contained in the FreeFEM code blocks. Documented Markdown examples in the distribution can be viewed on the documentation website. Markdown can also be previewed with Visual Studio Code, with the FreeFEM VS Code extension providing syntax highlighting for FreeFEM code blocks.

8.6 Reduced models for highly oscillating differential equations

In this part we are concerned with solving Vlasov equation involving several time scales. This work shows that finding and solving reduced models, though complicated to derive, turns out to be very useful for reducing the computational cost of an equation. In addition, using the parareal algorithm, a method performing parallel computing in time, allows to accelerate the computations.

Thus, in 19 we derive a reduced first-order model from a two-scale asymptotic expansion in a small parameter, in order to approximate the solution of a stiff differential equation. The problem of interest is a multi-scale Newton-Lorentz equation modeling the dynamics of a charged particle under the influence of a linear electric field and of a perturbed strong magnetic field. First, we show that in short times, the first-order model provides a much better approximation than the zero-order one. Then, we propose a volume-preserving method using a particular splitting technique to solve numerically the first-order model. Finally, it turns out that the first-order model does not systematically provide a satisfactory approximation in long times. To overcome this issue, we implement a recent strategy based on the Parareal algorithm, in which the first-order approximation is used for the coarse solver. This approach allows to perform efficient and accurate long-time simulations for any small parameter. Numerical results for two realistic Penning traps are provided to support these statements.

Second, we extended the previous idea to a more general equation, which models the motion of a charged particle in a tokamak-like magnetic field, dependent on space and time. Thus, by symbolic computations, we obtained first-order equations approximating stiff differential equations in the general case. This is a work in progress.

Finally, during the internship of Crédo Fanou (stage de M2) we studied theoretical results on the error estimation of reduced models, obtained by averaging, for generic stiff differential equations. We applied these results to Newton-Lorentz equation in order to obtain such error estimation. Elarif Djabir (stage de M1) performed, during his internship, numerical simulations of such stiff equations, in order to study the error bounds and their optimality. This is ongoing work.

8.7 Stable Trefftz approximation of Helmholtz solutions using evanescent plane waves

Superpositions of plane waves are known to approximate well the solutions of the Helmholtz equation. Their use in discretizations is typical of Trefftz methods for Helmholtz problems, aiming to achieve high accuracy with a small number of degrees of freedom. However, Trefftz methods lead to ill-conditioned linear systems, and it is often impossible to obtain the desired accuracy in floating-point arithmetic.

In 26 we extend the analysis of 12 from two to three space dimensions. We show that a judicious choice of plane waves can ensure high-accuracy solutions in a numerically stable way, in spite of having to solve such ill-conditioned systems. Numerical accuracy of plane wave methods is linked not only to the approximation space, but also to the size of the coefficients in the plane wave expansion. We show that the use of plane waves can lead to exponentially large coefficients, regardless of the orientations and the number of plane waves, and this causes numerical instability. We prove that all Helmholtz fields are continuous superposition of evanescent plane waves, i.e., plane waves with complex propagation vectors associated with exponential decay, and show that this leads to bounded representations. We provide a constructive scheme to select a set of real and complex-valued propagation vectors numerically. This results in an explicit selection of plane waves and an associated Trefftz method that achieves accuracy and stability. A non-trivial challenge in 3D was the parametrization of the complex direction set. Our approach involves defining a complex-valued reference direction and then consider its rigid-body rotations via Euler angles. Then, by generalizing the Jacobi–Anger identity to complex-valued directions, we prove that any solution of the Helmholtz equation on a three-dimensional ball can be written as a continuous superposition of evanescent plane waves in a stable way. Our numerical results showcase a great improvement of our proposed method on the achievable accuracy compared to standard method.

We are currently exploiting these new results to construct discontinuous Trefftz schemes based on plane-wave approximation spaces, which was the subject of the internship of Vincent Robert . In addition, and this is very novel with respect to the current literature, we are able to obtain conforming Galerkin discretization from plane-wave approximation spaces, which opens up to a whole new approach for Trefftz methods. Finally, we are working on new integral representations of any Helmholtz solution in a smooth convex domain, thereby generalyzing the above works on disk and balls.

8.8 Domain decomposition methods for random multi-scale Helmholtz problems arising in ultrasound imaging

The development of new quantitative ultrasound imaging algorithms, which aim to reconstruct a map of the local ultrasound speed in the medium, requires a validation process that can be achieved through numerical simulation. With this application in mind, we consider the scattering of plane waves by a tissue-mimicking medium where up to a hundred unresolved scatterers per wavelength are randomly distributed throughout the medium. The domains (about a hundred wavelengths in size) require billion degrees of freedom in a simulation, which corresponds to the state of the art in terms of direct numerical simulation capacity.

In the internship of Chris Teku , we investigated the efficiency and scalability of one-level and two-levels domain decomposition techniques to accurately solve the full scale model. The primary objective was to validate quantitative stochastic homogenization results obtained in 31, particularly the asymptotic expansions of the scattered field with respect to the size of the scatterers. This is still on-going work.

8.9 Parallel approximation of the exponential of Hermitian matrices

In 9, we consider a rational approximation of the exponential function to design an algorithm for computing matrix exponential in the Hermitian case. Using partial fraction decomposition, we obtain a parallelizable method, where the computation reduces to independent resolutions of linear systems. We analyze the effects of rounding errors on the accuracy of our algorithm. We complete this work with numerical tests showing the efficiency of our method and a comparison of its performances with Krylov algorithms.

8.10 A finite element toolbox for the Bogoliubov-de Gennes stability analysis of Bose-Einstein condensates

In 20, we present a finite element toolbox for the computation of Bogoliubov-de Gennes modes used to assess the linear stability of stationary solutions of the Gross-Pitaevskii (GP) equation. Applications concern one (single GP equation) or two-component (a system of coupled GP equations) Bose-Einstein condensates in one, two and three dimensions of space. An implementation using the free software FreeFem++ is distributed with this paper. For the computation of the GP stationary (complex or real) solutions we use a Newton algorithm coupled with a continuation method exploring the parameter space (the chemical potential or the interaction constant). Bogoliubov-de Gennes equations are then solved using dedicated libraries for the associated eigenvalue problem. Mesh adaptivity is proved to considerably reduce the computational time for cases implying complex vortex states. Programs are validated through comparisons with known theoretical results for simple cases and 97 numerical results reported in the literature.

8.11 Mathematics and Finite Element Discretizations of Incompressible Navier-Stokes Flows

With : C. Bernardi, V. Girault, F. Hecht, P.-A. Raviart, B. Riviere.

This book in final preparation (2024) is a revised, updated, and augmented version of the out-of-print classic book “Finite Element Methods for Navier–Stokes Equations" by Girault and Raviart published by Springer in 1986 [GR]. The incompressible Navier–Stokes equations model the flow of incom- pressible Newtonian fluids and are used in many practical applications, including computational fluid dynamics. In addition to the basic theoretical analysis, this book presents a fairly exhaustive treatment of the up-to-date finite element discretizations of incompressible Navier–Stokes equa- tions and a variety of numerical algorithms used in their computer implementation. It covers the cases of standard and non-standard boundary conditions and their numerical discretizations via the finite element methods. Both conforming and non-conforming finite elements are examined in detail, as well as their stability or instability. The topic of time-dependent Navier–Stokes equa- tions, which was missing from [GR], is now presented in several chapters. In the same spirit as [GR], we have tried as much as possible to make this book self-contained and therefore we have either proved or recalled all the theoretical results required. This book can be used as at textbook for advanced graduate students.

8.12 Randomized Householder QR

With Laura Grigori, Edouard Timsit

In 8, we introduce a randomized Householder QR factorization (RHQR). This factorization can be used to obtain a well conditioned basis of a set of vectors and thus can be employed in a variety of applications. We discuss in particular the usage of this randomized Householder factorization in the Arnoldi process. Numerical experiments show that RHQR produces a well conditioned basis and an accurate factorization. We observe that for some cases, it can be more stable than Randomized Gram-Schmidt (RGS) in both single and double precision.

8.13 Randomized Implicitly Restarted Arnoldi method for the non-symmetric eigenvalue problem

With Jean-Guillaume de Damas, Laura Grigori

In  25 we introduce a randomized algorithm for solving the non-symmetric eigenvalue problem, referred to as randomized Implicitly Restarted Arnoldi (rIRA). This method relies on using a sketchorthogonal basis during the Arnoldi process while maintaining the Arnoldi relation and exploiting a restarting scheme to focus on a specific part of the spectrum. We analyze this method and show that it retains useful properties of the Implicitly Restarted Arnoldi (IRA) method, such as restarting without adding errors to the Ritz pairs and implicitly applying polynomial filtering. Experiments are presented to validate the numerical efficiency of the proposed randomized eigenvalue solver.

8.14 PDE-constrained optimization within FreeFEM

With Gontran Lance, and Emmanuel Trélat

The book  22 is aimed at students and researchers who want to learn how to efficiently solve constrained optimization problems involving partial differential equations (PDE) using the FreeFEM software. PDE-constrained optimization problems are frequently encountered in many academic and industrial contexts. Readers should have a basic knowledge of the analysis and numerical solution of partial differential equations using finite element methods, optimization, algorithms and numerical implementation.

9 Bilateral contracts and grants with industry

9.1 Bilateral contracts with industry

Participants: F. Nataf, B. Antoine dit Urban.

  • Contract with IFPEN and ONERA, January 2025 - December 2027, that funds half of the Phd thesis of Brieuc Antoine dit Urban on "Méthodes d’apprentissage pour les solveurs linéaires et préconditionneurs pour matrices creuses" . Supervisor Frédéric Nataf , T. Faney and E. Martin.

10 Partnerships and cooperations

10.1 European initiatives

10.1.1 ERC Synergy

PSINumScat

Participants: Pierre-Henri Tournier, Frédéric Nataf.

  • Title:
    Phase-space-inspired Numerical Methods for High Frequency Wave Scattering
  • Duration:
    From May 1, 2025 to 2031
  • Partners:
    • UNIVERSITY OF BATH, UK
    • UNIVERSITY COLLEGE LONDON, UK
    • CENTRE NATIONAL DE LA RECHERCHE SCIENTIFIQUE CNRS (CNRS), France
  • Inria contact:
    Pierre-Henri Tournier (Alpines)
  • Coordinator:
  • Summary:
    Phase-space-inspired Numerical Methods for High Frequency Wave Scattering (PSINumScat) is an ERC Synergy funded project. Designing fast and reliable algorithms to numerically simulate the behaviour of high-frequency acoustic and electromagnetic waves is a longstanding open problem in computational mathematics. These waves underpin a plethora of communication and imaging technologies; therefore any progress towards solving this problem will have wide impact. By exploiting techniques from pure mathematics specifically designed to study high-frequency problems, PSINumScat aims to design, analyse, and implement in open-source software new methods for the numerical solution of high-frequency acoustic and electromagnetic wave scattering problems.

10.1.2 Inno4Scale

AMCG

Participants: Frédéric Nataf, Pierre-Henri Tournier, Lukas Spies.

  • Title:
    Asynchronous Modified Conjugate Gradient Method
  • Duration:
    From February 1, 2024 to March 28, 2025
  • Partners:
    • INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET AUTOMATIQUE (INRIA), France
    • CENTRE NATIONAL DE LA RECHERCHE SCIENTIFIQUE CNRS (CNRS), France
    • SORBONNE UNIVERSITE, France
    • WIKKI, Germany
  • Inria contact:
    Frédéric Nataf (Alpines)
  • Coordinator:
    Henrik Rusche (WIKKI)
  • Summary:

    This Innovation Study uses an asynchronous linear solver technology, focusing on a new approach for parallel computing. Unlike traditional methods that rely on global sums and synchronized pairwise communication, this study operates without these constraints. Asynchronous techniques offer advantages by eliminating synchronization barriers, thereby enhancing the efficiency of computer systems through reduced periods of inactivity. Furthermore, their inherent ability to tolerate communication latencies enhances fault tolerance, ensuring robustness in complex computing environments. In addition, asynchronous methods unlock the potential for leveraging power-saving and energy-efficient features offered by modern hardware architectures. This positions them as promising solutions for highly parallel and diverse hardware platforms.

    The experimental concept is based on the classic conjugate gradient (CG) method, adapted for partitioned data environments. This study however diverges significantly from conventional usage by treating each partition as an independent problem. Consequently, it advances a Krylov subspace for each partition autonomously, without requiring global synchronization. Global problem coupling and convergence are achieved by perturbing the local solution using information from remote partitions. Crucially, these perturbations are updated asynchronously, eliminating the need for complete partition progress synchronization. Through empirical validation and performance analysis, this study demonstrates the efficacy and potential of the asynchronous linear solver concept. By embracing the principles of asynchrony and partitioned computation, this technology is able to revolutionize scalable and fault-tolerant parallel computing paradigms, as well as enhance computational efficiency.

10.2 National initiatives

10.2.1 ANR

PEPR Numpex - PC Exama

Participants: Brieuc Antoine Dit Urban, Tom Caruso, Frédéric Hecht, Frédéric Nataf, Emile Parolin, Pierre-Henri Tournier.

Numpex is a French program dedicated to Exascale which is divided into five Research Projects. The first one is Exa-MA (Methods and Algorithms of Exascale) aimed at advancing scientific simulations and modeling capabilities to reach and surpass the exascale barrier. The project is organized into several work packages (WP) that focus on different aspects of the research objectives. Frédéric Nataf is in charge of the Work Package 3 (WP3): Numerical Kernels and Coupled Solvers Focuses on designing and implementing efficient and possibly provable numerical kernels and solvers for largescale problems. Tasks include domain decomposition methods, data sparsity techniques, multiple precision, adaptive solution strategies, and efficient coupling of multiphysics simulations.

Partners in WP3: LIP6 Sorbonne University, INRIA Bordeaux, CEA.

Muffin

Participants: Sever Hirstoaga, Pierre-Henri Tournier.

: ANR appel à projet générique 2019-2024.

S. Hirstoaga and P.-H. Tournier are members of the project MUFFIN, whose objective is to explore and optimize original computational scenarios for multi-scale and high dimensional transport codes, with priority applications in plasma physics. Several approximation methods are planed to be developed. It is at the frontier of computing and numerical analysis and intends to reduce the computational burden in the context of intensive calculation. Principal Investigator: B. Després (Sorbonne University).

11 Dissemination

11.1 Promoting scientific activities

11.1.1 Scientific events: organisation

Member of the organizing committees

Frédéric Nataf , Emile Parolin : SIAM Conference on Applied Linear Algebra (LA24)

11.1.2 Journal

Member of the editorial boards

Frédéric Nataf Associate Editor at Journal of Numerical Mathematics

Reviewer - reviewing activities

Emile Parolin reviewed 2 articles for Mathematics of Computation and IMA-Journal of Numerical Analysis journals.

11.1.3 Invited talks

Frédéric Nataf :

Indo-german Workshop on Hardware-aware Scientific Computing, Heidelberg, 2024.

Two-day course at University of Sevilla, Adaptive Coarse Space for Domain Decomposition Method, Feb. 2024

Emile Parolin was invited to the 10th International Conference on Computational Methods in Applied Mathematics (CMAM-10), June 2024.

Sever Hirstoaga was invited to SIAM Conference on Parallel Processing for Scientific Computing, Session MS56: Novel Algorithms and HPC Implementations for Exascale Particle-In-Cell Methods, Baltimore, March 2024.

11.2 Teaching - Supervision - Juries

11.2.1 Teaching

  • Licence 3: Sever Hirstoaga , Méthodes numériques, 38 HETD (cours + TP), EISE 3, Polytech Sorbonne.
  • Doctorat: Sever Hirstoaga , Parallel-in-time numerical methods and applications, 7 HETD (cours + TP), summer school "New trends in computing", IRMA, université de Strasbourg.
  • Master 2: Frédéric Nataf, Course on Domain Decomposition Methods, Sorbonne University, 22h.
  • Master 1: Emile Parolin , Numerical simulation and approximation of solutions to elliptic PDEs, Sorbonne Université, 40h.
  • Master 2: Emile Parolin , High Performance Computing, Sorbonne Université, 15h.

12 Scientific production

12.1 Major publications

12.2 Publications of the year

International journals

International peer-reviewed conferences

  • 21 inproceedingsS.Sahar Borzooei, P.-H.Pierre-Henri Tournier, V.Victorita Dolean and C.Claire Migliaccio. A SVM-based approach for detecting tendon Injury.AP-S/URSI 2024 - International Symposium on Antennas and Propagation and ITNC-USNC-URSI Radio ScienceFlorence, ItalySeptember 2024, pp. 1511-1512HALDOIback to text

Scientific books

Doctoral dissertations and habilitation theses

  • 23 thesisS.Siwar Badreddine. Leveraging symmetries and low-rank structure of matrices and tensors in high-dimensional quantum chemistry problems.Sorbonne UniversitéMarch 2024HAL

Reports & preprints

Other scientific publications

  • 28 inproceedingsS.Sahar Borzooei, P.-H.Pierre-Henri Tournier, V.Victoria Dolean, C.Christian Pichot, N.Nadine Joachimowicz, H.Helene Roussel and C.Claire Migliaccio. Modèle numérique d’un système d’imagerie microonde pour détecter les blessures de l’épaule par apprentissage machine.JNM 2024 - 23ième Journées Nationales MicroondesJuan les Pins, FranceJune 2024, 4HAL

12.3 Cited publications

  • 29 articleF.F Assous, M.M Kray, F.F Nataf and E.E Turkel. Time-reversed absorbing condition: application to inverse problems.Inverse Problems2762011, 065003URL: http://stacks.iop.org/0266-5611/27/i=6/a=065003back to text
  • 30 articleF.Franck Assous and F.Frédéric Nataf. Full-waveform redatuming via a TRAC approach: a first step towards target oriented inverse problem.Journal of Computational Physics4402021, 110377back to text
  • 31 unpublishedJ.Josselin Garnier, L.Laure Giovangigli, Q.Quentin Goepfert and P.Pierre Millien. Scattered wavefield in the stochastic homogenization regime.September 2023, working paper or preprintHALback to text
  • 32 articleP.-H.P.-H. Tournier, I.I. Aliferis, M.M. Bonazzoli, M.M. de Buhan, M.M. Darbas, V.V. Dolean, F.F. Hecht, P.P. Jolivet, I. E.I. El Kanfoud, C.C. Migliaccio, F.F. Nataf, C.Ch. Pichot and S.S. Semenov. Microwave tomographic imaging of cerebrovascular accidents by using high-performance computing.Parallel Computing852019, 88-97URL: https://www.sciencedirect.com/science/article/pii/S0167819118301042DOIback to text
  • 33 articleP.-H.Pierre-Henri Tournier, M.Marcella Bonazzoli, V.Victorita Dolean, F.Francesca Rapetti, F.Frederic Hecht, F.Frederic Nataf, I.Iannis Aliferis, I.Ibtissam El Kanfoud, C.Claire Migliaccio, M.Maya de Buhan, M.Marion Darbas, S.Serguei Semenov and C.Christian Pichot. Numerical Modeling and High-Speed Parallel Computing: New Perspectives on Tomographic Microwave Imaging for Brain Stroke Detection and Monitoring.IEEE Antennas and Propagation Magazine5952017, 98-110DOIback to text
  • 34 articleP.-H.Pierre-Henri Tournier, P.Pierre Jolivet, V.Victorita Dolean, H. S.Hossein S. Aghamiry, S.Stéphane Operto and S.Sebastian Riffo. 3D finite-difference and finite-element frequency-domain wave simulation with multilevel optimized additive Schwarz domain-decomposition preconditioner: A tool for full-waveform inversion of sparse node data sets.GEOPHYSICS8752022, T381-T402DOIback to text