# Structure-preserving reduced order model for parametric cross-diffusion systems

Jad Dabaghi, Virginie Ehrlacher

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

## To cite this version:

Jad Dabaghi, Virginie Ehrlacher. Structure-preserving reduced order model for parametric cross-diffusion systems. 2024. hal-03696025v2

HAL Id: hal-03696025

https://inria.hal.science/hal-03696025v2

Preprint submitted on 4 Mar 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# STRUCTURE-PRESERVING REDUCED ORDER MODEL FOR PARAMETRIC CROSS-DIFFUSION SYSTEMS

JAD DABAGHI<fnref n="1" />,* AND VIRGINIE EHRLACHER<fnref n="2" />

**Abstract.** In this work, we construct a structure-preserving Galerkin reduced-order model for the resolution of parametric cross-diffusion systems. Cross-diffusion systems are often used to model the evolution of the concentrations or volumic fractions of mixtures composed of different species, and can also be used in population dynamics (as for instance in the SKT system). These systems often read as nonlinear degenerated parabolic partial differential equations, the numerical resolutions of which are highly expensive from a computational point of view. We are interested here in cross-diffusion systems which exhibit a so-called entropic structure, in the sense that they can be formally written as gradient flows of a certain entropy functional which is actually a Lyapunov functional of the system. In this work, we propose a new reduced-order modelling method, based on a reduced basis paradigm, for the resolution of parameter-dependent cross-diffusion systems. Our method preserves, at the level of the reduced-order model, the main mathematical properties of the continuous solution, namely mass conservation, non-negativeness, preservation of the volume-filling property and entropy-entropy dissipation relationship. The theoretical advantages of our approach are illustrated by several numerical experiments.

**Mathematics Subject Classification.** 65M08, 65M60.

Received June 15, 2022. Accepted April 8, 2024.

## 1. INTRODUCTION

The study of cross-diffusion systems for the modeling of diffusive phenomena with multi-component mixtures has received a steadily growing amount of attention from mathematicians in the last decade. These systems arise in various application fields such as population dynamics [40] or growth of vascular tumors in medical biology [24]. In chemistry, such systems are used in order to model the evolution of densities, concentrations or volumic fractions of the various chemical species composing the mixture of interest [1].

### 1.1. Mathematical formulation of parametric cross-diffusion systems

Despite their practical relevance in a wide range of application fields, significant steps in their mathematical understanding have only been achieved recently. This is due to the fact that these systems usually read as coupled degenerate nonlinear systems, for which traditional analysis tools do not apply.

*Keywords and phrases.* Cross-diffusion systems, finite volumes, proper orthogonal decomposition.

<footnote n="1">¹ CERMICS, Ecole des Ponts ParisTech & INRIA, MATHERIALS project-team, 77455 Marne-la-Vallée Cedex 2, France.</footnote>

<footnote n="2">² Léonard de Vinci Pôle Universitaire, Research Center, 92 916 Paris La Défense, France.</footnote>

<footnote>*Corresponding author: jad.dabaghi@devinci.fr</footnote>
Structure-preserving reduced order model for parametric
cross-diffusion systems *

Jad Dabaghi †

Virginie Ehrlacher ‡

March 4, 2024

Abstract

In this work, we construct a structure-preserving Galerkin reduced-order model for the resolution of parametric cross-diffusion systems. Cross-diffusion systems are often used to model the evolution of the concentrations or volumic fractions of mixtures composed of different species, and can also be used in population dynamics (as for instance in the SKT system). These systems often read as nonlinear degenerated parabolic partial differential equations, the numerical resolutions of which are highly expensive from a computational point of view. We are interested here in cross-diffusion systems which exhibit a so-called entropic structure, in the sense that they can be formally written as gradient flows of a certain entropy functional which is actually a Lyapunov functional of the system. In this work, we propose a new reduced-order modelling method, based on a reduced basis paradigm, for the resolution of parameter-dependent cross-diffusion systems. Our method preserves, at the level of the reduced-order model, the main mathematical properties of the continuous solution, namely mass conservation, non-negativeness, preservation of the volume-filling property and entropy-entropy dissipation relationship. The theoretical advantages of our approach are illustrated by several numerical experiments.

**keywords:** cross-diffusion systems, finite volumes, Proper orthogonal decomposition.

# 1 Introduction

The study of cross-diffusion systems for the modeling of diffusive phenomena with multi-component mixtures has received a steadily growing amount of attention from mathematicians in the last decade. These systems arise in various application fields such as population dynamics [40] or growth of vascular tumors in medical biology [24]. In chemistry, such systems are used in order to model the evolution of densities, concentrations or volumic fractions of the various chemical species composing the mixture of interest [1].

## 1.1 Mathematical formulation of parametric cross-diffusion systems

Despite their practical relevance in a wide range of application fields, significant steps in their mathematical understanding have only been achieved recently. This is due to the fact that these systems usually read as coupled degenerate nonlinear systems, for which traditional analysis tools do not apply.

Let $\Omega \subset \mathbb{R}^d$, $d \ge 1$, be a polygonal domain, let $n$ denote its outward unit normal vector, and let $T > 0$ be the final simulation time. We consider a diffusive mixture of $N_s \ge 2$, chemical species, and assume that the diffusion laws of the different species within the mixture are parametrized by a vector of parameters $\mu \in \mathbb{R}^p$ for some $p \in \mathbb{N}^*$ and belonging to a set of parameter values $\mathcal{P} \subset \mathbb{R}^p$. For all $t \in [0, T]$, $x \in \Omega$ and $1 \le i \le N_s$, we denote by $u_{\mu,i}(t, x)$ the value of the local volumic fraction of species $1 \le i \le N_s$ at time $t \in [0, T]$, point $x \in \Omega$ and parameter value $\mu \in \mathcal{P}$. We also denote by $u_\mu := (u_{\mu,1}, \dots, u_{\mu,N_s})$ the set of all

<footnote>*This project has received funding from the ANR project "COMODO" (ANR-19-CE46-0002).</footnote>

<footnote>‡CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée cedex 2</footnote>

<footnote>†Léonard de Vinci Pôle Universitaire, Research Center, 92 916 Paris La Défense, France</footnote>
Let $\Omega \subset \mathbb{R}^d$, $d \ge 1$, be a polygonal domain, let $\mathbf{n}$ denote its outward unit normal vector, and let $T > 0$ be the final simulation time. We consider a diffusive mixture of $N_s \ge 2$, chemical species, and assume that the diffusion laws of the different species within the mixture are parametrized by a vector of parameters $\mu \in \mathbb{R}^p$ for some $p \in \mathbb{N}^*$ and belonging to a set of parameter values $\mathcal{P} \subset \mathbb{R}^p$. For all $t \in [0, T]$, $x \in \Omega$ and $1 \le i \le N_s$, we denote by $u_{\mu,i}(t,x)$ the value of the local volumic fraction of species $1 \le i \le N_s$ at time $t \in [0, T]$, point $x \in \Omega$ and parameter value $\mu \in \mathcal{P}$. We also denote by $\mathbf{u}_\mu := (u_{\mu,1}, \dots, u_{\mu,N_s})$ the set of all $N_s$ volumic fractions. We assume that the evolution of $\mathbf{u}_\mu$ is ruled by a parametrized cross-diffusion model with no-flux boundary conditions of the following form for all $\mu \in \mathcal{P}$:

$$ 
\begin{align} 
\partial_t u_{\mu,i} - \nabla \cdot \left( \sum_{j=1}^{N_s} A_{\mu,ij}(\mathbf{u}_\mu) \nabla u_{\mu,j} \right) &= 0, \quad \text{in } \Omega \times [0,T], \quad \forall i \in [1, N_s], \nonumber \\
\left( \sum_{j=1}^{N_s} A_{\mu,ij}(\mathbf{u}_\mu) \nabla u_{\mu,j} \right) \cdot \mathbf{n} &= 0, \quad \text{on } \partial\Omega \times [0,T], \quad \forall i \in [1, N_s], \nonumber \\
u_{\mu,i}(0, \mathbf{x}) &= u_i^0(\mathbf{x}), \quad \text{in } \Omega, \quad \forall i \in [1, N_s], \tag{1.1}
\end{align} 
$$

for some initial conditions $(u_1^0, \dots, u_{N_s}^0) \in [L^\infty(\Omega)]^{N_s}$, and where for all $1 \le i,j \le N_s$, $A_{\mu,ij}: \mathbb{R}^{N_s} \to \mathbb{R}$ is a smooth cross-diffusion coefficient.

It has been recently understood [8, 25] that various cross-diffusion systems can be formally seen as gradient flows of a certain entropy functional with respect to a metric which can be seen as an extension of the so-called Wasserstein metric. As a consequence, the corresponding entropy functional can be seen as a Lyapunov function for the cross-diffusion system. The rate of decay of the entropy functional can then give precious insight on the long-time behaviour of the solutions of the cross-diffusion system.

The development of numerical methods for the resolution of such systems has been a very active field of research in the recent years. An important challenge is to design numerical schemes which preserve the main mathematical properties of the continuous model at the discrete level such as non-negativity of the solutions and entropy-entropy dissipation relationship. Finite element schemes have been proposed for instance in [3], while structure-preserving finite-volume schemes have been studied in [10, 11, 27, 28]. The resolution of such schemes, which is usually based on implicit numerical schemes for nonlinear models, remains quite costly from a computational point of view especially when the number of species in the system or the number of spatial degrees of freedom is large.

## 1.2. Reduced-order model

The complexity of these computations becomes prohibitive when the cross-diffusion model depends on some parameters since it is necessary to compute the solution of these systems for many values of these parameters. The resolution of inverse problems so as to identify optimal values of these parameters which enable to reproduce experimental results is an example of such a context, where a significant parametric study has to be performed. As an alternative, model-order reduction techniques have been developed for many type of problems so as to alleviate the computational burden of parametric studies. Reduced order model techniques have been proved to yield computationally effective approaches to approximate the solution of parametrized partial differential equations encountered in many problems in science and engineering [22,38]. Several methods have been developed in the literature. Among them we mention the Proper Orthogonal Decomposition (POD) method [6, 16, 17, 19, 31,32], the Proper Generalized Decomposition (PGD) method [30,33,36] or reduced basis methods [22,34,35,38]. A vivid field of research consists in designing reduced-order models which preserve the mathematical properties of the high-fidelity model. We refer the reader to the nice overview [21] on reduced-order models for time-dependent problems, where these issues are tackled. Let us also mention for instance Hamiltonian systems [18,20], Hamiltonian port-systems [12,39], the list given here being far from exhaustive.
$N_s$ volumic fractions. We assume that the evolution of $u_\mu$ is ruled by a parametrized cross-diffusion model with no-flux boundary conditions of the following form for all $\mu \in \mathcal{P}$:

$$ \begin{gathered} \partial_t u_{\mu,i} - \nabla \cdot \left( \sum_{j=1}^{N_s} A_{\mu,ij}(u_\mu) \nabla u_{\mu,j} \right) = 0, \quad \text{in } \Omega \times [0,T], \quad \forall i \in [1, N_s], \\ \left( \sum_{j=1}^{N_s} A_{\mu,ij}(u_\mu) \nabla u_{\mu,j} \right) \cdot n = 0, \quad \text{on } \partial\Omega \times [0,T], \quad \forall i \in [1, N_s], \\ u_{\mu,i}(0, \boldsymbol{x}) = u_0^0(\boldsymbol{x}), \quad \text{in } \Omega, \quad \forall i \in [1, N_s], \end{gathered} \tag{1.1} $$

for some initial conditions $(u_1^0, \dots, u_{N_s}^0) \in [L^\infty(\Omega)]^{N_s}$, and where for all $1 \le i, j \le N_s$, $A_{\mu,ij}: \mathbb{R}^{N_s} \to \mathbb{R}$ is a smooth cross-diffusion coefficient.

It has been recently understood [8, 25] that various cross-diffusion systems can be formally seen as gradient flows of a certain entropy functional with respect to a metric which can be seen as an extension of the so-called Wasserstein metric. As a consequence, the corresponding entropy functional can be seen as a Lyapunov function for the cross-diffusion system. The rate of decay of the entropy functional can then give precious insight on the long-time behaviour of the solutions of the cross-diffusion system.

The development of numerical methods for the resolution of such systems has been a very active field of research in the recent years. An important challenge is to design numerical schemes which preserve the main mathematical properties of the continuous model at the discrete level such as non-negativity of the solutions and entropy-entropy dissipation relationship. Finite element schemes have been proposed for instance in [3], while structure-preserving finite-volume schemes have been studied in [10, 27, 11, 28]. The resolution of such schemes, which is usually based on implicit numerical schemes for nonlinear models, remains quite costly from a computational point of view especially when the number of species in the system or the number of spatial degrees of freedom is large.

## 1.2 Reduced-order model

The complexity of these computations becomes prohibitive when the cross-diffusion model depends on some parameters since it is necessary to compute the solution of these systems for many values of these parameters. The resolution of inverse problems so as to identify optimal values of these parameters which enable to reproduce experimental results is an example of such a context, where a significant parametric study has to be performed. As an alternative, model-order reduction techniques have been developed for many type of problems so as to alleviate the computational burden of parametric studies. Reduced order model techniques have been proved to yield computationally effective approaches to approximate the solution of parametrized partial differential equations encountered in many problems in science and engineering [22, 38]. Several methods have been developed in the literature. Among them we mention the Proper Orthogonal Decomposition (POD) method [17, 6, 31, 32, 19, 16], the Proper Generalized Decomposition (PGD) method [36, 33, 30] or reduced basis methods [35, 34, 22, 38]. A vivid field of research consists in designing reduced-order models which preserve the mathematical properties of the high-fidelity model. We refer the reader to the nice overview [21] on reduced-order models for time-dependent problems, where these issues are tackled. Let us also mention for instance Hamiltonian systems [18, 20], Hamiltonian port-systems [12, 39], the list given here being far from exhaustive.

The aim of the present work is to develop a new structure-preserving Galerkin reduced-order model, relying on a POD/Reduced Basis paradigm, so as to compute the solution of some parametrized cross-diffusion equations which exhibit a specific entropic structure as highlighted in [10, 11, 28]. This includes in particular the so-called Maxwell-Stefan model. Indeed, the reduced-order model we propose here preserves the main theoretical features of the continuous cross-diffusion model, namely

(P1) mass conservation,

(P2) non-negativeness of solutions,

(P3) volume-filling constraints,
The aim of the present work is to develop a new structure-preserving Galerkin reduced-order model, relying on a POD/Reduced Basis paradigm, so as to compute the solution of some parametrized cross-diffusion equations which exhibit a specific entropic structure as highlighted in [10, 11, 28]. This includes in particular the so-called Maxwell-Stefan model. Indeed, the reduced-order model we propose here preserves the main theoretical features of the continuous cross-diffusion model, namely

(P1) mass conservation,

(P2) non-negativeness of solutions,

(P3) volume-filling constraints,

(P4) entropic structure of the cross-diffusion system.

We emphasize on the fact that a standard POD reduced-order model does not satisfy the four properties listed above.

We explain in more details these four properties below. For the sake of simplicity, we will consider in the following one particular type of parametrized cross-diffusion system, namely the one highlighted in [10], which satisfies these properties, and explain the strategy we propose for building a structure preserving reduced order model on this particular system. However, we would like to emphasize the fact that the approach proposed here can be easily adapted to other types of cross-diffusion systems enjoying similar mathematical properties, such as the ones studied in [11, 28], including the so-called Maxwell-Stefan system.

Let us emphasize here that we focus in this work on the definition of a structure-preserving Galerkin reduced-order model for this type of equations. As these systems are nonlinear, a fully efficient reduced-order model would require the design of an appropriate hyper-reduction method, like Empirical Interpolation Method for instance [2], to speed up the computation of the nonlinear terms appearing in the equations. This is however out of the scope of the present article and we leave the question of designing an adequate structure-preserving hyperreduction method for future research. The present structure-preserving Galerkin ROM actually reads as a *transformed Galerkin ROM*, in the sense that it builds on the reduction of a transformed set of snapshots, using so-called *entropic variables* associated to the type of cross-diffusion systems we consider here. Let us mention here that there is a vast literature on transformed Galerkin ROMs for various types of problems [9, 13, 23, 37, 41, 42]. However, at least up to the best of our knowledge, this is the first work dealing with entropic cross-diffusion systems. Notice however that such entropic variables have already been considered for the reduction of Euler equations, where their use is quite standard [4, 7].

## 1.3. Outline of the article

Our paper is organized as follows. In Section 2.1, we present the cross-diffusion model studied in [10] and its mathematical properties as a particular example of cross-diffusion systems enjoying a favourable entropic structure for our analysis. We want to stress again on the fact that the method we propose in this work can be easily adapted to other cross-diffusion systems falling into the scope of the general setting studied in [28], including the Maxwell-Stefan model [11, 26]. However, we prefer to carry out the analysis here on this particular example for the sake of clarity. In particular, in Section 2.2 we present the two point flux approximation finite volume method preserving the structural properties of the solution proposed in [10], which is the high-fidelity solver we consider in this work. Section 3 is dedicated to the construction of the structure-preserving reduced model and its properties. Finally, in Section 4 we present numerical experiments showing the advantages of our approach.

## 2. A PARTICULAR CROSS-DIFFUSION SYSTEM

The aim of this section is to introduce the parametric cross-diffusion system considered in this work, together with the main mathematical properties of the solution which we wish to be preserved by any reduced-order model. We also describe here the discrete numerical scheme used in order to perform high-fidelity computations.
(P4) entropic structure of the cross-diffusion system.

We emphasize on the fact that a standard POD reduced-order model does not satisfy the four properties listed above.

We explain in more details these four properties below. For the sake of simplicity, we will consider in the following one particular type of parametrized cross-diffusion system, namely the one highlighted in [10], which satisfies these properties, and explain the strategy we propose for building a structure preserving reduced order model on this particular system. However, we would like to emphasize the fact that the approach proposed here can be easily adapted to other types of cross-diffusion systems enjoying similar mathematical properties, such as the ones studied in [11, 28], including the so-called Maxwell-Stefan system.

Let us emphasize here that we focus in this work on the definition of a structure-preserving Galerkin reduced-order model for this type of equations. As these systems are nonlinear, a fully efficient reduced-order model would require the design of an appropriate hyper-reduction method, like Empirical Interpolation Method for instance [2], to speed up the computation of the nonlinear terms appearing in the equations. This is however out of the scope of the present article and we leave the question of designing an adequate structure-preserving hyperreduction method for future research. The present structure-preserving Galerkin ROM actually reads as a transformed Galerkin ROM, in the sense that it builds on the reduction of a transformed set of snapshots, using so-called entropic variables associated to the type of cross-diffusion systems we consider here. Let us mention here that there is a vast literature on transformed Galerkin ROMs for various types of problems [41, 9, 42, 13, 23, 37]. However, at least up to the best of our knowledge, this is the first work dealing with entropic cross-diffusion systems. Notice however that such entropic variables have already been considered for the reduction of Euler equations, where their use is quite standard [7, 4].

## 1.3 Outline of the article

Our paper is organized as follows. In Section 2.1, we present the cross-diffusion model studied in [10] and its mathematical properties as a particular example of cross-diffusion systems enjoying a favourable entropic structure for our analysis. We want to stress again on the fact that the method we propose in this work can be easily adapted to other cross-diffusion systems falling into the scope of the general setting studied in [28], including the Maxwell-Stefan model [11, 26]. However, we prefer to carry out the analysis here on this particular example for the sake of clarity. In particular, in Section 2.2 we present the two point flux approximation finite volume method preserving the structural properties of the solution proposed in [10], which is the high-fidelity solver we consider in this work. Section 3 is dedicated to the construction of the structure-preserving reduced model and its properties. Finally, in Section 4 we present numerical experiments showing the advantages of our approach.

# 2 A particular cross-diffusion system

The aim of this section is to introduce the parametric cross-diffusion system considered in this work, together with the main mathematical properties of the solution which we wish to be preserved by any reduced-order model. We also describe here the discrete numerical scheme used in order to perform high-fidelity computations.

## 2.1 Continuous model and main mathematical properties

From now on, let us assume that

$$ \mathcal{P} := \{ \mu := (a_{ij})_{1 \le i \ne j \le N_s} \in (\mathbb{R}_+^*)^{N_s(N_s-1)}, a_{ij} = a_{ji} \quad \forall 1 \le i \ne j \le N_s \}, $$

and that for all $\mu := (a_{ij})_{1 \le i \ne j \le N_s} \in \mathcal{P}$, all $u := (u_1, \dots, u_{N_s}) \in \mathbb{R}^{N_s}$ and all $1 \le i \ne j \le N_s$,

$$ A_{\mu,ii}(\mathbf{u}) = \sum_{j=1, j \ne i}^{N_s} a_{ij} u_j, \quad \text{and} \quad A_{\mu,ij}(\mathbf{u}) = -a_{ij} u_i. $$
## 2.1. Continuous model and main mathematical properties

From now on, let us assume that

$$ \mathcal{P} := \{ \mu := (a_{ij})_{1 \le i \ne j \le N_s} \in (\mathbb{R}_+^*)^{N_s(N_s-1)}, a_{ij} = a_{ji} \quad \forall 1 \le i \ne j \le N_s \}, $$

and that for all $\mu := (a_{ij})_{1 \le i \ne j \le N_s} \in \mathcal{P}$, all $\mathbf{u} := (u_1, \dots, u_{N_s}) \in \mathbb{R}^{N_s}$ and all $1 \le i \ne j \le N_s$,

$$ A_{\mu,ii}(\mathbf{u}) = \sum_{j=1, j \neq i}^{N_s} a_{ij} u_j, \quad \text{and} \quad A_{\mu,ij}(\mathbf{u}) = -a_{ij} u_i. $$

Then, system (1.1) boils down to:

$$ \begin{gathered} \partial_t u_{\mu,i} - \nabla \cdot \left( \sum_{1 \le j \ne i \le N_s} a_{ij} (u_{\mu,j} \nabla u_{\mu,i} - u_{\mu,i} \nabla u_{\mu,j}) \right) = 0 \quad \text{in } \Omega \times [0, T], \quad \text{for } i \in [1, N_s], \\ \left( \sum_{1 \le j \ne i \le N_s} a_{ij} (u_{\mu,j} \nabla u_{\mu,i} - u_{\mu,i} \nabla u_{\mu,j}) \right) \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega \times [0, T], \quad \text{for } i \in [1, N_s], \\ u_{\mu,i}(\mathbf{x}, 0) = u_i^0(\mathbf{x}) \quad \text{in } \Omega, \quad \text{for } i \in [1, N_s]. \end{gathered} \quad (2.1) $$

This system has been originally introduced in [1, 43]. It is used in particular in [1] in order to model diffusion phenomena in the bulk of a thin-film layer solar cell during its production process by physical vapor deposition. Since for all $\mu \in \mathcal{P}$ and all $i \in [1, N_s]$, $u_{\mu,i}(t, x)$ represents the local volumic fraction of specie i at time t and point $x \in \Omega$, it is expected from a physical point of view that each function $u_{\mu,i}$ is non-negative and satisfy the so-called volumic constraint which reads as $\sum_{i=1}^{N_s} u_{\mu,i}(t, x) = 1$ for a.e. $(t, x) \in (0, T) \times \Omega$. This is the reason why we introduce the set $\mathcal{A} \subset \mathbb{R}^{N_s}$, which is defined by

$$ \mathcal{A} := \left\{ \mathbf{u} = (u_1, \dots, u_{N_s}) \in \mathbb{R}_+^{N_s}, \sum_{i=1}^{N_s} u_i = 1 \right\}. $$

From now on, let us assume that the initial condition $\mathbf{u}^0 := (u_1^0, \dots, u_{N_s}^0)$ satisfies $\mathbf{u}^0(x) \in \mathcal{A}$ for almost all $x \in \Omega$. Let us assume in addition that

$$ \int_{\Omega} u_i^0(x) \, dx > 0, \quad \forall i \in [1, N_s]. $$

Then, it is proved in [1] that there exists at least one weak solution $\mathbf{u}_{\mu} \in L^2((0, T), H^1(\Omega))^{N_s}$ to (2.1) such that $\partial_t \mathbf{u}_{\mu} \in L^2((0, T), H^1(\Omega))^{N_s}$. Moreover, such a weak solution can be proved to satisfy the following properties:

(P1) **mass conservation:** for all $t \in [0, T]$ and all $1 \le i \le N_s$, $\int_{\Omega} u_{\mu,i}(t, x) \, dx = \int_{\Omega} u_i^0(x) \, dx$.

(P2) **volume-filling constraint:** $\mathbf{u}_{\mu}(t, x) \in \mathcal{A}$ for almost all $(t, x) \in (0, T) \times \Omega$.

### 2.1.1. Gradient flow structure

System (2.1) has a formal gradient flow structure. Indeed, we introduce the entropy functional $E: L^{\infty}(\Omega, \mathcal{A}) \to \mathbb{R}$ defined by

$$ \forall \mathbf{u} := (u_i)_{i \in [1, N_s]} \in L^{\infty}(\Omega, \mathcal{A}), \quad E(\mathbf{u}) := \int_{\Omega} \left( \sum_{i=1}^{N_s} u_i(t, x) \ln(u_i(t, x)) \right) \, dx. \quad (2.2) $$
Then, system (1.1) boils down to:

