# Manifold reducibility for a Lagrangian finite element solver with remeshing

Max Beckermann, Massimiliano Cremonesi, Andrea Barbarulo

► To cite this version:

Max Beckermann, Massimiliano Cremonesi, Andrea Barbarulo. Manifold reducibility for a Lagrangian finite element solver with remeshing. Computer Methods in Applied Mechanics and Engineering, 2025, 446, Part A, pp.118212. 10.1016/j.cma.2025.118212. hal-05191181

HAL Id: hal-05191181

https://hal.science/hal-05191181v1

Submitted on 29 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Manifold reducibility for a Lagrangian finite element solver with remeshing

Max Beckermann aID a,b,*, Massimiliano Cremonesi bID, Andrea Barbarulo aID

## ARTICLE INFO

**Keywords:**

Particle finite element method

PFEM

Lagrangian fluid simulations

Proper orthogonal decomposition

Manifold reduction

Model order reduction

## ABSTRACT

Lagrangian solvers are well suited for large deformation fluid problems containing non-trivial mesh distortions. This includes problems with evolving free surfaces, breaking waves, fluid-structure interaction problems, or free-flowing fluids with complex rheology. We focus on the Particle Finite Element Method (PFEM). This Lagrangian solver combines the advantages of a particle description of the fluid with a finite element solver and efficient remeshing algorithms. It creates vast data structures of varying sizes as model outputs. As the mesh varies, regrouping the different fields of interest into a single solution manifold becomes non-trivial. This paper will discuss multiple strategies to reduce a *posteriori* the solution manifolds formed from PFEM simulations. Data reduction consists in a trade-off between preserving the original data precision and having smaller reduced data sets and is often an important first step towards surrogate models. Two novel mapping-based nonlinear reduction methods will be presented to reduce different high-fidelity PFEM simulations. The data reduction will be made on complete sets of time-dependent simulation snapshots and consists in the separation of space and time with the calculation of space and time modes. We will compare the precision - i.e. the $L^2$ error of the reconstructed snapshot - the reducibility - i.e. the number of required modes for a given precision in the reduction - and the implementation complexity of the proposed methods.

## 1. Introduction

The Particle Finite Element Method (PFEM) [1] is a powerful Lagrangian solver relying on a particle description of a medium paired with a finite element solver. The mesh displacement, the velocity, and pressure fields are computed every time step. To deal with the large displacement of the medium, it utilizes re-meshing algorithms that effectively ensure the non-distortion of the evolving finite element mesh. It is ideally suited to study free-surface fluid flows [2,3] as any interface is automatically recovered. It has been used for a wide array of academic and engineering applications: fluid-structure interaction problems [4–6], landslides [7,8], multi-phase flows [9], or nonlinear manufacturing processes [10–13]. Note that while we chose to focus on PFEM, the general scope of this paper is not limited to this solver. The main arguments laid out as well as some of the methods could be adapted to any Lagrangian solver for which the discretization evolves over time, such as Smoothed Particle Hydrodynamics [14] or Material Point Method [15].

PFEM is under active development to improve its speed, precision, and applicability [16–18] and is an expanding research field. Naturally, there is nowadays a growing interest in the development of surrogate models for PFEM. It is, however, in its early stages

<footnote>ª Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS - Laboratoire de Mécanique Paris-Saclay, 91190, Gif-Sur-Yvette, France</footnote>

<footnote>ᵇ Departement of Civil and Environmental Engineering, Politecnico di Milano, Piazza Leonardo da Vinci, 32, 20133, Milan, Italy</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: max.beckermann@centralesupelec.fr (M. Beckermann), massimiliano.cremonesi@polimi.it (M. Cremonesi), andrea.barbarulo@centralesupelec.fr (A. Barbarulo).</footnote>
and has a minimal presence in the existing literature. The key idea is to learn from a dataset of simulations and construct a cost-effective surrogate model. In this work, we focus on a beneficial preliminary step: the reduction of PFEM solution manifolds. This process minimizes storage needs while enabling the reconstruction of the original data with a certain accuracy. As a matter of fact, not only does it lie at the core of the learning phase in Model Order Reduction surrogate modeling [19,20], but it also greatly enhances storage efficiency, which is essential for handling the large datasets required in machine learning-based approaches [21].

The classical linear reduction technique, typically called Proper Orthogonal Decomposition (POD), Karhunen–Loève decomposition or Principal Component Analysis (PCA) [22–25], cannot be used with PFEM data for four major reasons:

* PFEM data is not defined on a fixed mesh but on successive different unstructured meshes. An extension of linear POD for arbitrary finite element discretization, relying on the construction of supermeshes [26], has been developed and used in Gräßle and Hinze [27], Little and Farhat [28], Ullmann et al. [29], Yano [30]; however, its use is mainly limited to structured and nested meshes due to cost.

* The problems that PFEM is efficient at solving, such as free surface fluid flows, are by nature hard to reduce with standard linear reduction techniques. The reducibility, characterized by their Kolmogorov N-width decay [31], does not decrease fast enough with the number of modes considered in the linear reduction. This is a well-identified problem in the bibliography [32,33], and nonlinear reduction techniques have been developed to tackle this issue. These techniques map the solution to a manifold that is then more reducible by construction. Lagrangian methods, nonlinear ansatz, transport mappings, and registration methods are all used interchangeably in model order reduction to refer to these nonlinear reduction methods. Numerous formulations have been proposed to construct the mappings via optimal transport [34,35], optimization-based strategies [36,37], velocity-driven registration [38–40]... Not only would the computation of such mapping be too expensive for PFEM, as every time step would require a map, but it would also require interpolation between the data points and the mapped meshes. Moreover, the interpolation errors are non-negligible compared to the reduction precision, making all these methods less appealing.

* To avoid interpolations, we choose to use the particle trajectories as mappings. PFEM updates its Lagrangian map at every time step according to the current velocity; however due to remeshing, it uses an updated Lagrangian procedure [1]. Therefore, there does not exist a single reference configuration to map our data back to in a bijective manner. At a given time in PFEM, the reference configuration does not coincide with the initial configuration at time $t=0$ but rather with the last remeshing instance. Note that this is intrinsic to PFEM as remeshing is how it handles large mesh distortions, a total Lagrangian formulation cannot be used in PFEM.

* Finally, the particle trajectories in PFEM are not continuous. As a means to prevent point concentration due to the Lagrangian nature of the solver, points are free to appear and disappear from one time step to another [1]. The individual particle trajectories must be identified if we wish to use the Lagrangian map provided by PFEM. These trajectories are considered gappy as the data is missing at points in time.

As a result, we decided to directly adapt the Lagrangian mesh computed with PFEM and perform a nonlinear reduction. Two methods have been developed and are presented in this paper:

* The first extends the updated Lagrangian map to a total Lagrangian map. If the bijectivity is conserved, the new maps go from every time step to time zero, and the data can be transported on the initial configuration. As the mapping is based on particle movement, all the nodes are superposed in the initial configuration and a supermesh can be computed efficiently. The reduction is then performed on the mapped supermesh.

* The second reduces along the gappy particle trajectories and deals directly with the missing data, ignoring the connectivity changes. As we know where and when the particles are added and removed, our first idea is to give mesh coherent extensions to the missing data points and perform a POD on the completed data. Another option is to directly use the Gappy POD technique [41–43], an iterative algorithm that constructs a reduction basis for gappy data.

The paper is organized as follows: Section 2 presents the Particle Finite Element Method and the reduction problem. Section 3 covers the tools for linear reduction for both the fixed mesh and arbitrary mesh settings, with time-dependent simulation data. Section 4 introduces the mapping approaches and our two nonlinear reduction methods based on the Lagrangian PFEM map. The two methods are illustrated with a simple 1D advection case. Finally, in Section 5, we compare our methods with four PFEM examples and conclude in Section 6.

## 2. The particle finite element method

To understand the unique PFEM data structure, let us start with an overview of the general PFEM algorithm given in Cremonesi et al. [1]:

1. Start the problem with a set of initial particles (Fig. 1a).

2. Generate a mesh using the set of particles with the Delaunay triangulation (Fig. 1b) and recover the free surface with the alpha shape method (Fig. 1c).

3. Solve the Lagrangian Finite Element equations with P1 linear finite elements for the main unknowns and the mesh displacement. The exact equations are problem-dependent; note that both implicit and explicit solvers have been developed for PFEM.

4. Update the mesh displacement (Fig. 1d).
<figure><img src="image_3.png" /><figcaption>Fig. 1. The PFEM solving steps, as seen in Cremonesi et al. [1].</figcaption></figure>

5. Proceed to the next time step; if mesh is not too distorted, then go to step 3. If it is, then go to step 2 (Fig. 1e-g). One might also add and remove nodes if necessary at this step. This is typically done if an element is too small or too large after the calculation of the new connectivity.

Importantly, the steps and the procedures presented in the paper are independent of the equations considered in step 3.

## 2.1. The resulting data structure

Let $\Omega \subset \mathbb{R}^d$ with $d \in \{1, 2, 3\}$ be an open and bounded space domain corresponding to the total bounding box of our PFEM simulation, $I = [0, T]$ with $0 < T < \infty$ the time domain corresponding to the simulation interval and $\Omega_i \subset \Omega$ be the current domain at time $t_i$, $u$ will designate the scalar field of interest computed with PFEM, which we will attempt to reduce. It can typically be a component of the velocity or displacement field, the pressure field, or any other data computed on our moving particles. Theoretically, we have $u(x, t) \in L^2(I, V)$ with $V \subset L^2(\Omega)$ as our functional space of choice. Practically, we much rather define the function on the current fluid domain and therefore introduce a time discretization $0 = t_1 < t_2 < \dots < t_n = T$ with $t_i - t_{i-1} = \Delta T = \frac{T}{n-1}$, $\forall i \in \llbracket 1, n \rrbracket$ and express our functional space at our current time $t_i$ with $V_i \subset L^2(\Omega_i)$. Our field of interest can therefore be written as

$$u(\cdot, t_i) := u(t_i) \subset V_i.$$

We deliberately will only work with equally spaced time discretization but all methods presented in this paper could be extended to arbitrary time discretization. The finite element space discretization is also distinct for every time step. $\forall i \in \llbracket 1, n \rrbracket$ we have a unique mesh connectivity, finite element discretization, and therefore discretized functional space $V_{h_i} \subset V_i \subset L^2(\Omega_{t_i})$. Note that $\dim(V_{h_i}) = N_i \in \mathbb{N}^*$ is not necessarily constant and can vary through time because of node additions and removals within the moving mesh. It corresponds to the number of active particles at a given time. For a given $u(t_i) \in V_i$ the space discretization can be expressed $\forall i \in \llbracket 1, n \rrbracket$, $\forall x \in \Omega_{t_i}$

$$u(x, t_i) \approx u^h(x, t_i) := u^{h_i}(x) = \mathcal{N}^{h_i}(x)^T u^{h_i},$$

with $u^{h_i} \in V_{h_i}$ the approximated field, $\mathcal{N}^{h_i}(x)$ the time-specific shape function arrays of size $N_i$, and $u^{h_i} \in \mathbb{R}^{N_i}$ the nodal values. To summarize, the PFEM simulation data for a given field $u$, which we will attempt to reduce, presents itself in the following Lagrangian finite element structure

$$D = \{ u^{h_i} \in \mathbb{R}^{N_i}, \mathcal{N}^{h_i} \in V_i^{N_i}, \forall i \in \llbracket 1, n \rrbracket \}. \qquad (1)$$

It corresponds to a set of $n$ snapshots of the field $u$ at different times, characterized by the nodal values $u^{h_i}$ and the shape functions $\mathcal{N}^{h_i}$ of the finite element discretization at time $t_i$. The shape functions $\mathcal{N}^{h_i}$ are each represented by their respective nodal positions and connectivity, of characteristic size $N_i$.

## 2.2. Reducing PFEM solution manifolds

The idea of data reduction is to construct an approximation space $V_N$ of low dimension $N \in \mathbb{N}$ for a so-called solution manifold $\mathcal{K} \subset V$ of a given parametric PDE. We have $V$ our Hilbert space, $|| \cdot ||$ its norm, and $\langle \cdot, \cdot \rangle$ its scalar product. In the case of our time-dependent PDE solving for the field $u$ that we defined above, the infinite-dimensional solution manifold can be written as

$$\mathcal{K} = \{ u(t_i) \mid i \in \llbracket 1, n \rrbracket \} \subset V \subset L^2(\Omega). \qquad (2)$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. The PFEM reduction diagram.</figcaption></figure>

$V_N$ is then built with $(\psi_i)_{1\le i \le N}$, an orthogonal family of $V$ computed with different techniques. They all revolve around finding a linear approximation space $V_N$ such that we get the closest possible to the theoretical optimum of

$$d_N(\mathcal{K}) = \inf_{V_N \subseteq V \atop \dim V_N \le N} \sup_{\psi \in \mathcal{K}} \inf_{w \in V_N} \|\psi - w\|, \quad (3)$$

called the Kolmogorov N-width of $\mathcal{K}$. This N-width is at the heart of linear reduction theory [31,32], and the decay of $d_N(\mathcal{K})$ with $N$ is fundamentally a measurement of the reducibility of the solution manifold $\mathcal{K}$. Finally, to compute the reduced field, the orthogonal projector $\Pi_N$ can also be introduced

$$V \to V_N \qquad (4) \\ u \to \Pi_N(u) = \sum_{i=1}^{N} \langle u, \psi_i \rangle \psi_i.$$

In practice we require $V$ to be finite-dimensional for almost every data reduction strategy. It is therefore replaced by a discrete approximation that is so highly resolved that the discretization error is negligible with regard to the model reduction error. Even if our space discretization Eq. (1) is indeed highly dimensional ($\forall i \in [1, n]$, $N_i >> N$), it is not the same throughout our parametric set. Therefore, we cannot easily build a discrete solution manifold. A new mapped manifold has to be constructed that will allow for one unique discretization of our data and the construction of the projector $\Pi_N$ Eq. (4).

For each time step $t_i$, we will construct a forward mapping $X_{t_i}$ and its associated backward mapping $\Xi_{t_i}$. We will then, with the individual $X_{t_i}$, map all our data to a common manifold $\mathcal{K}$, reduce the data, and construct the projector $\Pi_N$ Eq. (4). Finally, we can use the $\Xi_{t_i}$ to map our reduced data back to the original space. Combining the forward and backward mappings ($X_{t_i}$ and $\Xi_{t_i}$) with the projector ($\Pi_N$), we then have new reduction operators $\Pi_{i,N}$. The full reduction process can be visualized in the diagram in Fig. 2. Calling $\tilde{u}(t)$ the reduced form of $u(t)$, one can assess the total $L^2$ reduction error with respect to the original snapshots Eq. (1):

$$\epsilon_S^2(N) = \frac{\int_I \|u(t) - \tilde{u}(t)\|_{L^2(\Omega_i)}^2 dt}{\int_I \|u(t)\|_{L^2(\Omega_i)}^2 dt} \approx \frac{\sum_{i=1}^n \|u^{h_i} - \Pi_{i,N}(u^{h_i})\|_{M_{h_i}}^2}{\sum_{i=1}^n \|u^{h_i}\|_{M_{h_i}}^2}, \quad (5)$$

with $M_{h_i}$ the symmetric positive definite mass matrix associated with the finite element discretization such that $\{M_{h_i}\}_{j,k} = \langle N_{h_i}^j, N_{k}^{h_i} \rangle_{L^2(\Omega_i)}$ and with equal quadrature weights at all quadrature points.

To summarize, we will in this paper present two new ways to construct discretized manifolds $\mathcal{K} \subset V$ from our PFEM dataset $D$ Eq. (1), quantify their reducibility with Eq. (3), assess their total $L^2$ reduction error Eq. (5), and finally discuss their implementation complexity.

### 3. Computing a reduced basis with the proper orthogonal decomposition

We will now introduce the main tools that can be used to extract a linear approximation from time-dependent simulation data. These ideas are typically used in the offline phase of the Proper Orthogonal Decomposition (POD). We will only cover the tools used in our examples; we refer the reader to Cordier and Bergmann [22], Farhat [23], Volkwein [24,25] for more in-depth information about the subject. We will first look at the POD for data defined on a fixed mesh data and then for data defined on arbitrarily chosen meshes that can differ from one snapshot to another.

#### 3.1. Classical POD for time-dependent simulation snapshots on a fixed mesh

Let us start with some time-dependent Partial Differential Equation (PDE) simulation data following the exact time discretization described in Section 2.1 but with a fixed mesh. Discretizing space with a finite element space $V_h$ of dimension $m$, we can write the snapshots on the nodal values as:

$$\mathcal{K} = \{u^h(t_i) | i \in [1, n]\} \subset \mathbb{R}^m.$$

The discretized POD problem [24,25] consists in looking for $V_N = \text{Span}\{\{\psi_j\}_{j=1}^N\}$, with $\{\psi_j\}_{j=1}^N \subset \mathbb{R}^m$

$$\Psi^h = \{\psi_j\}_{j=1}^N = \arg \min_{\psi_1, \dots, \psi_N} \sum_{i=0}^n \|u^h(t_i) - \sum_{j=1}^N (u^h(t_j), \psi_j)_{M_h} \psi_j\|_{M_h}^2, \quad (6) \\ \text{s.t.} \langle \psi_i, \psi_j \rangle_{M_h} = \delta_{ij}, \forall i, j \in [1, r]$$
with $M_h$, the symmetric positive definite mass matrix associated with the finite element discretization such that $\{M_h\}_{i,j} = \langle N_i^h, N_j^h \rangle_{L^2}$. In case of non uniform time discretization, additional weights appear in the sum, but we will not consider this case in this paper. To solve problem Eq. (6), we need to solve for the first N eigenvectors $\{\psi_j\}_{j=1}^N \subset \mathbb{R}^m$ and eigenvalues $\{\lambda_j\}_{j=1}^N \subset \mathbb{R}$ of the eigenfunction problem defined as

$$ \sum_{i=1}^{n} \langle u^h(t_i), \psi \rangle_{M_h} u^h(t_i) = \lambda \psi. $$

It is known from the Eckart-Young theorem [22–25] that the minimizers also satisfy the following formula for the approximation error

$$ \sum_{i=0}^{n} \|u^h(t_i) - \sum_{j=1}^{N} \langle u^h(t_i), \psi_j \rangle_{M_h} \psi_j\|_{M_h}^2 = \sum_{i=r+1}^{n} \lambda_j. \quad (7) $$

This error quantifies the magnitude of the missing information, note that we are in a case for which $n < m$. In practice, we assemble the snapshots into a snapshot matrix $Y = [u^h(t_1), \dots, u^h(t_n)] \in \mathbb{R}^{m \times n}$ and compute these eigenfunctions via two possible methods:

* Solving Eq. (6) is equivalent to doing the Singular Value Decomposition (SVD) [22–25] of the weighted snapshots $\tilde{Y} = M_h^{\frac{1}{2}} Y$ and truncating at $r$ singular values. With $M_h^{\frac{1}{2}}$ originating from the Cholesky factorization of the positive-definite mass matrix $M_h$. We therefore get

$$ \tilde{Y} = U \Sigma \Phi^T, \quad (8) $$

where the left singular vectors - also called space modes - $U = [U_1, \dots, U_m] \in \mathbb{R}^{m \times m}$ and right singular vectors - or time modes - $\Phi = [\phi_1, \dots, \phi_n] \in \mathbb{R}^{n \times n}$ are both orthogonal matrices,

$$ \Sigma = \begin{pmatrix} D & 0 \\ 0 & 0 \end{pmatrix} \in \mathbb{R}^{m \times n}, $$

with $D = \text{diag}(\sigma_1, \dots, \sigma_d) \in \mathbb{R}^{d \times d}$ and $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_d > 0$ where $d = \text{rank}(Y)$. The equivalence of the two problems gives us that $\forall j \in [1, n]$, $\psi_j = M_h^{-\frac{1}{2}} U_j$ and $\lambda_j = \sigma_j^2$. Note that if $r$ is known or fixed beforehand, one can also directly compute the truncated SVD of $\tilde{Y}$.

* As we consistently have a larger space discretization than number of snapshots, Sirovich [44] has introduced to solve a cheap equivalent $n \times n$ eigenvalue problem that will compute the first $r$ time modes $\{\phi_j\}_{j=1}^N$ and $\{\lambda_j\}_{j=1}^N$ the eigenvalues of the original problem

$$ \tilde{Y}^T \tilde{Y} \phi_j = Y^T M_h Y \phi_j = \lambda_j \phi_j \quad \text{for } 1 \le j \le r. \quad (9) $$

The space modes $\{\psi_j\}_{j=1}^N$ minimizing the POD formulation can then be found with

$$ \psi_j = \frac{1}{\sqrt{\lambda_j}} Y \phi_j \quad \text{for } 1 \le j \le r. $$

Once the POD eigenfunctions are computed and knowing the approximation error presented in Eq. (7), we will quantify the neglected energy and the reducibility of our problem with

$$ 1 - \mathcal{E}_{\text{POD}}^2(N) = 1 - \frac{\sum_{j=1}^{N} \lambda_j}{\sum_{j=1}^{n} \lambda_j} = 1 - \frac{\sum_{j=1}^{N} \sigma_j^2}{\sum_{j=1}^{n} \|u^h(t_j)\|_{M_h}^2}, \quad (10) $$

where $\mathcal{E}_{\text{POD}}^2(N)$ represents the relative energy of the snapshots captured by the first N POD basis vectors. Note that the optimality of the POD approach ensures that the linear subspace constructed from the POD basis also corresponds to the optimal space required to compute the Kolmogorov N-width Eq. (3) and that, therefore, $1 - \mathcal{E}_{\text{POD}}^2(N)$ also corresponds to the square of the normalized Kolmogorov N-width of $\mathcal{K}$ denoted by $d_{\mathcal{K}}^2(\mathcal{K})$.