$$
\begin{equation} \tag{2.1}
\begin{aligned}
\partial_t u_{\mu,i} - \nabla \cdot \left( \sum_{1 \le j \ne i \le N_s} a_{ij} (u_{\mu,j} \nabla u_{\mu,i} - u_{\mu,i} \nabla u_{\mu,j}) \right) &= 0 \quad \text{in } \Omega \times [0, T], \text{ for } i \in [1, N_s], \\
\left( \sum_{1 \le j \ne i \le N_s} a_{ij} (u_{\mu,j} \nabla u_{\mu,i} - u_{\mu,i} \nabla u_{\mu,j}) \right) \cdot \mathbf{n} &= 0 \quad \text{on } \partial\Omega \times [0, T], \text{ for } i \in [1, N_s], \\
u_{\mu,i}(\mathbf{x}, 0) &= u_i^0(\mathbf{x}) \quad \text{in } \Omega, \text{ for } i \in [1, N_s].
\end{aligned}
\end{equation}
$$

This system has been originally introduced in [43, 1]. It is used in particular in [1] in order to model
diffusion phenomena in the bulk of a thin-film layer solar cell during its production process by physical
vapor deposition. Since for all $\mu \in \mathcal{P}$ and all $i \in [1, N_s]$, $u_{\mu,i}(t,x)$ represents the local volumic fraction
of specie $i$ at time $t$ and point $x \in \Omega$, it is expected from a physical point of view that each function
$u_{\mu,i}$ is non-negative and satisfy the so-called volumic constraint which reads as $\sum_{i=1}^{N_s} u_{\mu,i}(t,x) = 1$ for a.e. $(t,x) \in (0,T) \times \Omega$. This is the reason why we introduce the set $\mathcal{A} \subset \mathbb{R}^{N_s}$, which is defined by

$$
\mathcal{A} := \left\{ \mathbf{u} = (u_1, \dots, u_{N_s}) \in \mathbb{R}_+^{N_s}, \sum_{i=1}^{N_s} u_i = 1 \right\}.
$$

From now on, let us assume that the initial condition **u**⁰ := (u₀, ..., uₙ) satisfies **u**⁰(x) ∈ 𝒜 for almost all x ∈ Ω. Let us assume in addition that

$$
\int_{\Omega} u_{i}^{0}(x) \, dx > 0, \quad \forall i \in [1, N_{s}].
$$

Then, it is proved in [1] that there exists at least one weak solution **u**<sub>μ</sub> ∈ L<sup>2</sup>((0, T), H<sup>1</sup>(Ω))<sup>N<sub>s</sub></sup> to (2.1) such that ∂<sub>t</sub>**u**<sub>μ</sub> ∈ L<sup>2</sup>((0, T), H<sup>1</sup>(Ω))<sup>N<sub>s</sub></sup>. Moreover, such a weak solution can be proved to satisfy the following properties:

(P1) **mass conservation:** for all $t \in [0, T]$ and all $1 \le i \le N_s$, $\int_{\Omega} u_{\mu,i}(t,x) \, dx = \int_{\Omega} u_i^0(x) \, dx$.

(P2) volume-filling constraint: $u_\mu(t, x) \in \mathcal{A}$ for almost all $(t, x) \in (0, T) \times \Omega$.

2.1.1 Gradient flow structure

System (2.1) has a formal gradient flow structure. Indeed, we introduce the entropy functional $E$ :
$L^\infty(\Omega, \mathcal{A}) \to \mathbb{R}$ defined by

$$
\forall \mathbf{u} := (u_i)_{i \in [1, N_s]} \in L^\infty(\Omega, \mathcal{A}), \quad E(\mathbf{u}) := \int_{\Omega} \left( \sum_{i=1}^{N_s} u_i(t, x) \ln(u_i(t, x)) \right) dx. \quad (2.2)
$$

Following [1, Section 1.1.2] and [10, Proposition 1.3], it holds that system (2.1) enjoys a formal gradient flow structure and that E is a Lyapunov functional. More precisely, system (2.1) can be rewritten as

$$
\begin{align*}
\partial_t \boldsymbol{u} - \nabla \cdot (\mathcal{C}(\boldsymbol{u}) \nabla D E(\boldsymbol{u})) &= 0 && \text{in } \Omega \times [0, T] \\
(\mathcal{C}(\boldsymbol{u}) \nabla D E(\boldsymbol{u})) \cdot \boldsymbol{n} &= 0 && \text{on } \partial\Omega \times [0, T] \\
\boldsymbol{u}(\boldsymbol{x}, 0) &= \boldsymbol{u}^0(\boldsymbol{x}) && \text{in } \Omega
\end{align*}
$$

where $\mathcal{C}(u)$ is the mobility matrix and $DE$ is the entropy differential defined by $(DE(u))_i := D_{u_i} E(u) = 1 + \ln(u_i)$. Furthermore, the following entropy-entropy dissipation inequality holds for all $\mu := (a_{ij})_{1 \le i \ne j \le N_s} \in$
Following Section 1.1.2 of [1] and Proposition 1.3 of [10], it holds that system (2.1) enjoys a formal gradient flow structure and that E is a Lyapunov functional. More precisely, system (2.1) can be rewritten as

$$
\begin{align*}
\partial_t \boldsymbol{u} - \nabla \cdot (\mathbb{C}(\boldsymbol{u}) \nabla \mathrm{DE}(\boldsymbol{u})) &= 0 && \text{in } \Omega \times [0, T] \\
(\mathbb{C}(\boldsymbol{u}) \nabla \mathrm{DE}(\boldsymbol{u})) \cdot \boldsymbol{n} &= 0 && \text{on } \partial\Omega \times [0, T] \\
\boldsymbol{u}(\boldsymbol{x}, 0) &= \boldsymbol{u}^0(\boldsymbol{x}) && \text{in } \Omega
\end{align*}
$$

where $\mathbb{C}(u)$ is the mobility matrix and $\mathrm{DE}$ is the entropy differential defined by $(\mathrm{DE}(\boldsymbol{u}))_i := D_{u_i} E(\boldsymbol{u}) = 1 + \ln(u_i)$. Furthermore, the following entropy-entropy dissipation inequality holds for all $\mu := (a_{ij})_{1\le i \ne j \le N_s} \in \mathcal{P}$, for a.e. $t \in (0, T)$,

$$
\frac{d}{dt} E(\mathbf{u}_{\mu}(t, \cdot)) + \int_{\Omega} \left( \sum_{1 \le i < j \le N_s} a_{ij} u_{\mu,i}(t, x) u_{\mu,j}(t, x) |\nabla \ln(u_{\mu,i}(t, x)) - \nabla \ln(u_{\mu,j}(t, x))|^2 \right) dx = 0.
$$

**Proposition 2.1.** *Denoting by*

$$
a_{\mu}^{*} := \min_{1 \le i \ne j \le N_s} a_{ij}, \qquad (2.3)
$$

the quantity $\Gamma$ defined by

$$
\Gamma := \int_{\Omega} \left( \sum_{1 \le i < j \le N_s} a_{ij} u_{\mu,i}(t,x) u_{\mu,j}(t,x) |\nabla \ln(u_{\mu,i}(t,x)) - \nabla \ln(u_{\mu,j}(t,x))|^2 \right) dx \quad (2.4)
$$

satisfies

$$
\Gamma \geq a_{\mu}^{*} \sum_{i=1}^{N_s} \int_{\Omega} |\nabla \sqrt{u_{\mu,i}}(t,x)|^2 dx = a_{\mu}^{*} \sum_{i=1}^{N_s} \int_{\Omega} u_{\mu,i}(t,x) |\nabla \ln(u_{\mu,i})(t,x)|^2 dx. \quad (2.5)
$$

Proof. Since $a_{ij} = a_{ji}$ for all $1 \le i \ne j \le N_s$, we have

$$
\Gamma = \frac{1}{2} \int_{\Omega} \left( \sum_{1 \le i \le j \le N_s} a_{ij} u_{\mu,i}(t,x) u_{\mu,j}(t,x) |\nabla \ln(u_{\mu,i}(t,x)) - \nabla \ln(u_{\mu,j}(t,x))|^2 \right) dx.
$$

Next, employing the property $c\nabla(\ln(c)) = \nabla c$ and the fact that $\sum_{i=1}^{N_s} u_i = 1$, it yields

$$
\Gamma = \frac{1}{2} \int_{\Omega} \left( \sum_{1 \le i \le j \le N_s} a_{ij} \left( \frac{u_{\mu,j}(t,x)}{u_{\mu,i}(t,x)} |\nabla u_{\mu,i}(t,x)|^2 + \frac{u_{\mu,i}(t,x)}{u_{\mu,j}(t,x)} |\nabla u_{\mu,j}(t,x)|^2 \right) \right) dx. \quad (2.6)
$$

Furthermore the properties $2\sqrt{c}\nabla(\sqrt{c}) = \nabla c$ and $\sum_{i=1}^{N_s} u_{\mu,i}(t,x) = \sum_{j=1}^{N_s} u_{\mu,j}(t,x) = 1$ give

$$
\Gamma = 4 \int_{\Omega} \left( \sum_{1 \le i \le j \le N_s} a_{ij} \left| \nabla \left( \sqrt{u_{\mu,i}(t,x)} \right) \right|^2 \right) dx \ge a_{\mu}^{*} \int_{\Omega} \left( \sum_{i=1}^{N_s} \left| \nabla \left( \sqrt{u_{\mu,i}(t,x)} \right) \right|^2 \right) dx. \quad (2.7)
$$

In particular (2.5) implies that system (2.1) satisfies the following entropy-entropy dissipation relation.
P: for a.e. $t \in (0, T)$,

$$ \frac{d}{dt} E(\mathbf{u}_{\mu}(t, \cdot)) + \int_{\Omega} \left( \sum_{1 \le i < j \le N_s} a_{ij} u_{\mu,i}(t, x) u_{\mu,j}(t, x) |\nabla \ln(u_{\mu,i}(t, x)) - \nabla \ln(u_{\mu,j}(t, x))|^2 \right) dx = 0. $$

**Proposition 2.1.** Denoting by

$$ a_{\mu}^{*} := \min_{1 \le i \ne j \le N_s} a_{ij}, \qquad (2.3) $$

the quantity $\Gamma$ defined by

$$ \Gamma := \int_{\Omega} \left( \sum_{1 \le i < j \le N_s} a_{ij} u_{\mu,i}(t, x) u_{\mu,j}(t, x) |\nabla \ln(u_{\mu,i}(t, x)) - \nabla \ln(u_{\mu,j}(t, x))|^2 \right) dx \quad (2.4) $$

satisfies

$$ \Gamma \ge a_{\mu}^{*} \sum_{i=1}^{N_s} \int_{\Omega} |\nabla \sqrt{u_{\mu,i}}(t, x)|^2 dx = a_{\mu}^{*} \sum_{i=1}^{N_s} \int_{\Omega} u_{\mu,i}(t, x) |\nabla \ln(u_{\mu,i})(t, x)|^2 dx. \quad (2.5) $$

**Proof.** Since $a_{ij} = a_{ji}$ for all $1 \le i \ne j \le N_s$, we have

$$ \Gamma = \frac{1}{2} \int_{\Omega} \left( \sum_{1 \le i \le j \le N_s} a_{ij} u_{\mu,i}(t, x) u_{\mu,j}(t, x) |\nabla \ln(u_{\mu,i}(t, x)) - \nabla \ln(u_{\mu,j}(t, x))|^2 \right) dx. $$

Next, employing the property $c\nabla(\ln(c)) = \nabla c$ and the fact that $\sum_{i=1}^{N_s} u_i = 1$, it yields

$$ \Gamma = \frac{1}{2} \int_{\Omega} \left( \sum_{1 \le i \le j \le N_s} a_{ij} \left( \frac{u_{\mu,j}(t, x)}{u_{\mu,i}(t, x)} |\nabla u_{\mu,i}(t, x)|^2 + \frac{u_{\mu,i}(t, x)}{u_{\mu,j}(t, x)} |\nabla u_{\mu,j}(t, x)|^2 \right) \right) dx. \quad (2.6) $$

Furthermore the properties $2\sqrt{c}\nabla(\sqrt{c}) = \nabla c$ and $\sum_{i=1}^{N_s} u_{\mu,i}(t, x) = \sum_{j=1}^{N_s} u_{\mu,j}(t, x) = 1$ give

$$ \Gamma = 4 \int_{\Omega} \left( \sum_{1 \le i \le j \le N_s} a_{ij} |\nabla (\sqrt{u_{\mu,i}(t, x)})|^2 \right) dx \ge a_{\mu}^{*} \int_{\Omega} \left( \sum_{i=1}^{N_s} |\nabla (\sqrt{u_{\mu,i}(t, x)})|^2 \right) dx. \quad (2.7) $$

□

In particular (2.5) implies that system (2.1) satisfies the following entropy-entropy dissipation relation.

(P3) Entropy-entropy dissipation relation: for almost all $t \in (0, T)$,

$$ \frac{d}{dt} E(\mathbf{u}_{\mu}(t, \cdot)) + a_{\mu}^{*} \sum_{i=1}^{N_s} \int_{\Omega} u_{\mu,i}(t, x) |\nabla \ln(u_{\mu,i}(t, x))|^2 dx \le 0. $$

**Remark 2.2.** Other types of cross-diffusion systems enjoy similar mathematical properties as the ones highlighted here for system (2.1). In the most general setting, these can be written under the form highlighted in [28]. This includes the so-called Maxwell-Stefan [26, 11] cross-diffusion model, used for modeling gaseous mixtures, which is another example of system for which the reduced order modeling strategy we propose here can be easily adapted.

In the following section, we describe the finite volume numerical scheme we use here in order to compute high-fidelity solutions of system (2.1) and summarize its main mathematical properties. This finite volume scheme, which was introduced in [10], actually preserves analogous versions of properties (P1)-(P2)-(P3) on the discrete level.
(P3) **Entropy-entropy dissipation relation:** for almost all $t \in (0, T)$,

$$ \frac{d}{dt} E(\mathbf{u}_{\mu}(t, \cdot)) + a_{\mu}^{*} \sum_{i=1}^{N_b} \int_{\Omega} u_{\mu,i}(t, x) |\nabla \ln(u_{\mu,i}(t, x))|^2 dx \le 0. $$

**Remark 2.2.** Other types of cross-diffusion systems enjoy similar mathematical properties as the ones highlighted here for system (2.1). In the most general setting, these can be written under the form highlighted in [28]. This includes the so-called Maxwell-Stefan [11, 26] cross-diffusion model, used for modeling gaseous mixtures, which is another example of system for which the reduced order modeling strategy we propose here can be easily adapted.

In the following section, we describe the finite volume numerical scheme we use here in order to compute high-fidelity solutions of system (2.1) and summarize its main mathematical properties. This finite volume scheme, which was introduced in [10], actually preserves analogous versions of properties (P1)-(P2)-(P3) on the discrete level.

## 2.2. High-fidelity discretization: structure-preserving finite-volume scheme

### 2.2.1. Definition of the scheme and notation

For the time discretization, we introduce a division of the interval $[0, T]$ into subintervals $I_n := [t_{n-1}, t_n]$, $1 \le n \le N_T$, such that $0 = t_0 < t_1 < \dots < t_{N_T} = T$. The time steps are denoted by $\Delta t_n = t_n - t_{n-1}$, $n = 1, \dots, N_T$.

For the space discretization, we consider an admissible mesh of $\Omega$ in the sense of Definition 2.1 from [10] (see also [15]). We recall this definition here for the sake of completeness (see Fig. 1).

**Definition 2.3.** An admissible mesh of $\Omega$ is a triplet $(\mathcal{T}, \mathcal{E}, (x_K)_{K \in \mathcal{T}})$ such that the following conditions are fulfilled.

- Each control volume (or cell) $K \in \mathcal{T}$ is non-empty, open, polyhedral and convex. We assume that $K \cap L = \emptyset$ if $K, L \in \mathcal{T}_h$ with $K \neq L$, while $\bigcup_{K \in \mathcal{T}_h} \bar{K} = \bar{\Omega}$. We denote by $m_K$ the $d$-dimensional Lebesgue measure of the cell $K$.

- Each face $\sigma \in \mathcal{E}$ is closed and is contained in a hyperplane of $\mathbb{R}^d$, with positive $(d-1)$-dimensional Hausdorff (or Lebesgue) measure denoted by $m_\sigma := \mathcal{H}^{d-1}(\sigma) > 0$. We assume that $\mathcal{H}^{d-1}(\sigma \cap \sigma') = 0$ for $\sigma, \sigma' \in \mathcal{E}$ unless $\sigma' = \sigma$. For all $K \in \mathcal{T}$, we assume that there exists a subset $\mathcal{E}_K$ of $\mathcal{E}$ such that $\partial K = \bigcup_{\sigma \in \mathcal{E}_K} \sigma$. Moreover, we suppose that $\bigcup_{K \in \mathcal{T}} \mathcal{E}_K = \mathcal{E}$. Given two distinct cells $K, L \in \mathcal{T}$, the intersection $K \cap L$ either reduces to a single face $\sigma \in \mathcal{E}$ denoted by $K|L$, or its $(d-1)$-dimensional Hausdorff measure is 0.

- The cell-centers $(x_K)_{K \in \mathcal{T}}$ satisfy $x_K \in K$, and are such that, if $K, L \in \mathcal{T}$ share a face $K|L$, then the vector $x_L - x_K$ is orthogonal to $K|L$.

- For the boundary faces $\sigma \in \partial\Omega$ with $\sigma \in \mathcal{E}_K$ for some $K \in \mathcal{T}$, we assume additionally that there exists $x_\sigma \in \sigma$ such that $x_\sigma - x_K$ is orthogonal to $\sigma$. We denote by $\mathcal{E}^{\text{ext}}$ the set of boundary faces, i.e. the set of faces $\sigma \in \mathcal{E}$ such that $\sigma \subset \partial\Omega$ and by $\mathcal{E}^{\text{int}} := \mathcal{E} \setminus \mathcal{E}^{\text{ext}}$ the set of interior faces.

Denote by $h_K$ the diameter of the generic element $K \in \mathcal{T}$ and $h := \max_{K \in \mathcal{T}} h_K$. The number of cells in the mesh $\mathcal{T}$ is denoted by $N_h$. Given a vector $c := (c_K)_{K \in \mathcal{T}} \in \mathbb{R}^{N_h}$, we define for all $K \in \mathcal{T}$ and face $\sigma \in \mathcal{E}_K$, the mirror value $c_{K\sigma}$ of $c_K$ across $\sigma$ by setting:

$$ c_{K\sigma} := c_L \quad \text{if} \quad \sigma = K|L \in \mathcal{E}^{\text{int}} \quad \text{and} \quad c_{K\sigma} := c_K \quad \text{if} \quad \sigma \in \mathcal{E}^{\text{ext}}. \qquad (2.8) $$

We also define the oriented and absolute jumps of $c \in \mathbb{R}^{N_h}$ across any edge by

$$ D_{K\sigma}c := c_{K\sigma} - c_K, \quad \text{and} \quad D_{\sigma}c := |D_{K\sigma}c|, \quad \forall K \in \mathcal{T}, \forall \sigma \in \mathcal{E}_K. \qquad (2.9) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Illustration of an orthogonal mesh.</figcaption></figure>

## 2.2 High-fidelity discretization: structure-preserving finite-volume scheme

### 2.2.1 Definition of the scheme and notation

For the time discretization, we introduce a division of the interval $[0, T]$ into subintervals $I_n := [t_{n-1}, t_n]$, $1 \le n \le N_T$, such that $0 = t_0 < t_1 < \cdots < t_{N_T} = T$. The time steps are denoted by $\Delta t_n = t_n - t_{n-1}$, $n = 1, \cdots, N_T$.

For the space discretization, we consider an admissible mesh of $\Omega$ in the sense of [10, Definition 2.1] (see also [15]). We recall this definition here for the sake of completeness (see Figure 1).

**Definition 2.3.** An admissible mesh of $\Omega$ is a triplet $(\mathcal{T}, \mathcal{E}, (x_K)_{K \in \mathcal{T}})$ such that the following conditions are fulfilled.

* Each control volume (or cell) $K \in \mathcal{T}$ is non-empty, open, polyhedral and convex. We assume that $K \cap L = \emptyset$ if $K, L \in \mathcal{T}_h$ with $K \neq L$, while $\bigcup_{K \in \mathcal{T}_h} \bar{K} = \bar{\Omega}$. We denote by $m_K$ the $d$-dimensional Lebesgue measure of the cell $K$.

* Each face $\sigma \in \mathcal{E}$ is closed and is contained in a hyperplane of $\mathbb{R}^d$, with positive ($d-1$)-dimensional Hausdorff (or Lebesgue) measure denoted by $m_\sigma := \mathcal{H}^{d-1}(\sigma) > 0$. We assume that $\mathcal{H}^{d-1}(\sigma \cap \sigma') = 0$ for $\sigma, \sigma' \in \mathcal{E}$ unless $\sigma' = \sigma$. For all $K \in \mathcal{T}$, we assume that there exists a subset $\mathcal{E}_K$ of $\mathcal{E}$ such that $\partial K = \bigcup_{\sigma \in \mathcal{E}_K} \sigma$. Moreover, we suppose that $\bigcup_{K \in \mathcal{T}} \mathcal{E}_K = \mathcal{E}$. Given two distinct cells $K, L \in \mathcal{T}$, the intersection $K \cap L$ either reduces to a single face $\sigma \in \mathcal{E}$ denoted by $K|L$, or its ($d-1$)-dimensional Hausdorff measure is 0.

* The cell-centers $(x_K)_{K \in \mathcal{T}}$ satisfy $x_K \in K$, and are such that, if $K, L \in \mathcal{T}$ share a face $K|L$, then the vector $x_L - x_K$ is orthogonal to $K|L$.

* For the boundary faces $\sigma \in \partial\Omega$ with $\sigma \in \mathcal{E}_K$ for some $K \in \mathcal{T}$, we assume additionally that there exists $x_\sigma \in \sigma$ such that $x_\sigma - x_K$ is orthogonal to $\sigma$. We denote by $\mathcal{E}^{\text{ext}}$ the set of boundary faces, i.e. the set of faces $\sigma \in \mathcal{E}$ such that $\sigma \subset \partial\Omega$ and by $\mathcal{E}^{\text{int}} := \mathcal{E} \setminus \mathcal{E}^{\text{ext}}$ the set of interior faces.

Denote by $h_K$ the diameter of the generic element $K \in \mathcal{T}$ and $h := \max_{K \in \mathcal{T}} h_K$. The number of cells in the mesh $\mathcal{T}$ is denoted by $N_h$. Given a vector $c := (c_K)_{K \in \mathcal{T}} \in \mathbb{R}^{N_h}$, we define for all $K \in \mathcal{T}$ and face
<figure><img src="image_1.png" /><figcaption>FIGURE 1. Illustration of an orthogonal mesh.</figcaption></figure>

**Remark 2.4.** For all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}$ such that $\sigma = K|L$,

$$D_{L\sigma}c = c_{L\sigma} - c_L = c_K - c_L = -D_{K\sigma}c. \qquad (2.10)$$

A crucial property for meshes as described in Definition 2.3 is the discrete integration by parts formula which is provided further in Lemma 3.2.

Given $\sigma \in \mathcal{E}$, we also define by

$$d_{\sigma} := \begin{cases} |x_K - x_L| & \text{if } \sigma = K|L \in \mathcal{E}^{\text{int}}, \\ |x_K - x_{\sigma}| & \text{if } \sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{ext}}, \end{cases} \quad \text{and} \quad \tau_{\sigma} := \frac{m_{\sigma}}{d_{\sigma}}.$$

Let $\mu \in \mathcal{P}$. Using the two point flux approximation finite volume scheme introduced in [10], the unknown of the model (2.1) is discretized using a couple of constant values per cell: for all $1 \le n \le N_T$, we let

$$U_{\mu}^{n} := (u_{\mu,i,K}^{n})_{K \in \mathcal{T}, i \in [1, N_{s}]} \in \mathbb{R}^{N_{h} \times N_{s}},$$

where $u_{\mu,i,K}^n$ is a numerical approximation of $\frac{1}{m_K} \int_K u_{\mu,i}(t_n, x) dx$.

The finite volume scheme we use is then the following: for $U_{\mu}^0 := (u_{i,K}^0)_{K \in \mathcal{T}, i \in [1, N_s]} \in \mathbb{R}^{N_h \times N_s}$ given, where $u_{i,K}^0 = \frac{1}{m_K} \int_K u_i^0(x) dx$, for all $1 \le n \le N_T$, we find $U_{\mu}^n \in \mathbb{R}^{N_h \times N_s}$ satisfying