## 3.2. POD on arbitrary finite element discretization

An approach for snapshots on arbitrary meshes has been developed and used in Gräßle and Hinze [27], Little and Farhat [28]; it requires the evaluation of the inner products of the snapshots in a common Hilbert space. In practice, they consider directly the infinite-dimensional solution manifold

$$ \mathcal{K} = \{u^{h_i} | i \in [1, n]\} \subset L^2(\Omega), \quad (11) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. An supermesh construction example.</figcaption></figure>

The technique then looks to solve the space-continuous version of problem Eq. (6) but with our finite element snapshots. We will therefore look for $V_N = \text{Span}\{\{\psi_j\}_{j=1}^N\}$, with $\{\psi_j\}_{j=1}^N \subset L^2(\Omega)$ which solve

$$ \Psi = \{\psi_j\}_{j=1}^N = \arg \min_{\psi_1, \dots, \psi_N} \sum_{i=0}^n \|u^{h_i} - \sum_{j=1}^N \langle u^{h_i}, \psi_j \rangle_{L^2} \psi_j\|_{L^2}^2. \quad (12) $$

$$ \text{s.t.} \langle \psi_i, \psi_j \rangle_{L^2} = \delta_{ij}, \forall i, j \in [1, r] $$

To solve this problem, [27] proposes to assemble a matrix $W \in \mathbb{R}^{n \times n}$

$$ W_{ij} = \langle u^{h_i}, u^{h_j} \rangle_{L^2} \quad \forall (i, j) \in [1, n]^2. \quad (13) $$

The finite-dimensional time modes are the eigenvectors of $W$, which, along with the eigenvalues, are given by $\{(\lambda_k, \phi_k)\}_{k=1}^n$. If we sort them in order such that $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_n$, the optimal infinite dimensional basis functions $\{\psi_j\}_{j=1}^N \subset L^2(\Omega)$ solutions of the problem Eq. (12) are given by a linear combination of the snapshots weighted with the first $N$ time modes

$$ \psi_j = \frac{1}{\sqrt{\lambda_j}} \sum_{i=1}^{n} u^{h_i} \phi_i(t_i) \quad \forall j \in [1, N]. \quad (14) $$

The fundamental difficulty and cost of this method lies in the calculation of the $L^2$ scalar products Eq. (13), requiring a so-called supermesh for every pair of snapshots. A triangular mesh $T_+$ is called a supermesh of two triangular meshes $T_1$ and $T_2$, if and only if:

* Every node of both $T_1$ and $T_2$ is also a node of $T_+$.

* Every edge of $T_1$ and $T_2$ can be obtained by joining one or multiple edges of $T_+$.

* $T_+$ contains the minimum number of required edges to verify these first two requirements.

An example can be shown in Fig. 3; note that $T_+$ is not necessarily unique as the intersections are not all triangles. In practice, one needs to be able to detect all the intersecting polygons and then remesh the polyhedra that contain more than three edges with, for example, a Delaunay triangulation. Refer to Farrell and Maddison [26] for practical algorithms to construct the supermesh.

The only way to express the infinite-dimensional space modes $\psi_j$ without additional errors is to compute them on the total supermesh of all snapshots. The construction of this final supermesh is typically expensive. Therefore, this method has mainly been used on structured and nested meshes [27–30], for which you can drastically reduce the calculation costs.

## 4. Transporting the PFEM data and constructing a mapped manifold

To reduce our PFEM data, one could use the method introduced in Section 3.2 for arbitrary finite element discretization. The supermeshes will, however, be very cumbersome to calculate as they will contain all the points of the concerned snapshots and all the intersecting points of all the connectivities. Furthermore, PFEM data is particularly difficult to reduce with classical linear POD. Indeed, PFEM is commonly used for transport-dominated problems, the most common being free surface problems [1,7,12,13], which corresponds to the advection of an interface. It is well known that these problems do not have a good Kolmogorov-N width decay and require many modes to reduce [32]. In this section, we will introduce mappings to tackle both of these problems as well as two methods that directly exploit the PFEM Lagrangian displacement to construct the necessary mapped solution manifolds mentioned in Section 2.2.

### 4.1. The benefits of mapped manifolds

Let us first introduce a classical 1D example of slow Kolmogorov N-width decay that will illustrate our methods: the linear advection of a jump discontinuity.

$$ \begin{aligned} & \partial_t u(x, t) + \mu \cdot \partial_x u(x, t) = 0 \quad x, t \in [0, 1] \\ & u(x, 0) = 0, \quad u(0, t) = 1. \end{aligned} \quad (15) $$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Solution of Eq. (15) for $\mu = 1$.</figcaption></figure>

Let us consider the solution manifold, also depicted in Fig. 4:

$$ \mathcal{K} = \{u(x, t) | t \in [0, 1]\} \subset L^2([0, 1]) = \{1_{[-1,0]}(x - \mu t) | t \in [0, 1]\}, $$

where $1_{[-1,0]}$ is the indicator function of the interval $[-1, 0]$, equal to one within the interval and zero outside. We know from Ohlberger and Rave [32] that $d_N(\mathcal{K}) \propto N^{-1/2}$, which is a lot worse than the proven exponential decay for parametric elliptic problems. To improve the reducibility, we introduce the mappings $\forall t \in [0, 1]$:

$$ X_t := \mathbb{R} \to \mathbb{R} $$

$$ \xi \to \xi + \mu t, $$

and consider the mapped solution manifold:

$$ \bar{\mathcal{K}} = \{u(t) \circ X_t | t \in [0, 1]\} = \{1_{[-1,0]}\}. $$

Finally, $d_N(\bar{\mathcal{K}}) = 0$ as $\bar{\mathcal{K}}$ is of dimension 1. One can retrieve the original solutions by applying the reversed translations, also known as the inversed mappings.

Applying this in the context of our original PFEM problems, we are looking for $n$ reference configuration $\Omega_0^i \subseteq \Omega \subset \mathbb{R}^d$ and $n$ diffeomorphic mappings

$$ X_{t_i} := \Omega_0^i \to \Omega_{t_i} $$

$$ \xi \to x, $$

such that

$$ \mathcal{K} = \{u(t_i) \circ X_{t_i} | i \in \llbracket 1, n \rrbracket\} \subset L^2(\Omega), $$

(17)

has a better Kolmogorov $n$-width decay than the original problem. We then call $\Xi_{t_i}: \Omega_{t_i}^i \to \Omega_0^i$ the inverse mapping required to reconstruct the reduced snapshots.

## 4.2. Using the Lagrangian displacement as a reduction map

There exists a range of so-called velocity-driven registration methods [38–40] that exploit the fact that if you have a time-dependent field $v(t)$, sufficiently smooth such that $v \in L^1(I, C_0(\mathbb{R}^d, \mathbb{R}^d))$, then

$$ \left\{ \begin{array}{l} X^v(\xi; t = 0) = \xi \\ \partial_t X^v(t) = v(t) \circ X^v(t), \end{array} \right. $$

(18)

has one unique continuous solution $X^v$ called the flow of $v$ and $\forall t \in I$, $X^v(t)$ is a diffeomorphism and can be used as a map to construct a solution manifold. These methods will not be used here as they require expensive adjoint optimization problems to determine $v$. However, we can take inspiration from this and identify that at the heart of PFEM, we are already solving for a discrete version of equation Eq. (18) with the fluid velocity field and the Lagrangian displacement field, as we saw in Section 2.

Crucially, we only know the displacement along the particle trajectories, which are not continuous due to PFEM remeshing. Let us consider a particle $p \in \llbracket 1, N_e \rrbracket$ appearing at time $t_{a^p}$ and disappearing at time $t_{d^p}$, with $N_e$ the total number of particles in the simulation and $a^p, d^p \in \llbracket 1, n \rrbracket$. The PFEM data only gives us a partial particle displacement path

$$ X_{\text{Lag}} := \llbracket 1, N_e \rrbracket \times \llbracket a^p, d^p \rrbracket \to \Omega_{t_i} $$

(19)

$$ (p, i) \to x. $$

Identifying $X_{\text{Lag}}$ from pre-existing PFEM data can be tricky as the unique particle indexing is not necessarily tracked and would require a nonstandard data format in the original solver. Post-simulation identification had to be performed for the examples in our
<figure><img src="image_7.png" /><figcaption>Fig. 5. PFEM resolution of Eq. (15) for $\mu = 1$, before (left) and after (right) $X_{Lag}$ identification.</figcaption></figure>

<figure><img src="image_8.png" /></figure>

### Method 1: a supermesh manifold

<figure><img src="image_9.png" /></figure>

### Method 2: a particle manifold

<figure><img src="image_10.png" /></figure>

and the full connectivities

<figure><img src="image_11.png" /><figcaption>Fig. 6. The two reduction strategies, for a simplified case with 3 snapshots and one connectivity change. Method 1 requires 1 supermesh connectivity and 2 modes to fully rebuild the solution. Method 2 only requires 1 mode but the full original connectivities.</figcaption></figure>

result Section 5. An example of $X_{Lag}$ identification is given in Fig. 5 for a PFEM simulation solving equation Eq. (15). The color scheme is the one of Fig. 4, and unique particle paths are identified with the colored arrows. The left figure is the result of using the original PFEM nodal indexing, as any new particle is given the index of the last particle that was removed. The right figure corresponds to the correct identification.

We will now present two methods to construct solution manifolds using $X_{Lag}$:

- Method 1 will attempt to build n diffeomorphisms $X_{t_i}$ and perform the reduction on one supermesh, reducing the full functional space.

- Method 2 will take the unique particle indexing from $X_{Lag}$ and perform the reduction along the particle trajectories, not reducing the connectivities. Therefore, in case of addition and removal of nodes, we have missing data and get a gappy solution manifold.

The two methods are illustrated in Fig. 6, with an example with no addition and removal of nodes. Let us now look at the general implementation of both methods.

### 4.2.1. Method 1, reducing on a supermesh:

Let us first construct the required mappings $X_{t_i}$, $\forall i \in [1, n]$ Eq. (16) from $X_{Lag}$ Eq. (19). For that, we need to extend $X_{Lag}$ to time $t_1 = 0$ for every particle; we then call it the extended particle path $X_{Lag}^e$. This entails giving a displacement value to particles that have yet to exist and constructing an initial configuration at time zero. To do so, we will re-use the interpolation used by PFEM when a given particle $p$ appears at time $t_{ap} > 0$. The interpolation involves a set of local neighbor particles $A^p \subset [1, N_e]$ and
the local weights {$\alpha_j^p$}$_{j \in A^p} \subset [0, 1]$ s.t. $\sum_{j \in A^p} \alpha_j^p = 1$ and $X_{Lag}(p, t_{ap}) = \sum_{j \in A^p} \alpha_j^p X_{Lag}(j, t_{ap})$, we then define the extended mapping $\forall p \in \{p \in \llbracket 1, N_e \rrbracket | t_{ap} > 0\}, \forall i \in \llbracket 1, a^p - 1 \rrbracket$

$$X_{Lag}^{e}(p, i) = \sum_{j \in A^{p}} \alpha_{j}^{p} X_{Lag}(j, i). \qquad (20)$$

In practice, one has to start to extend the particles appearing the earliest and finish with the ones that appear last to prevent any dependencies. The inverse mappings can now be computed

$$\Xi_{t_i} := \Omega_{t_i} \rightarrow \Omega \\ x \rightarrow \sum_{p=1}^{N_e} \mathcal{N}^h(x, t_i)|_p^e X_{Lag}^e(p, 1), \qquad (21)$$

where $\mathcal{N}^h(x, t_i)|_p^e$ gives the shape function at $t_i$ of the particle with unique indexing $p$, the shape function can be zero if the particle $p$ is not present at time $t_i$. The reference configurations are now $\Omega_0^l = \mathrm{Im}(\Xi_{t_i})$, and the direct mappings $X_{t_i}$ are $\Xi_{t_i}^{-1}: \Omega_0^l \rightarrow \Omega_{t_i}$. The last and most important thing to verify is the bijectivity of all $X_{t_i}$, this can be done by verifying the positivity of the jacobian of all $\Xi_{t_i}$, we need $\forall i \in \llbracket 1, n \rrbracket, \forall x \in \Omega_{t_i}$:

$$\mathfrak{F}(\Xi_{t_i})(x) = \det(\nabla \Xi_{t_i})(x) > 0, \qquad (22)$$

this cannot be enforced and has to be done element by element. Referring back to the upper part of Fig. 6, a non-bijective case would be if any of the mesh elements flipped over after the application of $\Xi_{t_i}$. Assuming we verified that all the $X_{t_i}$ are indeed bijective, we can now map back our PFEM data into a single infinite dimensional manifold.

$$\mathcal{K} = \{u^{h_i} \circ X_{t_i} | i \in \llbracket 1, n \rrbracket\} \\ = \{\mathcal{N}^h(\Xi_{t_i}^{-1}(\xi), t_i)^T u^{h_i} | i \in \llbracket 1, n \rrbracket\} \subset L^2(\Omega). \qquad (23)$$

This now requires supermeshes to reduce; however, they should be relatively easy to compute as if a given particle is present at two different time steps, they now coincide in the common mapped manifold. Therefore, most nodes coincide and we mainly only observe connectivity changes. This simplifies the supermesh calculation significantly, we should expect the number of nodes in the total supermesh of all snapshots to be reasonable. The summary of Method 1 in cases where all the $X_{t_i}$ are bijective is given in Algorithm 1. We kept both the optionsof using or not using the supermesh for the reduction. In the latter, the supermesh is only used to interpolate the data and conserve the history of mesh changes. Generally, we recommend not using the mass matrix in the reduction only if you observe high point concentration in certain areas of the supermesh.

In 1D, the $X_{t_i}$ are always bijective because the extension Eq. (20) ensures no particle paths cross each other. We can apply this procedure to the PFEM example from Fig. 5; the results can be seen in Fig. 7. The resulting mapped manifold can be expressed on a single 1D supermesh containing all the mapped mesh positions (here, seven nodes). Interpolating the mapped solution manifold seen in Fig. 4 on this supermesh, one can observe two regimes: one for $t \le 0.7$ with light blue for $0 \le x \le 0.1$ a ramp for $0.1 \le x \le 0.35$ and dark blue for $0.35 \le x \le 1$; one for $t > 0.7$ with light blue for $0 \le x \le 0.1$ and a ramp for $0.1 \le x \le 1$. We would only require two modes to reconstruct the original solution fully. Note that if we do not use a mapped manifold, we require all fifteen snapshots as modes to reconstruct the solution exactly.

This method can only be used in bijective cases. An extension to this method for non-bijective cases relying on a piecewise linear approximation is presented in Appendix A.

### 4.2.2. Method 2, reducing along the particle trajectories:

The second idea is to reduce the data only along the particle trajectories without considering the connectivities. This means the original PFEM connectivities will not be reduced and will be stored to reconstruct the data. We can rearrange our particle data to form some gappy solution manifold

$$\mathcal{K} = \{u^e(t_i) | i \in \llbracket 1, n \rrbracket\} \subset \mathbb{R}^{N_e} \qquad (24)$$

<figure><img src="image_12.png" /><figcaption>Fig. 7. Method 1 on the PFEM resolution of Eq. (15), $X_{Lag}$ prolongation (left) and mapped solution manifold (right).</figcaption></figure>
**Algorithm 1 Reducing on a supermesh.**

**Input** - {$u^{h_i}, x^{h_i} \in \mathbb{R}^{N_i}, \mathcal{N}^{h_i} \in V^{N_i}, \forall i \in [1, n]$}, $N \in \mathbb{N}^*$ the required number of modes and ReduceWithMassMatrix a boolean indicating if the reduction should include the supermesh.

With $u^{h_i}$ the snapshots of the field of interest to reduce and $x^{h_i}$ the Lagrangian displacement, both computed with PFEM.

**Initialization** - Identification of $X_{Lag}$ (19), $\forall p \in P_a \{A^p, \{\alpha_j^p\}_{j \in A^p}\}$

This initialization can be skipped if this information has been saved during the PFEM solve.

Sort $P_a = \{p \in [1, N_e]\mid t_{ap} > 0\}$ by increasing $t_{ap}$

Compute $X_{Lag}^e(p, i) = \sum_{j \in A^p} \alpha_j^p X_{Lag}(j, i)$ for $p \in P$ and $0 < t < t_{ap}$ (20)

$M^+ : M_1, \mathcal{K} = \{u^{h_1}(x)\}$

**for** $i \in [2, n]$ **do**

Calculate $\Xi_{t_i}$ (A.1)

Verify that $\Xi_{t_i}$ is bijective (22)

Calculate mapped mesh $M$ with $\Xi_{t_i}$

Update the supermesh $M^+$, combining both $M$ and $M^+$, see Section 3.2

Add $u^{h_i} \circ \Xi_{t_i}^{-1}$ to $\mathcal{K}$

**end for**

Let $l \in \mathbb{N}^*$ be the number of nodes of $M^+$

Project without committing additional errors all $f \in \mathcal{K}$ onto $M^+$ and add them to $Y \in \mathbb{R}^{l \times n}$

**if** ReduceWithMassMatrix **then**

Compute $M^+ \in \mathbb{R}^{l \times l}$ the symmetric, positive definite mass matrix of $M^+$

Compute the eigenvalue decomposition $\{(\lambda_j, \phi_j)\}_{j=1}^{n_k} = \text{Eig}(Y^T M^+ Y)$

**else**

Compute the eigenvalue decomposition $\{(\lambda_j, \phi_j)\}_{j=1}^{n_k} = \text{Eig}(Y^T Y)$

**end if**

Compute the space modes with $\psi_j = \frac{1}{\sqrt{\lambda_j}}Y\phi_j$ for $1 \le j \le n$

Compute the reducibility with $1 - \varepsilon_{POD}^2(N) = 1 - \sum_{j=1}^N \frac{\lambda_j}{\sum_{j=1}^n \lambda_j}$

Reconstruct partial solutions and map them back to every configuration with $X_{t_i}$ and compute all $\varepsilon_S^2(N)$

<figure><img src="image_13.png" /><figcaption>**Fig. 8.** Concept of method 2 on the PFEM resolution of Eq. (15), Gappy expansion along the particle trajectories (left) and the manifold in $\mathbb{R}^{N_s}$ (right).</figcaption></figure>

where $\forall p \in [1, N_e]$, $\forall i \in [a^p, d^p]$ $u^e(t_i)|_p = u^{h_i}|_p$, with $\cdot|_p$ designating the component corresponding to particle $p$ in either format. We can illustrate our approach in Fig. 8 on the example used up to now, the idea is to extend the values for every particle trajectory for all time steps. We have unknown values in gray; if particles 3 and 4 are completed with dark blue values and particles 6 and 7 are completed with cyan values, only one mode is required to reconstruct the solution. This section will present two algorithms to systematically complete $u^e(t_i)$: an immediate geometrical interpolation and an iterative gappy POD approach. The second approach can be seen as a generalization of the first one, as one can use the geometrical interpolation as the initial guess for the iterative approach.

**A geometrical interpolation.** The first idea will be similar to the $X_{Lag}$ extension made in method 1 Eq. (20). Once again, we will look at $p \in P_a = \{p \in [1, N_e]\mid t_{ap} > 0\}$ and use the PFEM addition interpolation with weights $\{\alpha_j^p\}_{j \in A^p} \subset [0, 1]$ and after sorting $P_a$ from the earliest addition to the last, we can $\forall p \in P_a, \forall i \in [1, a^p - 1]$

$$u^e(t_i)|_p = \sum_{j \in A^p} \alpha_j^p u^{h_i}|_j. \qquad (25)$$

However, for $p \in P_d = \{p \in [1, N_e]\mid t_{dp} < n\}$, PFEM does not provide any exact interpolation, and more work will be required to determine the interpolation weights, see Fig. 9.
<figure><img src="image_14.png" /><figcaption>Fig. 9. Example of removal of particle p at time $t_{d^p}$ and the resulting mesh at time $t_{d^{p+1}}$, assuming the points do not move.</figcaption></figure>

To determine the set of neighbors $D_p \subset [1, N_e]$ and their associated weights $\{\delta_j^p\}_{j \in D^p} \subset [0, 1]$, we had multiple options:

* Looking at the configuration at time $t_{d^p}$, the set of neighbors $D_p = \{v_j\}_{j=1}^{|D_p|} | D_p$ that share an edge with p constitutes a polytope. A large range of generalized barycentric interpolation methods exist that exploit the exact geometry of the polytope [45,46]. We kept it simpler and took the equal weights $\delta_j^p = \frac{1}{|D^p|}, \forall j \in D^p$, approximating p with the barycenter $p_B$.

* Projecting the particle p one step in the future at time $t_{d^{p+1}}$, find the triangle in which p is situated, take as neighbors the three vertices of this triangle, and as weights the exact barycentric weights calculated from the relative areas.

To prevent any dependencies, we have to sort $P_d$ from the latest removal to the first, and we have $\forall p \in P_d, \forall i \in \llbracket d^p + 1, n \rrbracket$

$$ \mathbf{u}^e(t_i)|_p = \sum_{j \in D^p} \delta_j^p \mathbf{u}^{h_i}|_j. \qquad (26) $$

The full summary of Method 2.1 is given in Algorithm 2.

**Algorithm 2 Reducing on particle trajectories with a geometrical interpolation.**

**Input** - $(\mathbf{u}^{h_i}, \mathbf{x}^{h_i}) \in \mathbb{R}^{N_i}, \mathcal{N}^{h_i} \in V_i^{N_i}, \forall i \in [1, n]$ and $N \in \mathbb{N}^*$ the required number of modes, with $\mathbf{u}^{h_i}$ the snapshots of the field of interest to reduce and $\mathbf{x}^{h_i}$ the Lagrangian displacement, both computed with PFEM.
**Initialization** - Identification of $X_{Lag}$ (19), $\forall p \in P_a \{A^p, \{\alpha_j^p\}_{j \in A^p}\}$ and $\forall p \in P_d \{D^p, \{\delta_j^p\}_{j \in D^p}\}$
This initialization can be skipped if this information had been saved during the PFEM solve.

Let $\mathcal{K} = \{0 | i \in [1, n]\} \subset \mathbb{R}^{N_e}$
With $X_{Lag}$, add the values of $\mathbf{u}^{h_i}$ to $\mathcal{K}$ for $1 \le i \le n$

Sort $P_a = \{p \in [1, N_e] | t_{a^p} > 0\}$ by increasing $t_{a^p}$

Compute $\mathbf{u}^e(t_i)|_p = \sum_{j \in A^p} \alpha_j^p \mathbf{u}^{h_i}|_j$ for $p \in P_a$ (25)

Sort $P_d = \{p \in [1, N_e] | t_{d^p} < n\}$ by decreasing $t_{d^p}$

Compute $\mathbf{u}^e(t_i)|_p = \sum_{j \in D^p} \delta_j^p \mathbf{u}^{h_i}|_j$ for $p \in P_d$ (26)

Combine the content of $\mathcal{K}$ in a snapshot matrix $Y \in \mathbb{R}^{N_e \times n}$

Compute the eigenvalue decomposition $\{(\lambda_j, \phi_j)\}_{j=1}^n = \text{Eig}(Y^T Y)$

Compute the first $N$ space modes with $\psi_j = \frac{1}{\sqrt{\lambda_j}} Y \phi_j$ for $1 \le j \le N$

Compute the reducibility with $1 - \varepsilon_{POD}^2(N) = 1 - \sum_{j=1}^N \lambda_j / \sum_{j=1}^N \lambda_j$

Reconstruct partial solutions, express them in every configuration with $X_{Lag}$ and compute all $\varepsilon_S^2(N)$

An iterative gappy POD. We can go one step further and compare our technique to the state-of-the-art gappy POD approach developed by Everson and Sirovich [41]. It was developed to construct a POD basis with an incomplete—also known as gappy—snapshot set. At the heart of the method, there is an iterative procedure optimizing the missing values to improve the reducibility of the snapshots. It works as follows:

1. An initial guess is given for the unknown values in the snapshot set. Typically, the mean over all available data for a given data point in the snapshot set. For our examples, one can also use the previously mentioned geometrical interpolation as an initial guess.

2. A certain amount of POD modes $\{\psi_j\}_{j=1}^N$ are computed from this completed data set.

3. The time modes are then recomputed for every snapshot, minimizing this time the $L^2$ error over the available data. It is called the reconstruction phase and provides new values to the missing data.

4. The snapshots are updated with this new missing data, and we go back to step 2 until convergence.

We adapted this algorithm to our particle problem. It has the benefit of not requiring the identification of the interpolation weights, but it is iterative instead. The exact outline is given in Algorithm 3.

A few modifications have been made compared to the original algorithm found in Everson and Sirovich [41], heavily inspired by what was developed in [42,43]. First, the $L^2$ minimizations are done using $M_{h_i}$, representing the actual mesh configuration at time
$t_i$, This ensures that the reconstruction is optimized in relation to the full reduction error $e_N^2$ Eq. (5) and not just to the reducibility in $\mathbb{R}^{N_e}$. Furthermore, our gappy POD iterations progressively contain more and more modes. Performing the iterations directly with $N$ modes, as indicated in Everson and Sirovich [41], would be optimal only for $N$ modes and not for smaller values. Instead, with the progressive approach, we have near-optimal values of $e_r^2$ at any value $1 \le r \le N$; therefore, you can compare the full error decrease with other methods. This progressive approach is also found in Venturi and Karniadakis [42], Gunes et al. [43]. We fully acknowledge that Algorithm 3 is the most expensive version of the gappy POD approach, as simplifications exist in the literature [43] regarding the iterative process. However, it remains the most precise and will give us a satisfying lower bound for our reduction analysis.

**Algorithm 3 Reducing on particle trajectories with the gappy POD.**

**Input** - $(\{u^{h_i}, x^{h_i}\}_{i \in \mathbb{N}^N}, \mathcal{N}^{h_i} \in V_i^{N_i}, \forall i \in [1, n])$, $N \in \mathbb{N}^*$ the required number of modes and $\text{Tol} \in \mathbb{R}$ the Gappy POD tolerance. With $u^{h_i}$ the snapshots of the field of interest to reduce and $x^{h_i}$ the Lagrangian displacement, both computed with PFEM.

**Initialization** - Identification of $X_{Lag}$ (19)

This initialization can be skipped if this information had been saved during the PFEM solve.

Let $\mathcal{K} = \{0 | i \in [1, n]\} \subset \mathbb{R}^{N_e}$

With $X_{Lag}$, add the values of $u^{h_i}$ to $\mathcal{K}$ for $1 \le i \le n$

First mean value guess $u^e(t_i)|_p = \frac{1}{d^p-a^{p+1}} \sum_{j=a^p}^{d^p} u^e(t_j)|_p$ for $1 \le p \le N_e$ and for $1 \le i < a^p, d^p < i \le n$, for $k \in [1, R]$ do

Combine the content of $\mathcal{K}$ in a snapshot matrix $Y \in \mathbb{R}^{N_e \times n}$

Compute the eigenvalue decomposition $\{(\lambda_j, \phi_j)\}_{j=1}^n = \text{Eig}(Y^T Y)$

Compute the first $k$ space modes with $\psi_j = \frac{1}{\sqrt{\lambda_j}} Y \phi_j$ for $1 \le j \le k$

Let $\lambda_{Old} = \lambda_k$ and $\tau = 1$

while $\tau > \text{Tol}$ do

for $i \in [1, n]$ do

Let $\mathbf{M}_{h_i} \in \mathbb{R}^{N_i \times N_i}$ be the symmetric, positive mass matrix related to time $t_i$

Let $W \in \mathbb{R}^{k \times k}$ s.t. $W_{lm} = (\psi_l^{h_i})^T \mathbf{M}_{h_i} \psi_m^{h_i}$ and $F \in \mathbb{R}^k$ s.t. $F_i = (u^{h_i})^T \mathbf{M}_{h_i} \psi_l^{h_i}$

with $\psi^{h_i}$ being the mode $\psi$ expressed on only the nodes present at time $i$.

Solve for $\mathbf{\Phi} = \{\phi_j\}_{j=1}^k$ in $W \mathbf{\Phi} = F$

Reconstruction of the snapshot $u^e(t_i)|_p = \sum_{j=1}^k \phi_j \psi_j|_p$ for $1 \le p \le N_e$ not present at time $t_i$,

end for

Combine the content of $\mathcal{K}$ in a new snapshot matrix $Y \in \mathbb{R}^{N_e \times n}$

Compute the eigenvalue decomposition $\{(\lambda_j, \phi_j)\}_{j=1}^n = \text{Eig}(Y^T Y)$

Compute the first $k$ space modes with $\psi_j = \frac{1}{\sqrt{\lambda_j}} Y \phi_j$ for $1 \le j \le k$

Let $\tau = |\lambda_k - \lambda_{Old}|/|\lambda_{Old}|$ and $\lambda_{Old} = \lambda_k$

end while

Compute the near optimal reducibility with $1 - \mathcal{E}_{\text{POD}}^2(k) = 1 - \sum_{j=1}^k \lambda_j / \sum_{j=1}^n \lambda_j$

Reconstruct partial solutions, express them in every configuration with $X_{Lag}$

and compute near optimal $e_S^2(k)$

end for

## 5. PFEM reduction examples

We present the results obtained using the two methods introduced in Sections 4.2.1 and 4.2.2, demonstrated through a series of test cases. We begin with a 1D example based on the viscous Burgers' equation implemented in Python. Following this, we examine three 2D cases of increasing complexity derived from our in-house explicit PFEM solver written in FORTRAN.

### 5.1. Example 1: 1D resolution of the viscous Burgers' equation

Let us start with a demanding 1D example containing sharp discontinuities, for which we have an analytical solution. We will solve the viscous Burgers' equation with a Lagrangian solver, including a simple remeshing procedure.

$$ \partial_t u(x,t) + u(x,t) \partial_x u(x,t) = v \partial_{xx} u(x,t) \quad x \in [0, 1], \ 0 \le t < \infty, \qquad (27) $$

with the kinematic viscosity $v = 0.003$, the analytical solution is:

$$ u(x,t) = \frac{e^{-\omega_1(x)} + 0.5e^{-\omega_2(x,t)} + 0.1e^{-\omega_3(x,t)}}{e^{-\omega_1(x)} + e^{-\omega_2(x,t)} + e^{-\omega_3(x,t)}} \qquad (28) $$

$$ \text{with } \omega_1(x) = \frac{0.5}{v}(x - 0.375), \quad \omega_2(x,t) = \frac{0.25}{v}(x - 0.5 + 0.75t), \quad \omega_3(x,t) = \frac{0.05}{v}(x - 0.5 + 4.95t). $$
<figure><img src="image_15.png" /><figcaption>(a) The solution of equation (27) for 7 values of t</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Identified particle paths $X_{Lag}$ for only 20 particles out of the $N_e = 136$ particles used for the Lagrangian solve of equation (27)</figcaption></figure>

<caption>Fig. 10. The 1D Burgers' equation example.</caption>

This specific example was taken from Babu et al. [47]; the evolving solution can be seen in Fig. 10(a). In a Lagrangian setting, we lose the convection term and work with a material derivative instead,

$$ 
\begin{aligned}
D_t u(X(\xi, t), t) &= \nu \partial_{x_x} u(X(\xi, t), t) \quad X \in [0, 1], 0 \le t < \infty \\
\partial_t X(\xi, t) &= u(X(\xi, t), t),
\end{aligned} 
\qquad (29) $$

where the displacement $X(\xi, t)$ is described with a moving Lagrangian mesh that can be assimilated to particle trajectories. The goal is to mimic the PFEM procedure from Section 2 to solve our problem. Eq. (29) is therefore solved with standard P1 finite elements, an implicit time integration for both the speed and the mesh, and the following remeshing rules:

* If the size of an element is smaller than $0.5\bar{h}$, with $\bar{h}$ the mean mesh size, then remove one of the associated mesh points, taking care that none of the two boundary nodes are removed.

* If the size of an element is larger than $2\bar{h}$, add a node between the two and interpolate the solution accordingly.

These two procedures can happen independently; therefore, the number of nodes at a given time is not necessarily constant. We start with $N_0 = 75$ particles and study an interval of time $I = [0, 1]$ with $\Delta t = 0.01$, which gives us 100 snapshots to work with. Because of the remeshing, one can identify $N_e = 136$ unique particles in the whole simulation (see Fig. 10(b)). We can calculate an analytical error defined as:

$$ \epsilon_A^2 = \frac{\sum_{i=1}^{n} \|u^h(t_i) - u_a(t_i)\|_{L^2(\Omega_i)}^2}{\sum_{i=1}^{n} \|u_a(t_i)\|_{L^2(\Omega_i)}^2}. \qquad (30) $$

In this specific case, $\epsilon_A$ is of the order of $10^{-2}$.

The reduction performances of the proposed methods are compared in Fig. 11. All three graphs show monotonic decreasing error functions, with the reducibility graph being the smoothest. This is precisely what we expect as we directly reduce the manifold $\mathcal{K}$. The optimality of the POD process ensures both a smooth error decrease and a convergence to zero, which is reached if we take as many modes as there are snapshots ($r = 100$). Note that all the errors with respect to the snapshots also decrease to zero as we do not introduce any additional errors with the proposed methods, and we circumvent all potential interpolation errors. The analytical error does not decrease to zero as it converges to the numerical error of the original solver.

The only technique that attempts to use no mapping (the blue curve) performs the worst of all the techniques. This should come as no surprise as we are looking at a problem with transported sharp discontinuities that are, by nature, hard to reduce (see Section 4.1). Not only that, but it is also the least practical to use as the modes $\Psi$ are considered infinite-dimensional and would require to be interpolated onto a common mesh. This would, in turn, introduce additional errors in both the snapshot and analytical error graphs.
<figure><img src="image_17.png" /><figcaption>Fig. 11. Overview of the performance of our reduction methods on the 1D Burgers' equation example.</figcaption></figure>

For Method 1, we see two very different performances for the two options regarding the supermesh. This is due to the structure of the supermesh built for those examples: with the mesh prolongations Eq. (20), our supermesh $X_{super} : [1, N_e] \to [0, 1]$ is:

$$ 
\begin{aligned}
X_{super}(p) &= h(p-1) \quad \text{for } 1 \le p \le N_0 \\
X_{super}(p) &= h/2^{p-N_0} \quad \text{for } N_0 < p \le N_e,
\end{aligned} 
\qquad (31) $$

with $h = 1/(N_0 - 1)$, the mesh size of the original mesh. This means that all the nodes not present at time 0 get positioned between 0 and $h$ in the supermesh. The proximity of all these mesh points explains the poor performance of the option in green and why it is not beneficial to use the supermesh operator in the reduction. Indeed, an error decrease in the solution manifold does not systematically translate to a similar decrease in the other errors. Note that the other option for method 1 in orange remains competitive compared to all the other methods. Furthermore, it only requires one supermesh to reconstruct the entire solution without any additional errors.

In Method 2, we pay the price to keep all the connectivity changes in memory and only reduce along the gappy particle trajectory. We have $N_e = 136$ total particles in the simulation; on average, 70 are active each time step. This translates to a high percentage of missing information - the so-called gappyness factor - of 48%. We can observe a significant increase in reduction performance, going from very poor when adding zeros to excellent with the Gappy POD. The high gappyness factor can explain the impressive performance of the Gappy POD. However, the geometrical interpolation remains an immediate and straightforward construction, while the Gappy POD is a costly iterative algorithm.

In terms of computational costs, the Gappy POD is by far the most expensive method. Both the Arbitrary FE POD method and the method 1 are fairly cheap as supermeshes are very simple to construct in 1D and the geometrical interpolation is the cheapest of them all.

A typical rule of thumb in reduction theory is to keep enough modes to have an error $\epsilon < 1\%$. In the 1D example given here, it is the case with less than 10 modes for three of our methods, and all three are more than satisfactory.
<figure><img src="image_18.png" /><figcaption>Fig. 12. The sloshing configuration and the free surface at different times.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 13. Reducibility and total error with respect to the snapshots of the pressure and velocity fields for the no remeshing sloshing example.</figcaption></figure>

## 5.2. Example 2: 2D sloshing

Let us start the 2D examples with a simple sloshing with no remeshing; the configuration can be seen in Fig. 12.

We choose a length $L = 2.5$ m, a density of $\rho = 1000 \text{ kg m}^{-3}$, and a dynamic viscosity of $\mu = 1$ Pa·s. The walls are given slip-free boundary conditions (except the corners that remain fixed) to preserve the Lagrangian mappings' bijectivity in method 1. Indeed, we consistently observed that fully constrained boundaries are more likely to lead to non-bijective mappings, as elements for which you have a single fixed node can easily flip over and break the bijectivity. We chose to look at one oscillation period, simulated the first 1.8 s, and sampled 90 snapshots. We have decided to study the reducibility of three fields of interest: the pressure and both the vertical and horizontal velocity fields. The reducibility of the different manifolds and the total error made reconstructing the reduced solutions with the different methods can be found in Fig. 13.

The cost of the methods using supermeshes now become non negligible in 2D. Especially for the Arbitrary FE case in which you need to compute supermeshes of non matching meshes. For the mapped supermeshes of Method 1, one can implement techiques to compute them efficiently as most points match after the mapping.

The first thing that can be observed in the results is the difference in reducibility of the pressure vs the two velocity fields. The pressure field seems less affected positively by the mapping approaches in terms of reducibility. This can be explained by the fact that the pressure, contrary to the velocities, is at equilibrium on the free surface. Therefore, the geometric registration performed by
<figure><img src="image_20.png" /><figcaption>Fig. 14. The newtonian drop [13], its initial position and the free surface at different times.</figcaption></figure>

our Lagrangian mapping will not help its reducibility. As predicted, the arbitrary FE approach also remains the worst option for all fields of interest. Both of those observations will be consistent for all our surface fluid examples.

We start with 1610 nodes, and the computed supermesh for method 1 is of size 1892. In this case, for method 1, it is better to consider the supermesh mass matrix in the reduction as we get slightly better results. Method 2, in this no-remeshing example, is a direct POD of the non-gappy data and shows results similar to those of method 1. Both techniques only require 2, 5, and 6 modes to reconstruct the pressure, the horizontal and vertical velocity, respectively, with a precision of $\epsilon = 1\%$.

### 5.3. Example 3: 2D newtonian drop

The third example is the impact of a 2D drop on a rigid plane, a typical benchmark for free surface fluid flows. It was first introduced in Tomé et al. [48] and has recently been implemented in PFEM [13]. A representation of the initial and deformed configurations can be found in Fig. 14. We have $D = 0.01$ m, a drop height $H = 0.04$ m, an initial velocity of $U = -1$ m s$^{-1}$, a density of $\rho = 1000$ kg m$^{-3}$, a dynamic viscosity of $\mu = 4$ Pa s, and gravity is considered. The rigid plane is considered a no-slip boundary condition. The simulation lasted 0.1s, and we decided to sample 201 snapshots for the reduction procedure. Note that routine PFEM remeshing procedures [1] were used this time. The reducibility of the different manifolds can be found in Fig. 15, and the total error made reconstructing the reduced solutions with the different methods can be found in Fig. 16.

We first observe, once again, that the pressure field is much harder to reduce than the two velocity components. The drop impact generates a travelling discontinuity in both pressure and vertical velocity, which affects the reducibility negatively. We are working with an initial mesh of 1783 nodes and a total number of particles $N_e = 1849$, translating to a gappyness factor of 3.6%. This explains the tiny difference between adding zeros and the other options for method 2. Note that the gappy POD and geometrical interpolation results are quasi-identical and show that the expensive gappy POD is not required when dealing with small gappyness values.

For method 1, we had to resort to the piecewise linear decomposition described in Appendix A as the mappings were no longer bijective after some time, Eq. (22) not being verified with locally flipped triangles. A single mapped manifold was no longer sufficient, and we resorted to a decomposition with three bijective mappings instead. The three mapped manifolds lead to three supermeshes instead of one, corresponding to the configurations at times $t_1, t_{82}$, and $t_{98}$. These three time steps correspond to the time at which the bijectivity is broken. The free surfaces at these time instances can be seen in Fig. 14. A POD is then performed for each mapped manifold corresponding to the three time intervals/supermeshes, hence the groups of three curves for method 1 in Fig. 15. The results also show that performing a reduction with the supermesh mass matrices is more interesting. This is most likely due to the heterogeneity in element sizes within the supermeshes. The reconstruction of the complete solution is then made by taking modes from the manifold that have the most effect on the total error. This piecewise approach shows comparable performances to method 2 for the same number of total modes. Note that each supermesh contains a different number of nodes: 3624 nodes for the first interval, 2217 in the second, and 2548 in the third. While this means that the space modes are larger than the ones in method 2, for which we had $N_e = 1849$, we do not have to save the connectivities of each time step, but only the three supermesh connectivities.

Overall, both mapped methods are good strategies for this example, and have advantages and disadvantages. Aiming for $\epsilon < 1\%$ for both methods, we need roughly 15, 10 and, 5 modes for the pressure, horizontal, and vertical velocity, respectively, which seems reasonable.

### 5.4. Example 4: 2D viscous dam break

Finally, let us look at a typical demanding validation case for PFEM simulations, the dam break test. The setup is given in Fig. 17, with a block of liquid affected by gravity crashing into the walls of a box. We have a length $L = 1$ m, a density of $\rho = 1000$ kg m$^{-3}$, and a dynamic viscosity of $\mu = 5$ Pa s. The walls are considered as no-slip boundary conditions. The simulation lasted 3 seconds, we sampled 300 snapshots in time, and the evolution of the free surface for the first third of the sample time can also be observed in
<figure><img src="image_21.png" /><figcaption>**Fig. 15.** Reducibility of the pressure and velocity fields for the newtonian drop example.</figcaption></figure>

Fig. 17. The reconstructed total error can be found in Fig. 18 for both the velocity and pressure fields. First and foremost, in this very demanding example, method 1 cannot be used. Indeed, the large deformations combined with the non-slip boundary conditions on the borders make the mappings Eq. (21) non-bijective after almost every time step as we observe triangles flipping in the proximity of the boundaries. If we want to attempt the reduction on multiple supermeshes described in Appendix A, we would require a similar number of supermeshes as we have snapshots. This is not a viable solution and defeats the purpose of a mapped reduction approach.

On the other hand, Method 2 shows very similar results to previous cases, with the gappy POD results being the best. We are working with an initial mesh of 1123 nodes and a total number of particles $N_e = 1219$, which translates to a gappyness factor of 7.9%, hence the differences between the gappy and the other two methods. However, even in our best-case scenario, we require many modes to reconstruct the solution with appropriate precision. Indeed, with our usual error criteria, we need 150, 60, and 80 modes for the pressure, horizontal, and vertical speed, respectively. This is too many compared to the original snapshot number and shows that, even with the methods proposed in this article, some problems are too complex to be reduced efficiently. As a matter of fact, we recovered a reasonably reducible solution when we took only the initial 50 time steps, corresponding to a viscous slumping. It is the discontinuities generated by the impact of the fluid on the walls and the subsequent splashing that are too complex to be reduced efficiently.

In terms of costs, the Arbitrary FE case becomes too expensive to be worth computing with our current implementation. The Gappy method remains the most expensive of the ones we computed, this reiterates that it should only be used when the gappyness factor is high enough to justify its cost.

In practice, $\epsilon_S$ is not the only error engineers are interested in. Concentrating on less stringent criteria for the reduction, one could drastically reduce the number of modes required. For example, we observed that the overall shape of the free surface requires significantly less modes to reconstruct. Similarly, any integrated quantities on the boundaries could be obtained accurately with a much lower number of modes.
<figure><img src="image_22.png" /><figcaption>Fig. 16. Total error with respect to the snapshots of the pressure and velocity fields for the newtonian drop example.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 17. The dam break configuration and the free surface at different times.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 18. Total error with respect to the snapshots of the pressure and velocity fields for the dam break example.</figcaption></figure>

## 6. Conclusions

In this paper, we proposed a series of algorithms to reduce data produced by PFEM solvers efficiently. This means we want to store less data and still be able to reconstruct, with a given precision, the full PFEM simulation. We first reviewed the different numerical tools required to reduce data produced by typical finite element solvers. We then acknowledged that transport-dominated problems benefit immensely from mapping approaches to reduce the data efficiently and that, fittingly, PFEM produces such hard-to-reduce data by nature of the free surface fluid problems it is typically used to solve. Two methods are then proposed, both of which attempt to use the Lagrangian mapping computed with PFEM as a map to reduce the data. The first one constructs exact supermeshes and projects the data onto these supermeshes, committing no errors. The second one directly reduces the data along the gappy particle trajectories and preserves fully the original connectivity. Then, looking at the performance of both methods on four examples, we concluded the following. It is mainly the velocity fields that benefit from our geometrical registration using our Lagrangian map, the pressure field is not affected as much. We observed that both proposed methods showed very similar results regarding reductibility performance. The strength of method 1 is that everything is computed on a single supermesh, and therefore, we can imagine using this method for classical -single mesh- model order reduction in the future. However, the computation of the supermeshes remains fairly expensive and limited as it cannot be used for large deformations due to the bijectivity constraints. Method 2, conversely, can always be used but can be cumbersome as connectivities have to be conserved. A proposed gappy iterative POD algorithm consistently provides a reliable lower bound for reducibility. However, it remains the most computationally expensive approach. In cases with low gappyness factors, it shows comparable performance to the cheaper direct geometrical approach, which is preferable when remeshing is infrequent. Lastly, it is important to note that not every PFEM simulation is worth reducing with these methods. Complex free surface simulations, such as water dam breaks, remain too challenging to tackle and require too many modes to recreate accurately in their entirety. We showed that manifold reduction will rather shine with slower dynamics problems, high viscosity flows, or non-newtonian fluids.

The following steps will be to use this as a post-processing tool and integrate this reduction approach directly into our PFEM architecture. Indeed, we are often faced with huge outputs limiting our simulation capabilities, which this has the potential to solve. We also intend to use this work as a first step towards achieving model order reduction in PFEM for free surface fluid flows. Finally, method 2 could be adapted to other particle-based solvers, such as SPH and MPM, as it does not require any mesh connectivity.

## CRediT authorship contribution statement

**Max Beckermann:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization; **Massimiliano Cremonesi:** Writing – review & editing, Supervision, Methodology, Investigation, Funding acquisition, Conceptualization; **Andrea Barbarulo:** Writing – review & editing, Supervision, Methodology, Investigation, Funding acquisition, Conceptualization.

## Data availability

Data will be made available on request.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Reduction on multiple supermeshes

If, after a certain time, the bijectivity Eq. (22) is no longer verified in method 1 (see Section 4.2.1), the deformations are too large. To bypass this problem, we can build a piecewise linear approximation and build multiple manifolds and supermeshes to reconstruct our complete solution. We will divide our time interval into $N_k$ sections, let $f : [1, N_k + 1] \to [1, n]$ be a monotonic increasing function with $f(1) = 1$, $f(N_k + 1) = n + 1$ such that we can define $\forall k \in [1, N_k]$, $\forall i \in [f(k), f(k + 1) - 1]$

$$ \Xi_{t_i}^k := \Omega_{t_i} \to \Omega \\ x \to \sum_{p=1}^{N_k} \mathcal{N}^h(x, t_i)|_p^e X_{Lag}^e(p, f(k)), \qquad (A.1) $$

assuring this time that the $\Xi_{t_i}^k$ and therefore $X_{t_i}^k$ are bijective. We then end up with $N_k$ manifolds

$$ \begin{aligned} \mathcal{K}_k &= \{\mathbf{u}^{h_i} \circ X_{t_i}^k \mid i \in [f(k), f(k+1) - 1]\} \\ &= \{\mathcal{N}^h(\Xi_{t_i}^{k-1}(\xi), t_i)^T \mathbf{u}^{h_i} \mid i \in [f(k), f(k+1) - 1]\} \subset L^2(\Omega), \end{aligned} \qquad (A.2) $$

and therefore $N_k$ supermeshes and $N_k$ reduced spaces with $n_k \in \mathbb{N}^*$ eigenvalues each $\{\lambda_j^k\}_{j=1}^{n_k}$. A fully reconstructed reduced solution $\tilde{u}(t_i)$ with $r$ modes will depend on how many modes we take from each manifold. We will take $m_k^r$ modes from manifold $k$ with a $r$ mode reconstruction s.t. $\sum_k^{N_k} m_k^r = r$. A heuristic approach we developed takes the one contributing with the highest eigenvalue each time a new mode is required. This can be justified by looking at the total $L^2$ reduction error:

$$ \begin{aligned} \epsilon_S^2(r) &\approx \frac{\sum_{i=1}^n \|\mathbf{u}^{h_i} - \tilde{u}(t_i)\|_{M_{h_i}}^2}{\sum_{i=1}^n \|\mathbf{u}^{h_i}\|_{M_{h_i}}^2} = \frac{\sum_{k=1}^{N_k} \sum_{i=f(k)}^{f(k+1)-1} \|\mathbf{u}^{h_i} - \tilde{u}(t_i)\|_{M_{h_i}}^2}{\sum_{i=1}^n \|\mathbf{u}^{h_i}\|_{M_{h_i}}^2} \\ &\approx \frac{\sum_{k=1}^{N_k} \sum_{i=f(k)}^{f(k+1)-1} \|\mathbf{u}^{h_i}\|_{M_{h_i}}^2 - \sum_{j=1}^{m_k^r} \lambda_j^k}{\sum_{i=1}^n \|\mathbf{u}^{h_i}\|_{M_{h_i}}^2} = 1 - \frac{\sum_{k=1}^{N_k} \sum_{j=1}^{m_k^r} \lambda_j^k}{\sum_{i=1}^n \|\mathbf{u}^{h_i}\|_{M_{h_i}}^2}. \end{aligned} \qquad (A.3) $$

Note that this is only correct under the strong assumption that the deformation is small in each manifold and that each $\|\cdot\|_{M_{h_i}}^2$ can be approximated by its equivalent mapped version in the associated manifold.

The full summary of Method 1 with multiple supermeshes is given in Algorithm 4.
**Algorithm 4 Reducing on multpliple supermeshes.**

**Input** - {$\{u^{h_i}, x^{h_i}\} \in \mathbb{R}^{N_i}, \mathcal{N}^{h_i} \in V^{N_i}, \forall i \in [1, n]\} , N \in \mathbb{N}^*$ the required number of modes and ReduceWithMassMatrix a boolean indicating if the reduction should include the supermesh.

With $u^{h_i}$ the snapshots of the field of interest to reduce and $x^{h_i}$ the Lagrangian displacement, both computed with PFEM.

**Initialization** - Identification of $X_{Lag}$ (19), $\forall p \in P_a \{A^p, \{\alpha_j^p\}_{j \in A^p}\}$

This initialization can be skipped if this information has been saved during the PFEM solve.

Sort $P_a = \{p \in [1, N_e]\}|t_{a_p} > 0\}$ by increasing $t_{a_p}$