$$m_K \frac{u_{\mu,i,K}^n - u_{\mu,i,K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K} F_{\mu,i,K\sigma}(U_{\mu}^n) = 0, \quad \forall K \in \mathcal{T}, \forall i \in [1, N_s] \qquad (2.11)$$

where for all $U := (u_{i,K})_{i \in [1, N_s], K \in \mathcal{T}} \in \mathbb{R}^{N_h \times N_s}$, we will denote by

$$F_{\mu,i,K\sigma}(U) := -a_{\mu}^{*}\tau_{\sigma}D_{K\sigma}u_{i} - \tau_{\sigma}\left(\sum_{j=1, j \neq i}^{N_s}(a_{ij} - a_{\mu}^{*})(u_{j,\sigma}D_{K\sigma}u_{i} - u_{i,\sigma}D_{K\sigma}u_{j})\right) \quad \text{for } \sigma \in \mathcal{E}_{K}^{\text{int}} \qquad (2.12)$$

where

$$u_{i,\sigma} := \begin{cases} 0 & \text{if } \min(u_{i,K}, u_{i,K\sigma}) < 0, \\ \frac{u_{i,K} - u_{i,K\sigma}}{\ln(u_{i,K}) - \ln(u_{i,K\sigma})} & \text{if } u_{i,K} = u_{i,K\sigma} \ge 0, \\ \frac{u_{i,K} - u_{i,K\sigma}}{\ln(u_{i,K}) - \ln(u_{i,K\sigma})} & \text{otherwise.} \end{cases} \qquad (2.13)$$
$\sigma \in \mathcal{E}_K$, the mirror value $c_{K\sigma}$ of $c_K$ across $\sigma$ by setting:

$$c_{K\sigma} := c_L \quad \text{if} \quad \sigma = K|L \in \mathcal{E}^{\text{int}} \quad \text{and} \quad c_{K\sigma} := c_K \quad \text{if} \quad \sigma \in \mathcal{E}^{\text{ext}}. \qquad (2.8)$$

We also define the oriented and absolute jumps of $c \in \mathbb{R}^{N_h}$ across any edge by

$$D_{K\sigma}c := c_{K\sigma} - c_K, \quad \text{and} \quad D_{\sigma}c := |D_{K\sigma}c|, \quad \forall K \in \mathcal{T}, \forall \sigma \in \mathcal{E}_K. \qquad (2.9)$$

**Remark 2.4.** For all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}$ such that $\sigma = K|L$,

$$D_{L\sigma}c = c_{L\sigma} - c_L = c_K - c_L = -D_{K\sigma}c. \qquad (2.10)$$

A crucial property for meshes as described in Definition 2.3 is the discrete integration by parts formula which is provided further in Lemma 3.2.

Given $\sigma \in \mathcal{E}$, we also define by

$$d_{\sigma} := \begin{cases} |x_K - x_L| & \text{if } \sigma = K|L \in \mathcal{E}^{\text{int}}, \\ |x_K - x_{\sigma}| & \text{if } \sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{ext}}, \end{cases} \quad \text{and} \quad \tau_{\sigma} := \frac{m_{\sigma}}{d_{\sigma}}.$$

Let $\mu \in \mathcal{P}$. Using the two point flux approximation finite volume scheme introduced in [10], the unknown of the model (2.1) is discretized using a couple of constant values per cell: for all $1 \le n \le N_T$, we let

$$U_{\mu}^{n} := (u_{\mu,i,K}^{n})_{K \in \mathcal{T}, i \in [1, N_s]} \in \mathbb{R}^{N_h \times N_s},$$

where $u_{\mu,i,K}^n$ is a numerical approximation of $\frac{1}{m_K} \int_K u_{\mu,i}(t_n, x) dx$.

The finite volume scheme we use is then the following: for $U_{\mu}^0 := (u_{i,K}^0)_{K \in \mathcal{T}, i \in [1, N_s]} \in \mathbb{R}^{N_h \times N_s}$ given, where $u_{i,K}^0 = \frac{1}{m_K} \int_K u_{i}^0(x) dx$, for all $1 \le n \le N_T$, we find $U_{\mu}^n \in \mathbb{R}^{N_h \times N_s}$ satisfying

$$m_K \frac{u_{\mu,i,K}^n - u_{\mu,i,K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K} F_{\mu,i,K\sigma}(U_{\mu}^n) = 0, \quad \forall K \in \mathcal{T}, \forall i \in [1, N_s] \qquad (2.11)$$

where for all $U := (u_{i,K})_{i \in [1, N_s], K \in \mathcal{T}} \in \mathbb{R}^{N_h \times N_s}$, we will denote by

$$F_{\mu,i,K\sigma}(U) := -a_{\mu}^{*}\tau_{\sigma}D_{K\sigma}u_{i} - \tau_{\sigma}\left(\sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^{*})(u_{j,\sigma}D_{K\sigma}u_{i} - u_{i,\sigma}D_{K\sigma}u_{j})\right) \quad \text{for } \sigma \in \mathcal{E}_{K}^{\text{int}} \qquad (2.12)$$

where

$$u_{i,\sigma} := \begin{cases} 0 & \text{if } \min(u_{i,K}, u_{i,K\sigma}) < 0, \\ u_{i,K} & \text{if } u_{i,K} = u_{i,K\sigma} \ge 0, \\ \frac{u_{i,K} - u_{i,K\sigma}}{\ln(u_{i,K}) - \ln(u_{i,K\sigma})} & \text{otherwise.} \end{cases} \qquad (2.13)$$

**Remark 2.5.** Observe that the numerical flux defined in (2.12) is conservative in the sense that for $\sigma \in \mathcal{E}^{\text{int}}$, $\sigma = K|L$, and all $U \in \mathbb{R}^{N_h \times N_s}$,

$$F_{\mu,i,L\sigma}(U) = -F_{\mu,i,K\sigma}(U).$$

This implies in particular that, for all $1 \le i \le N_s$ and for all $U \in \mathbb{R}^{N_h \times N_s}$

$$\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu,i,K\sigma}(U) = 0. \qquad (2.14)$$
**Remark 2.5.** Observe that the numerical flux defined in (2.12) is conservative in the sense that for $\sigma \in \mathcal{E}^{\text{int}}$, $\sigma = K|L$, and all $U \in \mathbb{R}^{N_h \times N_s}$,

$$F_{\mu,i,L\sigma}(U) = -F_{\mu,i,K\sigma}(U).$$

This implies in particular that, for all $1 \le i \le N_s$ and for all $U \in \mathbb{R}^{N_h \times N_s}$

$$\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu,i,K\sigma}(U) = 0. \qquad (2.14)$$

### 2.2.2. Properties of the discrete solution

It is proved in [10] that the numerical scheme defined in the previous section is well-defined and preserves analogous properties as (P1)–(P2)–(P3) at the discrete level. These properties have already been demonstrated in Theorem 2.2 from [10]. We recall it in Proposition 2.6 below for the sake of clarity.

**Proposition 2.6** ([10], Thm. 2.2). Let $u^0 \in L^\infty(\Omega; \mathcal{A})$ so that $U^0 \in \mathbb{R}^{N_h \times N_s}$ is defined as $u_{i,K}^0 := \frac{1}{m_K} \int_K u_i^0(x) dx$. Then, for all $1 \le n \le N_T$, there exists at least one solution $U_\mu^n \in \mathbb{R}^{N_h \times N_s}$ to the numerical scheme (2.11) which satisfies the following properties:

(P1$^h$) **mass conservation:**

$$\sum_{K \in \mathcal{T}} m_K u_{\mu,i,K}^n = \sum_{K \in \mathcal{T}} m_K u_{\mu,i,K}^0 = \int_{\Omega} u_i^0(x) dx \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T]. \qquad (2.15)$$

(P2$^h$) **volume-filling constraint:**

$$u_{\mu,i,K}^n > 0 \quad \forall K \in \mathcal{T}, \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T], \qquad (2.16)$$

and

$$\sum_{i=1}^{N_s} u_{\mu,i,K}^{n} = 1 \quad \forall K \in \mathcal{T}, \quad \forall n \in [0, N_T]. \qquad (2.17)$$

(P3$^h$) **entropy-entropy dissipation inequality:** for all $U := (u_{i,K})_{K \in \mathcal{T}, [1, N_s]} \in \mathbb{R}_+^{N_h \times N_s}$, let us define

$$E_{\mathcal{T}}(U) := \sum_{K \in \mathcal{T}} \sum_{i=1}^{N_s} m_K u_{i,K} \ln(u_{i,K}). \qquad (2.18)$$

Then, it holds that

$$\frac{1}{\Delta t_n} (E_{\mathcal{T}}(U_{\mu}^n) - E_{\mathcal{T}}(U_{\mu}^{n-1})) + a_{\mu}^{*} \sum_{\sigma \in \mathcal{E}} \sum_{i=1}^{N_s} u_{\mu,i,\sigma}^{n} |D_{K\sigma} \ln(u_{\mu,i}^{n})|^2 \le 0, \quad \forall n \in [1, N_T]. \qquad (2.19)$$

### 2.2.3. Newton resolution

Let us define for all $\mu \in \mathcal{P}$, $\forall K \in \mathcal{T}$, $\forall i \in [1, N_s]$, $\forall 1 \le n \le N_T$, the nonlinear function $G_{\mu,i,K}^n: \mathbb{R}^{N_s \times N_h} \to \mathbb{R}$ defined as follows: for all $U = (u_{i,K})_{1 \le i \le N_s, K \in \mathcal{T}} \in \mathbb{R}^{N_h \times N_s}$,

$$G_{\mu,i,K}^{n}(U) := m_{K} \frac{u_{i,K} - u_{\mu,i,K}^{n-1}}{\Delta t_{n}} + \sum_{\sigma \in \mathcal{E}_{K} \cap \mathcal{E}^{\text{int}}} F_{\mu,i,K\sigma}(U), \qquad (2.20)$$

and define the application $G_{\mu}^n: \mathbb{R}^{N_s \times N_h} \to \mathbb{R}^{N_s \times N_h}$ as follows: for all $U \in \mathbb{R}^{N_s \times N_h}$,

$$G_{\mu}^{n}(U) := (G_{\mu,i,K}^{n}(U))_{K \in \mathcal{T}, [1, N_s]}.$$
### 2.2.2 Properties of the discrete solution

It is proved in [10] that the numerical scheme defined in the previous section is well-defined and preserves analogous properties as (P1)-(P2)-(P3) at the discrete level. These properties have already been demonstrated in [10, Theorem 2.2]. We recall it in Proposition 2.6 below for the sake of clarity.

**Proposition 2.6 (Theorem 2.2 of [10].)** Let $\mathbf{u}^0 \in L^\infty(\Omega; \mathcal{A})$ so that $\mathbf{U}^0 \in \mathbb{R}^{N_h \times N_s}$ is defined as $u_{i,K}^0 := \frac{1}{m_K} \int_K u_i^0(x) dx$. Then, for all $1 \le n \le N_T$, there exists at least one solution $\mathbf{U}_\mu^n \in \mathbb{R}^{N_h \times N_s}$ to the numerical scheme (2.11) which satisfies the following properties:

(P1$^h$) mass conservation:

$$ \sum_{K \in \mathcal{T}} m_K u_{\mu,i,K}^n = \sum_{K \in \mathcal{T}} m_K u_{\mu,i,K}^0 = \int_{\Omega} u_i^0(x) dx \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T]. \quad (2.15) $$

(P2$^h$) volume-filling constraint:

$$ u_{\mu,i,K}^n > 0 \quad \forall K \in \mathcal{T}, \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T], \quad (2.16) $$

and

$$ \sum_{i=1}^{N_s} u_{\mu,i,K}^n = 1 \quad \forall K \in \mathcal{T}, \quad \forall n \in [0, N_T]. \quad (2.17) $$

(P3$^h$) entropy-entropy dissipation inequality: For all $\mathbf{U} := (u_{i,K})_{K \in \mathcal{T}, [1, N_s]} \in (\mathbb{R}_+)^{N_h \times N_s}$, let us define

$$ E_{\mathcal{T}}(\mathbf{U}) := \sum_{K \in \mathcal{T}} \sum_{i=1}^{N_s} m_K u_{i,K} \ln(u_{i,K}). \quad (2.18) $$

Then, it holds that

$$ \frac{1}{\Delta t_n} (E_{\mathcal{T}}(\mathbf{U}_{\mu}^n) - E_{\mathcal{T}}(\mathbf{U}_{\mu}^{n-1})) + a_{\mu}^{*} \sum_{\sigma \in \mathcal{E}} \sum_{i=1}^{N_s} u_{\mu,i,\sigma}^{n} |D_{K\sigma} \ln(\mathbf{u}_{\mu,i}^{n})|^2 \le 0, \quad \forall n \in [1, N_T]. \quad (2.19) $$

### 2.2.3 Newton resolution

Let us define for all $\mu \in \mathcal{P}$, $\forall K \in \mathcal{T}$, $\forall i \in [1, N_s]$, $\forall 1 \le n \le N_T$, the nonlinear function $G_{\mu,i,K}^n: \mathbb{R}^{N_s \times N_h} \to \mathbb{R}$ defined as follows: for all $\mathbf{U} = (u_{i,K})_{1 \le i \le N_s, K \in \mathcal{T}} \in \mathbb{R}^{N_h \times N_s}$,