Compute $X_{Lag}^e(p, i) = \sum_{j \in A^p} \alpha_j^e X_{Lag}(j, i)$ for $p \in P$ and $0 < t < t_{a_p}$ (20)

Let $k = 1, f(1) = 1, i = 1, M_k^+ : M_1, \mathcal{K}_k = \{u^{h_1}(x)\}$

**for** $i \in [2, n]$ **do**

Calculate $\Xi_{t_i}^k$ (A.1)

if $\Xi_{t_i}^k$ is bijective (22) then

Calculate mapped mesh $M$ with $\Xi_{t_i}^k$

Update the supermesh $M_k^+$, combining both $M$ and $M_k^+$, see Section 3.2

Add $u^{h_i} \circ \Xi_{t_i}^{k-1}$ to $\mathcal{K}_k$

else

Let $k = k + 1, f(k) = i, M_k^+ : M_i, \mathcal{K}_k = \{u^{h_i}(x)\}$

end if

**end for**

Let $N_k = k, f(k+1) = n+1$

**for** $k \in [1, N_k]$ **do**

Let $n_k = f(k+1) - f(k) - 1$ be the number of snapshots and $l_k \in \mathbb{N}^*$ be the number of nodes of $M_k^+$

Project without committing additional errors all $f \in \mathcal{K}_k$ onto $M_k^+$ and add them to $Y \in \mathbb{R}^{l_k \times n_k}$

if ReduceWithMassMatrix then

Compute $M^+ \in \mathbb{R}^{l_k \times l_k}$ the symmetric, positive definite mass matrix from $M_k^+$

Compute the eigenvalue decomposition $\{(\lambda_j^k, \phi_j^k)\}_{j=1}^{n_k} = \text{Eig}(Y^T M^+ Y)$

else

Compute the eigenvalue decomposition $\{(\lambda_j^k, \phi_j^k)\}_{j=1}^{n_k} = \text{Eig}(Y^T Y)$

end if

Compute the space modes with $\psi_j^k = \frac{1}{\sqrt{\lambda_j^k}} Y \phi_j^k$ for $1 \le j \le n_k$

**end for**

Let $d = \{0 : \{0\}_{k=1}^{N_k}\}$ be the dictionary of the chosen number of modes for each manifold

**for** $r \in [1, R]$ **do**

$\{m_k^r\}_{k=1}^{N_k} = d[r-1]$

$i = \arg \max_{k \in [1, N_k]} \left( \frac{\lambda_k^r}{m_k^r + 1} \right)$

$m_i^r = m_i^r + 1$ and $d[r] = \{m_k^r\}_{k=1}^{N_k}$

**end for**

Reconstruct partial solutions following the order given in $d$, express them in every configuration with $X_{t_i}^k$ and compute $\epsilon_N^2$.