$$ G_{\mu,i,K}^n(\mathbf{U}) := m_K \frac{u_{i,K} - u_{\mu,i,K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu,i,K\sigma}(\mathbf{U}), \quad (2.20) $$

and define the application $G_{\mu}^n: \mathbb{R}^{N_s \times N_h} \to \mathbb{R}^{N_s \times N_h}$ as follows: for all $\mathbf{U} \in \mathbb{R}^{N_s \times N_h}$,

$$ G_{\mu}^{n}(\mathbf{U}) := (G_{\mu,i,K}^{n}(\mathbf{U}))_{K \in \mathcal{T}, [1, N_s]} $$

A solution $\mathbf{U}_\mu^n \in \mathbb{R}^{N_h \times N_s}$ to the scheme (2.11) can then be rewritten equivalently as a solution to the following nonlinear system

$$ G_{\mu}^{n}(\mathbf{U}_{\mu}^{n}) = 0, \quad (2.21) $$

which we solve in practice using a classical Newton procedure [14, 29, 5].

The resolution of problem (2.21) can be very expensive due to the number of mesh elements, the number of species, the number of time steps, and the number of cross diffusion matrix to be tested. It is then crucial to speed up the computation of the numerical solution. The aim of our work is to propose a new reduced-order model to efficiently compute approximations of $\mathbf{U}_\mu^n$ for a large number of values of $\mu \in \mathcal{P}$ while preserving analogous properties as those listed in Proposition 2.6.
A solution $U_{\mu}^n \in \mathbb{R}^{N_h \times N_s}$ to the scheme (2.11) can then be rewritten equivalently as a solution to the following nonlinear system

$$G_{\mu}^{n}(U_{\mu}^{n}) = 0, \qquad (2.21)$$

which we solve in practice using a classical Newton procedure [5, 14, 29].

The resolution of problem (2.21) can be very expensive due to the number of mesh elements, the number of species, the number of time steps, and the number of cross diffusion matrix to be tested. It is then crucial to speed up the computation of the numerical solution. The aim of our work is to propose a new reduced-order model to efficiently compute approximations of $U_{\mu}^n$ for a large number of values of $\mu \in \mathcal{P}$ while preserving analogous properties as those listed in Proposition 2.6.

### 3. STRUCTURE-PRESERVING REDUCED-ORDER MODEL

In Section 3.2, we expose the structure-preserving reduced-order model we propose in order to efficiently compute numerical approximations of $U_{\mu}^n$ for $1 \le n \le N_T$ and $\mu \in \mathcal{P}$. For the sake of comparison, we first describe how a standard reduced-order modeling method based on a POD approach would work in Section 3.1. We will compare the two approaches in our numerical experiments presented in Section 4.

In the following, we introduce $\mathcal{P}_{train} \subset \mathcal{P}$ a finite training set of parameter values, and assume that we have computed high-fidelity solutions $(U_{\mu}^n)_{1 \le n \le N_T}$ for all values of the parameter $\mu$ in the training set $\mathcal{P}_{train}$ according to the finite volume procedure described in Section 2.2. We will denote by $p \in \mathbb{N}^*$ the cardinality of the set $\mathcal{P}_{train}$.

#### 3.1. Standard POD reduced model

The aim of this section is to describe the classical POD reduced-order model (POD-ROM) we will use in our numerical experiments to compare with the structure-preserving reduced-order model (SP-ROM) we will introduce in Section 3.2.

Given the family of snapshot solutions $(U_{\mu}^n)_{1 \le n \le N_T, \mu \in \mathcal{P}_{train}} \in \mathbb{R}^{N_h \times N_s}$, for all $r \in \mathbb{N}^*$, we denote by $V^1, \dots, V^r \in \mathbb{R}^{N_h \times N_s}$ the $r$ first POD modes of this family. These first POD modes are computed as follows:

(1) We store all the snapshot solutions in the matrix $M \in \mathbb{R}^{N_h \times N_s, p \times N_T}$.

(2) We employ a SVD decomposition on the matrix $M$ as

$$M = VSW^T.$$

Here, the matrix $V \in \mathbb{R}^{N_s \times N_h, N_s \times N_h}$ and the matrix $W \in \mathbb{R}^{N_T \times p, N_T \times p}$ are orthogonal, and the diagonal matrix $S \in \mathbb{R}^{N \times N_h, N_T \times p}$ is rectangular. More precisely, $S_{ii} := \sqrt{\sigma_i}$ for $1 \le i \le \min(N_s \times N_h, N_T \times p)$ where $\sigma_i$ are the eigenvalues of the matrix $MM^T$ ranged in decreasing order.

(3) From the matrix $V$, we extract the matrix $V^r \in \mathbb{R}^{N_s \times N_h, r}$ which is the submatrix of $V$ composed by its first $r$ rows. Note that $r \ll \min(N_h \times N_s, N_T \times p)$. The precision of the integer $r \in \mathbb{N}^*$ is determined thanks to the criterion $\sum_{k \ge r+1} \sigma_k \le \varepsilon$ for $\varepsilon \ge 0$ a fixed tolerance.

In the classical POD-ROM, for all $\mu \in \mathcal{P}$ and $1 \le n \le N_T$, we compute a reduced-order model approximation $\tilde{U}_{\mu}^n \in \mathbb{R}^{N_s \times N_h}$ of $U_{\mu}^n$ as a linear combination of $V^1, \dots, V^r$

$$\tilde{U}_{\mu}^{n} := \sum_{k=1}^{r} \tilde{c}_{\mu}^{k,n} V^{k}, \qquad (3.1)$$

where the coefficients $\tilde{c}_{\mu}^n := (\tilde{c}_{\mu}^{k,n})_{1 \le k \le r} \in \mathbb{R}^r$ are computed as follows.
# 3 Structure-preserving reduced-order model

In Section 3.2, we expose the structure-preserving reduced-order model we propose in order to efficiently compute numerical approximations of $U_{\mu}^n$ for $1 \le n \le N_T$ and $\mu \in \mathcal{P}$. For the sake of comparison, we first describe how a standard reduced-order modeling method based on a POD approach would work in Section 3.1. We will compare the two approaches in our numerical experiments presented in Section 4.

In the following, we introduce $\mathcal{P}_{\text{train}} \subset \mathcal{P}$ a finite training set of parameter values, and assume that we have computed high-fidelity solutions $(U_{\mu}^n)_{1 \le n \le N_T}$ for all values of the parameter $\mu$ in the training set $\mathcal{P}_{\text{train}}$ according to the finite volume procedure described in Section 2.2. We will denote by $p \in \mathbb{N}^*$ the cardinality of the set $\mathcal{P}_{\text{train}}$.

## 3.1 Standard POD reduced model

The aim of this section is to describe the classical POD reduced-order model (POD-ROM) we will use in our numerical experiments to compare with the structure-preserving reduced-order model (SP-ROM) we will introduce in Section 3.2.

Given the family of snapshot solutions $(U_{\mu}^n)_{1 \le n \le N_T, \mu \in \mathcal{P}_{\text{train}}} \in \mathbb{R}^{N_h \times N_s}$, for all $r \in \mathbb{N}^*$, we denote by $\mathbf{V}^1, \dots, \mathbf{V}^r \in \mathbb{R}^{N_h \times N_s}$ the $r$ first POD modes of this family. These first POD modes are computed as follows:

1. We store all the snapshot solutions in the matrix $\mathbf{M} \in \mathbb{R}^{N_h \times N_s, p \times N_T}$.

2. We employ a SVD decomposition on the matrix $\mathbf{M}$ as

$$ \mathbf{M} = \mathbf{V}\mathbf{S}\mathbf{W}^T. $$

Here, the matrix $\mathbf{V} \in \mathbb{R}^{N_s \times N_h, N_s \times N_h}$ and the matrix $\mathbf{W} \in \mathbb{R}^{N_T \times p, N_T \times p}$ are orthogonal, and the diagonal matrix $\mathbf{S} \in \mathbb{R}^{N \times N_h, N_T \times p}$ is rectangular. More precisely, $\mathbf{S}_{ii} := \sqrt{\sigma_i}$ for $1 \le i \le \min(N_s \times N_h, N_T \times p)$ where $\sigma_i$ are the eigenvalues of the matrix $\mathbf{MM}^T$ ranged in decreasing order.

3. From the matrix $\mathbf{V}$, we extract the matrix $\mathbf{V}^r \in \mathbb{R}^{N_s \times N_h, r}$ which is the submatrix of $\mathbf{V}$ composed by its first $r$ rows. Note that $r \ll \min(N_h \times N_s, N_T \times p)$. The precision of the integer $r \in \mathbb{N}^*$ is determined thanks to the criterion $\sum_{k \ge r+1} \sigma_k \le \varepsilon$ for $\varepsilon \ge 0$ a fixed tolerance.

In the classical POD-ROM, for all $\mu \in \mathcal{P}$ and $1 \le n \le N_T$, we compute a reduced-order model approximation $\tilde{U}_{\mu}^n \in \mathbb{R}^{N_s \times N_h}$ of $U_{\mu}^n$ as a linear combination of $\mathbf{V}^1, \dots, \mathbf{V}^r$

$$ \tilde{U}_{\mu}^n := \sum_{k=1}^{r} \tilde{c}_{\mu}^{k,n} \mathbf{V}^k, \qquad (3.1) $$

where the coefficients $\tilde{c}_{\mu}^n := (\tilde{c}_{\mu}^{k,n})_{1 \le k \le r} \in \mathbb{R}^r$ are computed as follows.

Let us denote by $(V_{i,K}^k)_{K \in \mathcal{T}, 1 \le i \le N_s} \in \mathbb{R}^{N_h \times N_s}$ the components of $\mathbf{V}^k$ for all $1 \le k \le r$. For any vector $\mathbf{c} := (c^k)_{1 \le k \le r} \in \mathbb{R}^r$, let us denote by $\tilde{U}(\mathbf{c}) := (\tilde{u}_{i,K}(\mathbf{c}))_{K \in \mathcal{T}, [1, N_s]} \in \mathbb{R}^{N_h \times N_s}$ the vector defined so that

$$ \tilde{u}_{i,K}(\mathbf{c}) := \sum_{k=1}^{r} c^k V_{i,K}^k \quad \forall i \in [1, N_s], \quad \forall K \in \mathcal{T}. \qquad (3.2) $$

Let us now define for all $\mu \in \mathcal{P}$, $1 \le n \le N_T$, $[1, N_s]$ and $K \in \mathcal{T}$, the application $\tilde{G}_{\mu,i,K}^n : \mathbb{R}^r \to \mathbb{R}$ defined by

$$ \tilde{G}_{\mu,i,K}^n(\mathbf{c}) := m_K \frac{\tilde{u}_{i,K}(\mathbf{c}) - \tilde{u}_{\mu,i,K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K} F_{\mu,i,K\sigma}(\tilde{U}(\mathbf{c})), $$

where $\tilde{u}_{\mu,i,K}^0(\mathbf{c}) := u_{i,K}^0$, and by $\tilde{G}_{\mu}^n(\mathbf{c}) := (\tilde{G}_{\mu,i,K}^n(\mathbf{c}))_{K \in \mathcal{T}, [1, N_s]}$.
Let us denote by $(V_{i,K}^k)_{K \in \mathcal{T}, 1 \le i \le N_s} \in \mathbb{R}^{N_h \times N_s}$ the components of $\mathbf{V}^k$ for all $1 \le k \le r$. For any vector $\mathbf{c} := (c^k)_{1 \le k \le r} \in \mathbb{R}^r$, let us denote by $\tilde{\mathbf{U}}(\mathbf{c}) := (\tilde{u}_{i,K}(\mathbf{c}))_{K \in \mathcal{T}, [1, N_s]} \in \mathbb{R}^{N_h \times N_s}$ the vector defined so that

$$ \tilde{u}_{i,K}(\mathbf{c}) := \sum_{k=1}^{r} c^k V_{i,K}^k \quad \forall i \in [1, N_s], \quad \forall K \in \mathcal{T}. \qquad (3.2) $$

Let us now define for all $\mu \in \mathcal{P}$, $1 \le n \le N_T$, $[1, N_s]$ and $K \in \mathcal{T}$, the application $\tilde{G}_{\mu,i,K}^n : \mathbb{R}^r \to \mathbb{R}$ defined by

$$ \tilde{G}_{\mu,i,K}^n(\mathbf{c}) := m_K \frac{\tilde{u}_{i,K}(\mathbf{c}) - \tilde{u}_{\mu,i,K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K} F_{\mu,i,K\sigma}(\tilde{\mathbf{U}}(\mathbf{c})), $$

where $\tilde{u}_{\mu,i,K}^0(\mathbf{c}) := u_{i,K}^0$, and by $\tilde{G}_{\mu}^n(\mathbf{c}) := (\tilde{G}_{\mu,i,K}^n(\mathbf{c}))_{K \in \mathcal{T}, [1, N_s]}$.

Let us finally denote by $\tilde{H}_{\mu}^n : \mathbb{R}^r \to \mathbb{R}^r$ the application such that, for all $\mathbf{c} \in \mathbb{R}^r$, $\tilde{H}_{\mu}^n = (\tilde{H}_{\mu}^{n,l}(\mathbf{c}))_{1 \le l \le r}$ where for all $1 \le l \le r$

$$ \tilde{H}_{\mu}^{n,l}(\mathbf{c}) := \langle \mathbf{V}^l, \tilde{G}_{\mu}^n(\mathbf{c}) \rangle. $$

Here, $\langle \cdot, \cdot \rangle$ denotes the euclidean $\ell_2$ scalar product of $\mathbb{R}^{N_h \times N_s}$ meaning that for all $(\mathbf{d}, \mathbf{e}) \in \mathbb{R}^{N_h \times N_s} \times \mathbb{R}^{N_h \times N_s}$ $\langle \mathbf{d}, \mathbf{e} \rangle = \sum_{\alpha=1}^{N_h \times N_s} \mathbf{d}_{\alpha} \mathbf{e}_{\alpha}$. Then, the vector $\tilde{c}_{\mu}^n \in \mathbb{R}^r$ is defined as a solution to

$$ \tilde{H}_{\mu}^{n}(\tilde{c}_{\mu}^{n}) = 0, \qquad (3.3) $$

and the numerical approximation $\tilde{U}_{\mu}^n$ of $U_{\mu}^n$ given by the classical POD-ROM is given by (3.1). A solution to problem (3.3) is computed using a classical Newton algorithm.

Let us highlight here that the standard POD reduced-order model may not preserve the structural properties of the solutions of the cross-diffusion system we mentioned in the previous sections, such as non-negativeness, mass conservation, volumic constraint and decay of the entropy with respect to time evolution.

## 3.2. Structure-preserving (SP) reduced-order model

For all $\mu \in \mathcal{P}$ and $0 \le n \le N_T$, let us denote by

$$ z_{\mu,j,K}^n := \ln(u_{\mu,j,K}^n) - \ln(u_{\mu,N_s,K}^n) \quad \forall j \in [1, N_s - 1] \quad \forall K \in \mathcal{T}. \qquad (3.4) $$

We then denote by $\mathbf{V}^1, \dots, \mathbf{V}^r \in \mathbb{R}^{N_h \times (N_s-1)}$ the r first POD modes of the family $(\mathbf{Z}_{\mu}^n)_{\mu \in \mathcal{P}_{\text{train}}, 1 \le n \le N_T} \in \mathbb{R}^{N_h \times (N_s-1)}$ where $\mathbf{Z}_{\mu}^n := (z_{\mu,j,K}^n)_{K \in \mathcal{T}, j \in [1, N_s-1]}$. In addition, for all $j \in [1, N_s-1]$, we define

$$ \mathbf{V}^{r+j} := \mathbf{I}^j, $$

where $\mathbf{I}^j \in \mathbb{R}^{(N_s-1) \times N_h}$ is the vector defined by

$$ \mathbf{I}^j := (\mathbf{I}_{j',K}^j)_{K \in \mathcal{T}, j' \in [1, N_s-1]} \quad \text{with} \quad \mathbf{I}_{j',K}^j = \delta_{j',k}^j. \qquad (3.5) $$

In other words, $\mathbf{I}_{j',K}^j = 1$ if $j = j'$ and 0 otherwise. We also denote by $r_* := r + N_s - 1$. The principle of the SP-ROM method we present here is to construct for all $\mu \in \mathcal{P}$, a reduced-order model approximation $\tilde{\mathbf{Z}}_{\mu}^n \in \mathbb{R}^{N_h \times (N_s-1)}$ of $\mathbf{Z}_{\mu}^n$ that reads as a linear combination of $\mathbf{V}^1, \dots, \mathbf{V}^{r_*}$ given by

$$ \tilde{\mathbf{Z}}_{\mu}^{n} := \sum_{k=1}^{r_{*}} \tilde{c}_{\mu}^{k,n} \mathbf{V}^{k}, \qquad (3.6) $$
Let us finally denote by $\tilde{H}_{\mu}^{n}: \mathbb{R}^{r} \to \mathbb{R}^{r}$ the application such that, for all $c \in \mathbb{R}^{r}$, $\tilde{H}_{\mu}^{n} = (\tilde{H}_{\mu}^{n,l}(c))_{1 \le l \le r}$
where for all $1 \le l \le r$

$$
\tilde{H}_{\mu}^{n,l}(\mathbf{c}) := \langle \mathbf{V}^{l}, \tilde{G}_{\mu}^{n}(\mathbf{c}) \rangle.
$$

Here, $\langle \cdot, \cdot \rangle$ denotes the euclidean $\ell_2$ scalar product of $\mathbb{R}^{N_h \times N_s}$ meaning that for all $(d, e) \in \mathbb{R}^{N_h \times N_s} \times \mathbb{R}^{N_h \times N_s}$
$\langle d, e \rangle = \sum_{\alpha=1}^{N_h \times N_s} d_{\alpha} e_{\alpha}$. Then, the vector $\tilde{c}_{\mu}^{n} \in \mathbb{R}^{r}$ is defined as a solution to

$$
\tilde{H}_{\mu}^{n}(\tilde{c}_{\mu}^{n}) = 0, \tag{3.3}
$$

and the numerical approximation $\tilde{U}_{\mu}^{n}$ of $U_{\mu}^{n}$ given by the classical POD-ROM is given by (3.1). A solution
to problem (3.3) is computed using a classical Newton algorithm.

Let us highlight here that the standard POD reduced-order model may not preserve the structural
properties of the solutions of the cross-diffusion system we mentioned in the previous sections, such as
non-negativeness, mass conservation, volumic constraint and decay of the entropy with respect to time
evolution.

3.2 Structure-preserving (SP) reduced-order model

For all $\mu \in \mathcal{P}$ and $0 \le n \le N_T$, let us denote by

$$
z_{\mu,j,K}^{n} := \ln(u_{\mu,j,K}^{n}) - \ln(u_{\mu,N_s,K}^{n}) \quad \forall j \in [1, N_s - 1] \quad \forall K \in \mathcal{T}. \quad (3.4)
$$

We then denote by $V^1, \dots, V^r \in \mathbb{R}^{N_h \times (N_s-1)}$ the $r$ first POD modes of the family $(Z_\mu^n)_{\mu \in \mathcal{P}_{\text{train}}, 1 \le n \le N_T} \in \mathbb{R}^{N_h \times (N_s-1)}$ where $Z_\mu^n := (z_{\mu,j,K}^n)_{K \in \mathcal{T}, j \in [1, N_s-1]}$. In addition, for all $j \in [1, N_s-1]$, we define

$$
V^{r+j} := I^j,
$$

where $I^j \in \mathbb{R}^{(N_s-1) \times N_h}$ is the vector defined by

$$
I^j := (I_{j',K}^j)_{K \in \mathcal{T}, j' \in [1, N_s-1]} \quad \text{with} \quad I_{j',K}^j = \delta_{j',K}^j. \tag{3.5}
$$

In other words, $I_{j',K}^j = 1$ if $j = j'$ and $0$ otherwise. We also denote by $r_* := r + N_s - 1$. The principle of the SP-ROM method we present here is to construct for all $\mu \in \mathcal{P}$, a reduced-order model approximation $\bar{Z}_\mu^n \in \mathbb{R}^{N_h \times (N_s-1)}$ of $Z_\mu^n$ that reads as a linear combination of $V^1, \dots, V^{r_*}$ given by

$$
\bar{Z}_{\mu}^{n} := \sum_{k=1}^{r_{*}} \bar{c}_{\mu}^{k,n} V^{k}, \qquad (3.6)
$$

where the vector $\bar{c}_{\mu}^{n} := (\bar{c}_{\mu}^{k,n})_{1 \le k \le r_*} \in \mathbb{R}^{r_*}$ is computed as described below. Then, assuming that $\bar{Z}_{\mu}^{n} = (\bar{z}_{\mu,j,K}^{n})_{K \in \mathcal{T}, j \in [1,N_s-1]}$ is computed, a reduced-order model approximation $\bar{U}_{\mu}^{n} = (\bar{u}_{\mu,i,K}^{n})_{K \in \mathcal{T}, [1,N_s-1]}$ of $U_{\mu}^{n}$ is then computed as follows. Let the nonlinear operator $\mathcal{Q}: \mathbb{R}^{N_s-1} \to \mathcal{A}$ be defined as follows: for all $s = (s_i)_{i \in [1, N_s-1]} \in \mathbb{R}^{N_s-1}$, $\mathcal{Q}(s) = ([\mathcal{Q}(s)]_i)_{i \in [1, N_s]}$ where

$$
[\mathcal{Q}(s)]_i = \frac{e^{s_i}}{1 + \sum_{j=1}^{N_s-1} e^{s_j}} \quad \forall i \in [1, N_s - 1] \quad \text{and} \quad [\mathcal{Q}(s)]_{N_s} := \frac{1}{1 + \sum_{j=1}^{N_s-1} e^{s_j}}. \quad (3.7)
$$

We then define the reduced-order model approximation $\bar{U}_{\mu}^{n} := (\bar{u}_{\mu,i,K}^{n})_{K \in \mathcal{T}, i \in [1,N_s]}$ from $\bar{Z}_{\mu}^{n} := (\bar{z}_{\mu,i,K}^{n})_{K \in \mathcal{T}, i \in [1,N_s-1]}$
by

$$
\bar{u}_{\mu, i, K}^{n} := [\mathcal{Q}(\bar{z}_{\mu, K}^{n})]_{i}, \quad (3.8)
$$

where $\bar{z}_{\mu,K}^n = (\bar{z}_{\mu,j,K}^n)_{j \in [1,N_s-1]}$. Let us point out here that $\mathcal{Q}$ defines an injective map from $\mathbb{R}^{N_s-1}$ to $\mathcal{A}$.
Moreover, note that, by construction, a reduced-order model $\bar{U}_{\mu}^n$ necessarily satisfies
where the vector $\bar{c}_{\mu}^n := (\bar{c}_{\mu}^{k,n})_{1\le k\le r_*} \in \mathbb{R}^{r_*}$ is computed as described below. Then, assuming that $\bar{Z}_{\mu}^n = (\bar{z}_{\mu,j,K}^n)_{K \in \mathcal{T}, j \in [1,N_s-1]}$ is computed, a reduced-order model approximation $\bar{U}_{\mu}^n = (\bar{u}_{\mu,i,K}^n)_{K \in \mathcal{T}, j \in [1,N_s-1]}$ of $\bar{U}_{\mu}^n$ is then computed as follows. Let the nonlinear operator $\mathcal{Q} : \mathbb{R}^{N_s-1} \to \mathcal{A}$ be defined by: for all $\mathbf{s} = (s_i)_{i \in [1,N_s-1]} \in \mathbb{R}^{N_s-1}$, $\mathcal{Q}(\mathbf{s}) = ([\mathcal{Q}(\mathbf{s})]_i)_{i \in [1,N_s]}$ where

$$[\mathcal{Q}(\mathbf{s})]_i = \frac{e^{s_i}}{1 + \sum_{j=1}^{N_s-1} e^{s_j}} \quad \forall i \in [1, N_s-1] \quad \text{and} \quad [\mathcal{Q}(\mathbf{s})]_{N_s} := \frac{1}{1 + \sum_{j=1}^{N_s-1} e^{s_j}}. \qquad (3.7)$$

We then define the reduced-order model approximation $\bar{U}_{\mu}^n := (\bar{u}_{\mu,i,K}^n)_{K \in \mathcal{T}, i \in [1,N_s]}$ from $\bar{Z}_{\mu}^n := (\bar{z}_{\mu,i,K}^n)_{K \in \mathcal{T}, i \in [1,N_s-1]}$ by

$$\bar{u}_{\mu,i,K}^n := [\mathcal{Q}(\bar{z}_{\mu,K}^n)]_i, \qquad (3.8)$$

where $\bar{z}_{\mu,K}^n = (\bar{z}_{\mu,j,K}^n)_{j \in [1,N_s-1]}$. Let us point out here that $\mathcal{Q}$ defines an injective map from $\mathbb{R}^{N_s-1}$ to $\mathcal{A}$. Moreover, note that, by construction, a reduced-order model $\bar{U}_{\mu}^n$ necessarily satisfies

(P2<sup>red</sup>) volume-filling constraint:

$$\bar{u}_{\mu,i,K}^n > 0 \quad \forall K \in \mathcal{T}, \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T], \qquad (3.9)$$

and

$$\sum_{i=1}^{N_s} \bar{u}_{\mu,i,K}^n = 1 \quad \forall K \in \mathcal{T}, \quad \forall n \in [0, N_T], \qquad (3.10)$$

which is the analogous version of (P2<sup>h</sup>) for the SP-ROM.

Let us now explain in details how the vector $\bar{c}_{\mu}^n$ is computed. For any $\mathbf{c} := (c^k)_{1\le k\le r_*} \in \mathbb{R}^{r_*}$, let us introduce $\bar{Z}(\mathbf{c}) := (\bar{z}_{j,K}(\mathbf{c}))_{K \in \mathcal{T}, j \in [1,N_s-1]} \in \mathbb{R}^{N_h \times (N_s-1)}$ the vector defined by

$$\bar{Z}(\mathbf{c}) := \sum_{k=1}^{r_*} c^k \mathbf{V}^k.$$

Moreover, we define $\bar{U}(\mathbf{c}) := (\bar{u}_{i,K}(\mathbf{c}))_{K \in \mathcal{T}, i \in [1,N_s]} \in \mathbb{R}^{N_h \times N_s}$ by

$$\forall i \in [1, N_s], \forall K \in \mathcal{T}, \quad \bar{u}_{i,K}(\mathbf{c}) := [\mathcal{Q}(\bar{z}_K(\mathbf{c}))]_i, \qquad (3.11)$$

where $\bar{z}_K(\mathbf{c}) = (\bar{z}_{j,K}(\mathbf{c}))_{j \in [1,N_s-1]}$. Let us now define for all $\mu \in \mathcal{P}$, $1 \le n \le N_T$, $1 \le j \le N_s - 1$ and $K \in \mathcal{T}$, the application $\bar{G}_{\mu,i,K}^n : \mathbb{R}^{r_*} \to \mathbb{R}$ defined by

$$\bar{G}_{\mu,j,K}^n(\mathbf{c}) := m_K \frac{\bar{u}_{j,K}(\mathbf{c}) - \bar{u}_{\mu,j,K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K} F_{\mu,j,K\sigma}(\bar{U}(\mathbf{c})), \qquad (3.12)$$

where for all $i \in [1, N_s]$ and $K \in \mathcal{T}$, $\bar{u}_{\mu,i,K}^0(\mathbf{c}) := u_{i,K}^0$, and by $\bar{G}_{\mu}^n(\mathbf{c}) := (\bar{G}_{\mu,j,K}^n(\mathbf{c}))_{K \in \mathcal{T}, j \in [1, N_s - 1]}$. Let us finally denote by $\bar{H}_{\mu} : \mathbb{R}^{r_*} \to \mathbb{R}^{r_*}$ the application such that, for all $\mathbf{c} \in \mathbb{R}^{r_*}$, $\bar{H}_{\mu}^n = (\bar{H}_{\mu}^{n,l}(\mathbf{c}))_{1\le l\le r_*}$ where for all $1 \le l \le r_*$

$$\bar{H}_{\mu}^{n,l}(\mathbf{c}) := \langle \mathbf{V}^l, \bar{G}_{\mu}^n(\mathbf{c}) \rangle, \qquad (3.13)$$

where $\langle \cdot, \cdot \rangle$ denotes the euclidean scalar product of $\mathbb{R}^{N_h \times (N_s-1)}$. Then, the vector $\bar{c}_{\mu}^n \in \mathbb{R}^{r_*}$ is defined as a solution to the nonlinear problem

$$\bar{H}_{\mu}^{n}(\bar{c}_{\mu}^{n}) = 0, \qquad (3.14)$$

and the numerical approximation $\bar{U}_{\mu}^n$ of $\bar{U}_{\mu}^n$ given by the SP-ROM is then equal to $\bar{U}(\bar{c}_{\mu}^n)$, which is exactly (3.8) with $\bar{Z}_{\mu}^n$ defined by (3.6). A solution to problem (3.14) is computed using a classical Newton algorithm.

The following proposition states that the SP-ROM indeed produces an approximation $\bar{U}_{\mu}^n$ of $\bar{U}_{\mu}^n$ which enables to preserve analogous properties to (P1<sup>h</sup>)-(P2<sup>h</sup>)-(P3<sup>h</sup>) on the level of the reduced-order model.
(P2<sup>red</sup>) volume-filling constraint:

$$
\bar{u}_{\mu,i,K}^{n} > 0 \quad \forall K \in \mathcal{T}, \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T], \qquad (3.9)
$$

and

$$
\sum_{i=1}^{N_s} \bar{u}_{\mu,i,K}^{n} = 1 \quad \forall K \in \mathcal{T}, \quad \forall n \in [0, N_T], \qquad (3.10)
$$

which is the analogous version of (P2<sup>h</sup>) for the SP-ROM.

Let us now explain in details how the vector $\bar{c}_{\mu}^n$ is computed. For any $c := (c^k)_{1\le k\le r_*} \in \mathbb{R}^{r_*}$, let us introduce $\bar{Z}(c) := (\bar{z}_{j,K}(c))_{K\in\mathcal{T}, j\in[1,N_s-1]} \in \mathbb{R}^{N_h \times (N_s-1)}$ the vector defined by

$$
\bar{Z}(c) := \sum_{k=1}^{r_*} c^k \mathbf{V}^k.
$$

Moreover, we define $\bar{U}(c) := (\bar{u}_{i,K}(c))_{K \in \mathcal{T}, i \in [1, N_s]} \in \mathbb{R}^{N_h \times N_s}$ by

$$
\forall i \in [1, N_s], \forall K \in \mathcal{T}, \quad \bar{u}_{i,K}(c) := [\mathcal{Q}(\bar{z}_K(c))]_i, \qquad (3.11)
$$

where $\bar{z}_K(c) = (\bar{z}_{j,K}(c))_{j \in [1, N_s - 1]}$. Let us now define for all $\mu \in \mathcal{P}$, $1 \le n \le N_T$, $1 \le j \le N_s - 1$ and $K \in \mathcal{T}$, the application $\bar{G}_{\mu, i, K}^n : \mathbb{R}^{r_*} \to \mathbb{R}$ defined by

$$
\bar{G}_{\mu,j,K}^{n}(\mathbf{c}) := m_K \frac{\bar{u}_{j,K}(\mathbf{c}) - \bar{u}_{\mu,j,K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K} F_{\mu,j,K\sigma} (\bar{U}(\mathbf{c})), \quad (3.12)
$$

where for all $i \in [1, N_s]$ and $K \in \mathcal{T}$, $\bar{u}_{\mu,i,K}^0(c) := u_{i,K}^0$, and by $\bar{G}_{\mu}^n(c) := (\bar{G}_{\mu,j,K}^n(c))_{K \in \mathcal{T}}$, $j \in [1, N_s - 1]$.
Let us finally denote by $\bar{H}_{\mu}^n : \mathbb{R}^{r_*} \to \mathbb{R}^{r_*}$ the application such that, for all $c \in \mathbb{R}^{r_*}$, $\bar{H}_{\mu}^n = (\bar{H}_{\mu}^{n,l}(c))_{1 \le l \le r_*}$
where for all $1 \le l \le r_*$

$$
\bar{H}_{\mu}^{n,l}(c) := \langle \mathbf{V}^l, \bar{G}_{\mu}^{n}(c) \rangle, \qquad (3.13)
$$

where $\langle \cdot, \cdot \rangle$ denotes the euclidean scalar product of $\mathbb{R}^{N_h \times (N_s-1)}$. Then, the vector $\bar{c}_{\mu}^{n} \in \mathbb{R}^{r_{*}}$ is defined as a
solution **to the nonlinear problem**

$$
\bar{H}_{\mu}^{n}(\bar{c}_{\mu}^{n}) = 0, \qquad (3.14)
$$

and the numerical approximation $\bar{U}_{\mu}^{n}$ of $U_{\mu}^{n}$ given by the SP-ROM is then equal to $\bar{U}(\bar{c}_{\mu}^{n})$, which is exactly
(3.8) with $\bar{Z}_{\mu}^{n}$ defined by (3.6). A solution to problem (3.14) is computed using a classical Newton algorithm.

The following proposition states that the SP-ROM indeed produces an approximation $\bar{U}_{\mu}^n$ of $U_{\mu}^n$ which enables to preserve analogous properties to (P1$^h$)-(P2$^h$)-(P3$^h$) on the level of the reduced-order model.

**Proposition 3.1.** Let $\mu \in \mathcal{P}$ and let $(\bar{U}_{\mu}^n)_{1 \le n \le N_T}$ be the solution of the SP-ROM described above. Then, denoting by $\bar{U}_{\mu} := (\bar{u}_{\mu,i,K}^n)_{K \in \mathcal{T}, 1 \le i \le N_s}$ for all $1 \le n \le N_T$, the following properties hold:

(P1<sup>red</sup>) mass conservation:

$$
\sum_{K \in \mathcal{T}} m_K \bar{u}_{\mu, i, K}^{n} = \sum_{K \in \mathcal{T}} m_K u_{\mu, i, K}^{0} = \int_{\Omega} u_{i}^{0}(x) \, dx \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T]. \quad (3.15)
$$

(P2<sup>red</sup>) volume-filling constraint:

$$
\bar{u}_{\mu, i, K}^{n} > 0 \quad \forall K \in \mathcal{T}, \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T], \qquad (3.16)
$$

and

$$
\sum_{i=1}^{N_s} \bar{u}_{\mu,i,K}^{n} = 1 \quad \forall K \in \mathcal{T}, \quad \forall n \in [0, N_T]. \qquad (3.17)
$$
**Proposition 3.1.** Let $\mu \in \mathcal{P}$ and let $(\overline{U}_{\mu}^n)_{1 \le n \le N_T}$ be the solution of the SP-ROM described above. Then, denoting by $\overline{U}_{\mu}^n := (\overline{u}_{\mu,i,K}^n)_{K \in \mathcal{T}, 1 \le i \le N_s}$ for all $1 \le n \le N_T$, the following properties hold:

(P1<sup>red</sup>) **mass conservation:**

$$ \sum_{K \in \mathcal{T}} m_K \overline{u}_{\mu,i,K}^n = \sum_{K \in \mathcal{T}} m_K u_{\mu,i,K}^0 = \int_{\Omega} u_i^0(x) \, dx \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T]. \quad (3.15) $$

(P2<sup>red</sup>) **volume-filling constraint:**

$$ \overline{u}_{\mu,i,K}^n > 0 \quad \forall K \in \mathcal{T}, \quad \forall i \in [1, N_s], \quad \forall n \in [0, N_T], \quad (3.16) $$

and

$$ \sum_{i=1}^{N_s} \overline{u}_{\mu,i,K}^n = 1 \quad \forall K \in \mathcal{T}, \quad \forall n \in [0, N_T]. \quad (3.17) $$

(P3<sup>red</sup>) **entropy-entropy dissipation inequality:** For all $1 \le n \le N_T$, it holds that

$$ \frac{1}{\Delta t_n} \left( E_{\mathcal{T}}\left(\overline{U}_{\mu}^n\right) - E_{\mathcal{T}}\left(\overline{U}_{\mu}^{n-1}\right) \right) + a_{\mu}^{*} \sum_{\sigma \in \mathcal{E}_h} \sum_{i=1}^{N_s} \overline{u}_{\mu,i,\sigma}^n |D_{K\sigma} \ln(\overline{u}_{\mu,i}^n)|^2 \le 0, \quad \forall n \in [1, N_T]. \quad (3.18) $$

The next section is devoted to the proof of Proposition 3.1.

### 3.3. Proof of Proposition 3.1

Before proving Proposition 3.1 we introduce the following intermediate results.

**Lemma 3.2.** Let $e = (e_K)_{K \in \mathcal{T}}$, $v = (v_K)_{K \in \mathcal{T}} \in \mathbb{R}^{N_h}$. Then,

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} e_K D_{K\sigma} v = - \sum_{\sigma \in \mathcal{E}^{\text{int}}, \sigma=K|L} D_{K\sigma} e D_{K\sigma} v. \quad (3.19) $$

*Proof.* We have

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} e_K D_{K\sigma} v = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} e_K (v_{K\sigma} - v_K). $$

Next, observe that for $\sigma \in \mathcal{E}^{\text{int}}$ such that $\sigma = K|L$, we have,

$$ e_K(v_{K\sigma} - v_K) + e_L(v_{L\sigma} - v_L) = -(e_K - e_L)(v_K - v_L) = -D_{K\sigma}eD_{K\sigma}v. $$

We sum on all internal edges to get the desired result. $\square$

**Lemma 3.3.** For all $U := (u_{i,K})_{K \in \mathcal{T}, i \in [1, N_s]} \in \mathbb{R}^{N_h \times N_s}$ such that $\sum_{i=1}^{N_s} u_{i,K} = 1$ for all $K \in \mathcal{T}$, it holds that for all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}$, we have

$$ \sum_{i=1}^{N_s} F_{\mu,i,K\sigma}(U) = 0. \quad (3.20) $$

*Proof.* Let $U := (u_{i,K})_{K \in \mathcal{T}, i \in [1, N_s]} \in \mathbb{R}^{N_h \times N_s}$ such that $\sum_{i=1}^{N_s} u_{i,K} = 1$ for all $K \in \mathcal{T}$. Then, it holds that

$$ \sum_{i=1}^{N_s} F_{\mu,i,K\sigma}(U) = \sum_{i=1}^{N_s} -a_{\mu}^{*}\tau_{\sigma}D_{K\sigma}\mathbf{u}_i - \tau_{\sigma}\sum_{i=1}^{N_s}\left(\sum_{j=1,j \neq i}^{N_s}(a_{ij} - a_{\mu}^{*})(u_{j,\sigma}D_{K\sigma}\mathbf{u}_i - u_{i,\sigma}D_{K\sigma}\mathbf{u}_j)\right) $$
(P3red) entropy-entropy dissipation inequality: For all $1 \le n \le N_T$, it holds that

$$ \frac{1}{\Delta t_n} (E_T(\bar{u}_\mu^n) - E_T(\bar{u}_\mu^{n-1})) + a_\mu^* \sum_{\sigma \in \mathcal{E}_h} \sum_{i=1}^{N_s} \bar{u}_{\mu,i,\sigma}^n |D_{K\sigma} \ln(\bar{u}_{\mu,i}^n)|^2 \le 0, \quad \forall n \in [1, N_T]. \quad (3.18) $$

The next section is devoted to the proof of Proposition 3.1.

## 3.3 Proof of Proposition 3.1

Before proving Proposition 3.1 we introduce the following intermediate results.

**Lemma 3.2.** Let $e = (e_K)_{K \in \mathcal{T}}$, $v = (v_K)_{K \in \mathcal{T}} \in \mathbb{R}^{N_h}$. Then,

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} e_K D_{K\sigma} v = - \sum_{\sigma \in \mathcal{E}^{\text{int}}, \sigma=K|L} D_{K\sigma} e D_{K\sigma} v. \quad (3.19) $$

*Proof.* We have

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} e_K D_{K\sigma} v = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} e_K (v_{K\sigma} - v_K). $$

Next, observe that for $\sigma \in \mathcal{E}^{\text{int}}$ such that $\sigma = K|L$, we have,

$$ e_K(v_{K\sigma} - v_K) + e_L(v_{L\sigma} - v_L) = -(e_K - e_L)(v_K - v_L) = -D_{K\sigma}eD_{K\sigma}v. $$

We sum on all internal edges to get the desired result.

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

**Lemma 3.3.** For all $U := (u_{i,K})_{K \in \mathcal{T}, i \in [1, N_s]} \in \mathbb{R}^{N_h \times N_s}$ such that $\sum_{i=1}^{N_s} u_{i,K} = 1$ for all $K \in \mathcal{T}$, it holds that for all $K \in \mathcal{T}$ and all $\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}$, we have

$$ \sum_{i=1}^{N_s} F_{\mu,i,K\sigma}(U) = 0. \quad (3.20) $$

*Proof.* Let $U := (u_{i,K})_{K \in \mathcal{T}, i \in [1, N_s]} \in \mathbb{R}^{N_h \times N_s}$ such that $\sum_{i=1}^{N_s} u_{i,K} = 1$ for all $K \in \mathcal{T}$. Then, it holds that

$$ \begin{aligned} \sum_{i=1}^{N_s} F_{\mu,i,K\sigma}(U) &= \sum_{i=1}^{N_s} -a_{\mu}^{*} \tau_{\sigma} D_{K\sigma} u_i - \tau_{\sigma} \sum_{i=1}^{N_s} \left( \sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^{*}) (u_{j,\sigma} D_{K\sigma} u_i - u_{i,\sigma} D_{K\sigma} u_j) \right) \\ &= -a_{\mu}^{*} \tau_{\sigma} D_{K\sigma} \left( \sum_{i=1}^{N_s} u_i \right) - \tau_{\sigma} \left( \sum_{i=1}^{N_s} \sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^{*}) (u_{j,\sigma} D_{K\sigma} u_i - u_{i,\sigma} D_{K\sigma} u_j) \right). \end{aligned} $$

We first have $D_{K\sigma}(\sum_{i=1}^{N_s} u_i) = 0$ since $\sum_{i=1}^{N_s} u_{i,K} = 1$ for all $K \in \mathcal{T}$. Furthermore,

$$ \sum_{i=1}^{N_s} \sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^{*}) (u_{j,\sigma} D_{K\sigma} u_i - u_{i,\sigma} D_{K\sigma} u_j) = 0, $$

since $a_{ij} = a_{ji}$ for all $1 \le i, j \le N_s$. Hence the result.

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

We are now in position to prove Proposition 3.1.

*Proof of Proposition 3.1.* First, we observe that the positivity of the SP-ROM model is easily checked following the construction (3.7). Furthermore, $\forall \mu \in \mathcal{P} \forall K \in \mathcal{T}_h$

$$ \sum_{i=1}^{N_s-1} \bar{u}_{\mu,i,K}^n + \bar{u}_{\mu,N_s,K}^n = \frac{1}{1 + \sum_{j=1}^{N_s-1} e^{\bar{z}_{\mu,j,K}^n}} \left( 1 + \sum_{i=1}^{N_s-1} e^{\bar{z}_{\mu,i,K}^n} \right) $$
$$= -a_{\mu}^{*} \tau_{\sigma} D_{K\sigma} \left( \sum_{i=1}^{N_s} \mathbf{u}_{i} \right) - \tau_{\sigma} \left( \sum_{i=1}^{N_s} \sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^{*}) (u_{j,\sigma} D_{K\sigma} \mathbf{u}_{i} - u_{i,\sigma} D_{K\sigma} \mathbf{u}_{j}) \right).$$

We first have $D_{K\sigma}(\sum_{i=1}^{N_s} \mathbf{u}_i) = 0$ since $\sum_{i=1}^{N_s} u_{i,K} = 1$ for all $K \in \mathcal{T}$. Furthermore,

$$\sum_{i=1}^{N_s} \sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^{*}) (u_{j,\sigma} D_{K\sigma} \mathbf{u}_{i} - u_{i,\sigma} D_{K\sigma} \mathbf{u}_{j}) = 0,$$

since $a_{ij} = a_{ji}$ for all $1 \le i, j \le N_s$. Hence the result.

We are now in position to prove Proposition 3.1.

*Proof of Proposition 3.1.* First, we observe that the positivity of the SP-ROM model is easily checked following the construction (3.7). Furthermore, $\forall \mu \in \mathcal{P} \forall K \in \mathcal{T}_h$

$$\sum_{i=1}^{N_s-1} \bar{u}_{\mu, i, K}^{n} + \bar{u}_{\mu, N_s, K}^{n} = \frac{1}{1 + \sum_{j=1}^{N_s-1} e^{\bar{z}_{\mu, j, K}^{n}}} \left( 1 + \sum_{i=1}^{N_s-1} e^{\bar{z}_{\mu, i, K}^{n}} \right)$$

so that $\sum_{i=1}^{N_s} \bar{u}_{\mu, i, K}^{n} = 1$ which shows that (P2<sup>red</sup>) is checked. We are thus left to prove (P1<sup>red</sup>) and (P3<sup>red</sup>). On the one hand, we have for all $j \in [1, N_s - 1]$, and $1 \le n \le N_T$ since $\bar{U}_{\mu}^{n} = \bar{U}(\bar{c}_{\mu}^{n})$,

$$\langle \mathbf{V}^{r+j}, \bar{G}_{\mu}^{n}(\bar{U}_{\mu}^{n}) \rangle = \sum_{K \in \mathcal{T}} m_K \frac{\bar{u}_{\mu, j, K}^{n} - \bar{u}_{\mu, j, K}^{n-1}}{\Delta t_n} + \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}_{\text{int}}} F_{\mu, j, K\sigma}(\bar{U}_{\mu}^{n}) = 0.$$

Let us prove here that

$$\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}_{\text{int}}} F_{\mu, j, K\sigma}(\bar{U}_{\mu}^{n}) = 0. \qquad (3.21)$$

Indeed, it holds that

$$\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}_{\text{int}}} F_{\mu, j, K\sigma}(\bar{U}_{\mu}^{n}) = \sum_{\sigma \in \mathcal{E}_{\text{int}}, \sigma = K|L} (F_{\mu, j, K\sigma}(\bar{U}_{\mu}^{n}) + F_{\mu, j, L\sigma}(\bar{U}_{\mu}^{n})).$$

Furthermore, for all $\sigma = K|L \in \mathcal{E}_{\text{int}}$, it holds that $F_{\mu, j, K\sigma}(\bar{U}_{\mu}^{n}) + F_{\mu, j, L\sigma}(\bar{U}_{\mu}^{n}) = 0$ from Remark 2.5. As a consequence, we obtain that

$$\sum_{K \in \mathcal{T}} m_K \bar{u}_{\mu, j, K}^{n} = \sum_{K \in \mathcal{T}} m_K \bar{u}_{\mu, j, K}^{n-1}. \qquad (3.22)$$

Reasoning by induction, we then obtain that

$$\sum_{K \in \mathcal{T}} m_K \bar{u}_{\mu, j, K}^{n} = \sum_{K \in \mathcal{T}} m_K \bar{u}_{\mu, j, K}^{0} = \sum_{K \in \mathcal{T}} m_K u_{\mu, j, K}^{0} = \int_{\Omega} u_{j}^{0}(x) dx$$

which yields (P1<sup>red</sup>) since the relationship naturally holds for $j = N_s$ too.

Let us now prove (P3<sup>red</sup>). It holds that

$$E_{\mathcal{T}}(\bar{U}_{\mu}^{n}) - E_{\mathcal{T}}(\bar{U}_{\mu}^{n-1}) = \sum_{K \in \mathcal{T}} \sum_{i=1}^{N_s} m_K (\bar{u}_{\mu, i, K}^{n} \ln(\bar{u}_{\mu, i, K}^{n}) - \bar{u}_{\mu, i, K}^{n-1} \ln(\bar{u}_{\mu, i, K}^{n-1})). \qquad (3.23)$$
so that $\sum_{i=1}^{N_s} \bar{u}_{\mu,i,K}^n = 1$ which shows that (P2<sup>red</sup>) is checked. We are thus left to prove (P1<sup>red</sup>) and (P3<sup>red</sup>). On the one hand, we have for all $j \in [1, N_s - 1]$, and $1 \le n \le N_T$ since $\bar{U}_{\mu}^n = \bar{U}(\bar{c}_{\mu}^n)$,

$$ \langle \mathbf{V}^{r+j}, \bar{G}_{\mu}^{n}(\bar{U}_{\mu}^{n}) \rangle = \sum_{K \in \mathcal{T}} m_K \frac{\bar{u}_{\mu,j,K}^{n} - \bar{u}_{\mu,j,K}^{n-1}}{\Delta t_n} + \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu,j,K\sigma}(\bar{U}_{\mu}^{n}) = 0. $$

Let us prove here that

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu,j,K\sigma}(\bar{U}_{\mu}^n) = 0. \qquad (3.21) $$

Indeed, it holds that

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu,j,K\sigma}(\bar{U}_{\mu}^n) = \sum_{\sigma \in \mathcal{E}^{\text{int}}, \sigma=K|L} \left( F_{\mu,j,K\sigma}(\bar{U}_{\mu}^n) + F_{\mu,j,L\sigma}(\bar{U}_{\mu}^n) \right). $$

Furthermore, for all $\sigma = K|L \in \mathcal{E}^{\text{int}}$, it holds that $F_{\mu,j,K\sigma}(\bar{U}_{\mu}^n) + F_{\mu,j,L\sigma}(\bar{U}_{\mu}^n) = 0$ from Remark 2.5. As a consequence, we obtain that

$$ \sum_{K \in \mathcal{T}} m_K \bar{u}_{\mu,j,K}^n = \sum_{K \in \mathcal{T}} m_K \bar{u}_{\mu,j,K}^{n-1}. \qquad (3.22) $$

Reasoning by induction, we then obtain that

$$ \sum_{K \in \mathcal{T}} m_K \bar{u}_{\mu,j,K}^n = \sum_{K \in \mathcal{T}} m_K \bar{u}_{\mu,j,K}^0 = \sum_{K \in \mathcal{T}} m_K u_{\mu,j,K}^0 = \int_{\Omega} u_j^0(x) \, dx $$

which yields (P1<sup>red</sup>) since the relationship naturally holds for $j = N_s$ too.

Let us now prove (P3<sup>red</sup>). It holds that

$$ E_{\mathcal{T}}(\bar{U}_{\mu}^{n}) - E_{\mathcal{T}}(\bar{U}_{\mu}^{n-1}) = \sum_{K \in \mathcal{T}} \sum_{i=1}^{N_s} m_K \left( \bar{u}_{\mu,i,K}^{n} \ln(\bar{u}_{\mu,i,K}^{n}) - \bar{u}_{\mu,i,K}^{n-1} \ln(\bar{u}_{\mu,i,K}^{n-1}) \right). \qquad (3.23) $$

Besides, the coefficients $\bar{c}_{\mu}^n$ are solution to $\bar{H}_{\mu}^n(\bar{c}_{\mu}^n) = 0$. As a consequence of (3.13), we have for all $1 \le l \le r_*$, $\langle \mathbf{V}^l, \bar{G}_{\mu}^n(\bar{U}_{\mu}^n) \rangle = 0$ and then $(\sum_{l=1}^{r_*} \bar{c}_{\mu}^{l,n} \mathbf{V}^l)^T \bar{G}_{\mu}^n(\bar{U}_{\mu}^n) = 0$, which equivalently reads as

$$ \langle \bar{Z}_{\mu}^{n}, \bar{G}_{\mu}^{n}(\bar{U}_{\mu}^{n}) \rangle = \sum_{i=1}^{N_s-1} \sum_{K \in \mathcal{T}} \bar{z}_{\mu,i,K}^{n} G_{i,K}^{n}(\bar{U}_{\mu}^{n}) = 0. $$

We also introduce the notation $\bar{\rho}_{\mu,K}^n := 1 + \sum_{j=1}^{N_s-1} \exp(\bar{z}_{\mu,j,K}^n)$. It follows immediately that

$$ \ln(\bar{\rho}_{\mu,K}^n) + \ln(\bar{u}_{\mu,j,K}^n) = \bar{z}_{\mu,j,K}^n \quad \forall j \in [1, N_s - 1], \\ \ln(\bar{\rho}_{\mu,K}^n) + \ln(\bar{u}_{\mu,N_s,K}^n) = 0. \qquad (3.24) $$

We then get $A + B = 0$, with

$$ A := \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu,i,K}^n) \left( m_K \frac{\bar{u}_{\mu,i,K}^n - \bar{u}_{\mu,i,K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu,i,K\sigma}(\bar{U}_{\mu}^n) \right), \qquad (3.25) $$

and

$$ B := \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{\rho}_{\mu,K}^n) \left( m_K \frac{\bar{u}_{\mu,i,K}^n - \bar{u}_{\mu,i,K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu,i,K\sigma}(\bar{U}_{\mu}^n) \right). $$
Besides, the coefficients $\bar{c}_{\mu}^n$ are solution to $\overline{H}_{\mu}^n(\bar{c}_{\mu}^n) = 0$. As a consequence of (3.13), we have for all $1 \le l \le r_*$, $\langle \mathbf{V}^l, \overline{G}_{\mu}^n(\overline{U}_{\mu}^n) \rangle = 0$ and then $(\sum_{l=1}^{r_*} \bar{c}_{\mu}^{l,n} \mathbf{V}^l)^T \overline{G}_{\mu}^n(\overline{U}_{\mu}^n) = 0$, which equivalently reads as

$$ \langle \overline{Z}_{\mu}^{n}, \overline{G}_{\mu}^{n}(\overline{U}_{\mu}^{n}) \rangle = \sum_{i=1}^{N_s-1} \sum_{K \in \mathcal{T}} \bar{z}_{\mu, i, K}^{n} G_{i, K}^{n}(\overline{U}_{\mu}^{n}) = 0. $$

We also introduce the notation $\bar{\rho}_{\mu,K}^n := 1 + \sum_{j=1}^{N_s-1} \exp(\bar{z}_{\mu,j,K}^n)$. It follows immediately that

$$ \begin{aligned} \ln(\bar{\rho}_{\mu,K}^n) + \ln(\bar{u}_{\mu,j,K}^n) &= \bar{z}_{\mu,j,K}^n \quad \forall j \in [1, N_s - 1], \\ \ln(\bar{\rho}_{\mu,K}^n) + \ln(\bar{u}_{\mu,N_s,K}^n) &= 0. \end{aligned} \quad (3.24) $$

We then get $A + B = 0$, with

$$ A := \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu, i, K}^n) \left( m_K \frac{\bar{u}_{\mu, i, K}^n - \bar{u}_{\mu, i, K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu, i, K \sigma} (\bar{U}_{\mu}^n) \right), \quad (3.25) $$

and

$$ B := \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{\rho}_{\mu,K}^n) \left( m_K \frac{\bar{u}_{\mu, i, K}^n - \bar{u}_{\mu, i, K}^{n-1}}{\Delta t_n} + \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu, i, K \sigma} (\bar{U}_{\mu}^n) \right). $$

Besides, we observe that, using (3.22) and Lemma 3.3, we have $B = 0$ so that $A = 0$. Therefore, it holds that

$$ \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu, i, K}^n) \left( m_K \frac{\bar{u}_{\mu, i, K}^n - \bar{u}_{\mu, i, K}^{n-1}}{\Delta t_n} \right) + \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu, i, K}^n) \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu, i, K \sigma} (\bar{U}_{\mu}^n) = 0. $$

Furthermore, we have

$$ F_{\mu,i,K\sigma}(\bar{U}_{\mu}^n) := -a_{\mu}^* \tau_{\sigma} \bar{u}_{i,\sigma}^n D_{K\sigma}(\ln(\bar{u}_{\mu,i}^n)) \\ - \tau_{\sigma} \left( \sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^*) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n (D_{K\sigma}(\ln(\bar{u}_{\mu,i}^n)) - D_{K\sigma}(\ln(\bar{u}_{\mu,j}^n))) \right). $$

Thus we get,

$$ A = \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu, i, K}^n) \left( m_K \frac{\bar{u}_{\mu, i, K}^n - \bar{u}_{\mu, i, K}^{n-1}}{\Delta t_n} \right) + A_1 + A_2 \quad (3.26) $$

with

$$ A_1 := \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu, i, K}^n) \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} (-a_{\mu}^* \tau_{\sigma} \bar{u}_{\mu, i, \sigma}^n D_{K\sigma}(\ln(\bar{u}_{\mu, i}^n))) $$

and

$$ A_2 := \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu, i, K}^n) \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} \tau_{\sigma} \left( \sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^*) \bar{u}_{\mu, j, \sigma}^n \bar{u}_{\mu, i, \sigma}^n (D_{K\sigma}(\ln(\bar{u}_{\mu, i}^n)) - D_{K\sigma}(\ln(\bar{u}_{\mu, j}^n))) \right). $$

Besides,

$$ A_2 = \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} \tau_{\sigma} \left( \sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^*) \bar{u}_{\mu, j, \sigma}^n \bar{u}_{\mu, i, \sigma}^n \ln(\bar{u}_{\mu, i, K}^n) (D_{K\sigma}(\ln(\bar{u}_{\mu, i}^n)) - \ln(\bar{u}_{\mu, j}^n)) \right). \quad (3.27) $$
Besides, we observe that, using (3.22) and Lemma 3.3, we have $B = 0$ so that $A = 0$. Therefore, it holds that

$$ \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu,i,K}^n) \left( m_K \frac{\bar{u}_{\mu,i,K}^n - \bar{u}_{\mu,i,K}^{n-1}}{\Delta t_n} \right) + \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu,i,K}^n) \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} F_{\mu,i,K\sigma}(\bar{U}_{\mu}^n) = 0. $$

Furthermore, we have

$$ F_{\mu,i,K\sigma}(\bar{U}_{\mu}^n) := -a_{\mu}^* \tau_{\sigma} \bar{u}_{i,\sigma}^n D_{K\sigma} \left( \ln(\bar{u}_{\mu,i}^n) \right) - \tau_{\sigma} \left( \sum_{j=1, j \neq i}^{N_s} \left( a_{ij} - a_{\mu}^* \right) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n \left( D_{K\sigma} \left( \ln(\bar{u}_{\mu,i}^n) \right) - D_{K\sigma} \left( \ln(\bar{u}_{\mu,j}^n) \right) \right) \right). $$

Thus we get,

$$ A = \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu,i,K}^n) \left( m_K \frac{\bar{u}_{\mu,i,K}^n - \bar{u}_{\mu,i,K}^{n-1}}{\Delta t_n} \right) + A_1 + A_2 \quad (3.26) $$

with

$$ A_1 := \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu,i,K}^n) \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} (-a_{\mu}^* \tau_{\sigma} \bar{u}_{\mu,i,\sigma}^n D_{K\sigma} \left( \ln(\bar{u}_{\mu,i}^n) \right)) $$

and

$$ A_2 := \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \ln(\bar{u}_{\mu,i,K}^n) \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} \tau_{\sigma} \left( \sum_{j=1, j \neq i}^{N_s} \left( a_{ij} - a_{\mu}^* \right) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n \left( D_{K\sigma} \left( \ln(\bar{u}_{\mu,i}^n) \right) - D_{K\sigma} \left( \ln(\bar{u}_{\mu,j}^n) \right) \right) \right). $$

Besides,

$$ A_2 = \sum_{i=1}^{N_s} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K \cap \mathcal{E}^{\text{int}}} \tau_{\sigma} \left( \sum_{j=1, j \neq i}^{N_s} \left( a_{ij} - a_{\mu}^* \right) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n \ln(\bar{u}_{\mu,i}^n) D_{K\sigma} \left( \ln(\bar{u}_{\mu,i}^n) - \ln(\bar{u}_{\mu,j}^n) \right) \right). \quad (3.27) $$

Using Lemma 3.2, we get

$$ A_2 = - \sum_{i=1}^{N_s} \sum_{\sigma \in \mathcal{E}^{\text{int}}} \tau_{\sigma} \left( \sum_{j=1, j \neq i}^{N_s} \left( a_{ij} - a_{\mu}^* \right) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n D_{K\sigma} \left( \ln(\bar{u}_{\mu,i}^n) \right) D_{K\sigma} \left( \ln(\bar{u}_{\mu,j}^n) - \ln(\bar{u}_{\mu,i}^n) \right) \right), $$

and

$$ A_1 = \sum_{i=1}^{N_s} \sum_{\sigma \in \mathcal{E}^{\text{int}}} a_{\mu}^* \tau_{\sigma} \bar{u}_{\mu,i,\sigma}^n D_{K\sigma} \left( \ln(\bar{u}_{\mu,i}^n) \right) D_{K\sigma} \left( \ln(\bar{u}_{\mu,j}^n) \right) = \sum_{i=1}^{N_s} \sum_{\sigma \in \mathcal{E}^{\text{int}}} a_{\mu}^* \tau_{\sigma} \bar{u}_{\mu,i,\sigma}^n \left( D_{K\sigma} \ln(\bar{u}_{\mu,i}^n) \right)^2. $$

Next, we develop the term $D_{K\sigma}(\ln(\bar{u}_{\mu,i}^n)) D_{K\sigma}(\ln(\bar{u}_{\mu,j}^n) - \ln(\bar{u}_{\mu,i}^n))$ to get

$$ \begin{aligned} A_2 &= - \sum_{\sigma \in \mathcal{E}^{\text{int}}} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^*) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n (D_{K\sigma}(\ln(\bar{u}_{\mu,i}^n))^2 \\ &\quad + \sum_{\sigma \in \mathcal{E}^{\text{int}}} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^*) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n D_{K\sigma}(\ln(\bar{u}_{\mu,i}^n)) D_{K\sigma}(\ln(\bar{u}_{\mu,j}^n)) \\ &= - \sum_{\sigma \in \mathcal{E}^{\text{int}}} \frac{1}{2} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^*) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n (D_{K\sigma}(\ln(\bar{u}_{\mu,i}^n))^2 \\ &\quad + \sum_{\sigma \in \mathcal{E}^{\text{int}}} \frac{1}{2} \times 2 \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^*) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n D_{K\sigma}(\ln(\bar{u}_{\mu,i}^n)) D_{K\sigma}(\ln(\bar{u}_{\mu,j}^n)) \\ &\quad - \sum_{\sigma \in \mathcal{E}^{\text{int}}} \frac{1}{2} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{j,i} - a_{\mu}^*) \bar{u}_{\mu,j,\sigma}^n \bar{u}_{\mu,i,\sigma}^n (D_{K\sigma} \ln(\bar{u}_{\mu,j}^n))^2. \end{aligned} $$
Using Lemma 3.2, we get