**References**

[1] M. Cremonesi, A. Franci, S. Idelsohn, E. Oñate, A state of the art review of the particle finite element method (PFEM), Arch. Comput. Methods Eng. 27 (5) (2020) 1709–1735. https://doi.org/10.1007/s11831-020-09468-4

[2] E. Oñate, S.R. Idelsohn, F. Del Pin, R. Aubry, et al., The particle finite element method—An overview, Int. J. Comput. Methods 01 (02) (2004) 267–307. https://doi.org/10.1142/S0219876204000204

[3] S.r. Idelsohn, E. Oñate, F. Del Pin, et al., The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves, Int. J. Numer. Methods Eng. 61 (7) (2004) 964–989. https://doi.org/10.1002/nme.1096

[4] M. Cremonesi, U. Perego, A lagrangian finite element approach for the analysis of fluid-structure interaction problems, Int. J. Numer. Methods Eng. 84 (2010) 610–630. https://doi.org/10.1002/nme.2911

[5] M.L. Ceruglia, D. Thomas, R. Boman, V. Terrapon, J.P. Ponthot, et al., A fully partitioned Lagrangian framework for FSI problems characterized by free surfaces, large solid deformations and displacements, and strong added-mass effects, Comput. Methods Appl. Mech. Eng. 348 (2019) 409–442. https://doi.org/10.1016/j.cma.2019.01.021