$$A_2 = - \sum_{i=1}^{N_s} \sum_{\sigma \in \mathcal{E}^{\text{int}}} \tau_{\sigma} \left( \sum_{j=1, j \neq i}^{N_s} (a_{ij} - a_{\mu}^{*}) \bar{u}_{\mu,j,\sigma}^{n} \bar{u}_{\mu,i,\sigma}^{n} D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n})) D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n}) - \ln(\bar{u}_{\mu,j}^{n})) \right),$$

and

$$A_1 = \sum_{i=1}^{N_s} \sum_{\sigma \in \mathcal{E}^{\text{int}}} a_{\mu}^{*} \tau_{\sigma} \bar{u}_{\mu,i,\sigma}^{n} D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n})) D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n})) = \sum_{i=1}^{N_s} \sum_{\sigma \in \mathcal{E}^{\text{int}}} a_{\mu}^{*} \tau_{\sigma} \bar{u}_{\mu,i,\sigma}^{n} (D_{K\sigma} \ln(\bar{u}_{\mu,i}^{n}))^2.$$

Next, we develop the term $D_{K\sigma}(\ln(\bar{u}_{\mu,i}^n))D_{K\sigma}(\ln(\bar{u}_{\mu,i}^n) - \ln(\bar{u}_{\mu,j}^n))$ to get

$$\begin{aligned}
A_2 &= - \sum_{\sigma \in \mathcal{E}^{\text{int}}} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^{*}) \bar{u}_{\mu,j,\sigma}^{n} \bar{u}_{\mu,i,\sigma}^{n} (D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n}))^2) \\
&\quad + \sum_{\sigma \in \mathcal{E}^{\text{int}}} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^{*}) \bar{u}_{\mu,j,\sigma}^{n} \bar{u}_{\mu,i,\sigma}^{n} D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n})) D_{K\sigma}(\ln(\bar{u}_{\mu,j}^{n})) \\
&= - \sum_{\sigma \in \mathcal{E}^{\text{int}}} \frac{1}{2} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^{*}) \bar{u}_{\mu,j,\sigma}^{n} \bar{u}_{\mu,i,\sigma}^{n} (D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n}))^2) \\
&\quad + \sum_{\sigma \in \mathcal{E}^{\text{int}}} \frac{1}{2} \times 2 \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^{*}) \bar{u}_{\mu,j,\sigma}^{n} \bar{u}_{\mu,i,\sigma}^{n} D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n})) D_{K\sigma}(\ln(\bar{u}_{\mu,j}^{n})) \\
&\quad - \sum_{\sigma \in \mathcal{E}^{\text{int}}} \frac{1}{2} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{j,i} - a_{\mu}^{*}) \bar{u}_{\mu,j,\sigma}^{n} \bar{u}_{\mu,i,\sigma}^{n} (D_{K\sigma} \ln(\bar{u}_{\mu,j}^{n}))^2.
\end{aligned}$$

We thus obtain

$$A_2 = - \sum_{\sigma \in \mathcal{E}^{\text{int}}} \frac{1}{2} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^{*}) \bar{u}_{\mu,j,\sigma}^{n} \bar{u}_{\mu,i,\sigma}^{n} (D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n})) - D_{K\sigma}(\ln(\bar{u}_{\mu,j}^{n})))^2.$$

Now, let us consider again the equality (3.23). It follows from the convexity of the function $h: \mathbb{R}_+^* \ni x \mapsto x \ln(x)$ that

$$\bar{u}_{\mu,i,K}^{n} - \bar{u}_{\mu,i,K}^{n-1} + \bar{u}_{\mu,i,K}^{n} \ln(\bar{u}_{\mu,i,K}^{n}) - \bar{u}_{\mu,i,K}^{n-1} \ln(\bar{u}_{\mu,i,K}^{n}) \geq \bar{u}_{\mu,i,K}^{n} \ln(\bar{u}_{\mu,i,K}^{n}) - \bar{u}_{\mu,i,K}^{n-1} \ln(\bar{u}_{\mu,i,K}^{n-1}). \quad (3.28)$$

Thus, employing (3.28) and the mass conservation property (3.22), we obtain that

$$E_{T}(\bar{U}_{\mu}^{n}) - E_{T}(\bar{U}_{\mu}^{n-1}) \leq \Delta t_{n} \sum_{K \in \mathcal{T}} \sum_{i=1}^{N_s} m_{K} \ln(\bar{u}_{\mu,i,K}^{n}) \left( \frac{\bar{u}_{\mu,i,K}^{n} - \bar{u}_{\mu,i,K}^{n-1}}{\Delta t} \right) = \Delta t(A_2 - A_1).$$

Using (3.26) yields

$$E_T(\bar{U}_\mu^n) - E_T(\bar{U}_\mu^{n-1}) + \Delta t_n a_\mu^* \sum_{\sigma \in \mathcal{E}^{\text{int}}} \sum_{i=1}^{N_s} \tau_\sigma \bar{u}_{\mu,i,\sigma}^n (D_{K\sigma} \ln(\bar{u}_{\mu,i}^n))^2 \le \Delta t_n A_2 \le 0.$$

Hence (P3$^{\text{red}}$) and the desired result. □
We thus obtain

$$A_2 = - \sum_{\sigma \in \mathcal{E}^{\text{int}}} \frac{1}{2} \sum_{1 \le i \ne j \le N_s} \tau_{\sigma} (a_{ij} - a_{\mu}^{*}) \bar{u}_{\mu,j,\sigma}^{n} \bar{u}_{\mu,i,\sigma}^{n} (D_{K\sigma}(\ln(\bar{u}_{\mu,i}^{n})) - D_{K\sigma}(\ln(\bar{u}_{\mu,j}^{n})))^2.$$

Now, let us consider again the equality (3.23). It follows from the convexity of the function $h: \mathbb{R}_+^* \ni x \mapsto x \ln(x)$ that

$$\bar{u}_{\mu,i,K}^{n} - \bar{u}_{\mu,i,K}^{n-1} + \bar{u}_{\mu,i,K}^{n} \ln(\bar{u}_{\mu,i,K}^{n}) - \bar{u}_{\mu,i,K}^{n-1} \ln(\bar{u}_{\mu,i,K}^{n}) \geq \bar{u}_{\mu,i,K}^{n} \ln(\bar{u}_{\mu,i,K}^{n-1}) - \bar{u}_{\mu,i,K}^{n-1} \ln(\bar{u}_{\mu,i,K}^{n-1}). \quad (3.28)$$

Thus, employing (3.28) and the mass conservation property (3.22), we obtain that

$$E_{\tau}(\bar{U}_{\mu}^{n}) - E_{\tau}(\bar{U}_{\mu}^{n-1}) \leq \Delta t_n \sum_{K \in \mathcal{T}} \sum_{i=1}^{N_s} m_K \ln(\bar{u}_{\mu,i,K}^{n}) \left( \frac{\bar{u}_{\mu,i,K}^{n} - \bar{u}_{\mu,i,K}^{n-1}}{\Delta t} \right) = \Delta t(A_2 - A_1).$$

Using (3.26) yields

$$E_{\tau}(\bar{U}_{\mu}^{n}) - E_{\tau}(\bar{U}_{\mu}^{n-1}) + \Delta t_n a_{\mu}^{*} \sum_{\sigma \in \mathcal{E}^{\text{int}}} \sum_{i=1}^{N_s} \tau_{\sigma} \bar{u}_{\mu,i,\sigma}^{n} (D_{K\sigma} \ln(\bar{u}_{\mu,i}^{n}))^2 \leq \Delta t_n A_2 \leq 0.$$

Hence (P3red) and the desired result. □

# 4 Numerical experiments

The aim of this section is to give a numerical illustration of our theoretical developments on simple one-dimensional test cases. In particular, the purpose of these numerical illustrations here is to assess the validity of Proposition 3.1. To get a fully efficient reduced-order model, it would be necessary to consider in addition adequate hyper-reduction techniques such as the Empirical Interpolation Method for instance. However, developing a structure-preserving hyperreduction approach does not seem obvious, and we leave this issue for future work. This is the reason why we consider here only simplistic one-dimensional test cases and do not provide estimates in terms of gains in computational times.

All the tests have been performed using a Python code and can be found at url https://jdabaghi.github.io/. We consider two different test cases with respectively $N_s = 3$ and $N_s = 4$ that are reported respectively in Section 4.1 and Section 4.2.

We consider a one-dimensional domain $\Omega = (0, 1)$. We consider a uniform spatial discretization grid of step $\Delta x = 1, 4 \times 10^{-2}$ so that the total number of cells is equal to 70. We use the Newton solver described in Section 2.2.3 with the following stopping criterion: $\sup |U_{\mu}^{n,k} - U_{\mu}^{n,k-1}| < \varepsilon_{\text{lin}}$ with $\varepsilon_{\text{lin}} = 10^{-12}$, where $U_{\mu}^{n,k} \in \mathbb{R}^{N_h \times N_s}$ is the approximate solution provided by the Newton solver at a given step $k \geq 1$. In both test cases, we compare the POD reduced-order model described in Section 3.1 and the SP reduced-order model we propose in Section 3.2. In particular, we numerically illustrate the fact that the SP reduced model satisfies the structural properties of the discrete solution listed in Lemma 3.1. Furthermore, the final simulation time is $T := 0.5$ and we use for the high-fidelity resolution a constant time step $\Delta t = \Delta t_n := 2, 5 \times 10^{-4}$ so that $N_T = 2000$. Therefore, for each parameter $\mu$ we need to solve $N_T$ nonlinear systems yielding a huge computational cost and thus justifying the importance of using a reduced-order-model technique. The value of the parameter $a_{\mu}^{*} > 0$ raised in (2.12) is chosen as follows

$$a_{\mu}^{*} := \min \left\{ \max_{i \neq j} a_{ij}, \max \left\{ \min_{i \neq j} a_{ij}, \frac{1}{2} \frac{\Delta x^2}{\Delta t_n} \right\} \right\}, \quad (4.1)$$

as suggested in [10]. We refer to [10] for a complete discussion.
## 4. NUMERICAL EXPERIMENTS

The aim of this section is to give a numerical illustration of our theoretical developments on simple one-dimensional test cases. In particular, the purpose of these numerical illustrations here is to assess the validity of Proposition 3.1. To get a fully efficient reduced-order model, it would be necessary to consider in addition adequate hyper-reduction techniques such as the Empirical Interpolation Method for instance. However, developing a structure-preserving hyperreduction approach does not seem obvious, and we leave this issue for future work. This is the reason why we consider here only simplistic one-dimensional test cases and do not provide estimates in terms of gains in computational times.

All the tests have been performed using a Python code and can be found at url https://jdabaghi.github.io/. We consider two different test cases with respectively $N_s = 3$ and $N_s = 4$ that are reported respectively in Sections 4.1 and 4.2.

We consider a one-dimensional domain $\Omega = (0, 1)$. We consider a uniform spatial discretization grid of step $\Delta x = 1, 4 \times 10^{-2}$ so that the total number of cells is equal to 70. We use the Newton solver described in Section 2.2.3 with the following stopping criterion: $\sup |U_{\mu}^{n,k} - U_{\mu}^{n,k-1}| < \varepsilon_{\text{lin}}$ with $\varepsilon_{\text{lin}} = 10^{-12}$, where $U_{\mu}^{n,k} \in \mathbb{R}^{N_h \times N_s}$ is the approximate solution provided by the Newton solver at a given step $k \ge 1$. In both test cases, we compare the POD reduced-order model described in Section 3.1 and the SP reduced-order model we propose in Section 3.2. In particular, we numerically illustrate the fact that the SP reduced model satisfies the structural properties of the discrete solution listed in Lemma 3.1. Furthermore, the final simulation time is $T := 0.5$ and we use for the high-fidelity resolution a constant time step $\Delta t = \Delta t_n := 2, 5 \times 10^{-4}$ so that $N_T = 2000$. Therefore, for each parameter $\mu$ we need to solve $N_T$ nonlinear systems yielding a huge computational cost and thus justifying the importance of using a reduced-order-model technique. The value of the parameter $a_{\mu}^{*} > 0$ raised in (2.12) is chosen as follows

$$a_{\mu}^{*} := \min \left\{ \max_{i \neq j} a_{ij}, \max \left\{ \min_{i \neq j} a_{ij}, \frac{1}{2} \frac{\Delta x^2}{\Delta t_n} \right\} \right\}, \quad (4.1)$$

as suggested in [10]. We refer to [10] for a complete discussion.

### 4.1. Test case 1: three species

In this first test case, the initial guess is chosen as

$$u_1^0(x) := \begin{cases} 1 - 2\delta & \text{if } x \in [\frac{3}{8}, \frac{5}{8}] \\ \delta & \text{else,} \end{cases} \qquad u_2^0(x) := \begin{cases} 1 - 2\delta & \text{if } x \in [\frac{1}{8}, \frac{3}{8}] \cap [\frac{5}{8}, \frac{7}{8}] \\ \delta & \text{else} \end{cases} \quad (4.2)$$

and

$$u_3^0(x) := \begin{cases} 1 - 2\delta & \text{if } x \in [0, \frac{1}{8}] \cap [\frac{7}{8}, 1] \\ \delta & \text{else.} \end{cases} \quad (4.3)$$

Note that the parameter $\delta$ is chosen equal to $\delta = 0.1$. Furthermore, the initial solution $u^0$ satisfies the volume filling constraint property [(P3)] of Section 1.2. Within the offline stage, we compute 20 snapshots of solutions, i.e. we compute the collection $(U_{\mu}^n)_{1 \le n \le N_T}$ for all $\mu$ belonging to a subset $\mathcal{P}_{\text{train}} \subset \mathcal{P}$ so that $\text{Card}(\mathcal{P}_{\text{train}}) = 20$. For the sake of clarity, the elements of $\mathcal{P}_{\text{train}}$ which are in fact $N_s \times N_s$ matrices, are determined by selecting random numbers in the set $[0, 1]$. More precisely, as the matrices of $\mathcal{P}_{\text{train}}$ are symmetric with diagonal off coefficients $\frac{N_s \times (N_s - 1)}{2}$ real values belonging to the interval $[0, 1]$ are sampled. We denote in the following by $\mu^0, \dots, \mu^{19}$ the elements of $\mathcal{P}_{\text{train}}$. An illustration of the initial conditions $u_1^0, u_2^0$, and $u_3^0$ is provided in Figure 2.

#### 4.1.1. The high-fidelity problem

In Figure 3, we represent the behavior of the numerical solution at several time steps for one selected parameter $\mu^0 \in \mathcal{P}_{\text{train}}$ and when the Newton solver has converged. In fact, it corresponds to high-fidelity
<figure><img src="image_6.png" /><figcaption>Figure 2: Initial condition for test case 1.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3: $u_1(t, x)$, $u_2(t, x)$, and $u_3(t, x)$ at different values of the time $t$. $t = 0.005$ (left), $t = 0.05$ (middle), and $t = 0.2$ (right).</figcaption></figure>

## 4.1 Test case 1: three species

In this first test case, the initial guess is chosen as

$$ u_1^0(x) := \begin{cases} 1 - 2\delta & \text{if } x \in [\frac{3}{8}, \frac{5}{8}] \\ \delta & \text{else,} \end{cases} \quad u_2^0(x) := \begin{cases} 1 - 2\delta & \text{if } x \in [\frac{1}{8}, \frac{3}{8}] \cap [\frac{5}{8}, \frac{7}{8}] \\ \delta & \text{else} \end{cases} \quad (4.2) $$

and

$$ u_3^0(x) := \begin{cases} 1 - 2\delta & \text{if } x \in [0, \frac{1}{8}] \cap [\frac{7}{8}, 1] \\ \delta & \text{else.} \end{cases} \quad (4.3) $$

Note that the parameter $\delta$ is chosen equal to $\delta = 0.1$. Furthermore, the initial solution $u^0$ satisfies the volume filling constraint property [(P3)] of Section 1.2. Within the offline stage, we compute 20 snapshots of solutions, i.e. we compute the collection $(U_{\mu}^n)_{1 \le n \le N_T}$ for all $\mu$ belonging to a subset $\mathcal{P}_{\text{train}} \subset \mathcal{P}$ so that $\text{Card}(\mathcal{P}_{\text{train}}) = 20$. For the sake of clarity, the elements of $\mathcal{P}_{\text{train}}$ which are in fact $N_s \times N_s$ matrices, are determined by selecting random numbers in the set $[0, 1]$. More precisely, as the matrices of $\mathcal{P}_{\text{train}}$ are symmetric with diagonal off coefficients $\frac{N_s \times (N_s - 1)}{2}$ real values belonging to the interval $[0, 1]$ are sampled. We denote in the following by $\mu^0, \dots, \mu^{19}$ the elements of $\mathcal{P}_{\text{train}}$. An illustration of the initial conditions $u_1^0, u_2^0$, and $u_3^0$ is provided in Figure 2.

### 4.1.1 The high-fidelity problem

In Figure 3, we represent the behavior of the numerical solution at several time steps for one selected parameter $\mu^0 \in \mathcal{P}_{\text{train}}$ and when the Newton solver has converged. In fact, it corresponds to high-fidelity numerical solutions. In this case, the parameter $\mu^0 = (a_{ij}^0)_{1 \le i,j \le N_s}$ is given by

$$ a_{12}^0 = 0.75, \quad a_{13}^0 = 0.73, \quad a_{23}^0 = 0.84. \quad (4.4) $$

We observe that the local volumic fractions evolve over time to reach constant profiles (around 0.44 for $u_2$, and around 0.26 for $u_1$ and $u_3$) in the long time limit. This behavior is typical of the solutions of the particular cross-diffusion system we consider here.
<figure><img src="image_2.png" /><figcaption>FIGURE 2. Initial condition for test case 1.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 3. $u_1(t, x)$, $u_2(t, x)$, and $u_3(t, x)$ at different values of the time $t$. $t = 0.005$ (left), $t = 0.05$ (middle), and $t = 0.2$ (right).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 4. Properties of the fine numerical solution. Left: non-negativity of the solution, middle: preservation of the volume filling constraint, right: conservation of mass.</figcaption></figure>

numerical solutions. In this case, the parameter $\mu^0 = (a_{ij}^0)_{1 \le i,j \le N_s}$ is given by

$$a_{12}^0 = 0.75, \quad a_{13}^0 = 0.73, \quad a_{23}^0 = 0.84. \qquad (4.4)$$

We observe that the local volumic fractions evolve over time to reach constant profiles (around 0.44 for $u_2$, and around 0.26 for $u_1$ and $u_3$) in the long time limit. This behavior is typical of the solutions of the particular cross-diffusion system we consider here.

In Figure 4 we show that the numerical solution obtained by the finite volume resolution (see Sect. 2.2) preserves the structural properties listed in Section 2.2.2. We have simulated the cross-diffusion model (2.1) for all $\mu \in \mathcal{P}_{\text{train}}$. In the left figure, we have represented for each time step $n \in [1, N_T]$, the following quantity:

$$\mathcal{P}_U(t_n) := \inf_{\mu \in \mathcal{P}_{\text{train}}[1, N_s]} \inf_{K \in \mathcal{T}} \inf_{\mu, i, K} u_{\mu, i, K}^n. \qquad (4.5)$$

We observe that the functional $\mathcal{P}_U$ reaches its minimum around $0.065 > 0$ so the numerical solution is always positive. In the middle figure is displayed the volume filling constraint error $|1 - S_U(t_n)|$ obtained at each time
<figure><img src="image_8.png" /><figcaption>Figure 4: Properties of the fine numerical solution. Left: non-negativity of the solution, middle: preservation of the volume filling constraint, right: conservation of mass.</figcaption></figure>

In Figure 4 we show that the numerical solution obtained by the finite volume resolution (see Section 2.2) preserves the structural properties listed in Section 2.2.2. We have simulated the cross-diffusion model (2.1) for all $\mu \in \mathcal{P}_{\text{train}}$. In the left figure, we have represented for each time step $n \in [1, N_T]$, the following quantity:

$$ \mathcal{P}_U(t_n) := \inf_{\mu \in \mathcal{P}_{\text{train}}} \inf_{K \in \mathcal{T}} \inf_{n \in \mathcal{N}_s} u_{\mu, i, K}^n. \quad (4.5) $$

We observe that the functional $\mathcal{P}_U$ reaches its minimum around $0.065 > 0$ so the numerical solution is always positive. In the middle figure is displayed the volume filling constraint error $|1 - S_U(t_n)|$ obtained at each time step $n \in [1, N_T]$ where $S_U(t_n)$ is defined by

$$ S_U(t_n) := \inf_{\mu \in \mathcal{P}_{\text{train}}} \inf_{K \in \mathcal{T}} \sum_{i=1}^{N_s} u_{\mu, i, K}^n. \quad (4.6) $$

In particular we observe that this quantity is very close to the optimal value of 0 which shows that the finite volume procedure described in Section 2.2 preserves the volume filling constraint property. The right-hand-side figure illustrates the mass conservation property. We have displayed for each parameter $\mu \in \mathcal{P}_{\text{train}}$ the maximal deviation of the mass from the initial mass. This maximal deviation is given by

$$ \mathcal{M}_U(\mu) := \max_{i \in [1, N_s]} \max_{n \in [1, N_T]} \left| \sum_{K \in \mathcal{T}} m_K u_{\mu, i, K}^n - \int_{\Omega} u_i^0(x) \, dx \right|. \quad (4.7) $$

We observe that this deviation of the mass is of order $10^{-12}$ which is very close to computational accuracy. It proves the consistency of the finite volume procedure described in Section 2.2.

Finally, Figure 5 is a complement to Figure 4 and shows the exponential decay of the entropy to a constant value. Recall that the entropy at time $t_n$ is given by $E_T(U^n) := \sum_{K \in \mathcal{T}} \sum_{i=1}^{N_s} m_K u_{i, K}^n \ln(u_{i, K}^n)$.

Note that the left-hand side figure corresponds to the cross-diffusion coefficients $\mu^0$ defined by (4.4) whereas the middle and right figures corresponds to other sets of cross-diffusion coefficients: $\mu^9 = (a_{ij}^9)_{1 \le i \ne j \le N_s}$ and $\mu^{17} = (a_{ij}^{17})_{1 \le i \ne j \le N_s}$, with

$$ a_{12}^9 = 0.93, a_{13}^9 = 0.71, a_{23}^9 = 0.44 \quad \text{and} \quad a_{12}^{17} = 0.37, a_{13}^{17} = 0.004, a_{23}^{17} = 0.72. \quad (4.8) $$

### 4.1.2 The POD reduced order model

In this section, we show that the first POD reduced-order model does not preserve the structural properties of the solution.

In Figure 6 is displayed the violation of the structural properties of the solution for the first POD reduced model (see Section 3.1). In these figures, the cardinality of the reduced basis is $r = 2$. On the left figure, we considered the quantity $\mathcal{P}_{\bar{U}}(t_n)$ (see (4.5)) which is plotted as a function of $n \in [1, N_T]$. We observe that the positivity constraint is violated at many time steps. Next, in the middle figure, we have represented
<figure><img src="image_5.png" /><figcaption>FIGURE 5. Entropy functional for three values of the parameters $\mu \in \mathcal{P}_{\text{train}}$.</figcaption></figure>

step $n \in [1, N_T]$ where $S_U(t_n)$ is defined by

$$S_U(t_n) := \inf_{\mu \in \mathcal{P}_{\text{train}}} \inf_{K \in \mathcal{T}} \sum_{i=1}^{N_s} u_{\mu, i, K}^n. \qquad (4.6)$$

In particular we observe that this quantity is very close to the optimal value of 0 which shows that the finite volume procedure described in Section 2.2 preserves the volume filling constraint property. The right-hand-side figure illustrates the mass conservation property. We have displayed for each parameter $\mu \in \mathcal{P}_{\text{train}}$ the maximal deviation of the mass from the initial mass. This maximal deviation is given by

$$\mathcal{M}_U(\mu) := \max_{i \in [1, N_s]} \max_{n \in [1, N_T]} \left| \sum_{K \in \mathcal{T}} m_K u_{\mu, i, K}^n - \int_{\Omega} u_i^0(x) dx \right|. \qquad (4.7)$$

We observe that this deviation of the mass is of order $10^{-12}$ which is very close to computational accuracy. It proves the consistency of the finite volume procedure described is Section 2.2.

Finally, Figure 5 is a complement to Figure 4 and shows the exponential decay of the entropy to a constant value. Recall that the entropy at time $t_n$ is given by $E_T(U^n) := \sum_{K \in \mathcal{T}} \sum_{i=1}^{N_s} m_K u_{i, K}^n \ln(u_{i, K}^n)$.

Note that the left-hand side figure corresponds to the cross-diffusion coefficients $\mu^0$ defined by (4.4) whereas the middle and right figures corresponds to other sets of cross-diffusion coefficients: $\mu^9 = (a_{ij}^9)_{1 \le i \ne j \le N_s}$ and $\mu^{17} = (a_{ij}^{17})_{1 \le i \ne j \le N_s}$, with

$$a_{12}^9 = 0.93, a_{13}^9 = 0.71, a_{23}^9 = 0.44 \quad \text{and} \quad a_{12}^{17} = 0.37, a_{13}^{17} = 0.004, a_{23}^{17} = 0.72. \qquad (4.8)$$

### 4.1.2. The POD reduced order model

In this section, we show that the first POD reduced-order model does not preserve the structural properties of the solution.

In Figure 6 is displayed the violation of the structural properties of the solution for the first POD reduced model (see Sect. 3.1). In these figures, the cardinality of the reduced basis is $r = 2$. On the left figure, we considered the quantity $\mathcal{P}_{\tilde{U}}(t_n)$ (see (4.5)) which is plotted as a function of $n \in [1, N_T]$. We observe that the positivity constraint is violated at many time steps. Next, in the middle figure, we have represented the quantity $\mathcal{S}_{\tilde{U}}(t_n)$ (see (4.6)) as a function of $n \in [1, N_T]$. We also observe that $\mathcal{S}_{\tilde{U}}(t_n) \ne 1$ for many time steps so that the volume filling constraint is violated. On the right figure, we show the deviation of the mass $\mathcal{M}_{\tilde{U}}(\mu)$ (see (4.7)) of the reduced problem from the initial mass $\sum_{K \in \mathcal{T}} m_K \tilde{u}_{\mu, i, K}^0$. This deviation is important for several values of $\mu$ so that the mass conservation property is likewise violated.

Finally, in Figure 7 we show (left figure and middle figure) when $r = 2$ that the discrete entropy functional defined by (2.18) (where $U$ is replaced by $\tilde{U}_\mu^n$) is not a decreasing function. Note that here, the image in the middle coincides with the one in the left with a zoom on the time steps beyond 500. The right part of Figure 7 shows
<figure><img src="image_9.png" /><figcaption>Figure 5: Entropy functional for three values of the parameters $\mu \in \mathcal{P}_{\text{train}}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 6: Violation of the physical properties of the reduced solution. Left: violation of the positivity, middle: violation of the volume filling constraint, right: violation of the conservation of mass.</figcaption></figure>

the quantity $S_{\tilde{U}}(t_n)$ (see (4.6)) as a function of $n \in [1, N_T]$. We also observe that $S_{\tilde{U}}(t_n) \neq 1$ for many time steps so that the volume filling constraint is violated. On the right figure, we show the deviation of the mass $M_{\tilde{U}}(\mu)$ (see (4.7)) of the reduced problem from the initial mass $\sum_{K \in \mathcal{T}} m_K \tilde{u}_{\mu, i, K}^0$. This deviation is important for several values of $\mu$ so that the mass conservation property is likewise violated.

Finally, in Figure 7 we show (left figure and middle figure) when $r=2$ that the discrete entropy functional defined by (2.18) (where $U$ is replaced by $\tilde{U}_{\mu}^n$) is not a decreasing function. Note that here, the image in the middle coincides with the one in the left with a zoom on the time steps beyond 500. The right part of Figure 7 shows that the $L^{\infty}(\mathcal{P}_{\text{train}}, L^2(\Omega), L^2([0, T]))$ error between the POD-ROM solution and the fine solution decreases when the dimension of the reduced basis increases. This behavior is consistent showing a decrease to the optimal value of 0. Here, $u_{\mu}^i$ respectively $\tilde{u}_{\mu}^i$ is the functional representation of $U_{\mu}^i := (u_{\mu, i, K}^n)_{\mu \in \mathcal{P}_{\text{train}}, K \in \mathcal{T}, n \in [1, N_T]}$ respectively $\tilde{U}_{\mu}^i := (\tilde{u}_{\mu, i, K}^n)_{\mu \in \mathcal{P}_{\text{train}}, K \in \mathcal{T}, n \in [1, N_T]}$. Furthermore, the error is defined by

$$ \max_{i \in [1, N_s]} \|u_{\mu}^i - \tilde{u}_{\mu}^i\|_{L^{\infty}(\mathcal{P}_{\text{train}}, L^2(\Omega), L^2([0, T]))} := \max_{i \in [1, N_s]} \max_{\mu \in \mathcal{P}_{\text{train}}} \left( \int_0^T \|u_{\mu}^i - \tilde{u}_{\mu}^i\|_{L^2(\Omega)}^2 (t) \, dt \right)^{\frac{1}{2}}. \quad (4.9) $$

<figure><img src="image_11.png" /><figcaption>Figure 7: Violation of the decay of entropy and reduced model error.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 6. Violation of the physical properties of the reduced solution. Left: violation of the positivity, middle: violation of the volume filling constraint, right: violation of the conservation of mass.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 7. Violation of the decay of entropy and reduced model error.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 8. Violation of the physical properties of the reduced solution as a function of the dimension of the reduced problem.</figcaption></figure>

that the $L^\infty(\mathcal{P}_{\text{train}}, L^2(\Omega), L^2([0, T]))$ error between the POD-ROM solution and the fine solution decreases when the dimension of the reduced basis increases. This behavior is consistent showing a decrease to the optimal value of 0. Here, $u_\mu^i$ respectively $\tilde{u}_\mu^i$ is the functional representation of $U_\mu^i := (u_{\mu,i,K}^n)_{\mu \in \mathcal{P}_{\text{train}}, K \in \mathcal{T}, n \in [1,N_T]}$ respectively $\tilde{U}_\mu^i := (\tilde{u}_{\mu,i,K}^n)_{\mu \in \mathcal{P}_{\text{train}}, K \in \mathcal{T}, n \in [1,N_T]}$. Furthermore, the error is defined by

$$ \max_{i \in [1, N_s]} \|u_{\mu}^{i} - \tilde{u}_{\mu}^{i}\|_{L^{\infty}(\mathcal{P}_{\text{train}}, L^{2}(\Omega), L^{2}([0, T]))} := \max_{i \in [1, N_s]} \max_{\mu \in \mathcal{P}_{\text{train}}} \left( \int_{0}^{T} \|u_{\mu}^{i} - \tilde{u}_{\mu}^{i}\|_{L^{2}(\Omega)}^{2} (t) \, dt \right)^{\frac{1}{2}}. \quad (4.9) $$

In Figure 8 we have represented some inconsistencies of the reduced solution when the cardinality of the reduced basis increases. In the left figure is displayed for each value of $r$ the number of negative values that could have the numerical solution. In fact, it corresponds to the number of negative elements from the whole snapshot matrix $M \in \mathbb{R}^{N_h \times N_s \times p \times N_T}$. In the right part of Figure 8, we depicted the maximal deviation of the reduced mass
<figure><img src="image_12.png" /><figcaption>Figure 8: Violation of the physical properties of the reduced solution as a function of the dimension of the reduced problem</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 9: Preservation of the structural properties of the reduced solution. Left: positivity property, middle: volume filling constraint, right: conservation of mass.</figcaption></figure>

In Figure 8 we have represented some inconsistencies of the reduced solution when the cardinality of the reduced basis increases. In the left Figure is displayed for each value of $r$ the number of negative values that could have the numerical solution. In fact, it corresponds to the number of negative elements from the whole snapshot matrix $M \in \mathbb{R}^{N_h \times N_s \times p \times N_T}$. In the right part of Figure 8, we depicted the maximal deviation of the reduced mass from the initial one when the dimension of the reduced basis increases. More precisely, it corresponds to the quantity $\max_{\mu \in \mathcal{P}_{\text{train}}} M_{\bar{U}}(\mu)$. We observe that this quantity decreases fastly when the dimension of the reduced basis increases. This result is coherent with the fact that for smallest values of $r$ the POD-ROM is not accurate whereas for big values the constraints are no more violated.

### 4.1.3 SP reduced-order model

In this section, we propose numerical experiments for the structure preserving reduced order model described in Section 3.2.

Figure 9 is the analog of Figure 6 where we observe that the first three structural properties listed in Lemma 3.1 are satisfied. In these figures $r = 2$ and we employ (4.5), (4.6) and (4.7) with $\overline{U}_{\mu}^n$ instead of $\tilde{U}_{\mu}^n$.

Next, we represent in Figure 10 the evolution of the entropy functional for the cross-diffusion coefficient

<figure><img src="image_14.png" /><figcaption>Figure 10: Preservation of the entropy property: left and middle. Error of the reduced model: right.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIGURE 9. Preservation of the structural properties of the reduced solution. *Left*: positivity property, *middle*: volume filling constraint, *right*: conservation of mass.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 10. Preservation of the entropy property: *left* and *middle*. Error of the reduced model: *right*.</figcaption></figure>

from the initial one when the dimension of the reduced basis increases. More precisely, it corresponds to the quantity $\max_{\mu \in \mathcal{P}_{\text{train}}} M_{\bar{U}}(\mu)$. We observe that this quantity decreases fastly when the dimension of the reduced basis increases. This result is coherent with the fact that for smallest values of $r$ the POD-ROM is not accurate whereas for big values the constraints are no more violated.

### 4.1.3. SP reduced-order model

In this section, we propose numerical experiments for the structure preserving reduced order model described in Section 3.2.

Figure 9 is the analog of Figure 6 where we observe that the first three structural properties listed in Lemma 3.1 are satisfied. In these figures $r = 2$ and we employ (4.5), (4.6) and (4.7) with $\overline{U}_{\mu}^n$ instead of $\tilde{U}_{\mu}^n$.

Next, we represent in Figure 10 the evolution of the entropy functional for the cross-diffusion coefficient $\mu^0$ defined by (4.4) where the reduced solution is provided by the SP-reduced model of Section 3.2. In particular, we observe that the entropy decreases exponentially before reaching a constant profile, as for the high-fidelity model. Note that in the middle figure a zoom is performed on the time step beyond $n = 600$. The right figure compares the behavior of the two reduced-order model error. The blue curve corresponds to the error provided by the POD-ROM and the violet curve corresponds to the error given by the SP-ROM (this time we use (4.9) but with $\overline{v}_{\mu}^i$ instead of $\tilde{v}_{\mu}^i$). We see that the behavior of the two curves are similar. For small sizes of reduced basis (below 15), the error given by the SP-ROM is better than the error given by the POD-ROM. However, the reverse is observed for larger size of reduced basis. Let us mention here that this type of behaviour (loss of accuracy of structure-preserving reduced-order models) is also observed for other types of models like Port-Hamiltonian systems [12, 39]. Nevertheless, the SP-ROM is an interesting alternative since even for few POD modes it preserves the structural properties of the solution and has a suitable error behavior, which is of crucial importance in context where multiphysics simulations are performed.
<figure><img src="image_15.png" /><figcaption>Figure 11: Profile of the numerical solution at $n = 20$ (left), $n = 200$ (middle), and $n = 1700$ (right) for the cross diffusion coefficients $\mu^{17}$.</figcaption></figure>

$\mu^0$ defined by (4.4) where the reduced solution is provided by the SP-reduced model of Section 3.2. In particular, we observe that the entropy decreases exponentially before reaching a constant profile, as for the high-fidelity model. Note that in the middle figure a zoom is performed on the time step beyond $n = 600$. The right figure compares the behavior of the two reduced-order model error. The blue curve corresponds to the error provided by the POD-ROM and the violet curve corresponds to the error given by the SP-ROM (this time we use (4.9) but with $\bar{u}_{\mu}^i$ instead of $\tilde{u}_{\mu}^i$). We see that the behavior of the two curves are similar. For small sizes of reduced basis (below 15), the error given by the SP-ROM is better than the error given by the POD-ROM. However, the reverse is observed for larger size of reduced basis. Let us mention here that this type of behaviour (loss of accuracy of structure-preserving reduced-order models) is also observed for other types of models like Port-Hamiltonian systems [12, 39]. Nevertheless, the SP-ROM is an interesting alternative since even for few POD modes it preserves the structural properties of the solution and has a suitable error behavior, which is of crucial importance in context where multiphysics simulations are performed.

## 4.2 Test case 2: four species

In this test case, we consider an example inspired from applications related to the PVD fabrication process of thin film solar cells as proposed in [1] with $N_s = 4$ species (referring to Copper, Indium, Gallium and Selenium in photovoltaics applications). The numerical resolution of this particular parametric cross-diffusion system is particularly expensive and justifies the use of reduced-order-model techniques. We also compare in this test case the behaviors of the POD-ROM and SP-ROM. The final simulation time is $T = 0.5$ and we use a constant time step $\Delta t = \Delta t_n = 2.5 \times 10^{-4}$ so that $N_T = 2 \times 10^3$.

We again use the Newton solver described in Section 2.2.3. In this Section, we compute 20 snapshots of solutions so that Card($\mathcal{P}_{train}$) = 20. The definition of the parameter $a_{\mu}^*$ follows (4.1). Concerning the initial conditions we take

$$w_1^0(x) := e^{-25(x-0.5)^2}, \quad w_2^0(x) := x^2 + \varepsilon, \quad w_3^0(x) := 1 - e^{-25(x-0.5)^2}, \quad w_4^0(x) := |\sin(\pi x)| \quad (4.10)$$

where $\varepsilon = 10^{-6}$, and define

$$u_i^0(x) = \frac{w_i^0(x)}{\sum_{l=1}^{N_s} w_l^0(x)} \quad (4.11)$$

for all $1 \le i \le 4$ and $x \in (0, 1)$.

In Figure 11 we represent the shape of the high-fidelity numerical solution at three time steps at $n = 20$, $n = 200$, and $n = 1700$. Here, the cross-diffusion coefficients $\mu^{17} = (a_{ij}^{17})_{1 \le i \ne j \le 4}$ are given by

$$a_{12}^{17} = 0.64, \quad a_{13}^{17} = 0.31, \quad a_{14}^{17} = 0.53, \quad a_{23}^{17} = 0.99, \quad a_{24}^{17} = 0.84, \quad a_{34}^{17} = 0.99. \quad (4.12)$$

In Figure 12 are displayed the structural properties of the numerical solution, namely the positivity (left figure), the volume filling constraint (middle figure), and the mass conservation (right figure). As in Section 4.1.1, we observe that the three properties are preserved along the simulation.
<figure><img src="image_11.png" /><figcaption>FIGURE 11. Profile of the numerical solution at $n = 20$ (left), $n = 200$ (middle), and $n = 1700$ (right) for the cross diffusion coefficients $\mu^{17}$.</figcaption></figure>

## 4.2. Test case 2: four species

In this test case, we consider an example inspired from applications related to the PVD fabrication process of thin film solar cells as proposed in [1] with $N_s = 4$ species (referring to Copper, Indium, Gallium and Selenium in photovoltaics applications). The numerical resolution of this particular parametric cross-diffusion system is particularly expensive and justifies the use of reduced-order-model techniques. We also compare in this test case the behaviors of the POD-ROM and SP-ROM. The final simulation time is $T = 0.5$ and we use a constant time step $\Delta t = \Delta t_n = 2.5 \times 10^{-4}$ so that $N_T = 2 \times 10^3$.

We again use the Newton solver described in Section 2.2.3. In this section, we compute 20 snapshots of solutions so that $\text{Card}(\mathcal{P}_{\text{train}}) = 20$. The definition of the parameter $a_{\mu}^{*}$ follows (4.1). Concerning the initial conditions we take

$$w_1^0(x) := e^{-25(x-0.5)^2}, \quad w_2^0(x) := x^2 + \varepsilon, \quad w_3^0(x) := 1 - e^{-25(x-0.5)^2}, \quad w_4^0(x) := |\sin(\pi x)| \quad (4.10)$$

where $\varepsilon = 10^{-6}$, and define

$$u_i^0(x) = \frac{w_i^0(x)}{\sum_{l=1}^{N_s} w_l^0(x)} \quad (4.11)$$

for all $1 \le i \le 4$ and $x \in (0, 1)$.

In Figure 11 we represent the shape of the high-fidelity numerical solution at three time steps at $n = 20$, $n = 200$, and $n = 1700$. Here, the cross-diffusion coefficients $\mu^{17} = (a_{ij}^{17})_{1 \le i \ne j \le 4}$ are given by

$$a_{12}^{17} = 0.64, a_{13}^{17} = 0.31, a_{14}^{17} = 0.53, a_{23}^{17} = 0.99, a_{24}^{17} = 0.84, a_{34}^{17} = 0.99. \quad (4.12)$$

In Figure 12 are displayed the structural properties of the numerical solution, namely the positivity (left figure), the volume filling constraint (middle figure), and the mass conservation (right figure). As in Section 4.1.1, we observe that the three properties are preserved along the simulation.

In Figure 13 we represent the behavior of the entropy functional for three sets of cross-diffusion coefficients $\mu^0 := (a_{ij}^0)_{1 \le i \ne j \le 4}$, $\mu^8 := (a_{ij}^8)_{1 \le i \ne j \le 4}$, where

$$a_{12}^0 = 0.36, a_{13}^0 = 0.19, a_{14}^0 = 0.64, a_{23}^0 = 0.07, a_{24}^0 = 0.61, a_{34}^0 = 0.51,$$

$$a_{12}^8 = 0.69, a_{13}^8 = 0.30, a_{14}^8 = 0.16, a_{23}^8 = 0.37, a_{24}^8 = 0.95, a_{34}^8 = 0.38,$$

and $\mu^{17}$ defined in (4.12).

We observe again the fact that the entropy of the system converges exponentially fast to some limit value.
<figure><img src="image_16.png" /><figcaption>Figure 12: Properties of the fine numerical solution. Left: positivity of the solution, middle: preservation of the volume filling constraint, right: conservation of mass.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 13: Entropy functional for three parameters.</figcaption></figure>

In Figure 13 we represent the behavior of the entropy functional for three sets of cross-diffusion coefficients $\mu^0 := (a_{ij}^0)_{1 \le i \ne j \le 4}$, $\mu^8 := (a_{ij}^8)_{1 \le i \ne j \le 4}$, where

$$a_{12}^0 = 0.36, a_{13}^0 = 0.19, a_{14}^0 = 0.64, a_{23}^0 = 0.07, a_{24}^0 = 0.61, a_{34}^0 = 0.51,$$

$$a_{12}^8 = 0.69, a_{13}^8 = 0.30, a_{14}^8 = 0.16, a_{23}^8 = 0.37, a_{24}^8 = 0.95, a_{34}^8 = 0.38,$$

and $\mu^{17}$ defined in (4.12).

We observe again the fact that the entropy of the system converges exponentially fast to some limit value.

### 4.2.1 The POD reduced-order model

In Figure 14, we depict the violation of the structural properties of the numerical solution given by the POD-ROM. In the left Figure, we represent the number of cells where the solution could be negative as a function of the dimension of the reduced basis. More precisely, for each $r > 0$, is associated 20 set of $N_T$

<figure><img src="image_18.png" /><figcaption>Figure 14: Violation of the physical properties of the POD-ROM. Left: violation of the positivity, middle: violation of the volume filling constraint, right: violation of the mass conservation.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIGURE 12. Properties of the fine numerical solution. *Left*: positivity of the solution, *middle*: preservation of the volume filling constraint, *right*: conservation of mass.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 13. Entropy functional for three parameters.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 14. Violation of the physical properties of the POD-ROM. *Left*: violation of the positivity, *middle*: violation of the volume filling constraint, *right*: violation of the mass conservation.</figcaption></figure>

### 4.2.1. The POD reduced-order model

In Figure 14, we depict the violation of the structural properties of the numerical solution given by the POD-ROM. In the left figure, we represent the number of cells where the solution could be negative as a function of the dimension of the reduced basis. More precisely, for each $r > 0$, is associated 20 set of $N_T$ solutions. The results obtained in the left figure correspond to the number of negative values that appear in the whole set of $20 \times N_T$ possible vector of solutions.

We observe that for many configurations, we can find several negative elements. Note that in the middle figure is displayed the quantity $\inf_{n \in \llbracket 1, N_T \rrbracket} \inf_{\mu \in \mathcal{P}_{\text{train}}} \inf_{K \in \mathcal{T}} \sum_{i=1}^{N_s} \tilde{u}_{\mu, i, K}^n$ for several values of $r$. Concerning the
<figure><img src="image_19.png" /><figcaption>Figure 15: Violation of the decay of entropy for $r = 1$ (left), $r = 3$ (right) and error for POD-ROM.</figcaption></figure>

solutions. The results obtained in the left figure correspond to the number of negative values that appear in the whole set of $20 \times N_T$ possible vector of solutions.

We observe that for many configurations, we can find several negative elements. Note that in the middle figure is displayed the quantity $\inf_{n \in [1, N_T]} \inf_{\mu \in \mathcal{P}_{\text{train}}} \inf_{K \in \mathcal{T}} \sum_{i=1}^{N_s} \bar{u}_{\mu, i, K}^n$ for several values of $r$. Concerning the mass deviation, we have represented the quantity $\max_{\mu \in \mathcal{P}_{\text{train}}} \mathcal{M}_{\bar{U}}(\mu)$ as a function of $r$ where $\mathcal{M}_{\bar{U}}(\mu)$ is defined by (4.7). Thus, to obtain a reduced model respecting the first three properties of Lemma 3.1 we need to consider $r = 22$ basis vectors which is not appropriate. Figure 15 completes the results obtained from Figure 14. The discrete entropy functional is represented as a function of the time steps. We observe for $r = 1$ and $r = 3$ that the entropy increases beyond some time step which constitutes a physical violation. Note that for the middle figure we zoomed on the time steps where we observe the growth of the entropy. Finally, the figure on the right shows that the maximal deviation of the error defined by (4.9) is important for small values of $r$ and becomes very small when $r$ increases. This behavior is coherent and characteristic of POD reduced-order model.

### 4.2.2 The SP reduced-order model

This section is devoted to numerical experiments for the SP-reduced order model. We represent in Figure 16 the structural properties of the reduced solution. On the left Figure, we represented for several values of $r > 0$ the minimum value that can takes the numerical solution. More precisely we have represented the quantity

$$ \inf_{n \in [1, N_T]} \inf_{\mu \in \mathcal{P}_{\text{train}}} \inf_{K \in \mathcal{T}} \inf_{i \in [1, N_s]} \bar{u}_{\mu, i, K}^n $$

for several values of $r > 0$. We observe that the solution is always nonnegative which shows that the positivity constraint is satisfied for this reduced model. In the middle Figure, is verified the volume filling constraint property. We displayed the quantity

$$ \inf_{n \in [1, N_T]} \inf_{\mu \in \mathcal{P}_{\text{train}}} \inf_{K \in \mathcal{T}} \sum_{i=1}^{N_s} \bar{u}_{\mu, i, K}^n $$

for several values of $r > 0$. We observe that this quantity is always equal to the optimal value of 1. Finally, in the right Figure we have represented the error of the two reduced models as a function of the cardinality of the reduced basis. Recall that these error are defined by (4.9). We observe that the error decreases when the dimension of the reduced basis $r$ increases, at a similar rate for both approaches when $r \le 25$. For larger values, we also observed decreasing error with a small roughly constant shift between the two curves.

Then, in Figure 17 we have illustrated the behavior of the discrete reduced entropy for the cross-diffusion coefficients $\mu^{17}$ defined by (4.12) and for $r = 1, 2, 4$. Contrary to Figure 15 we can see that the entropy is a decreasing function with respect to the time steps. Thus, our SP-ROM indeed preserves the structural properties of the numerical solution.

### 4.2.3 Validation stage

In this section, we check the validity of the various reduced-order models for cross-diffusion coefficients which do not belong to the training set $\mathcal{P}_{\text{train}}$. To this aim, we select values of parameters $\mu$ in a subset $\mathcal{P}_{\text{valid}}$ of $\mathcal{P}$, which is different from $\mathcal{P}_{\text{train}}$ and compare the error between the high-fidelity model and the
<figure><img src="image_15.png" /><figcaption>FIGURE 15. Violation of the decay of entropy for $r = 1$ (left), $r = 3$ (right) and error for POD-ROM.</figcaption></figure>

mass deviation, we have represented the quantity $\max_{\mu \in \mathcal{P}_{\text{train}}} \mathcal{M}_{\bar{U}}(\mu)$ as a function of $r$ where $\mathcal{M}_{\bar{U}}(\mu)$ is defined by (4.7). Thus, to obtain a reduced model respecting the first three properties of Lemma 3.1 we need to consider $r = 22$ basis vectors which is not appropriate.

Figure 15 completes the results obtained from Figure 14. The discrete entropy functional is represented as a function of the time steps. We observe for $r = 1$ and $r = 3$ that the entropy increases beyond some time step which constitutes a physical violation. Note that for the middle figure we zoomed on the time steps where we observe the growth of the entropy. Finally, the figure on the right shows that the maximal deviation of the error defined by (4.9) is important for small values of $r$ and becomes very small when $r$ increases. This behavior is coherent and characteristic of POD reduced-order model.

### 4.2.2. The SP reduced-order model

This section is devoted to numerical experiments for the SP-reduced order model. We represent in Figure 16 the structural properties of the reduced solution. On the left figure, we represented for several values of $r > 0$ the minimum value that can takes the numerical solution. More precisely we have represented the quantity