[6] E. Oñate, A. Cornejo, F. Zárate, K. Kashiyama, A. Franci, et al., Combination of the finite element method and particle-based methods for predicting the failure of reinforced concrete structures under extreme water forces, Eng. Struct. 251 (2022) 113510. https://doi.org/10.1016/j.engstruct.2021.113510

[7] A. Franci, M. Cremonesi, U. Perego, G. Crosta, E. Oñate, 3D simulation of vajont disaster. Part 1: numerical formulation and validation. Eng. Geol. 279 (2020) 105854.

[8] X. Zhang, K. Krabbenhoft, D. Sheng, W. Li, et al., Numerical simulation of a flow-like landslide using the particle finite element method, Comput. Mech. 55 (1) (2015) 167–177. https://doi.org/10.1007/s00466-014-1088-z
[9] S. Idelsohn, M. Mier-Torrecilla, E. Oñate, et al., Multi-fluid flows with the particle finite element method, Comput. Methods Appl. Mech. Eng. 198 (33) (2009) 2750-2767. https://doi.org/10.1016/j.cma.2009.04.002

[10] J.M. Rodriguez, J.M. Carbonell, J.C. Cante, J. Oliver, et al., The particle finite element method (PFEM) in thermo-mechanical problems, Int. J. Numer. Methods Eng. 107 (9) (2016) 733-785. https://doi.org/10.1002/nme.5186

[11] P.B. Ryzhakov, J. García, E. Oñate, et al., Lagrangian finite element model for the 3D simulation of glass forming processes, Comput. Struct. 177 (2016) 126-140. https://doi.org/10.1016/j.compstruc.2016.09.007

[12] G. Rizzieri, L. Ferrara, M. Cremonesi, et al., Numerical simulation of the extrusion and layer deposition processes in 3D concrete printing with the particle finite element method, Comput. Mech. 73 (2) (2024) 277-295. https://doi.org/10.1007/s00466-023-02367-y

[13] G. Rizzieri, L. Ferrara, M. Cremonesi, et al., Simulation of viscoelastic free-surface flows with the particle finite element method, Comput. Part. Mech. (2024). https://doi.org/10.1007/s40571-024-00730-1

[14] T. Ye, D. Pan, C. Huang, M. Liu, Smoothed particle hydrodynamics (SPH) for complex fluid flows: recent developments in methodology and applications, Phys. Fluids 31 (1) (2019) 011301. https://doi.org/10.1063/1.5068697

[15] W. Solowski, B. Martin, W. Coombs, J. Guilkey, M. Matthias, Q. Tran, A. Tito, S. Seyedmohammadjavad, T. Roel, S. Kenichi, Material point method: overview and challenges ahead, Adv. Appl. Mech. 54 (2021) 113-204. https://doi.org/10.1016/bs.aams.2020.12.002

[16] C. Fu, M. Cremonesi, U. Perego, B. Hudobivnik, P. Wriggers, et al., Particle virtual element method (PVEM): an agglomeration technique for mesh optimization in explicit Lagrangian free-surface fluid modelling, Comput. Methods Appl. Mech. Eng. 433 (2025) 117461. https://doi.org/10.1016/j.cma.2024.117461

[17] T. Leyssens, M. Henry, J. Lambrechts, J.F. Remacle, et al., A Delaunay Refinement Algorithm for the Particle Finite Element Method Applied to Free Surface Flows, 2024.

[18] R. Zorrilla, A. Franci, Cut-PFEM: a particle finite element method using unfitted boundary meshes, Eng. Comput. 40 (5) (2024) 2739-2760. https://doi.org/10.1007/s00366-024-01956-6

[19] P. Benner, M. Ohlberger, A. Cohen, K. Willcox, Model Reduction and Approximation, Society for Industrial and Applied Mathematics, 2017. https://doi.org/10.1137/1.9781611974829

[20] Volume 2: Snapshot-Based Methods and Algorithms, De Gruyter, 2021. https://doi.org/10.1515/9783110671490. https://doi.org/10.1515/9783110671490

[21] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, P.W. Battaglia, et al., Learning to simulate complex physics with graph networks. International Conference on Machine Learning (2020) 8459-8468.

[22] L. Cordier, M. Bergmann, Réduction de dynamique par Décomposition Orthogonale aux Valeurs Propres POD, 2006. https://apisemanticscholar.org/CorpusID:53585466.

[23] C. Farhat, Lecture notes in Model reduction - AA216/CME345, 2017.

[24] S. Volkwein, Lecture Notes in "Model Reduction Using Proper Orthogonal Decomposition" (2011).

[25] S. Volkwein, Lecture Notes in "Proper Orthogonal Decomposition: Theory and Reduced-Order Modelling" (2013).

[26] P.E. Farrell, J.R. Maddison, Conservative interpolation between volume meshes by local Galerkin projection, Comput. Methods Appl. Mech. Eng. 200 (1) (2011) 89-100. https://doi.org/10.1016/j.cma.2010.07.015

[27] C. Gräßle, M. Hinze, POD reduced-order modeling for evolution equations utilizing arbitrary finite element discretizations, Adv. Comput. Math. 44 (6) (2018) 1941-1978. https://doi.org/10.1007/s10444-018-9620-x

[28] C. Little, C. Farhat, Projection-based dimensional reduction of adaptively refined nonlinear models, Commun. Appl. Math. Comput. (2023). https://doi.org/10.1007/s42967-023-00308-3

[29] S. Ullmann, M. Rotkvic, J. Lang, et al., POD-Galerkin reduced-order modeling with adaptive finite element snapshots, J. Comput. Phys. 325 (2016) 244-258. 1603.00668https://doi.org/10.1016/j.jcp.2016.08.018

[30] M. Yano, A minimum-residual mixed reduced basis method: exact residual certification and simultaneous finite-element reduced-basis refinement, ESAIM: Math. Model. Numer. Anal. 50 (1) (2016) 163-185. https://doi.org/10.1051/m2an/2015039

[31] A. Kolmogoroff, Uber die beste annaherung von funktionen einer gegebenen funktionenklasse, Ann. Math. 37 (1) (1936) 107. 1968691https://doi.org/10.2307/1968691

[32] M. Ohlberger, S. Rave, Reduced Basis Methods: Success, Limitations and Future Challenges, 2016.

[33] V. Ehrlacher, A. Iollo, D. Lombardi, T. Taddei, et al., Non-Linear Interpolations - Summer School Inria 2024 Lecture Notes, 2024.

[34] A. Iollo, D. Lombardi, Advection modes by optimal mass transfer, Phys. Rev. E 89 (2) (2014) 022923. https://doi.org/10.1103/PhysRevE.89.022923

[35] A. Iollo, T. Taddei, Mapping of coherent structures in parameterized flows by learning optimal transportation with Gaussian models, J. Comput. Phys. 471 (2022) 111671. https://doi.org/10.1016/j.jcp.2022.111671

[36] T. Taddei, A registration method for model order reduction: data compression and geometry reduction, SIAM J. Sci. Comput. 42 (2) (2020) A997-A1027. https://doi.org/10.1137/19M1271270

[37] T. Taddei, L. Zhang, Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs, ESAIM: Math. Model. Numer. Anal. 55 (1) (2021) 99-130. https://doi.org/10.1051/m2an/2020073

[38] M.F. Beg, M.I. Miller, A. Trouvé, L. Younes, et al., Computing large deformation metric mappings via geodesic flows of diffeomorphisms, Int. J. Comput. Vis. 61 (2) (2005) 139-157. https://doi.org/10.1023/B:VISI.000043755.93987 aa

[39] M. de Buhan, C. Dapogny, P. Frey, C. Nardoni, et al., An optimization method for elastic shape matching, C. R. Math. 354 (8) (2016) 783-787. https://doi.org/10.1016/j.cma.2016.05.007

[40] A. Kabalan, F. Casenave, F. Bordeu, V. Ehrlacher, A. Ern, et al., Elasticity-Based Mesh Morphing Technique with Application to Reduced-Order Modeling, 2024, https://doi.org/10.48550/arXiv.2407.02433

[41] R. Everson, L. Sirovich, Karhunen-Loève procedure for gappy data, JOSA A 12 (8) (1995) 1657-1664. https://doi.org/10.1364/JOSAA.12.001657

[42] D. Venturi, G.E. Karniadakis, Gappy data and reconstruction procedures for flow past a cylinder, J. Fluid Mech. 519 (2004) 315-336. https://doi.org/10.1017/S002211004001338

[43] H. Gunes, S. Sirisup, G.E. Karniadakis, et al., Gappy data: to krig or not to krig?, J. Comput. Phys. 212 (1) (2006) 358-382. https://doi.org/10.1016/j.jcp.2005.06.023

[44] L. Sirovich, Turbulence and the dynamics of coherent structures. I - Coherent structures. II - Symmetries and transformations. III - dynamics and scaling, Q. Appl. Math. 45 (1987). https://doi.org/10.1090/qam/910463

[45] E.L. Wachspress, A Rational Finite Element Basis, number 114 in Mathematics in Science and Engineering, Academic Press, New York, 1975.

[46] M.S. Floater, Mean value coordinates, Comput. Aided Geom. Des. 20 (1) (2003) 19-27. https://doi.org/10.1016/S0167-8396(03)00002-5

[47] A. Babu, B. Han, N. Asharaf, et al., Numerical solution of the viscous Burgers' equation using localized differential quadrature method, Part. Differ. Equ. Appl. Math. 4 (2021) 100044. https://doi.org/10.1016/j.padiff.2021.100044

[48] M.F. Tomé, R. Merejolli, G.S. Paulo, S. McKee, et al., An assessment of the PTT model on the impacting drop problem, J. Braz. Soc. Mech. Sci. Eng. 40 (9) (2018) 468. https://doi.org/10.1007/s40430-018-1388-x