$$ \inf_{n \in [1, N_T]} \inf_{\mu \in \mathcal{P}_{\text{train}}} \inf_{K \in \mathcal{T}} \inf_{i \in [1, N_s]} \overline{u}_{\mu, i, K}^n $$

for several values of $r > 0$. We observe that the solution is always nonnegative which shows that the positivity constraint is satisfied for this reduced model. In the middle figure, is verified the volume filling constraint property. We displayed the quantity

$$ \inf_{n \in [1, N_T]} \inf_{\mu \in \mathcal{P}_{\text{train}}} \inf_{K \in \mathcal{T}} \sum_{i=1}^{N_s} \overline{u}_{\mu, i, K}^n $$

for several values of $r > 0$. We observe that this quantity is always equal to the optimal value of 1. Finally, in the right figure we have represented the error of the two reduced models as a function of the cardinality of the reduced basis. Recall that these error are defined by (4.9). We observe that the error decreases when the dimension of the reduced basis $r$ increases, at a similar rate for both approaches when $r \le 25$. For larger values, we also observed decreasing error with a small roughly constant shift between the two curves.

Then, in Figure 17 we have illustrated the behavior of the discrete reduced entropy for the cross-diffusion coefficients $\mu^{17}$ defined by (4.12) and for $r = 1, 2, 4$. Contrary to Figure 15 we can see that the entropy is a decreasing function with respect to the time steps. Thus, our SP-ROM indeed preserves the structural properties of the numerical solution.

### 4.2.3. Validation stage

In this section, we check the validity of the various reduced-order models for cross-diffusion coefficients which do not belong to the training set $\mathcal{P}_{\text{train}}$. To this aim, we select values of parameters $\mu$ in a subset $\mathcal{P}_{\text{valid}}$ of
<figure><img src="image_20.png" /><figcaption>Figure 16: Preservation of the structural properties of the SP-ROM. Left : positivity property, middle : volume filling constraint, right : reduced model error.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 17: Entropy property of the SP-ROM. Left : r = 1, middle : r = 2, right : r = 4.</figcaption></figure>

ROMs. Here, Card($P_{valid}$) = 20 and to sample the coefficients of the cross-diffusion matrices we proceed as for the matrices of $P_{train}$. The results reported in Figure 18 show that the solutions given by a fine resolution evolve to constant profiles which is the typical scenario of cross-diffusion system simulations.

Figure 19 is the analog of Figure 17 and Figure 16. We observe that the positivity constraint, the volume filling constraint, and the mass conservation property are satisfied for each reduced model SP-ROM. Furthermore, the decay of entropy is also found back and we have displayed the particular case where the cross-diffusion coefficients $\mu^{32} = (a_{ij}^{32})_{1 \le i \ne j \le 4}$ are defined by

$$a_{12}^{32} = 0.22, a_{13}^{32} = 0.37, a_{14}^{32} = 0.17, a_{23}^{32} = 0.53, a_{24}^{32} = 0.97, a_{34}^{32} = 0.82.$$

Finally, the results that we obtained for other parameters $\mu \in P_{valid}$ are similar to the results obtained in the previous section with $\mu \in P_{train}$.

# 5 Conclusion

In this work, we proposed a structure-preserving reduced-order model for a cross-diffusion system, which

<figure><img src="image_22.png" /><figcaption>Figure 18: Profile of the numerical solution at n = 20 (left), n = 420 (middle), and n = 1960 (right) for the cross-diffusion matrix $A_{27}$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIGURE 16. Preservation of the structural properties of the SP-ROM. Left: positivity property, middle: volume filling constraint, right: reduced model error.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 17. Entropy property of the SP-ROM. Left: $r = 1$, middle: $r = 2$, right: $r = 4$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 18. Profile of the numerical solution at $n = 20$ (left), $n = 420$ (middle), and $n = 1960$ (right) for the cross-diffusion matrix $A_{27}$.</figcaption></figure>

$P$, which is different from $P_{train}$ and compare the error between the high-fidelity model and the ROMs. Here, $\text{Card}(P_{\text{valid}}) = 20$ and to sample the coefficients of the cross-diffusion matrices we proceed as for the matrices of $P_{train}$. The results reported in Figure 18 show that the solutions given by a fine resolution evolve to constant profiles which is the typical scenario of cross-diffusion system simulations.

Figure 19 is the analog of Figures 17 and 16. We observe that the positivity constraint, the volume filling constraint, and the mass conservation property are satisfied for each reduced model SP-ROM. Furthermore, the decay of entropy is also found back and we have displayed the particular case where the cross-diffusion coefficients $\mu^{32} = (a_{ij}^{32})_{1 \le i \ne j \le 4}$ are defined by

$$a_{12}^{32} = 0.22, a_{13}^{32} = 0.37, a_{14}^{32} = 0.17, a_{23}^{32} = 0.53, a_{24}^{32} = 0.97, a_{34}^{32} = 0.82.$$
<figure><img src="image_23.png" /><figcaption>Figure 19: Validation stage for SP-ROM. Top left: positivity constraint, top middle : volume filling constraint property, top right : mass deviation, down left : entropy property for $r = 1$, down middle : entropy property for $r = 10$.</figcaption></figure>

yields comparable accuracy with respect to a standard POD-ROM while preserving the main features of the original model. The numerical experiments confirmed the theoretical properties of the proposed approach. Several open questions remain to be tackled: first, for the SP-ROM to be online efficient, it will be needed to use it in conjunction with an appropriate hyperreduction method, like for instance an Empirical Interpolation Method. It is not clear however how to design such an empirical interpolation method which would ensure that the resulting ROM would still preserve the desired properties such as non-negativeness of the solutions or entropy decay. Moreover, the design of appropriate a posteriori estimators would make possible the use of a greedy algorithm for the selection of the reduced basis, which would be cheaper than the POD selection which is done here. We leave these questions for future work.

## References

[1] A. BAKHTA AND V. EHRLACHER, *Cross-diffusion systems with non-zero flux and moving boundary conditions*, ESAIM Math. Model. Numer. Anal., 52 (2018), pp. 1385–1415.

[2] M. BARRAULT, Y. MADAY, N. C. NGUYEN, AND A. T. PATERA, *An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations*, Comptes Rendus Mathematique, 339 (2004), pp. 667–672.

[3] J. W. BARRETT AND J. F. BLOWEY, *Finite element approximation of a nonlinear cross-diffusion population model*, Numerische Mathematik, 98 (2004), pp. 195–221.

[4] T. J. BARTH, *Numerical methods for gasdynamic systems on unstructured meshes*, in An Introduction to Recent Developments in Theory and Numerics for Conservation Laws: Proceedings of the International School on Theory and Numerics for Conservation Laws, Freiburg/Littenweiler, October 20–24, 1997, Springer, 1999, pp. 195–285.

[5] I. BEN GHARBIA, J. DABAGHI, V. MARTIN, AND M. VOHRALÍK, *A posteriori error estimates for a compositional two-phase flow with nonlinear complementarity constraints*, Comput. Geosci., 24 (2020), pp. 1031–1055.

[6] G. BERKOOZ, P. HOLMES, AND J. L. LUMLEY, *The proper orthogonal decomposition in the analysis of turbulent flows*, in Annual review of fluid mechanics, Vol. 25, Annual Reviews, Palo Alto, CA, 1993, pp. 539–575.
<figure><img src="image_19.png" /><figcaption>FIGURE 19. Validation stage for SP-ROM. *Top left:* positivity constraint, *top middle:* volume filling constraint property, *top right:* mass deviation, *down left:* entropy property for $r = 1$, *down middle:* entropy property for $r = 10$.</figcaption></figure>

Finally, the results that we obtained for other parameters $\mu \in P_{\text{valid}}$ are similar to the results obtained in the previous section with $\mu \in P_{\text{train}}$.

## 5. CONCLUSION

In this work, we proposed a structure-preserving reduced-order model for a cross-diffusion system, which yields comparable accuracy with respect to a standard POD-ROM while preserving the main features of the original model. The numerical experiments confirmed the theoretical properties of the proposed approach. Several open questions remain to be tackled: first, for the SP-ROM to be online efficient, it will be needed to use it in conjunction with an appropriate hyperreduction method, like for instance an Empirical Interpolation Method. It is not clear however how to design such an empirical interpolation method which would ensure that the resulting ROM would still preserve the desired properties such as non-negativeness of the solutions or entropy decay. Moreover, the design of appropriate *a posteriori* estimators would make possible the use of a greedy algorithm for the selection of the reduced basis, which would be cheaper than the POD selection which is done here. We leave these questions for future work.

## ACKNOWLEDGEMENTS

This project was also funded by the European Union (ERC, HighLEAP, 101077204). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This project has received funding from the ANR project "COMODO" (ANR-19-CE46-0002).

## REFERENCES

[1] A. Bakhta and V. Ehrlacher, Cross-diffusion systems with non-zero flux and moving boundary conditions. *ESAIM Math. Model. Numer. Anal.* **52** (2018) 1385–1415.

[2] M. Barrault, Y. Maday, N.C. Nguyen and A.T. Patera, An “empirical interpolation” method: application to efficient reduced-basis discretization of partial differential equations. *C. R. Math. 339* (2004) 667–672.
[7] N. BOTTA AND M. PANDOLFI, *Upwind formulations for the euler equations in steady supersonic flows*, AIAA journal, 27 (1989), pp. 293-298.

[8] M. BURGER, M. DI FRANCESCO, J.-F. PIETSCHMANN, AND B. SCHLAKE, *Nonlinear cross-diffusion with size exclusion*, SIAM J. Math. Anal., 42 (2010), pp. 2842-2871.

[9] N. CAGNIART, Y. MADAY, AND B. STAMM, *Model order reduction for problems with large convection effects*, Contributions to partial differential equations and applications, (2019), pp. 131-150.

[10] C. CANCÈS AND B. GAUDEUL, *A convergent entropy diminishing finite volume scheme for a cross-diffusion system*, SIAM J. Numer. Anal., 58 (2020), pp. 2684-2710.

[11] C. CANCÈS, V. EHRLACHER, AND L. MONASSE, *Finite volumes for the Stefan-Maxwell cross-diffusion system*, IMA Journal of Numerical Analysis, (2023), p. drad032.

[12] S. CHATURANTABUT, C. BEATTIE, AND S. GUGERCIN, *Structure-preserving model reduction for non-linear port-hamiltonian systems*, SIAM Journal on Scientific Computing, 38 (2016), pp. B837-B865.

[13] V. EHRLACHER, D. LOMBARDI, O. MULA, AND F.-X. VIALARD, *Nonlinear model reduction on metric spaces. application to one-dimensional conservative pdes in wasserstein spaces*, ESAIM: Mathematical Modelling and Numerical Analysis, 54 (2020), pp. 2159-2197.

[14] S. C. EISENSTAT AND H. F. WALKER, *Globally convergent inexact Newton methods*, SIAM J. Optim., 4 (1994), pp. 393-422.

[15] R. EYMARD, T. GALLOUËT, AND R. HERBIN, *Finite volume methods*, in Handbook of numerical analysis, Vol. VII, Handb. Numer. Anal., VII, North-Holland, Amsterdam, 2000, pp. 713-1020.

[16] L. FICK, Y. MADAY, A. T. PATERA, AND T. TADDEI, *A stabilized POD model for turbulent flows over a range of Reynolds numbers: optimal parameter sampling and constrained projection*, J. Comput. Phys., 371 (2018), pp. 214-243.

[17] R. L. FOX AND H. MIURA, *An approximate analysis technique for design calculations*, AIAA Journal, 9 (1971), pp. 177-179.

[18] Y. GONG, Q. WANG, AND Z. WANG, *Structure-preserving galerkin pod reduced-order modeling of hamiltonian systems*, Computer Methods in Applied Mechanics and Engineering, 315 (2017), pp. 780-798.

[19] M. GUBISCH AND S. VOLKWEIN, *Proper orthogonal decomposition for linear-quadratic optimal control*, in Model reduction and approximation, vol. 15 of Comput. Sci. Eng., SIAM, Philadelphia, PA, 2017, pp. 3-63.

[20] J. S. HESTHAVEN, C. PAGLIANTINI, AND N. RIPAMONTI, *Structure-preserving model order reduction of hamiltonian systems*, arXiv preprint arXiv:2109.12367, (2021).

[21] J. S. HESTHAVEN, C. PAGLIANTINI, AND G. ROZZA, *Reduced basis methods for time-dependent problems*, Acta Numerica, 31 (2022), pp. 265-345.

[22] J. S. HESTHAVEN, G. ROZZA, AND B. STAMM, *Certified reduced basis methods for parametrized partial differential equations*, SpringerBriefs in Mathematics, Springer, Cham; BCAM Basque Center for Applied Mathematics, Bilbao, 2016. BCAM SpringerBriefs.

[23] C. HUANG, C. R. WENTLAND, K. DURAISAMY, AND C. MERKLE, *Model reduction for multi-scale transport problems using model-form preserving least-squares projections with variable transformation*, Journal of Computational Physics, 448 (2022), p. 110742.

[24] T. L. JACKSON AND H. M. BYRNE, *A mechanical model of tumor encapsulation and transcapsular spread*, vol. 180, 2002, pp. 307-328. John A. Jacquez memorial volume.
[3] J.W. Barrett and J.F. Blowey, Finite element approximation of a nonlinear cross-diffusion population model. *Numer. Math.* **98** (2004) 195–221.

[4] T.J. Barth, Numerical methods for gasdynamic systems on unstructured meshes, in *An Introduction to Recent Developments in Theory and Numerics for Conservation Laws: Proceedings of the International School on Theory and Numerics for Conservation Laws*, Freiburg/Littenweiler, October 20–24, 1997. Springer (1999) 195–285.

[5] I. Ben Gharbia, J. Dabaghi, V. Martin and M. Vohralík, A posteriori error estimates for a compositional two-phase flow with nonlinear complementarity constraints. *Comput. Geosci.* **24** (2020) 1031–1055.

[6] G. Berkooz, P. Holmes and J.L. Lumley, The proper orthogonal decomposition in the analysis of turbulent flows, in *Annual Review of Fluid Mechanics*. Vol. 25. Annual Reviews, Palo Alto, CA (1993) 539–575.

[7] N. Botta and M. Pandolfi, Upwind formulations for the euler equations in steady supersonic flows. *AIAA J.* **27** (1989) 293–298.

[8] M. Burger, M. Di Francesco, J.-F. Pietschmann and B. Schlake, Nonlinear cross-diffusion with size exclusion. *SIAM J. Math. Anal.* **42** (2010) 2842–2871.

[9] N. Cagniart, Y. Maday and B. Stamm, Model order reduction for problems with large convection effects, in *Contributions to Partial Differential Equations and Applications*. Springer (2019) 131–150.

[10] C. Cancès and B. Gaudeul, A convergent entropy diminishing finite volume scheme for a cross-diffusion system. *SIAM J. Numer. Anal.* **58** (2020) 2684–2710.

[11] C. Cancès, V. Ehrlacher and L. Monasse, Finite volumes for the Stefan–Maxwell cross-diffusion system. *IMA J. Numer. Anal.* **44** (2024) 1029–1060.

[12] S. Chaturantabut, C. Beattie and S. Gugercin, Structure-preserving model reduction for nonlinear port-Hamiltonian systems. *SIAM J. Sci. Comput.* **38** (2016) B837–B865.

[13] V. Ehrlacher, D. Lombardi, O. Mula and F.-X. Vialard, Nonlinear model reduction on metric spaces. Application to one-dimensional conservative PDEs in Wasserstein spaces. *ESAIM: Math. Modell. Numer. Anal.* **54** (2020) 2159–2197.

[14] S.C. Eisenstat and H.F. Walker, Globally convergent inexact Newton methods. *SIAM J. Optim.* **4** (1994) 393–422.

[15] R. Eymard, T. Gallouët and R. Herbin, Finite volume methods, in *Handbook of Numerical Analysis. Handb. Numer. Anal., VII*. Vol. VII. North-Holland, Amsterdam (2000) 713–1020.

[16] L. Fick, Y. Maday, A.T. Patera and T. Taddei, A stabilized POD model for turbulent flows over a range of Reynolds numbers: optimal parameter sampling and constrained projection. *J. Comput. Phys.* **371** (2018) 214–243.

[17] R.L. Fox and H. Miura, An approximate analysis technique for design calculations. *AIAA J.* **9** (1971) 177–179.

[18] Y. Gong, Q. Wang and Z. Wang, Structure-preserving Galerkin POD reduced-order modeling of Hamiltonian systems. *Comput. Methods Appl. Mech. Eng.* **315** (2017) 780–798.

[19] M. Gubisch and S. Volkwein, Proper orthogonal decomposition for linear-quadratic optimal control, in *Model Reduction and Approximation*. Vol. 15 of *Comput. Sci. Eng*. SIAM, Philadelphia, PA (2017) 3–63.

[20] J.S. Hesthaven, C. Pagliantini and N. Ripamonti, Structure-preserving model order reduction of Hamiltonian systems. Preprint arXiv:2109.12367 (2021).

[21] J.S. Hesthaven, C. Pagliantini and G. Rozza, Reduced basis methods for time-dependent problems. *Acta Numer.* **31** (2022) 265–345.

[22] J.S. Hesthaven, G. Rozza and B. Stamm, Certified Reduced Basis Methods for Parametrized Partial Differential Equations. SpringerBriefs in Mathematics, Springer, Cham; BCAM Basque Center for Applied Mathematics, Bilbao, BCAM SpringerBriefs (2016).

[23] C. Huang, C.R. Wentland, K. Duraisamy and C. Merkle, Model reduction for multi-scale transport problems using model-form preserving least-squares projections with variable transformation. *J. Comput. Phys.* **448** (2022) 110742.

[24] T.L. Jackson and H.M. Byrne, A mechanical model of tumor encapsulation and transcapsular spread. *Math. Biosci.* **180** (2002) 307–328.

[25] A. Jüngel, The boundedness-by-entropy method for cross-diffusion systems. *Nonlinearity* **28** (2015) 1963–2001.

[26] A. Jungel and I.V. Stelzer, Existence analysis of Maxwell–Stefan systems for multicomponent mixtures. *SIAM J. Math. Anal.* **45** (2013) 2421–2440.

[27] A. Jüngel and A. Zurek, A finite-volume scheme for a cross-diffusion model arising from interacting many-particle population systems, in *International Conference on Finite Volumes for Complex Applications*. Springer (2020) 223–231.

[28] A. Jüngel and A. Zurek, A discrete boundedness-by-entropy method for finite-volume approximations of cross-diffusion systems. *IMA J. Numer. Anal.* **43** (2023) 560–589.
[25] A. JÜNGEL, *The boundedness-by-entropy method for cross-diffusion systems*, Nonlinearity, 28 (2015), pp. 1963–2001.

[26] A. JUNGEL AND I. V. STELZER, *Existence analysis of maxwell–stefan systems for multicomponent mixtures*, SIAM Journal on Mathematical Analysis, 45 (2013), pp. 2421–2440.

[27] A. JÜNGEL AND A. ZUREK, *A finite-volume scheme for a cross-diffusion model arising from interacting many-particle population systems*, in International Conference on Finite Volumes for Complex Applications, Springer, 2020, pp. 223–231.

[28] ——, *A discrete boundedness-by-entropy method for finite-volume approximations of cross-diffusion systems*, IMA Journal of Numerical Analysis, 43 (2023), pp. 560–589.

[29] C. T. KELLEY, *Solving nonlinear equations with Newton's method*, vol. 1 of Fundamentals of Algorithms, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2003.

[30] K. KERGRENE, L. CHAMOIN, M. LAFOREST, AND S. PRUDHOMME, *On a goal-oriented version of the proper generalized decomposition method*, J. Sci. Comput., 81 (2019), pp. 92–111.

[31] K. KUNISCH AND S. VOLKWEIN, *Control of the Burgers equation by a reduced-order approach using proper orthogonal decomposition*, J. Optim. Theory Appl., 102 (1999), pp. 345–371.

[32] ——, *Galerkin proper orthogonal decomposition methods for parabolic problems*, Numer. Math., 90 (2001), pp. 117–148.

[33] P. LADEVÈZE AND L. CHAMOIN, *On the verification of model reduction methods based on the proper generalized decomposition*, Comput. Methods Appl. Mech. Engrg., 200 (2011), pp. 2032–2047.

[34] Y. MADAY, *Reduced basis method for the rapid and reliable solution of partial differential equations*, in International Congress of Mathematicians. Vol. III, Eur. Math. Soc., Zürich, 2006, pp. 1255–1270.

[35] A. K. NOOR AND J. M. PETERS, *Reduced basis technique for nonlinear analysis of structures*, AIAA Journal, 18 (1980), pp. 455–462.

[36] A. NOUY, *A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations*, Comput. Methods Appl. Mech. Engrg., 199 (2010), pp. 1603–1626.

[37] E. J. PARISH AND F. RIZZI, *On the impact of dimensionally-consistent and physics-based inner products for pod-galerkin and least-squares model reduction of compressible flows*, Journal of Computational Physics, 491 (2023), p. 112387.

[38] A. QUARTERONI, A. MANZONI, AND F. NEGRI, *Reduced basis methods for partial differential equations*, vol. 92 of Unitext, Springer, Cham, 2016. An introduction, La Matematica per il 3+2.

[39] P. SCHWERDTNER, T. MOSER, V. MEHRMANN, AND M. VOIGT, *Structure-preserving model order reduction for index one port-hamiltonian descriptor systems*, arXiv preprint arXiv:2206.01608, (2022).

[40] N. SHIGESADA, K. KAWASAKI, AND E. TERAMOTO, *Spatial segregation of interacting species*, J. Theoret. Biol., 79 (1979), pp. 83–99.

[41] G. WELPER, *Interpolation of functions with parameter dependent jumps by transformed snapshots*, SIAM Journal on Scientific Computing, 39 (2017), pp. A1225–A1250.

[42] ——, *Transformed snapshot interpolation with high resolution transforms*, SIAM Journal on Scientific Computing, 42 (2020), pp. A2037–A2061.

[43] N. ZAMPONI AND A. JÜNGEL, *Analysis of degenerate cross-diffusion population models with volume filling*, in *Annales de l'Institut Henri Poincaré C, Analyse non linéaire*, vol. 34, Elsevier, 2017, pp. 1–29.
[29] C.T. Kelley, Solving Nonlinear Equations with Newton'S Method. Vol. 1 of *Fundamentals of Algorithms*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA (2003).

[30] K. Kergrene, L. Chamoin, M. Laforest and S. Prudhomme, On a goal-oriented version of the proper generalized decomposition method. *J. Sci. Comput.* **81** (2019) 92–111.

[31] K. Kunisch and S. Volkwein, Control of the Burgers equation by a reduced-order approach using proper orthogonal decomposition. *J. Optim. Theory Appl.* **102** (1999) 345–371.

[32] K. Kunisch and S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems. *Numer. Math.* **90** (2001) 117–148.

[33] P. Ladevèze and L. Chamoin, On the verification of model reduction methods based on the proper generalized decomposition. *Comput. Methods Appl. Mech. Eng.* **200** (2011) 2032–2047.

[34] Y. Maday, Reduced basis method for the rapid and reliable solution of partial differential equations, in International Congress of Mathematicians. Vol. III. Eur. Math. Soc., Zürich (2006) 1255–1270.

[35] A.K. Noor and J.M. Peters, Reduced basis technique for nonlinear analysis of structures. *AIAA J.* **18** (1980) 455–462.

[36] A. Nouy, A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations. *Comput. Methods Appl. Mech. Eng.* **199** (2010) 1603–1626.

[37] E.J. Parish and F. Rizzi, On the impact of dimensionally-consistent and physics-based inner products for POD-Galerkin and least-squares model reduction of compressible flows. *J. Comput. Phys.* **491** (2023) 112387.

[38] A. Quarteroni, A. Manzoni and F. Negri, Reduced Basis Methods for Partial Differential Equations: An Introduction. Vol. 92 of *Unitext*. La Matematica per il 3+2. Springer, Cham (2016).

[39] P. Schwerdtner, T. Moser, V. Mehrmann and M. Voigt, Structure-preserving model order reduction for index one port-Hamiltonian descriptor systems. Preprint arXiv:2206.01608 (2022).

[40] N. Shigesada, K. Kawasaki and E. Teramoto, Spatial segregation of interacting species. *J. Theor. Biol.* **79** (1979) 83–99.

[41] G. Welper, Interpolation of functions with parameter dependent jumps by transformed snapshots. *SIAM J. Sci. Comput.* **39** (2017) A1225–A1250.

[42] G. Welper, Transformed snapshot interpolation with high resolution transforms. *SIAM J. Sci. Comput.* **42** (2020) A2037–A2061.

[43] N. Zamponi and A. Jüngel, Analysis of degenerate cross-diffusion population models with volume filling, in Annales de l'Institut Henri Poincaré C, Analyse non linéaire, Vol. 34. Elsevier (2017) 1–29.

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

Please help to maintain this journal in open access!

This journal is currently published in open access under the *Subscribe to Open model (S2O)*. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.