Friedrichs' systems discretized with the Discontinuous Galerkin method:
  domain decomposable model order reduction and Graph Neural Networks
  approximating vanishing viscosity solutions
Friedrichs' systems (FS) are symmetric positive linear systems of first-order
partial differential equations (PDEs), which provide a unified framework for
describing various elliptic, parabolic and hyperbolic semi-linear PDEs such as
the linearized Euler equations of gas dynamics, the equations of compressible
linear elasticity and the Dirac-Klein-Gordon system. FS were studied to
approximate PDEs of mixed elliptic and hyperbolic type in the same domain. For
this and other reasons, the versatility of the discontinuous Galerkin method
(DGM) represents the best approximation space for FS. We implement a
distributed memory solver for stationary FS in deal.II. Our focus is model
order reduction. Since FS model hyperbolic PDEs, they often suffer from a slow
Kolmogorov n-width decay. We develop two approaches to tackle this problem. The
first is domain decomposable reduced-order models (DD-ROMs). We will show that
the DGM offers a natural formulation of DD-ROMs, in particular regarding
interface penalties, compared to the continuous finite element method. We also
develop new repartitioning strategies to obtain more efficient local
approximations of the solution manifold. The second approach involves graph
neural networks used to infer the limit of a succession of projection-based
linear ROMs corresponding to lower viscosity constants: the heuristic behind is
to develop a multi-fidelity super-resolution paradigm to mimic the mathematical
convergence to vanishing viscosity solutions while exploiting to the most
interpretable and certified projection-based ROMs.
Francesco Romor, Davide Torlo, Gianluigi Rozza
Friedrichs' systems discretized with the Discontinuous Galerkin method: domain decomposable model order reduction and Graph Neural Networks approximating vanishing viscosity solutions

###### Abstract

Friedrichs' systems (FS) are symmetric positive linear systems of first-order partial differential equations (PDEs), which provide a unified framework for describing various elliptic, parabolic and hyperbolic semi-linear PDEs such as the linearized Euler equations of gas dynamics, the equations of compressible linear elasticity and the Dirac-Klein-Gordon system. FS were studied to approximate PDEs of mixed elliptic and hyperbolic type in the same domain. For this and other reasons, the versatility of the discontinuous Galerkin method (DGM) represents the best approximation space for FS. We implement a distributed memory solver for stationary FS in deal.II. Our focus is model order reduction. Since FS model hyperbolic PDEs, they often suffer from a slow Kolmogorov n-width decay. We develop two approaches to tackle this problem. The first is domain decomposable reduced-order models (DD-ROMs). We will show that the DGM offers a natural formulation of DD-ROMs, in particular regarding interface penalties, compared to the continuous finite element method. We also develop new repartitioning strategies to obtain more efficient local approximations of the solution manifold. The second approach involves graph neural networks used to infer the limit of a succession of projection-based linear ROMs corresponding to lower viscosity constants: the heuristic behind is to develop a multi-fidelity super-resolution paradigm to mimic the mathematical convergence to vanishing viscosity solutions while exploiting to the most interpretable and certified projection-based ROMs.

## 1 Introduction

Friedrichs' systems (FS) are a class of symmetric positive linear systems of first-order partial derivative equations (PDEs). They were introduced by Friedrichs [34] as a tool to study hyperbolic and elliptic phenomena in different parts of the domain within a unifying framework. The main ideas that allow recasting many models into the FS frameworks are the introduction of extra variables to lower the order of the higher derivatives and the linearization of nonlinear problems. FS are characterized by linear and positive operators and (non-uniquely defined) boundary operators that allow them to impose classical boundary conditions (BCs). Various works proved uniqueness, existence and well-posedness of the FS in their strong, weak and ultraweak formulation and the necessary conditions to properly define the boundary operators [34, 73, 74, 28, 31, 3, 24].

In the last decades, different numerical discretizations of the FS have been proposed to approximate the analytical solutions. The strategies vary among finite volume [79] and discontinuous Galerkin (DG) formulations [45, 52, 28, 29, 31, 30, 12, 17]. Along with the DG discretization, also error estimation analysis that provide, according to the type of edge penalization, optimal or sub-optimal estimates, have been carried out [28, 29, 24]. We focus on the DG method since it is more versatile to approximate both elliptic and hyperbolic PDEs and it fits naturally in the framework of domain decomposable ROMs (DD-ROMs).

In the context of parametric PDEs, for multi-query tasks or real-time simulations, fast and reliable simulations of the same problem for different parameters are often needed. This is especially true when the full-order models (FOMs) are based on expensive and high-order DG discretizations. Reduced order models (ROMs) decrease the computational costs looking for the solutions of unseen parametric instances on low-dimensional discretization spaces, called reduced basis spaces. This is possible because the new solutions to be predicted are expected to be highly correlated with the database of training DG solutions used to build the reduced spaces. ROMs have been proven to be a powerful tool for many applications
[69, 62, 44, 77]. In particular for linear problems, classical Galerkin and Petrov-Galerkin projection methods are very easy to set up and extremely convenient in terms of computational costs. FS are perfectly suited for such algorithms due to their linearity. This is a preliminary step needed to reduce parametric nonlinear PDEs whose linearization results in FS. In the most simple formulation, we will apply singular value decomposition (SVD) to compress a database of snapshots and provide a reliable reduced order model (ROM), with standard _a posteriori_ error estimators.

In the context of model order reduction, FS are particularly beneficial as theoretical frameworks for many reasons. They represent a new form of structure-preserving ROMs: the positive symmetric properties of FS are in fact easily inherited by the reduced numerical formulations. This advocates for the employment of FS for reduced order modelling whenever a PDE can be reformulated in the FS framework. This is the case for the Euler equations of gas dynamics, when they are written in terms of entropy variables [79, 66]. The same rationale is behind structure-preserving symplectic or Hamiltonian ROMs [43] and port-Hamiltonian ROMs [84, 10]. Moreover, since FS are often studied in their ultraweak formulation, they are good candidates for optimally stable error estimates [13] at the full-order level [12], also in a hybridized DG implementation in [17], and at the reduced order level, similarly to what has been achieved in the works [11, 39, 42]. Finally, from the point of view of software design, the possibility to implement in a unique maintainable and generic manner the realization of ROMs for PDEs ascribed to the class of FS is a convenient feature to search for.

Though being linear, FS are hyperbolic systems and often show an advection dominated character, which is not easily approximable through a simple proper orthogonal decomposition (POD). This leads to a slow Kolmogorov \(n\)-width (KnW) decay that results in very inefficient approximations of the reduced models. Several approaches have been studied to overcome this difficulty [80, 49, 67, 75, 16, 15, 2, 59, 20, 82, 51].

A strategy that has been developed to reduce PDEs solved numerically with domain decomposition approaches, like fluid-structure interaction systems, are domain decomposable ROMS (DD-ROMS). The initial formulations [62, 63, 48, 27, 48] involved continuous finite elements discretizations for which new ways to couple the solutions restricted to different subdomains needed to be devised, especially to enforce continuity at the interfaces. We show that the DGM imposes naturally flux interface penalties from the full-order discretizations and it is, thus, amenable for straightforward implementations of DD-ROMs. From the point of view of solution manifold approximability and so KnW decay, DD-ROMs are based on local linear approximants that are employed to reach a higher accuracy for unseen solutions. This is useful when the computational domain is divided in subdomains that are independently affected by the parametric instances. The typical case in which this may happen are parametric models for which discontinuous values of the parameters over fixed subdomains cause non correlated responses on their respective subdomains. Similar cases will be studied in sections 5.3.1 and 5.3.2. Another example is represented by parametric fluid-structure interaction systems in which the parameters cause complex interdependencies between the structure and fluid components in favor of partitioned linear solution manifold approximations (SVD is performed separately for the fluid, for the structure and for the interface) rather than monolithic ones. In our implementation of DD-ROMs, we exploit the partitions obtained from the distributed memory solver in deal.II. Since these domain decompositions typically satisfy constraints related to the computational efficiency, we devise some strategies to repartition the domain responding to solution manifold approximability concerns instead. Another work that implements this is [87] where the Reynolds stress tensor is employed, among others, as indicator for partitioning the computational domain. Similarly, we develop new indicators.

Another way to approach the problem of a slow KnW decay is exploiting the mathematical proofs of existence of vanishing viscosity solutions [65, 57, 25, 38]. In fact, solutions of hyperbolic problems can be obtained as a limit process of solutions associated to viscosity terms approaching zero. The crucial point is that ROMs associated to larger viscosity values may not suffer from a slow Kolmogorov \(n\)-width decay. Hence, we can set up classical projection based ROMs for the high viscosity solutions, and use graph neural networks (GNNs) [81] only to infer the vanishing viscosity solution in a very efficient manner. This procedure can be applied also to more general hyperbolic problems, not necessarily FS. The key features of this new methodology are the following: the employment of computationally heavy graph neural networks is reduced to a minimum and, at the same time, interpretable certified projection ROMs are exploited as much as possible in their regime of accurate linear approximability. In fact, GNNs, used generally to perform non-intrusive MOR, have high training computational costs and they are employed mainly for small academic benchmarks in terms of number of degrees of freedom, up to now. We avoid these high computational efforts with our multi-fidelity formulation: the GNNs are employed only to infer the vanishing viscosity solutions from the previous higher viscosity level, not to approximate and perform dimension reduction of the entire solution manifold. The overhead is the collection of additional full-order snapshots corresponding to high viscosity values, but this can be performed on coarser meshes as it will be done in section 6. Moreover, the support of our GNNs is the DG discretization space, so, we can enrich the typical machine learning framework of GNNs with data structure and operators from numerical analysis. We validate the use of data augmentation with numerical filters (discretized Laplacian, gradients), as proposed in [81].

In brief, we summarize our contributions with the present work:

* structure-preserving model order reduction for Friedrichs' systems. We synthetically describe the realization of ROMs for FS and the definition of standard _a posteriori_ error estimators. Hints towards the implementation of optimally stable ROMs are highlighted.
* domain decomposable reduced-order models for full-order models discretized with the discontinuous Galerkin method. We introduce DD-ROMs for DG discretizations and introduce novel indicators to repartition the computational domain with the aim of obtaining more efficient local solution manifold approximants.
* surrogate modelling of vanishing viscosity solutions with graph neural networks. We propose a new framework for the MOR of parametric hyperbolic PDEs with a slow Kolmogorov n-width decay.

The topics addressed in this work are presented as follows. In Section 2, we introduce the definition of FS and well-posedness results and we will provide several examples of models that fall into this framework: the Maxwell equations in stationary regime, the equations of linear compressible elasticity and the advection diffusion reaction equations. Then, we provide a DG discretization of the FS following [24] with related error estimates in Section 3. In Section 4, we introduce the projection-based MOR technique and some error bounds that can be effectively used. In Section 5, we will discuss a new implementation of domain decomposable ROMs for FOMs discretized with the DGM and we will test the approach on three parametric models. In Section 6, we introduce the concept of vanishing viscosity solutions and how graph neural network are exploited to overcome the problem of a slow Kolmogorov \(n\)-width decay. We will provide some numerical tests to show the effectiveness of the proposed approach. Finally, in Section 7 we summarize our results and we suggest further directions of research.

## 2 Friedrichs' systems

In this section, we will provide a summary of FS theory: their definition, existence, uniqueness and well-posedness results, their weak and ultraweak forms and many PDEs which can be rewritten into FS. The following discussion collects many results from [34, 73, 47, 74, 46, 52, 79, 52, 29, 31, 30, 3, 12, 24], but we will follow the notation in [24]. Let us represent with \(d\) the ambient space dimension and with \(m\geq 1\) the number of equations of the FS. We consider a connected Lipschitz domain \(\Omega\subset\mathbb{R}^{d}\), with boundary \(\partial\Omega\) and outward unit normal \(\mathbf{n}:\partial\Omega\to\mathbb{R}^{d}\).

A FS is defined through \((d+1)\) matrix-valued fields \(A^{0},A^{1}\ldots,A^{d}\in[L^{\infty}(\Omega)]^{m\times m}\) and the following differential operators \(\mathcal{X},A,\tilde{A}:\Omega\to\mathbb{R}^{m\times m}\). We suppose that \(\mathcal{X}\in[L^{\infty}(\Omega)]^{m\times m}\) and define

\[\mathcal{X}=\sum_{k=1}^{d}\partial_{k}A^{k}\,\qquad A=A^{0}+\sum_{i=1}^{d}A^{ i}\partial_{i}\,\qquad\tilde{A}=\left(A^{0}\right)^{t}-\mathcal{X}-\sum_{i=1}^{d}A^{i} \partial_{i}\, \tag{1}\]

assuming that

\[A^{k}=(A^{k})^{T}\,\text{a.e. in }\Omega, \text{for }k=1,\ldots,d, \text{(symmetry property)} \tag{2a}\] \[A^{0}+(A^{0})^{T}-\mathcal{X}\text{ is u.p.d. \text{ a.e. in }\Omega, \text{(positivity property)} \tag{2b}\]

thus, the name **symmetric positive operators** or **Friedrichs operators**, which is used to refer to (\(A\), \(\tilde{A}\)). We recall that the operator in (2b) is uniformly positive definite (u.p.d) if and only if

\[\exists\mu>0:A^{0}+(A^{0})^{T}-\mathcal{X}>2\mu_{0}\mathbb{I}\quad\text{a.e. in }\Omega. \tag{3}\]

If this property is not satisfied, it can be sometimes recovered as shown in Appendix A. A weaker condition can be required for two-field systems [30].

The boundary conditions are expressed through two boundary operators \(\mathcal{D}:\partial\Omega\to\mathbb{R}^{m\times m}\) with

\[\mathcal{D}=\sum_{k=1}^{d}n_{k}A^{k},\qquad\text{a.e. in }\partial\Omega \tag{4}\]

and \(\mathcal{M}:\partial\Omega\to\mathbb{R}^{m\times m}\) satisfying the following **admissible boundary conditions**

\[\mathcal{M}\quad\text{is nonnegative \ a.e. on }\quad\partial\Omega, \text{(monotonicity property)} \tag{5a}\] \[\ker(\mathcal{D}-\mathcal{M})+\ker(\mathcal{D}+\mathcal{M})= \mathbb{R}^{m}\quad\text{a.e. on }\quad\partial\Omega. \text{(strict adjointness property)} \tag{5b}\]

_Remark 1_ (Strict adjointness).: The term strict adjointness property comes from Jensen [52, Theorem 31]. The strict adjointness property is needed for the solution of the ultra-weak formulation of the FS to uniquely satisfy the boundary conditions: in a slightly different framework from the one presented here, see [52, Theorem 29] and [12, proof of Lemma 2.4].
**Theorem 1** (Friedrichs' system strong solution [34]).: _Let \(f\in[L^{2}(\Omega)]^{m}\), the strong solution \(z\in[C^{1}(\overline{\Omega})]^{m}\) to Friedrichs' system_

\[\begin{cases}Az=f,&\text{in }\Omega,\\ (\mathcal{D}-\mathcal{M})z=0,&\text{on }\partial\Omega.\end{cases} \tag{6}\]

_is unique. Moreover, there exists a solution of the ultra-weak formulation_

\[(z,\tilde{A}y)_{L^{2}}=(f,y)_{L^{2}},\qquad\forall y\in[C^{1}(\overline{\Omega })]^{m}\ s.t.\ (\mathcal{D}+\mathcal{M}^{t})y=0. \tag{7}\]

Let \(L=[L^{2}(\Omega)]^{m}\). We define the weak formulation on the graph space \(V=\{z\in L:Az\in L\}\), which amounts to differentiability in the characteristics directions: \(A\in\mathcal{L}(V,L)\) and \(\tilde{A}\in\mathcal{L}(V^{\prime},L)\). The boundary operator \(\mathcal{D}\) is translated into the abstract operator \(D\in\mathcal{L}(V,V^{\prime})\):

\[\left\langle Dz,y\right\rangle_{V,V^{\prime}}=(Az,y)_{L}-(z,\tilde{A}y)_{L}, \quad\forall z,y\in V. \tag{8}\]

When \(z\) is smooth, it can be seen as the integration by parts formula [52, 12]:

\[\left\langle Dz,y\right\rangle_{V,V^{\prime}}=\left\langle\mathcal{D}z,y \right\rangle_{H^{\frac{1}{2}}(\partial\Omega),H^{-\frac{1}{2}}(\partial \Omega)},\quad\forall z\in H^{1}(\Omega),\ y\in H^{1}(\Omega). \tag{9}\]

A sufficient condition for well-posedness of the weak formulation is provided by the cone formalism [3, 24] that poses the existence of two linear subdomains \((V_{0},V_{0}^{*})\) of \(V\):

\[V_{0}\text{ maximal in }C^{+},\quad V_{0}^{*}\text{ maximal in }C^{-} \tag{10a}\] \[V_{0}=D(V_{0}^{*})^{\perp},\quad V_{0}^{*}=D(V_{0})^{\perp}, \tag{10b}\]

such that \(A:V_{0}\to L\) and \(\tilde{A}:V_{0}^{*}\to L\) are isomorphism, where \(C^{\pm}=\{w\in V|\pm\left\langle Dw,w\right\rangle_{V,V^{\prime}}\geq 0\}\).

Provided that \(V_{0}+V_{0}^{*}\subset V\) is closed [3], the conditions in (10) are equivalent to the existence of the boundary operator \(M\in\mathcal{L}(V,V^{\prime})\) that satisfies **admissible boundary conditions** analogue to the ones in (5):

\[M\quad\text{is monotone}, \text{(monotonicity property)} \tag{11a}\] \[\ker(D-M)+\ker(D+M)=V, \text{(strict adjointness property)} \tag{11b}\]

identifying \(V_{0}=\ker(D-M)\) and \(V_{0}^{*}=\ker(D+M^{*})\).

**Theorem 2** (Friedrichs' System weak form [28, 29, 24, 12]).: _Let us assume that the boundary operator \(M\in\mathcal{L}(V,V^{\prime})\) satisfies the monotonicity and strict adjointness properties (11). Let us define for \(z,z^{*}\in V\) the bilinear forms_

\[a(z,y) =(Az,y)_{L}+\tfrac{1}{2}((D-M)z,y)_{V^{\prime},V},\quad\forall y \in V, \tag{12a}\] \[a^{*}(z^{*},y) =\left(\tilde{Az}^{*},y\right)_{L}+\tfrac{1}{2}((D+M^{*})z^{*},y )_{V^{\prime},V},\quad\forall y\in V. \tag{12b}\]

_Then, Friedrichs' operators \(A:V_{0}\to L\) and \(\tilde{A}:V_{0}^{*}\to L\) are isomorphisms: for all \(f\in L\) and \(g\in V\) there exists unique \(z,z^{*}\in V\) s.t._

\[a(z,y) =(f,y)_{L}+\langle(D-M)g,y\rangle_{V^{\prime},V}\quad\forall y \in V, \tag{13a}\] \[a^{*}(z^{*},y) =(f,y)_{L}+\langle(D+M^{*})g,y\rangle_{V^{\prime},V}\quad\forall y \in V, \tag{13b}\]

_that is_

\[\begin{cases}Az=f,&\text{in }L,\\ (M-D)(z-g)=0&\text{in }V^{\prime},\end{cases}\quad\begin{cases}\tilde{Az}^{*}=f,& \text{in }L,\\ (M^{*}+D)(z^{*}-g)=0&\text{in }V^{\prime}.\end{cases} \tag{14}\]

### A unifying framework

The theory of Friedrichs' systems provides a unified framework to study different classes of PDEs [52]: first-order uniformly hyperbolic, second-order uniformly hyperbolic, elliptic and parabolic partial differential equations. Originally, Friedrichs' aim was to study equations of mixed type (hyperbolic, parabolic, elliptic) inside the same domain such as the Tricomi equation [34] (or more generally the Frankl equation [52]) inspired by models from compressible gas dynamics for which the domain is subdivided in a hyperbolic supersonic and an elliptic subsonic part.

Some examples of FS can be found in the literature:

\[(x_{2}\partial_{1}^{2}\bullet+\partial_{1}^{2}\bullet)u=0,\quad\text{(Tricomi \@@cite[cite]{[\@@bibref{}{Tricomi}{}{}]})} \tag{15a}\]
\[\begin{bmatrix}-\partial_{1}\bullet&\partial_{2}\bullet\\ -\partial_{2}\bullet&\partial_{1}\bullet\end{bmatrix}\begin{pmatrix}u_{1}\\ u_{2}\end{pmatrix} =0,&\text{(Cauchy-Riemann \@@cite[cite]{[\@@bibref{}{Cauchy-Riemann \@@citephrase{}{}{}{}]}})}\\ (A(x_{2})\partial_{1}^{2}+\partial_{1}^{2})u =0,&\text{(Frankl \@@cite[cite]{[\@@bibref{}{Frankl \@@citephrase{}{}{}]}})}\end{bmatrix}\] (15c) \[\begin{bmatrix}\mathbb{I}_{3}&-\lambda^{-1}\mathbb{I}_{3}( \nabla\cdot\bullet)-\frac{(\nabla\bullet+(\nabla\bullet)^{t})}{2}\\ -\frac{1}{2}\nabla\cdot(\bullet+\bullet^{t})&\alpha\mathbb{I}_{3}\end{bmatrix} \begin{pmatrix}\sigma\\ \mathbf{u}\end{bmatrix} =0,&\text{(Compressible linear elasticity \@@cite[cite]{[\@@bibref{}{Compressible-linear elasticity \@@citephrase{}{}{}]}})}\end{bmatrix}\] (15d) \[\begin{bmatrix}\mu\mathbb{I}_{3}&\nabla\times\bullet\\ -\nabla\times\bullet&\sigma\mathbb{I}_{3}\end{bmatrix}\begin{pmatrix}\mathbf{ H}\\ \mathbf{E}\end{bmatrix} =0,&\text{(Maxwell eq. in stationary regime \@@cite[cite]{[\@@bibref{}{Compressible-linear elasticity \@@citephrase{}{}{}]}})}\] (15e) \[(-\nabla\cdot(\kappa\nabla\bullet)+\beta\cdot\nabla\cdot\epsilon+ \mu\bullet)\mathbf{u} =0,&\text{(Diffusion advection reaction \@@cite[cite]{[\@@bibref{}{Compressible-linear elasticity \@@citephrase{}{}{}]}})}\] (15f) \[(A_{0}\partial_{t}\bullet+\Sigma_{i=1}^{3}\tilde{A}_{i}\partial_{i }\bullet)\mathbf{V} =0,\text{(Linearized symmetric Euler \@@cite[cite]{[\@@bibref{}{Cauchy-Riemann \@@citephrase{}{}{}]}},\text{\@@cite[cite]{[\@@bibref{}{Compressible-linear elasticity \@@citephrase{}{}{}]}})}\] (15g) \[(a\gamma^{0}\partial_{t}\bullet+\gamma^{1}\partial_{1}\bullet+ \gamma^{2}\partial_{2}\bullet+\gamma^{3}\partial_{3}\bullet+B)\mathbf{\psi} =0,\text{(Dirac system \@@cite[cite]{[\@@bibref{}{Dirac-Klein-Gordon system \@@cite[cite]{[\@@bibref{}{Dirac-Klein-Gordon system \@@cite[cite]{[\@@bibref{}{Dirac-Klein-Gordon system \@@cite[cite]{[\@@bibref{Dirac-Klein-Gordon system \
with \(\mathcal{R}^{k}_{ij}=\epsilon_{ikj}\) being the Levi-Civita tensor. The graph space is \(V=H(\mathrm{curl},\Omega)\times H(\mathrm{curl},\Omega)\). The boundary operator is

\[\mathcal{D}=\sum_{k=1}^{d}n_{k}A^{k} =\begin{bmatrix}0_{d,d}&\mathcal{T}\\ \mathcal{T}^{T}&0_{d,d}\end{bmatrix},\qquad\text{with }\mathcal{T}\xi:=\mathbf{n} \times\mathbf{\xi}, \tag{16}\] \[\left\langle D(\mathbf{H},\mathbf{E}),(\mathbf{h},\mathbf{e}) \right\rangle_{V^{{}^{\prime}},V} =(\mathbf{n}\times\mathbf{E},\mathbf{e})_{L^{2}(\partial\Omega)}-( \mathbf{n}\times\mathbf{H},\mathbf{h})_{L^{2}(\partial\Omega)}. \tag{17}\]

We impose homogeneous Dirichlet boundary conditions tangential to the electric field \((\mathbf{n}\times\mathbf{E})_{|\partial\Omega}=0\) through

\[\mathcal{M}=\begin{bmatrix}0_{d,d}&-\mathcal{T}\\ \mathcal{T}^{T}&0_{d,d}\end{bmatrix},\qquad\left\langle M(\mathbf{H},\mathbf{ E}),(\mathbf{h},\mathbf{e})\right\rangle_{V^{{}^{\prime}},V}=-(\mathbf{n}\times \mathbf{E},\mathbf{e})_{L^{2}(\partial\Omega)}-(\mathbf{n}\times\mathbf{H}, \mathbf{h})_{L^{2}(\partial\Omega)}. \tag{18}\]

#### 2.1.2 Compressible linear elasticity

We consider the parametric compressible linear elasticity system in \(\mathbb{R}^{d}=\mathbb{R}^{3}\), where \(\mathbf{\sigma}\in\mathbb{R}^{d\times d}\) is the stress tensor and \(\mathbf{u}\in\mathbb{R}^{d}\) is the displacement vector. The system can be written as

\[\left(\begin{array}{c}\mathbf{\sigma}-\mu_{1}(\nabla\cdot\mathbf{u})\mathbb{I}_ {3,3}-2\mu_{2}\frac{(\nabla\mathbf{u}+(\nabla\mathbf{u})^{t})}{2}\\ -\frac{1}{2}\nabla\cdot(\mathbf{\sigma}+\mathbf{\sigma}^{t})+\mu_{3}\mathbf{u}\end{array} \right)=\left(\begin{array}{c}0\\ \mathbf{r}\end{array}\right),\quad\forall x\in\Omega, \tag{19}\]

where \(\mathbf{r}\in\mathbb{R}^{3}\), and \(\mu_{1},\ \mu_{2}>0\) are the Lame constants. Rescaling the displacement \(\mathbf{u}\) by \(2\mu_{2}\), we obtain

\[\left(\begin{array}{c}\mathbf{\sigma}-\frac{\mu_{1}}{2\mu_{2}+3\mu_{1}}\mathrm{ tr}(\mathbf{\sigma})\mathbb{I}_{3,3}-\frac{(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})}{2} \\ -\frac{1}{2}\nabla\cdot\left(\mathbf{\sigma}+\mathbf{\sigma}^{T}\right)+\frac{\mu_{3}}{2 \mu_{2}}\mathbf{u}\end{array}\right)=\left(\begin{array}{c}\mathbf{0}\\ \mathbf{r}\end{array}\right),\quad\forall x\in\Omega. \tag{20}\]

In this case, we consider the graph space

\[V=H_{\mathbf{\sigma}}\times[H^{1}(\Omega)]^{d},\quad H_{\mathbf{\sigma}}=\{\mathbf{\sigma} \in[L^{2}(\Omega)]^{d\times d}\,|\,\nabla\cdot(\mathbf{\sigma}+\mathbf{\sigma}^{t})\in [L^{2}(\Omega)]^{d}\}. \tag{21}\]

If we reorder the coefficients of \(\mathbf{\sigma}\) into a vector, we can define \(\mathbf{z}=\begin{pmatrix}\mathbf{\sigma}\\ \mathbf{u}\end{pmatrix}\) and have

\[A^{0}=\begin{bmatrix}\mathbb{I}_{d^{2},d^{2}}-\frac{\mu_{1}}{2\mu_{2}+3\mu_{1} }\mathcal{Z}&0_{d^{2},d}\\ 0_{d,d^{2}}&\frac{\mu_{3}}{2\mu_{2}}\mathbb{I}_{d,d}\end{bmatrix},\qquad A^{k} =\begin{bmatrix}0_{d^{2},d^{2}}&\mathcal{E}^{k}\\ (\mathcal{E}^{k})^{T}&0_{d,d}\end{bmatrix},\qquad f=\begin{bmatrix}\mathbf{0}_ {d^{2}}\\ \mathbf{r}\end{bmatrix}, \tag{22}\]

with \(\mathcal{Z}_{[ij],[kl]}=\delta_{ij}\delta_{kl}\) and \(\mathcal{E}^{k}_{[ij],l}=-\frac{1}{2}\left(\delta_{ik}\delta_{jl}+\delta_{il} \delta_{jk}\right)\). This leads to the definition of the boundary operator

\[\mathcal{D}=\sum_{k=1}^{d}n_{k}A^{k}=\begin{bmatrix}0_{d^{2},d^{2 }}&\mathcal{N}\\ \mathcal{N}^{T}&0_{d,d}\end{bmatrix}\text{ with }\mathcal{N}\mathbf{\xi}:=-\frac{1}{2}( \mathbf{n}\otimes\mathbf{\xi}+\mathbf{\xi}\otimes\mathbf{n}), \tag{23a}\] \[\left\langle D(\mathbf{\sigma},\mathbf{u}),(\mathbf{\tau},\mathbf{v}) \right\rangle_{V^{{}^{\prime}},V}=-\langle\tfrac{1}{2}(\mathbf{\sigma}+\mathbf{\sigma}^ {t})\cdot\mathbf{n},\mathbf{v}\rangle_{-\frac{1}{2},\frac{1}{2}}-\langle \tfrac{1}{2}(\mathbf{\tau}+\mathbf{\tau}^{t})\cdot\mathbf{n},\mathbf{u}\rangle_{-\frac {1}{2},\frac{1}{2}}. \tag{23b}\]

Mixed boundary conditions \(\mathbf{u}_{|\Gamma_{D}}=0\) and \((\mathbf{\sigma}\cdot\mathbf{n})_{|\Gamma_{N}}=0\) can be applied through the following boundary operator on the Dirichlet boundary \(\Gamma_{D}\) and on the Neumann boundary \(\Gamma_{N}\):

\[\begin{split}\left\langle M(\mathbf{\sigma},\mathbf{u}),(\mathbf{\tau}, \mathbf{v})\right\rangle_{V^{{}^{\prime}},V}=&-\langle\tfrac{1}{2}(\mathbf{ \sigma}+\mathbf{\sigma}^{t})\cdot\mathbf{n},\mathbf{v}\rangle_{-\frac{1}{2}, \frac{1}{2},\Gamma_{D}}+\langle\tfrac{1}{2}(\mathbf{\tau}+\mathbf{\tau}^{t})\cdot \mathbf{n},\mathbf{u}\rangle_{-\frac{1}{2},\frac{1}{2},\Gamma_{D}}\\ &+\langle\tfrac{1}{2}(\mathbf{\sigma}+\mathbf{\sigma}^{t})\cdot\mathbf{n}, \mathbf{v}\rangle_{-\frac{1}{2},\frac{1}{2},\Gamma_{N}}-\langle\tfrac{1}{2}(\bm {\tau}+\mathbf{\tau}^{t})\cdot\mathbf{n},\mathbf{u}\rangle_{-\frac{1}{2},\frac{1}{2}, \Gamma_{N}},\end{split} \tag{24}\]

the constructive procedure employed to define the boundary operator \(M\in\mathcal{L}(V,V^{{}^{\prime}})\) is reported in the Appendix B.

#### 2.1.3 Grad-div problem: advection-diffusion-reaction equations

Another example is the advection-diffusion-reaction equation

\[-\nabla\cdot(\kappa\nabla u)+\mathbf{\beta}\cdot\nabla u+\mu u=\mathbf{r}, \tag{25}\]

with \(\kappa\in[L^{\infty}(\Omega)]^{d\times d}\) and \(\mathbf{\beta}\in[W^{1,\infty}(\Omega)]^{d}\), under the hypothesis that \(\kappa\in[L^{\infty}(\Omega)]^{d\times d}\) and \(\mu-\nabla\cdot\mathbf{\beta}\in L^{\infty}(\Omega)\) are uniformly bounded from below to satisfy the positivity property (2b). Let us write the equation in the mixed form with \(\mathbf{\sigma}=-\kappa\nabla u\) and \(\mathbf{z}=\begin{pmatrix}\mathbf{\sigma}\\ u\end{pmatrix}\). Then, (25) can be rewritten as (6) with

\[A^{0}=\begin{bmatrix}\kappa^{-1}&0_{d,1}\\ 0_{1,d}&\mu\end{bmatrix},\qquad A^{k}=\begin{bmatrix}0_{d,d}&\mathbf{e}_{k}\\ (\mathbf{e}_{k})^{T}&\beta_{k}\end{bmatrix},\qquad\mathbf{f}=\begin{pmatrix}0\\ \mathbf{r}\end{pmatrix}. \tag{26}\]
Here, \(0_{m,\ell}\in\mathbb{R}^{m\times\ell}\) is a matrix of zeros and \(\mathbf{e}_{k}\) is the unitary vector with the \(k\)-th entry equal to \(1\). The graph space is \(V=H(\mathrm{div},\Omega)\times H^{2}(\Omega)\). The boundary operator \(D\) becomes

\[D=\sum_{k=1}^{d}n_{k}A^{k}=\begin{bmatrix}0_{d,d}&\mathbf{n}\\ \mathbf{n}^{t}&\beta\cdot\mathbf{n}\end{bmatrix},\qquad\langle D(\boldsymbol{ \sigma},u),(\boldsymbol{\tau},v)\rangle_{V^{\prime},V}=\langle\boldsymbol{ \sigma}\cdot\mathbf{n},v\rangle_{-\frac{1}{2},\frac{1}{2}}-\langle\boldsymbol{ \tau}\cdot\mathbf{n},u\rangle_{-\frac{1}{2},\frac{1}{2}}. \tag{27}\]

Homogeneous Dirichlet boundary conditions \(u_{|\partial\Omega}=0\) can be imposed with

\[\mathcal{M}=\begin{bmatrix}0_{d,d}&-\mathbf{n}\\ \mathbf{n}^{t}&0\end{bmatrix}, \tag{28}\]

while Robin/Neumann boundary conditions of the type \(\boldsymbol{\sigma}\cdot\mathbf{n}=\gamma u\) are imposed with

\[\mathcal{M}=\begin{bmatrix}0_{d,d}&\mathbf{n}\\ -\mathbf{n}^{t}&2\gamma+\boldsymbol{\beta}\cdot\mathbf{n}\end{bmatrix}. \tag{29}\]

For our test case in Section 6, we will consider as advection field \(\boldsymbol{\beta}:\Omega\to\mathbb{R}^{d}\) an incompressible velocity field from the solution of the 2d incompressible Navier-Stokes equations as described later. Similarly to the linear compressible elasticity mixed boundary conditions in Section 2.1.2, we want to impose \(u_{|\Gamma_{D}}=\mathbf{g}\in[L^{\frac{1}{2}}(\Gamma_{D})]^{d}\) and \(\left(\boldsymbol{\sigma}\cdot\mathbf{n}\right)_{|\Gamma_{N}}=0\). This is possible with

\[\langle M(\boldsymbol{\sigma},u),(\boldsymbol{\tau},v)\rangle_{V^{\prime},V}= \langle\boldsymbol{\sigma}\cdot\mathbf{n},v\rangle_{-\frac{1}{2},\frac{1}{2}, \Gamma_{D}}+\langle\boldsymbol{\tau}\cdot\mathbf{n},u\rangle_{-\frac{1}{2}, \frac{1}{2},\Gamma_{D}}-\langle\boldsymbol{\sigma}\cdot\mathbf{n},v\rangle_{- \frac{1}{2},\frac{1}{2},\Gamma_{N}}-\langle\boldsymbol{\tau}\cdot\mathbf{n},u \rangle_{-\frac{1}{2},\frac{1}{2},\Gamma_{N}}, \tag{30}\]

the proof is similar to the one reported in Appendix B.

## 3 Discontinuous Galerkin discretization

In the literature, a few discretization approaches for FS are presented, e.g. finite volume method [79] or discontinuous Galerkin (DG) method [28, 24, 12]. More recently, a hp-adaptive hybridizable DG formulation was introduced in [17]. In this work, we perform a DG discretization following the notation reported in [28, 24]. Consider a shape-regular tessellation \(\mathcal{T}_{h}\) of the domain \(\Omega\) and take a piecewise polynomial space \(V_{h}\) over \(\mathcal{T}_{h}\), defined by \(V_{h}=\{z\in V:z|_{T}\in\mathbb{P}_{h}^{d}(T),\forall T\in\mathcal{T}_{h}\}\), where \(k\) is the polynomial degree. We assume that there is a partition \(P_{\Omega}=\{\Omega_{i}\}_{1\leq i\leq N_{\Omega}}\) of \(\Omega\) into disjoint polyhedra such that the exact solution \(z\) belongs to \(V^{*}=V\cap[H^{1}(P_{\Omega})]^{m}\). We define the discrete bilinear form \(\forall y_{h}\in V_{h},\ z\in V^{*}\)

\[a_{h}^{cf}(z,y_{h}) =\sum_{T\in\mathcal{T}_{h}}(Az,y_{h})_{L^{2}(T)}+\tfrac{1}{2}\sum _{F\in\mathcal{F}_{h}^{b}}\left((\mathcal{M}-\mathcal{D})z,y_{h}\right)_{L^{2} (F)}-\sum_{F\in\mathcal{F}_{h}^{i}}\left(\mathcal{D}_{F}[\![z]\!],\{\!\{y_{h}\} \!\}\right)_{L^{2}(F)} \tag{31a}\] \[=\sum_{T\in\mathcal{T}_{h}}(z,\tilde{A}y_{h})_{L^{2}(T)}+\tfrac{1}{ 2}\sum_{F\in\mathcal{F}_{h}^{b}}\left((\mathcal{M}+\mathcal{D})z,y_{h}\right)_{L ^{2}(F)}+\sum_{F\in\mathcal{F}_{h}^{i}}\left(\mathcal{D}_{F}\{\!\{z\}\!\},[\![y _{h}]\!]_{L^{2}(F)}, \tag{31b}\]

where the first two terms are the piece-wise discontinuous discretization of the bilinear form (12a) and the last term penalizes the jump across neighboring cells and stabilizes the method. Here, \(\mathcal{F}_{h}^{k}\) is the collection of the faces of the triangulation \(\mathcal{T}_{h}\) belonging to the boundary of \(\Omega_{h}\), while \(\mathcal{F}_{h}^{i}\) is the collection of internal faces. The jump and the average of a function on a face \(F\) shared by two elements \(T_{1}\) and \(T_{2}\) are defined as \([\![u]\!]=u|_{T_{1}}-u|_{T_{2}}\) and \(\{[\![u]\!]=\tfrac{1}{2}(u|_{T_{1}}+u|_{T_{2}})\), respectively. The boundary operator \(\mathcal{D}:\partial\Omega\to\mathbb{R}^{m\times m}\) can be extended also on the internal faces \(F\in\mathcal{F}_{h}^{i}\) as \(\mathcal{D}_{F}=\sum_{k=1}^{d}n_{k}^{F}A^{k}\), where \(n^{F}\) is a normal to the face \(F\) and it is well-defined.

In order to obtain quasi-optimal error estimates, extra stabilization terms are needed. We additionally impose that \(A^{i}\in[C^{0,\frac{1}{2}}(\overline{\Omega}_{j})]^{m\times m},\ \forall T\in \mathcal{T}_{h},\ i=1,\ldots,d,\ \forall\Omega_{j}\in P_{\Omega}\). A possibility is given by the following stabilization term

\[s_{h}(z,y_{h})=\sum_{F\in\mathcal{F}_{h}^{k}}(S_{F}^{b}z,y_{h})_{L^{2}(F)}+\sum_ {F\in\mathcal{F}_{h}^{i}}(S_{h}^{i}[\![z]\!],[\![y_{h}]\!])_{L^{2}(F)}, \tag{32}\]

where the operators \(S_{h}^{i}\) and \(S_{F}^{b}\) have to satisfy the following constraints for some \(\alpha_{j}>0\) for \(j=1,\ldots,5\):

\[S_{F}^{b}z=0\quad\forall F\in\mathcal{F}_{h}^{b},\qquad S_{F}^{i}[ \![z]\!]=0\quad\forall F\in\mathcal{F}_{h}^{i},\qquad\text{with $z$ the exact solution,} \tag{33a}\] \[S_{F}^{b}\text{ and }S_{F}^{i}\text{ are symmetric and nonnegative,}\] (33b) \[S_{F}^{b}\leq\alpha_{1}\mathbb{I}_{m,m},\qquad\alpha_{2}|D_{F}| \leq S_{F}^{i}\leq\alpha_{3}\mathbb{I}_{m,m},\] (33c) \[|((M-D)y,z)_{L^{2}(F)}|\leq\alpha_{4}((S_{F}^{b}+M)y,y)_{L^{2}(F)}^ {(12)}\|z\|_{L^{2}(F)},\] (33d) \[|((M+D)y,z)_{L^{2}(F)}|\leq\alpha_{5}((S_{F}^{b}+M)z,z)_{L^{2}(F)}^ {1/2}\|y\|_{L^{2}(F)}. \tag{33e}\]
Specific definitions of these operators for our test cases are presented in [28, 24], properly declined for our mixed boundary conditions in the compressible linear elasticity and advection-diffusion-reaction test cases, see Sections 2.1.2 and 2.1.3, respectively. Finally, we can define the bilinear form and the right-hand side

\[a_{h}(z,y_{h})=a_{h}^{cf}(z,y_{h})+s_{h}(z,y_{h}),\quad l_{h}(y_{h})=\sum_{T\in \mathcal{T}_{h}}(f,y_{h})_{L^{2}(T)}+\tfrac{1}{2}\sum_{F\in\mathcal{F}_{h}^{b} }\left((M-D)g,y_{h}\right)_{L^{2}(F)} \tag{34}\]

that lead to the definition of the discrete problem.

**Definition 1** (DG Friedrichs' System).: _Given \(f\in L\) and \(g\in V_{h}\), the DG approximation of the FS constitute in finding a \(z_{h}\in V_{h}\) such that_

\[a_{h}(z_{h},y_{h})=l_{h}(y_{h}),\qquad\forall y_{h}\in V_{h}. \tag{35}\]

To prove the accuracy of the discrete problem, it is necessary to have the following conditions:

* Consistency, i.e., \(a_{h}(z,y_{h})=a(z,y_{h})\) for \(z\in V^{*}\);
* \(L^{2}\)-coercivity, i.e., \(a_{h}(y_{h},y_{h})\geq\rho_{0}\|y_{h}\|_{L}^{2}+\tfrac{1}{2}|y_{h}|_{M}^{2}\), with \(|y_{h}|_{M}^{2}=\int_{\partial\Omega}y^{t}My\);
* Inf-sup stability \[|||z_{h}|||\lesssim\sup_{y_{h}\neq 0}\frac{a_{h}(z_{h},y_{h})}{|||y_{h}|||}\] (36) with \(|||y|||^{2}=||y||_{L^{2}}^{2}+|y|_{M}^{2}+|y|_{S}^{2}+\sum_{T\in\mathcal{T}_ {h}}h_{T}\left\|A^{b}\partial_{k}y\right\|_{L^{2}(T)}^{2}\) and \(|y|_{S}^{2}=s_{h}(y,y)\);
* Boundedness \(a_{h}(w,y_{h})\lesssim|||w|||_{*}|||y_{h}|||\) with \[|||y|||_{*}^{2}=|||y|||^{2}+\sum_{T\in\mathcal{T}_{h}}\left(h_{T}^{-1}||y||_{ L^{2}(T)}^{2}+||y||_{L^{2}(\partial T)}^{2}\right).\] (37)

**Theorem 3** (Error estimate from [28, 24]).: _Let \(z\in V^{*}\) be the solution of the weak problem (13a) and \(z\in V_{h}\) be the solution of the discrete DG problem (35). Then, the consistency and inf-sup stability of the discrete system (35) imply_

\[|||z-z_{h}|||\lesssim\inf_{y_{h}\in V_{h}}|||z-y_{h}|||_{*}, \tag{38}\]

_in particular, if \(z\in[H^{k+1}(\Omega)]^{m}\) the following convergence rate holds_

\[|||z-z_{h}|||\lesssim h^{k+\tfrac{1}{2}}\|z\|_{[H^{k+1}(\Omega)]^{m}}. \tag{39}\]

## 4 Projection-based model order reduction

The computation of discrete solutions of parametrized PDEs can require a not negligible computational time. In particular, in multi-query context, when many evaluations for different parameters are required, the computations may become unbearable. In this section, we introduce a reduced order model (ROM) for the FS in case of parameter dependent problems [44, 77], in order to drastically reduce the computational costs. To do so, we exploit two aspects of the above presented FS: the linearity of the problems and the affine dependence of the operators on the physical parameters.

As we have seen in Section 2.1, all the problems are depending on some parameters \(\boldsymbol{\rho}\in\mathcal{P}\subset\mathbb{R}^{N_{\text{par}}}\) and the dependence is affine. This means that it is possible to find \(N_{\text{aff}}\) terms independent on the parameters for each form, such that they can be affinely combined with some parameter dependent functions to obtain the original operator, i.e.,

\[a_{h}(z,y_{h};\boldsymbol{\rho})=\sum_{\ell=1}^{N_{\text{aff}}}\theta_{\ell}^{ \varepsilon}(\boldsymbol{\rho})a_{\ell,h}(z,y_{h}),\qquad l_{h}(y_{h})=\sum_{ \ell=1}^{N_{\text{aff}}}\theta_{\ell}^{f}(\boldsymbol{\rho})l_{h,l}(y_{h}). \tag{40}\]

Then, we select a reduced space \(V_{r}\subset V_{h}\) provided by a compression algorithm, e.g. SVD/POD/PCA [53, 58, 83] or Greedy algorithm [70, 69, 44, 20]. We suppose that the reduced dimension \(r\) is much smaller than the dimension \(N_{h}\) of the full order model space \(V_{h}\). We look as ansatz for a reduced solution \(z_{\text{RB}}\in V_{r}\) a linear combination of the bases \(\{\psi_{j}^{\text{RB}}\}_{j=1}^{r}\) of \(V_{r}\), i.e.,

\[z_{\text{RB}}=\sum_{j=1}^{r}z_{\text{RB}}^{j}\psi_{j}^{\text{RB}}, \tag{41}\]

then, performing a standard Galerkin projection, we obtain the following RB problem.
**Definition 2** (Reduced Basis Problem).: _Find \(z_{RB}\in V_{r}\), given by the coefficients \(z_{RB}^{j}\), such that_

\[\sum_{i=1}^{r}z_{RB}^{j}(\boldsymbol{\rho})\sum_{\ell=1}^{N_{\text{ pdf}}}\theta_{\ell}^{\text{g}}(\boldsymbol{\rho})a_{\ell,h}(\psi_{j}^{RB},\psi_{i}^{ RB})=\sum_{\ell=1}^{N_{\text{ pdf}}}\theta_{\ell}^{f}(\boldsymbol{\rho})(f_{\ell},\psi_{i}^{RB}),\qquad\text{ for all }i=1,\ldots,r. \tag{42}\]

The obtained problem scales depend on the dimension \(r\) and \(N_{\text{aff}}\) in its assembly and only on \(r\) in its solution, and it is completely independent on \(N_{h}\). To obtain computational advantages for the parametric problem, we split the tasks into an expensive _offline phase_ and a cheap _online phase_. In the _offline phase_, we find the reduced space \(V_{r}\) and we assemble the reduced matrices and right hand sides

\[A_{\ell}:=\{a_{\ell,h}(\psi_{j}^{\text{RB}},\psi_{i}^{\text{RB}})\}_{i,j}, \qquad b_{\ell}:=\{(f_{\ell},\psi_{i}^{\text{RB}})\}_{i}. \tag{43}\]

In the _online phase_, we can simply evaluate the coefficients \(\theta_{\ell}^{\text{g}}(\boldsymbol{\rho})\) and \(\theta_{\ell}^{f}(\boldsymbol{\rho})\) and obtain the reduced linear system

\[A(\boldsymbol{\rho})z_{\text{RB}}=b(\boldsymbol{\rho}),\qquad\text{with }A( \boldsymbol{\rho}):=\sum_{\ell}\theta_{\ell}^{\text{g}}(\boldsymbol{\rho})A_{ \ell}\text{ and }b(\boldsymbol{\rho}):=\sum_{\ell}\theta_{\ell}^{f}( \boldsymbol{\rho})b_{\ell}. \tag{44}\]

This gives a great speed up in computational times.

### Reduced basis a posteriori error estimate

We derive two error estimators for the energy norm and the \(L^{2}\) norm of the reduced basis error \(e_{h}=z_{h}-z_{RB}\in V_{h}\) following the procedure in [44]. Exploiting the equality in (31), we obtain the following lower bound

\[a_{h}(e_{h},e_{h})= a_{h}^{cf}(e_{h},e_{h})+s_{h}(e_{h},e_{h}) \tag{45a}\] \[= \sum_{T\in\mathcal{T}_{h}}(Ae_{h},e_{h})_{L^{2}(T)}+\tfrac{1}{2} \sum_{F\in\mathcal{F}^{\text{k}}}\left((\mathcal{M}-\mathcal{D})e_{h},e_{h} \right)_{L^{2}(F)}-\sum_{F\in\mathcal{F}_{h}^{\text{k}}}(\mathcal{D}_{F}[\![ e_{h}]\!],\{\!(e_{h}\}\!])_{L^{2}(F)}+\] (45b) \[\sum_{F\in\mathcal{F}_{h}^{\text{k}}}(S_{F}^{b}e_{h},e_{h})_{L^{2 }(F)}+\sum_{F\in\mathcal{F}_{h}^{\text{k}}}(S_{h}^{i}[\![e_{h}]\!],[\![e_{h}] \!])_{L^{2}(F)}\] (45c) \[= \sum_{T\in\mathcal{T}_{h}}\left((A^{0}-\tfrac{1}{2}\mathcal{X})e_{ h},e_{h}\right)_{L^{2}(T)}+\tfrac{1}{2}\sum_{F\in\mathcal{F}^{\text{k}}} \left(\mathcal{M}e_{h},e_{h}\right)_{L^{2}(F)}+\sum_{F\in\mathcal{F}_{h}^{ \text{k}}}(S_{F}^{b}e_{h},e_{h})_{L^{2}(F)}+\sum_{F\in\mathcal{F}_{h}^{\text{ k}}}(S_{h}^{i}[\![e_{h}]\!],[\![e_{h}]\!])_{L^{2}(F)}\] (45d) \[\geq \mu_{0}\|e_{h}\|_{L}^{2}+\tfrac{1}{2}|e_{h}|_{M}^{2}+|e_{h}|_{S}^ {2}, \tag{45e}\]

where we have defined \(|\cdot|_{M}^{2}=\tfrac{1}{2}\sum_{F\in\mathcal{F}^{\text{k}}}\left(\mathcal{M} \cdot,\cdot\right)_{L^{2}(F)}\) and \(|\cdot|_{S}^{2}=s_{h}(\cdot,\cdot)\). We define the \(R\)-norm

\[||y_{h}||_{R}^{2}=\mu_{0}\|y_{h}\|_{L}^{2}+\tfrac{1}{2}|y_{h}|_{M}^{2}+|y_{h}|_ {S}^{2},\quad\forall y_{h}\in V_{h}, \tag{46}\]

that may depend on \(\rho\) only through \(\mu_{0}\) and is generated by the scalar product

\[\langle u_{h},v_{h}\rangle_{R}=\mu_{0}\sum_{T\in\mathcal{T}_{h}}(u_{h},v_{h})_ {L^{2}(T)}+\tfrac{1}{2}\sum_{F\in\mathcal{F}^{\text{k}}}\left(\mathcal{M}^{ \text{sym}}u_{h},v_{h}\right)_{L^{2}(F)}+\sum_{F^{\text{k}}}(S^{b}u_{h},v_{h})+ \sum_{F^{\text{i}}}(S^{i}u_{h},v_{h}). \tag{47}\]

The boundary operators we will employ in our benchmarks are all skew-symmetric so \(\mathcal{M}^{\text{sym}}=\frac{\mathcal{M}+\mathcal{M}^{\text{t}}}{2}\) is the null matrix and \(|e_{h}|_{M}=0\). Now, we can proceed to provide an _a posteriori_ error estimate for the \(R\)-norm and energy norm. Hence, let us define \(r_{RB}(y_{h})=l_{h}(y_{h};\rho)-a_{h}(z_{RB},y_{h};\rho)\) and its \(R\) and \(L\)-Riesz representations as \(\hat{r}_{R}\) and \(\hat{r}_{L}\) such that

\[r_{RB}(u_{h})=\langle\hat{r}_{R},u_{h}\rangle_{R},\quad X_{R}\hat{\mathbf{r}}_{ R}=\mathbf{L}_{h}-A_{h}\mathbf{z}_{RB},\quad r_{RB}(u_{h})=\langle\hat{r}_{L},u_{h} \rangle_{L},\quad X_{L}\hat{\mathbf{r}}_{L}=\mathbf{L}_{h}-A_{h}\mathbf{z}_{RB}, \tag{48}\]

where \(X_{R}\) and \(X_{L}\) are the \(R\)-norm and \(L\)-norm mass matrices, and \(L_{h}\), \(A_{h}\) and \(\mathbf{z}_{RB}\) are the representations of \(l_{h}(\cdot;\rho)\), \(a_{h}(\cdot,\cdot;\rho)\) and \(z_{RB}\) in the DG basis of \(V_{h}\). The \(\hat{r}_{L}\) representation can be computed cheaply when the parametric model is affinely decomposable with respect to the parameters, while \(\hat{r}_{R}\) requires the inversion of a possibly parametric dependent matrix \(X_{R}\).

Now, consider the energy norm of the error \(\|e_{h}\|_{nrg}^{2}=a_{h}(e_{h},e_{h})\) and the coercivity constant \(\|e_{h}\|_{nrg}^{2}\geq\mu_{0}\|e_{h}\|_{L}^{2}\) derived in (45e), we have the following _a posteriori_ error estimates

\[\frac{\|e_{h}\|_{nrg}}{\|z_{h}\|_{nrg}}\leq\frac{\|\hat{r}_{R}\|_{R}}{\|z_{h}\|_{ nrg}},\qquad\frac{\|e_{h}\|_{R}}{\|z_{h}\|_{R}}\leq\frac{\|\hat{r}_{R}\|_{R}}{\|z_{h}\|_{R}}, \qquad\frac{\|e_{h}\|_{nrg}}{\|z_{h}\|_{nrg}}\leq\frac{\|\hat{r}_{L}\|_{L}}{\|z_{h} \|_{nrg}},\qquad\frac{\|e_{h}\|_{lrg}}{\|z_{h}\|_{lrg}}\leq\frac{\|\hat{r}_{L}\|_{L}}{ \|z_{h}\|_{lrg}}, \tag{49}\]

namely, the relative energy error with the corresponding _a posteriori_\(R\)-norm energy estimate, the relative \(R\)-norm error with the corresponding _a posteriori_\(L\)-norm energy estimate, and the relative \(L\)-norm error with the corresponding _a posteriori_\(L\)-norm
### Optimally stable error estimates for the ultraweak Petrov-Galerkin formulation

In this section, we show that Friedrichs' systems are a desirable unifying formulation to consider when performing model order reduction also due to the possibility to achieve an optimally stable formulation. This can further simplify the error estimator analysis reaching the equality between the error and the residual norm. This is not the first case in which optimally stable formulations are introduced also at the reduced level, see [11, 39, 42]. In the following, we describe how to achieve this ultraweak formulation and we delineate the path one should follow to use such formulation. Nevertheless, we will not use this formulation in our numerical tests and we leave the implementation to future studies.

We introduce the following Discontinuous Petrov-Galerkin (DPG) formulation from [12]. To do so, we first define \(V(\mathcal{T}_{h})\) the broken graph space with norm \(\|\bullet\|_{V(\mathcal{T}_{h})}^{2}=\|\bullet\|_{L}^{2}+\sum_{T\in\mathcal{T }_{h}}\|A\bullet\|_{L^{2}(T)}^{2}\) and \(\tilde{V}=V/Q(\Omega)\) the quotient of the graph space \(V\) with

\[\begin{split} Q(\Omega)=&\left\{z\in V\,\bigg{|} \sum_{T\in\mathcal{T}_{h}}\langle Dz,y\rangle_{\tilde{V}(T),V(T)}+\tfrac{1}{2} ((M-D)z,y\rangle_{\tilde{V}(\Omega),V(\Omega)}=0,\quad\forall y\in V(\mathcal{ T}_{h})\right\}\\ =&\left\{z\in V\,\Big{|}\,a(z,y)=0,\quad\forall y\in V (\mathcal{T}_{h})\right\}.\end{split} \tag{50}\]

The DPG formulation reads: find \((z,q)\in L\times\tilde{V}\) such that, for all \(y\in V(\mathcal{T}_{h})\),

\[\sum_{T\in\mathcal{T}_{h}}(z,\tilde{A}y)_{L^{2}(T)}+\sum_{T\in\mathcal{T}_{h} }\langle Dq,y\rangle_{\tilde{V}(T),V(T)}+\tfrac{1}{2}((M-D)q,y\rangle_{\tilde {V},V}=\sum_{T\in\mathcal{T}_{h}}\,(f,v)_{\tilde{V}(T),V(T)}+\tfrac{1}{2} \langle(M-D)g,y\rangle_{\tilde{V},V}. \tag{51}\]

The introduction of the hybrid face variables \(q\in\tilde{V}\) is necessary since \(z\in L\) does not satify (8). In practice, assuming that the traces of \(y\in V(\mathcal{T}_{h})\) are well-defined and belong to a space \(X(\mathcal{F}_{i,b})\), we can formulate (51) as follows: find \((z,q)\in L\times X(\mathcal{F}_{i,b})\) such that, for all \(y\in V(\mathcal{T}_{h})\),

\[\sum_{T\in\mathcal{T}_{h}}(z,\tilde{A}y)_{L^{2}(T)}+\sum_{F\in\mathcal{F}_{i}} (\mathcal{D}q,[\![y]\!])_{X(F)}+\tfrac{1}{2}\sum_{F\in\mathcal{F}_{h}}\,(( \mathcal{M}-\mathcal{D})q,y)_{X(F)}=\sum_{T\in\mathcal{T}_{h}}\,(f,v)_{\tilde {V}(T),V(T)}+\tfrac{1}{2}\sum_{F\in\mathcal{F}_{h}}\,((\mathcal{M}-\mathcal{D })g,y)_{X(F)}, \tag{52}\]

where \(X(\mathcal{F}_{i,b})\) is, for example, \([H^{-\frac{1}{2}}(\mathcal{F}_{i,b})]^{d}\times[H^{\frac{1}{2}}(\mathcal{F}_{ i,b})]^{d}\) for compressible linear elasticity, \(H^{-\frac{1}{2}}(\mathcal{F}_{i,b})\times H^{\frac{1}{2}}(\mathcal{F}_{i,b})\) for the scalar advection-diffusion-reaction and \(L^{2}_{\mathcal{T}}(\mathcal{F}_{i,b})\times L^{2}_{\mathcal{T}}(\mathcal{F}_{ i,b})\) for the Maxwell equations in stationary regime, with \(L^{2}_{\mathcal{T}}(\mathcal{F}_{i,b})\) being the space of fields in \(H(\operatorname{curl},\mathcal{T}_{h})\) whose tangential component belongs to \([L^{2}(\mathcal{F}_{i,b})]^{3}\).

The problem (51) above is well-posed and consistent [12, Lemma 2.4] with the previous formulation in (13a). We consider the optimal norms,

\[\|(z,q)\|_{\mathcal{U}}^{2}=\sum_{T\in\mathcal{T}_{h}}\|z\|_{L(T)}^{2}+\|q\|_{ \tilde{V}}^{2},\qquad\|y\|_{\mathcal{Y}}^{2}=\sum_{T\in\mathcal{T}_{h}}\| \tilde{A}y\|_{L(T)}^{2}+\|[\![y]\!]\|_{\partial\Omega_{h}}^{2},\quad\text{with} \quad\|[\![y]\!]\|_{\partial\Omega_{h}}=\sup_{q\in\tilde{V}}\frac{a(q,y)}{\| q\|_{\tilde{V}}}, \tag{53}\]

or formally, considering (52)

\[\|(z,q)\|_{\mathcal{U}}^{2}=\sum_{T\in\mathcal{T}_{h}}\|z\|_{L(T)}^{2}+\sum_{F \in\mathcal{F}_{i,b}}\|q\|_{X(F)}^{2},\qquad\|y\|_{\mathcal{Y}}^{2}=\sum_{T\in \mathcal{T}_{h}}\|\tilde{A}y\|_{L(T)}^{2}+\sum_{F\in\mathcal{F}_{i}}\| \mathcal{D}_{F}[\![y]\!]\|_{X(F)}^{2}+\sum_{F\in\mathcal{F}_{h}}\|(\mathcal{ M}^{t}-\mathcal{D})y\|_{X(F)}^{2}. \tag{54}\]

With these optimal norms for the trial and test spaces we have the following result [13, Theorem 2.6].

**Theorem 4** (Optimally stable formulation).: _The bilinear form \(b:(L\times\tilde{V},\|\bullet\|_{\mathcal{U}})\to(V(\mathcal{T}_{h}),\| \bullet\|_{\mathcal{Y}})\) defined as_

\[b(u,y)=\sum_{T\in\mathcal{T}_{h}}(z,\tilde{A}y)_{L^{2}(T)}+\sum_{T\in\mathcal{T }_{h}}\langle Dq,y\rangle_{\tilde{V}(T),V(T)}+\tfrac{1}{2}\langle(M-D)q,y \rangle_{\tilde{V},V} \tag{55}\]

_with \(u=(z,q)\), is an isometry between \(L\times\tilde{V}\) and \(V^{\prime}(\mathcal{T}_{h})\): we have that \(\gamma=\beta=\beta^{*}=1,\) where_

\[\gamma:=\sup_{u\in\mathcal{U}}\sup_{y\in\mathcal{Y}}\frac{b(u,y)}{\|u\|_{ \mathcal{U}}\|y\|_{\mathcal{Y}}},\quad\beta:=\inf_{u\in\mathcal{U}}\sup_{y\in \mathcal{Y}}\frac{b(u,y)}{\|u\|_{\mathcal{U}}\|y\|_{\mathcal{Y}}},\quad\beta^{* }:=\inf_{y\in\mathcal{Y}}\sup_{u\in\mathcal{U}}\frac{b(u,y)}{\|u\|_{\mathcal{U }}\|y\|_{\mathcal{Y}}} \tag{56}\]

_with \(\mathcal{U}=(L\times\tilde{V},\|\bullet\|_{\mathcal{U}})\)._

This property is inherited at the discrete level as long as fixed \(U_{h}^{N_{h}}\subset\mathcal{U}\) a discretization of the trial space with \(\dim\,Z_{h}=N_{h}\), the discrete test space \(Y_{h}^{N_{h}}\subset V(\mathcal{T}_{h})\) is the set of supremizers

\[Y_{h}^{N_{h}}=\operatorname{span}\left\{y_{u_{h}}\in V(\mathcal{T}_{h})\,\bigg{|} \,y=\operatorname*{argmax}_{y\in V(\mathcal{T}_{h})}\frac{b(u_{h},y)}{\|y\|_{ \mathcal{Y}}},\quad u_{h}\in Z_{h}^{N_{h}}\right\} \tag{57}\]
and \(\dim\ U_{h}=\dim\ Y_{h}\), see [13, Lemma 2.8]. In particular, for every \(u_{h}=(z_{h},q_{h})\in U_{h}^{N_{h}}\), we have the optimal _a posteriori_ error estimate

\[\|u-u_{h}\|_{\mathcal{U}}=\sup_{y\in V(\mathcal{T}_{h})}\frac{b(u-u_{h},y)}{\|y \|_{Y}}=\|r_{h}(u_{h})\|_{(V(\mathcal{T}_{h}))^{\prime}},\quad\langle r_{h}(u_ {h}),v\rangle_{(V(\mathcal{T}_{h}))^{\prime},V(\mathcal{T}_{h})}=\langle f,v \rangle_{(V(\mathcal{T}_{h}))^{\prime},V(\mathcal{T}_{h})}-b(u_{h},v). \tag{58}\]

The same reasoning can be iterated another time to perform model order reduction with the choice \(V_{n}=\{\psi_{j}^{\rm RB}\}_{j=1}^{r}\subset U_{h}^{N_{h}}\subset\mathcal{U}\), and

\[Y^{RB}=\text{span}\left\{y_{u_{RB}}\in V(\mathcal{T}_{h})\,\Big{|}\,y=\mathop {\rm argmin}_{y\in V(\mathcal{T}_{h})}\frac{b(u_{RB},y)}{\|y\|_{\mathcal{Y}}},\quad u_{RB}\in V_{n}\right\}, \tag{59}\]

such that for \(u_{RB}\in V_{n}\),

\[\|u-u_{RB}\|_{\mathcal{U}}=\|r_{RB}(u_{RB})\|_{(V(\mathcal{T}_{h}))^{\prime}},\quad\langle r_{RB}(u_{RB}),v\rangle_{(V(\mathcal{T}_{h}))^{\prime},V( \mathcal{T}_{h})}=\langle f,v\rangle_{(V(\mathcal{T}_{h}))^{\prime},V( \mathcal{T}_{h})}-b(u_{RB},v). \tag{60}\]

The main difficulty is the evaluation of the trial spaces \(Y_{h}^{N_{h}}\) and \(Y^{RB}\) since the bilinear form \(b\) may depend on the parameters \(\boldsymbol{\rho}\). If the parameters affect only the source terms, the boundary conditions or the initial conditions for time-dependent FS, this problem is avoided. The evaluation of \(Y_{h}^{N_{h}}\) can be performed locally for each element \(T\in\mathcal{T}_{h}\), differently from \(Y^{RB}\). An example of the evaluation of the basis of \(Y_{h}^{N_{h}}\) is presented in [13, Equations 24, 25] for linear scalar hyperbolic equations that can be interpreted as FS.

## 5 Domain decomposable Discontinuous Galerkin ROMs

Extreme-scale parametric models are unfeasible to reduce with standard approaches due to the high computational costs of the offline stage. Parametric multi-physics simulations, such as fluid-structure interaction problems, are reduced inefficiently with a global reduced basis, depending on the complexity of the interactions between the physical models considered and the parametric dependency. In some cases, only a part of a decomposable system is reducible with a ROM, thus a possible solution is to implement a ROM-FOM coupling through an interface. In presence of moving shocks [7] affected by the parametrization, one may want to isolate these difficult features to approximate and apply different dimension reduction methodologies depending on the subdomain. These are the main reasons to develop domain decomposable or partitioned ROMs (DD-ROMs).

Some approaches from the literature are the reduced basis element methods [62, 63], the static condensation reduced basis element method [48, 27], non-intrusive methods based on local regressions in each subdomain [86, 87], overlapping Schwarz methods [23, 50], optimization-based MOR approaches [71, 72] and hyper-reduced ROMs [60]. In this last case, local approximations are useful because the local reduced dimensions are smaller and therefore more accurate local regressions can be designed to perform non-intrusive surrogate modelling. Little has been developed for the DG method, even though its formulation imposes naturally flux and solution interface penalties at the internal boundaries of the subdomains, in perspective of performing model order reduction. In our case, the linear systems associated to the parametric models are algebraically partitioned in disjoint subdomains coupled with the standard penalties from the weak DG formulations, without the need to devise additional operators to perform the coupling as long as the interfaces' cuts fall on the cell boundaries.

Another less explored feature of DD-ROMs is the possibility to repartition the computational domain, while keeping the data structures relative to each subdomain local in memory, with the aim of obtaining more efficient or accurate ROMs. In fact, one additional reason to subdivide the computational domain is to partition the solution manifold into local solution manifolds that have a faster decay of the Kolmogorov n-width. The repartition of the computational domain can be performed with _ad hoc_ domain decomposition strategies. To our knowledge, the only case found in the literature is introduced in [87], where the degrees of freedom are split in each subdomain minimising the communication and activity between them and balancing the computational load across them. Relevant is the choice of weights to assign to each degree of freedom: uniform weights, nodal values of Reynolds stresses for the turbulent Navier-Stokes equations or the largest singular value of the discarded local POD modes. In particular, the last option results in a balance of energy in \(L^{2}\) norm retained in each subdomain. We explore a different approach.

It must be remarked that in any case, fixed the value of the reconstruction error of the training dataset in Frobenious norm \(\|\cdot\|_{F}\), there must be at least a local reduced basis dimension greater or equal to the global reduced basis space dimension. In fact, for the Eckhart-Young theorem, if \(X\in\mathbb{R}^{d\times n}\) is the snapshots matrix ordered by columns, we have that the projection into the first \(k\) modes \(\{v_{k}\}_{i=1}^{k}\) achieves the best approximation error in the Frobenious norm in the space of matrices of rank \(k\):

\[P_{k}=\mathop{\rm argmin}_{P\in\mathbb{R}^{m\times n\times s.t.\ r(P)=k}}\|X- PX\|_{F},\quad P_{k}=\sum_{i=1}^{k}v_{i}\otimes v_{i}, \tag{61}\]
where \(r(\cdot)\) is the matrix rank. So, in general, it is not possible to achieve a better training approximation error in the Frobenious norm \(\|\cdot\|_{F}\) employing a number of local reduced basis smaller than what would be needed to achieve the same accuracy with a global reduced basis. So, differently from [87], instead of balancing the local reduced basis dimension among subdomains, we repartition the computational domain in regions whose restricted solution manifold is easily approximable by linear subspaces and regions for which more modes are needed. Anyway, for truly decomposable systems we expect that the reconstruction error on the test set is lower when considering local reduced basis instead of global ones, as will be shown for the Maxwell equation in stationary regime test case with discontinuous piecewise constant parameters, see Figure 7.

### Implementation of Domain Decomposable ROMs

Let us assume that the full-order model is implemented in parallel with distributed memory parallelism with \(K>1\) cores, i.e., each \(i\)-th core owns locally the data structures relevant only to its assigned subdomain \(\Omega_{i}\) of the whole computational domain \(\cup_{i=1}^{K}\Omega_{i}=\Omega_{h}\subset\mathbb{R}^{d}\), for \(i=1,\ldots,K\). We will employ the deal.II library [6] to discretize the FS with the DG method, assemble the associated linear systems and solve them in parallel [9]. In particular, we employ p4est[14] to decompose the computational domain, PETSc[8] to assemble the linear system and solve it at the full-order level and petsc4py[22] to assemble and solve the reduced order system. At the offline and online stages the computations are performed in a distributed memory setting in which each core assembles its own affine decomposition, so that the evaluation of the reduced basis and of the projected local operators is always performed in parallel.

The weak formulation (35) is easily decomposable thanks to the additive properties of the integrals. We recall the definition of the weak formulation \(\forall y_{h}\in V_{h},\ z_{h}\in V^{*}\)

\[\begin{split} a_{h}^{cf}(z_{h},y_{h})+s_{h}(z_{h},y_{h})=\sum_{T \in\mathcal{T}_{h}}(z_{h},\tilde{A}y_{h})_{L^{2}(T)}+\frac{1}{2}\sum_{F\in \mathcal{F}_{h}^{b}}\left((\mathcal{M}+\mathcal{D})z_{h},y_{h}\right)_{L^{2}(F )}+\sum_{F\in\mathcal{F}_{h}^{s}}\left(\mathcal{D}_{F}\{\!\{z_{h}\}\!\},\!\{y_ {h}\!\}\right)_{L^{2}(F)}+\\ \sum_{F\in\mathcal{F}_{h}^{s}}(S_{F}^{b}z_{h},y_{h})_{L^{2}(F)}+ \sum_{F\in\mathcal{F}_{h,i}^{s}}(S_{h}^{i}\!\!\{z_{h}\}\!\},\!\{\llbracket y_ {h}\!\rrbracket\})_{L^{2}(F)},\end{split} \tag{62}\]

and we decompose it into the \(K\) subdomains as

\[\begin{split} a_{h}^{cf}(z_{h},y_{h})+s_{h}(z_{h},y_{h})=& \sum_{i=1}^{K}\left(\sum_{T\in\mathcal{T}_{h,i}}(z_{h},\tilde{A}y_{h})_{L^{2}( T)}+\frac{1}{2}\sum_{F\in\mathcal{F}_{h,i}^{s}}\left((\mathcal{M}+\mathcal{D})z_{ h},y_{h}\right)_{L^{2}(F)}+\sum_{F\in\mathcal{F}_{h,i}^{s}}\left(\mathcal{D}_{F} \{\!\{z_{h}\}\!\},\!\{y_{h}\!\}\right)_{L^{2}(F)}+\\ &\sum_{F\in\mathcal{F}_{h,i}^{s}}(S_{F}^{b}z_{h},y_{h})_{L^{2}(F )}+\sum_{F\in\mathcal{F}_{h,i}^{s}}\left(S_{h}^{i}\!\!\{z_{h}\},\!\{y_{h}\!\} \right)_{L^{2}(F)}\right)+\sum_{j=i}^{K}\left(\sum_{F\in\mathcal{F}_{h,i,j}^{s} }\left(\mathcal{D}_{F}\{\!\{z_{h}\}\!\},\!\{y_{h}\!\}\right)_{L^{2}(F)}+\\ &\sum_{F\in\mathcal{F}_{h,i,j}^{s}}\left(S_{h}^{i}\!\!\{z_{h}\}, \!\{y_{h}\!\}\right)_{L^{2}(F)}\right),\end{split} \tag{63}\]

\[l_{h}(y_{h})=\sum_{T\in\mathcal{T}_{h}}(f,y_{h})_{L^{2}(T)}+\frac{1}{2}\sum_{F \in\mathcal{F}_{h}^{b}}\left((M-D)g,y_{h}\right)_{L^{2}(F)}=\sum_{i=1}^{K} \left(\sum_{T\in\mathcal{T}_{h,i}}(f,y_{h})_{L^{2}(T)}+\frac{1}{2}\sum_{F\in \mathcal{F}_{h,i}^{s}}\left((M-D)g,y_{h}\right)_{L^{2}(F)}\right), \tag{64}\]

where we have defined the internal subsets \(\mathcal{T}_{h,i}=\mathcal{T}_{h}\cap\Omega_{i}\), \(\mathcal{F}_{h,i}^{i}=\mathcal{F}_{h}^{i}\cap\bar{\Omega}_{i}\) and \(\mathcal{F}_{h,i}^{b}=\mathcal{F}_{h}^{b}\cap\bar{\Omega}_{i}\), \(\forall i=1,\ldots,K\) and the interfaces subsets \(\mathcal{F}_{h,i,j}^{i}=\mathcal{F}_{h}^{i}\cap\overline{\Omega}_{i}\cap \overline{\Omega}_{j}\) and \(\mathcal{F}_{h,i,j}^{b}=\mathcal{F}_{h}^{b}\cap\overline{\Omega}_{i}\cap \overline{\Omega}_{j}\), \(\forall i=1,\ldots,K\). We remark that the computational domain is always decomposed such that the cuts of the subdomains \(\{\partial\Omega_{i}\}_{i=1}^{K}\) fall on the interfaces of the triangulation \(\mathcal{F}_{h}^{i}\cup\mathcal{F}_{h}^{b}\).

We define the bilinear and linear operators in \(V_{h}^{*}\),

\[\begin{split}\mathcal{A}_{ii}=\sum_{T\in\mathcal{T}_{h,i}}( \bullet,\tilde{A}\bullet)_{L^{2}(T)}+\frac{1}{2}\sum_{F\in\mathcal{F}_{h,i}^{b} }\left((\mathcal{M}+\mathcal{D})\bullet,\bullet\right)_{L^{2}(F)}+\sum_{F\in \mathcal{F}_{h,i}^{i}}\left(\mathcal{D}_{F}\{\!\{\bullet\}\!\},\!\{\bullet \}\right)_{L^{2}(F)}+\\ \sum_{F\in\mathcal{F}_{h,i}^{s}}(S_{F}^{b}\bullet,\bullet)_{L^{2}( F)}+\sum_{F\in\mathcal{F}_{h,i}^{s}}(S_{h}^{i}\!\!\{\bullet\},\!\{\bullet \})_{L^{2}(F)},\quad\forall i=1,\ldots,K,\end{split} \tag{65}\] \[\mathcal{A}_{ij}=\mathcal{A}_{ji}=\sum_{F\in\mathcal{F}_{h,i,j}^{s} }\left(\mathcal{D}_{F}\{\!\{\bullet\}\!\},\!\{\bullet\}\right)_{L^{2}(F)}+\sum_{F \in\mathcal{F}_{h,i,j}^{s}}(S_{h}^{i}\!\!\{\bullet\},\!\{\bullet\},\!\{ \bullet\})_{L^{2}(F)},\quad\forall j,i=1,\ldots,K,\ i\neq j,\] (66) \[\mathcal{F}_{i}=\sum_{T\in\mathcal{T}_{h,i}}(f,\bullet)_{L^{2}(T)}+ \frac{1}{2}\sum_{F\in\mathcal{F}_{h,i}^{b}}\left((M-D)g,\bullet\right)_{L^{2}(F)},\quad\forall i=1,\ldots,K, \tag{67}\]
and their matrix representation in the discontinuous Galerkin basis of \(V_{h}\),

\[\left.\left(\mathcal{A}_{ii}\right)\right|_{V_{h}}=A_{ii},\qquad\mathcal{F}_{i} \big{|}_{V_{h}}=F_{i},\quad\forall i=1,\ldots,K,\qquad\left(\mathcal{A}_{ij} \right)\big{|}_{V_{h}}=\left.\left(\mathcal{A}_{ji}\right)\right|_{V_{h}}=A_{ ij}=A_{ji},\quad\forall j,i=1,\ldots,K,\ i\neq j, \tag{68}\]

and in the local reduced basis \(V_{i}=\{\psi_{j,i}^{\rm RB}\}_{j=1}^{r}\subset V_{h}(\Omega_{i}),\ i=1,\ldots,K\),

\[\left.\left(\mathcal{A}_{ii}\right)\right|_{V_{RB}}=B_{ii},\qquad\mathcal{F}_{ i}\big{|}_{V_{RB}}=L_{i},\quad\forall i=1,\ldots,K,\qquad\left(\mathcal{A}_{ ij}\right)\right|_{V_{RB}}=\left.\left(\mathcal{A}_{ji}\right)\right|_{V_{RB}}=B_{ ij}=B_{ji},\quad\forall j,i=1,\ldots,K,\ i\neq j. \tag{69}\]

As anticipated, in our test cases the subdomains interface penalties are naturally included inside \(\{\mathcal{A}_{ij}\}_{i,j=1,\ldots,K}\). In practice, additional penalty terms could be implemented:

\[\mathcal{S}_{ij}=\sum_{F\in\mathcal{F}_{h,i,j}^{k}}(S[\bullet],[\bullet])_{L^ {2}(F)},\qquad\mathcal{S}_{ij}\big{|}_{V_{h}}=S_{ij}\qquad\left.\left( \mathcal{A}_{ij}+\mathcal{S}_{ij}\right)\right|_{V_{RB}}=B_{ij},\quad\forall j,i=1,\ldots,K,\ i\neq j. \tag{70}\]

A matrix representation of the projection of the full-order block matrix \((A_{ij})_{i,j=1}^{K}\in\mathbb{R}^{d\times d}\) into the reduced order block matrix \((B_{ij})_{i,j=1}^{K}\in\mathbb{R}^{Kr\times Kr}\) is shown in Figure 1 for \(K=4\). We remark that, differently from continuous Galerkin formulations, the DG penalization on jumps across the interfaces is already enough to couple the subdomains and there is no need of further stabilization, as shown in Figure 1. Nonetheless, additional interface penalties terms can be easily introduced, taking also into account DG numerical fluxes. The reduced dimension is the number of subdomains \(K\) times the local reduced basis dimensions \(\{r_{i}\}_{i=1}^{K}\), here supposed equal \(r=r_{i},\ i=1,\ldots,K\), but in general can be different.

### Repartitioning strategy

A great number of subdomains can pollute the efficiency of the developed DD-ROMs at the online stage since the reduced dimension would be \(\sum_{i=1}^{K}r_{i}\) that scales linearly with the number of cores if the local reduced dimensions \(r_{i}\) are equal. In order to keep the computational savings in the assembly of the affine decomposition at the offline stage, we may want to preserve the distributed property of our ROM. One possible solution is to fix a reduced number of subdomains \(k\ll K\) such that \(\sum_{i=1}^{k}r_{i}\) is small enough to achieve a significant speedup with respect to the FOM. The additional cost with respect to a monodomain ROM is associated to the evaluation of the \(k\) local reduced basis with SVD and the assembly of the affine decomposition operators. The new \(k\) reduced subdomains do not need to be agglomerations of the FOM subdomains, hence, different strategies to assemble the new \(k\) reduced subdomains can be investigated.

The number of subdomains \(K\) was kept the same as the FOM since it is necessary to collect the snapshots efficiently at the full-order level through _p4est_. However, if we decide to repartition our computational domain, we can develop decomposition strategies that reduce \(\sum_{i=1}^{K}r_{i}\). Ideally, having in mind the Eckhart-Young theorem, a possible strategy is to lump together all the dofs of the cells that have a fast decaying Kolmogorov n-width, and focus on the remaining ones. We test this procedure in the practical case \(k=2,\ K=4\) to perform numerical experiments in section 5.3.

To solve the classification problem of partitioning the elements of the mesh into \(k\) subdomains, we describe here two scalar indicators that will be used as metrics. For \(k=2\) subdomains, it will be sufficient to choose the percentage of cells \(P_{l}\) corresponding to the lowest values of the chosen scalar indicator. Other strategies for \(k>2\) may also involve clustering

Figure 1: Assembly of the reduced block matrix \(\{B_{i,j}\}_{i,j=1}^{4}\) through the projection onto the local reduced basis \(\{V_{i}\}_{i=1}^{4}\) of the full-order partitioned matrix \(A=\{A_{i,j}\}_{i,j=1}^{4}\) when considering 4 subdomains. The natural DG penalty terms are included in the matrix \(A\) without the need for additional penalty terms \(\{S_{i,j}\}_{i\neq j,\,i,j=1}^{4}\) to impose stability at the reduced level.
algorithms and techniques to impose connectedness of the clusters, as done for local dimension reduction in parameter spaces in [76]. A first crude and cheap indicator to repartition the computational domain is the cellwise variance of the training snapshots, as it measures how well, in mean squared error, the training snapshots are approximated by their mean, \(\forall T\in\mathcal{T}_{h}\).

**Definition 3** (Cellwise variance indicator).: _We define the cellwise variace indicator \(I_{var}:\mathcal{T}_{h}\to\mathbb{R}^{+}\),_

\[I_{var}(T)=\int_{T}\lVert\text{Var}(\{\mathbf{z}(\boldsymbol{\rho}_{i})\}_{i=1 }^{n})\rVert_{L^{2}(\mathbb{R}^{m})}\ d\mathbf{x},\qquad(\text{Var}(\{\mathbf{ z}(\boldsymbol{\rho}_{i})\}_{i=1}^{n}))_{l}=\tfrac{1}{n}\sum_{i=1}^{n}\left| \mathbf{z}_{i}(\boldsymbol{\rho}_{i})-\tfrac{1}{n}\sum_{j=1}^{n}\mathbf{z}_{l} (\boldsymbol{\rho}_{j})\right|^{2},\quad l=1,\ldots,m, \tag{71}\]

_where \(n>0\) is the number of training DG solutions \(\{\mathbf{z}(\boldsymbol{\rho}_{i})\}_{i=1}^{n}\) with \(\mathbf{z}(\boldsymbol{\rho}_{i}):\Omega\subset\mathbb{R}^{d}\to\mathbb{R}^{m},\ \forall i\in\{1,\ldots,n\}\)._

Note that the indicator is a scalar function on the set of elements of the triangulation \(\mathcal{T}_{h}\). This is possible thanks to the assumption that boundaries of the subdomains belong to the interfaces of the elements of \(\mathcal{T}_{h}\). When this hypothesis is not fulfilled, we would need to evaluate additional operators to impose penalties at the algebraical interfaces between subdomains that are not included in the set \(\mathcal{F}_{h}^{i}\cup\mathcal{T}_{h}^{b}\), not to degrade the accuracy.

The cellwise variance indicator is effective for all the test cases for which there is a relatively large region that is not sensitive to the parametric instances, as in our advection diffusion reaction test case in Section 5.3.3. Common examples are all the CFD numerical simulations that have a far field with fixed boundary conditions. However, the variance indicator may be blind to regions in which the snapshots can be spanned by a one or higher dimensional linear subspace and are not well approximated by a constant field, as in the compressible linear elasticity test case in Section 5.3.2.

In these cases, a valid choice is represented by a cellwise Grassmannian dimension indicator. We denote with \(D_{\text{T}}\) the number of degrees of freedom associated to each element \(T\), assumed constant in our test cases.

**Definition 4** (Cellwise Grassmannian dimension indicator).: _Fixed \(1\leq r_{T}\in\mathbb{N}\), and \(1\leq n_{neig}\in\mathbb{N}\), we define the cellwise Grassmannian dimension indicator \(I_{G}:\mathcal{T}_{h}\to\mathbb{R}^{+}\),_

\[I_{G}(T)=\lVert X_{T}-U_{T}U_{T}^{T}X_{T}\rVert_{F}, \tag{72}\]

_where \(X_{T}\in\mathbb{R}^{n_{neigh}D_{T}\times n}\) is the snapshots matrix restricted to the cell \(T\) and its \(n_{neig}\) nearest neighbours, and \(U_{T}\in\mathbb{R}^{n_{neigh}D_{T}\times\tau_{T}}\) are the modes of the truncated SVD of \(X_{T}\) with dimension \(r_{T}\)._

The cellwise Grassmannian dimension indicator \(I_{G}\) is a measure of how well the training snapshots restricted to a neighbour of each cell are approximated by a \(r_{T}\) dimensional linear subspace. Employing this indicator, we recover an effective repartitioning of the computational subdomain of the compressible linear elasticity test case, see Section 5.3.2. The Grassmannian indicator has two hyper-parameters that we fix for each test case in section 5.3: the number of nearest-neighbour cells is \(n_{neigh}=3\) and the number of reduced local dimension used to evaluate the \(L^{2}\) reconstruction error is \(r_{T}=1\). The number of nearest-neighbour is chosen to deal with critical cases at the boundaries and the closest neighbouring cells are chosen based on the distance of barycenters. The reduced local dimension \(r_{T}=1\) is chosen very small as the computations must be done on very few cells.

We remark that both indicators do not guarantee that the obtained subdomains belong to the same connected components and, though this might be a problem in terms of connectivity and computational costs for the FOM, at the reduced level this does not affect the online computational costs. Nevertheless, in the tests we perform, the obtained subdomains are connected.

Now, the assembly of the affine decomposition proceeds as explained in Section 5.1 with the difference that at least one local reduced basis and reduced operator is split between at least 2 subdomains/cores. A schematic block matrix representation of the procedure is shown in Figure 2.

### Numerical experiments

In this section, we test the presented methodology for different linear parametric partial differential equations: the Maxwell equations in stationary regime in section 5.3.1 (**MS**), the compressible linear elasticity equations in section 5.3.2 (**CLE**) and the advection diffusion reaction equations in section 5.3.3 (**ADR**). We study two different parametrizations for the test cases **MS** and **CLE**: one with parameters that affect the whole domain **MS1** and **CLE1**, and one with parameters that affect independently different subdomains **MS2** and **CLE2**. We show a case in which DD-ROMs work effectively **MS2** and a case **CL2** in which the performance is analogous to single domain ROMs, even if the parameters have a local influence.

We test the effectiveness of the _a posteriori_ error estimates introduced in section 4.1, the accuracy of DD-ROMs for \(K=4\) and the results of repartitioning strategies with \(k=2\) subdomains. When performing a repartition of the computational domain \(\Omega\) in subdomains \(\{\Omega_{i}\}_{i=1}^{k}\) with reduced dimensions \(\{r_{\Omega_{i}}\}_{i=1}^{k}\), we call the subdomains with lower values of the variance indicator \(I_{\text{var}}\), see definition 3, _low variance_ regions and with lower values of the Grassmannian indicator \(I_{G}\), see definition 4,
low Grassmannian reconstruction error_. The complementary subdomains are the _high variance_ and _high Grassmannian reconstruction error_ regions, respectively. We show a case (**CLE1**) in which the Grassmannian indicator detects a better partition in terms of local reconstruction error with respect to the variance indicator.

We will observe that the relative errors in \(R\)-norm and energy norm and the \(L^{2}\) relative error estimator and \(L^{2}\) relative energy norm estimator are the most affected by the domain partitions.

The open-source software library employed for the implementation of the full-order Friedrichs' systems discontinuous Galerkin solvers is deal.II[6] and we have used piecewise \(\mathbb{P}^{2}\) basis functions in all simulations. The partition of the computational domain is performed in deal.II through the open-source p4est package [14]. The distributed affine decomposition data structures are collected in the offline stage and exported in the sparse NumPy format [41]. The reduced order models and the repartition of the computational domains are implemented in Python with MPI-based parallel distributed computing mpi4py[21] and petsc4py[8] for solving the linear full-order systems through MUMPS[1], a sparse direct solver.

#### 5.3.1 Maxwell equations in stationary regime (MS)

We consider the parametric Maxwell equations in the stationary regime in \(d=3\) spatial dimensions, with \(m=6\) equations, on a torus \(\Omega\subset\mathbb{R}^{3}\) with inner radius \(r=0.5\) and outer radius \(R=2\) centered in \(\mathbf{0}\) and lying along the \((x,z)\) plane:

\[\left(\begin{array}{c}\mu\mathbf{H}+\nabla\times\mathbf{E}\\ \sigma\mathbf{E}-\nabla\times\mathbf{H}\end{array}\right)=\left(\begin{array}[ ]{c}\mathbf{g}\\ \mathbf{f}\end{array}\right),\quad\forall\mathbf{x}\in\Omega, \tag{73}\]

the tangential homogeneous boundary conditions \(\mathbf{n}\times\mathbf{E}=\mathbf{0}\) are applied with the boundary operator (18). We vary the parameters in the interval \(\boldsymbol{\rho}=(\mu,\sigma)\in[0.5,2]\times[0.5,3]\subset\mathbb{R}^{2}\), leading to \(\mu_{0}=\min(\mu,\sigma)\).

We consider the exact solutions

\[\mathbf{H}_{\text{exact}}(\mathbf{x}) =-\frac{1}{\mu}\left(\frac{2xy}{\sqrt{x^{2}+z^{2}}},-\frac{4y^{2} \sqrt{x^{2}+z^{2}}+\sqrt{x^{2}+z^{2}}(-12(x^{2}+z^{2})-15)+32(x^{2}+z^{2})}{4 (x^{2}+z^{2})},\frac{2xy}{\sqrt{x^{2}+z^{2}}}\right),\] \[\mathbf{E}_{\text{exact}}(\mathbf{x}) =\left(\frac{z}{\sqrt{x^{2}+z^{2}}},0,-\frac{x}{\sqrt{x^{2}+z^{2} }}\right)\cdot\left(r^{2}-y^{2}-\left(R-\sqrt{x^{2}+z^{2}}\right)^{2}\right).\]

We remark that the exact solutions can be approximated with a linear reduced subspace of dimension \(1\), if we obtained the reduced basis with a partitioned SVD on the fields \((\mathbf{H},\mathbf{E})\) separately. We do not choose this approach and perform a monolithic SVD to test the convergence of the approximation with a DD-ROMs with respect to the local reduced dimensions. The source terms are defined consequently as

\[\mathbf{g}(\mathbf{x})=0,\qquad\mathbf{f}(\mathbf{x})=\sigma\mathbf{E}_{ \text{exact}}-\nabla\times\mathbf{H}_{\text{exact}}. \tag{74}\]

We consider two parametric spaces:

\[\boldsymbol{\rho}=(\mu,\sigma)\in[0.5,2]\times[0.5,3] =\mathcal{P}_{1}\subset\mathbb{R}^{2},\qquad(\mathbf{MS1}) \tag{75a}\] \[\boldsymbol{\rho}=(\mu_{1},\sigma_{1},\mu_{2},\sigma_{2})\in[0.5,2]\times[0.5,3]\times[0.5,2]\times[0.5,3]=\mathcal{P}_{2}\subset\mathbb{R}^ {4},\qquad(\mathbf{MS2}) \tag{75b}\]

where in the second case, the parameters \(\mu\) and \(\sigma\) are now piecewise constant:

\[\mu=\begin{cases}\mu_{1},&x<0,\\ \mu_{2},&x\geq 0,\end{cases}\qquad\sigma=\begin{cases}\sigma_{1},&x<0,\\ \sigma_{2},&x\geq 0,\end{cases} \tag{76}\]

Figure 2: Repartitioning of the reduced block matrix shown in Figure 1 from \(K=4\) subdomains to \(k=2\) repartitioned subdomains. The projection from the full-order matrix \(A=\{A_{i,j}\}_{i,j=1}^{4}\) to the reduced matrix \(B=\{B_{i,j}\}_{i,j=1}^{2}\) is sketched. It is performed locally in a distributed memory setting, the re-ordering shown by the arrows is reported only to visually see the which block structure of the full-order matrix \(A\) would correspond to the blocks of the reduced matrix \(B\).
where \(\mathbf{x}=(x,y,z)\in\Omega\subset\mathbb{R}^{3}\). In Figure 3, we show solutions for \(\mu=\sigma=1\) and for discontinuous values of the parameters: \(\mu_{1}=\sigma_{1}=1\) in \(\{x<0\}\cap\Omega\) and \(\mu_{2}=\sigma_{2}=2\) in \(\{x\geq 0\}\cap\Omega\). The FOM partitioned and DD-ROM repartitioned subdomains are shown in Figure 4. For **MS1**, we choose the variance indicator to repartition the computational subdomain in two subsets: \(P_{l}=20\%\) of the cells for the _low variance_ part and \(80\%\) for the _high variance_ part. For **MS2**, we split the computational domain in two parts with the Grassmannian indicator and \(P_{l}=50\%\).

At the end of this subsection a comparison of the effectiveness of DD-ROMs with and without discontinuous parameters will be performed, the associated error plots are reported in Figure 6 and Figure 7. We will see that, for this simple test case **MS2**, there is an appreciable improvement of the accuracy when the computational domain subdivisions match the regions \(\{x<0\}\cap\Omega\) and \(\{x\geq 0\}\cap\Omega\) in which \(\mu\) and \(\sigma\) are constant. Such subdivision is detected by the Grassmannian indicator with \(P_{l}=50\%\), as shown in Figure 4 on the right. This is the archetypal case in which DD-ROMs are employed successfully, in comparison with **MS1** for which there is no significant improvement with respect to classical global linear reduced basis.

In Figure 5, we show how the different thresholds applied to the two indicators can affect the reconstruction error on a reduced space with \(r_{\Omega_{i}}=3\). All the lines plot the local relative error computed on different subdomains (either one of the \(k\) DD-ROM subdomains or on the whole domain). On the \(x\)-axis it is shown the percentage of cells that are grouped into the

Figure 4: **MS. Left**: FOM computational domain partitioned in \(K=4\) subdomains inside deal.II. **Center**: **MS1**, DD-ROM repartition of the computational subdomain \(k=2\) with the cellwise variance indicator \(I_{\text{var}}\), definition 3: \(20\%\) of the cells belong to the _low variance_ part, represented in blue inside the torus, and the other \(80\%\) belong to the _high variance_ part, represented in red. **Right**: **MS2**, DD-ROM repartition with variance indicator \(P_{l}=50\%\). The computational domain is exactly split at the interfaces that separate the subdomains \(\{x<0\}\cap\Omega\) and \(\mu=\sigma=2\) in which the parameters \(\mu\) and \(\sigma\) are constant.

Figure 3: **MS. Electric and magnetic fields of the Maxwell equations in stationary regime with Dirichlet homogeneous boundary conditions \(\mathbf{n}\times\mathbf{E}=\mathbf{0}\). The vectors of the magnetic and electric fields are scaled by \(0.5\) and \(2\) of their magnitude respectively. **Left**: **MS1**, \(\mu=\sigma=1\), test case errors shown in Figure 6. **Right**: **MS2**, \(\mu=\sigma=1\) in \(\{x<0\}\cap\Omega\) and \(\mu=\sigma=2\) in \(\{x\geq 0\}\cap\Omega\), test case errors shown in Figure 7.
low variance or low Grassmannian DD-ROM subdomain. We observe that the cellwise variance indicator is a good choice for the purpose of repartitioning the subdomain from \(K=4\) to \(k=2\). Indeed, it is possible to build a low variance subdomain (value of the abscissa \(20\%\) in Figure 5) with a low local relative reconstruction error (\(5\cdot 10^{-4}\)) with respect to the global one (\(8.6\cdot 10^{-4}\)). This means that choosing the threshold \(P_{l}=20\%\) for the low variance subdomain, we should be able to use less reduced basis functions for that subdomain without affecting too much the global error.

Test case **MS1**. We evaluate \(n_{\text{train}}=20\) training full-order solutions and \(n_{\text{test}}=80\) test full-order solutions, corresponding to a uniform independent sampling from the parametric domain \(\mathcal{P}_{1}\subset\mathbb{R}^{2}\). Figure 6 shows the result relative to the relative \(L^{2}\)-error and relative errors in energy norm, with associated _a posteriori_ estimators. The number abscissae \(0,5,10,\ldots,95\) represents the train parameters \(n_{\text{train}}=20\) while the others \(n_{\text{test}}=80\) parameters are the test set. For these studies, we have fixed the local reduced dimensions to \(r_{\Omega_{i}}=3,\;i=1,\ldots,K\) for \(K=4\), \(r_{\Omega}=3\) for the whole computational domain and \(r_{\Omega_{1}}=2,\;r_{\Omega_{2}}=3\) for the DD-ROM repartitioned case with \(k=2\). This choice of repartitioning with the \(20\%\) of _low variance_ cells and local reduced dimension \(r_{\Omega_{1}}=2\) does not deteriorate significantly the accuracy and the errors almost coincide for all approaches. However, unless the parameters \(\sigma,\mu\) assume different discontinuous values in the computational domain \(\Omega\), DD-ROMs are not advisable for this test case if the objective is improving the predictions' accuracy.

Test case **MS2**. Similarly to the previous case, we evaluate \(n_{\text{train}}=20\) training full-order solutions and \(n_{\text{test}}=80\) test full-order solutions, corresponding to a uniform independent sampling from the parametric domain \(\mathcal{P}_{2}\subset\mathbb{R}^{4}\). As mentioned above, if we vary the parameters \(\boldsymbol{\rho}=(\mu,\sigma)\) discontinuously on the subdomains \(\{x\geq 0\}\cap\Omega\) and \(\mathbf{x}\in\{x<0\}\cap\Omega\), we obtain the results shown in Figure 7. It can be seen that repartitioning \(\Omega\) in \(k=2\) DD-ROM subdomains with the local Grassmannian indicator \(I_{G}\) and \(P_{l}=50\%\) produces effective DD-ROMs compared to the case of a single reduced solution manifold for the whole computational domain and for the DD-ROM with \(k=4\) for which the subdomains do not match \(\{x<0\}\cap\Omega\) and \(\{x\geq 0\}\cap\Omega\). In this case, we kept the local dimension of DD-ROM repartitioned case with \(k=2\) equal \(r_{\Omega_{1}}=r_{\Omega_{2}}=3\). For this simple test case, there is an appreciable improvement for some test parameters in the accuracy for \(k=2\) instead of \(K=4\) or a classical global linear basis ROM.

In Table 1, we list the computational times and speedups for a simulation with the different methods. For an error convergence analysis with respect to the size of the reduced space, we refer to Appendix C.

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|c|} \hline \multicolumn{3}{|c|}{FOM} & \multicolumn{3}{|c|}{ROM} & \multicolumn{3}{|c|}{DD-ROM} \\ \hline \(N_{h}\) & time & \(r\) & time & speedup & \(r_{i}\) & time & speedup \\ \hline \hline
6480 & 254.851 [ms] & 3 & 51.436 [\(\mu\)s] & \(\sim\) 495 & [3, 3, 3, 3] & 62.680 [\(\mu\)s] & \(\sim\) 406 \\ \hline \end{tabular}
\end{table}
Table 1: **MS1**. Average computational times and speedups for ROM and DD-ROM approaches for Maxwell equations. The speedup is computed as the FOM computational time over the ROM one. The FOM runs in parallel with 4 cores, so “FOM time” refers to wallclock time.

Figure 5: **MS1**. Local relative \(L^{2}\)-reconstruction errors of the snapshots matrix restricted to the two subdomains of the repartitioning performed with the indicator \(I_{\text{var}}\) (in red and light-blue), Definition 3, and \(I_{\text{G}}\) (in orange and blue), Definition 4. The relative \(L^{2}\)-reconstruction error attained on the whole domain is shown in black for the indicator \(I_{\text{var}}\) and in brown for the indicator \(I_{\text{G}}\). The local reduced dimensions used to evaluate the local reconstruction errors is \(r_{\Omega_{i}}=3,\;i=1,2\).
Figure 6: **MS1. Errors and estimators for Maxwell equations corresponding to the \(n_{\text{train}}=20\) uniformly sampled training snapshots corresponding to the abscissae \(0,5,10,\ldots,95\), and \(n_{\text{test}}=80\) uniformly sampled test snapshots, corresponding to the other abscissae. The reduced dimensions of the ROMs are \(\{r_{\Omega_{i}}\}_{i=1}^{K}=[3,3,3,3]\) for \(K=4\) partitions, \(r_{\Omega}=3\) for \(k=1\) partition, and \(\{r_{\Omega_{i}}\}_{i=1}^{k}=[2,3]\) for \(k=2\) partitions. For the case \(k=2\) we employed the cellwise variance indicator \(I_{G}\), Definition 3, with \(P_{l}=20\%\). It can be seen that even reducing the local dimension from \(3\) to \(2\) of one of the \(k=2\) repartitioned subdomains, the accuracy of the predictions does not decrease sensibly.**
Figure 7: **MS2**. Errors and estimators for Maxwell equations with discontinuous \(\mu\) and \(\sigma\) corresponding to the \(n_{\text{train}}=20\) uniformly sampled training snapshots corresponding to the abscissae \(0,5,10,\ldots,95\), and \(n_{\text{test}}=80\) uniformly sampled test snapshots, corresponding to the other abscissae. The reduced dimensions of the ROMs are \(\{r_{\Omega_{i}}\}_{i=1}^{K}=[3,3,3,3]\) for \(K=4\) partitions, \(r_{\Omega}=3\) for \(k=1\) partition, and \(\{r_{\Omega_{i}}\}_{i=1}^{k}=[3,3]\) for \(k=2\) partitions. For the case \(k=2\) we employed the cellwise local Grassmannian dimension indicator \(I_{G}\), Definition 3, with \(P_{l}=50\%\). The subdivisions detected exactly match the subdomains \(\{x<0\}\cap\Omega\) and \(\{x\geq 0\}\cap\Omega\) in which the parameters are constant. An improvement of the predictions can be appreciated for some test parameters when employing \(k=2\) repartitions.
#### 5.3.2 Compressible linear elasticity (CLE)

Next, we consider the parametric compressible linear elasticity system in \(d=3\) physical dimensions with a cylindrical shell along the z-axis as domain: the inner radius is \(1\), outer radius \(3\) and height \(10\), and the base centered in \(\mathbf{0}\). The \(m=12\) equations of the FS are

\[\left(\begin{array}{c}\boldsymbol{\sigma}-\mu_{1}(\nabla\cdot\mathbf{u}) \mathbb{I}_{3}-2\mu_{2}\frac{\left(\nabla\mathbf{u}+(\nabla\mathbf{u})^{t} \right)}{2}\\ -\frac{1}{2}\nabla\cdot(\boldsymbol{\sigma}+\boldsymbol{\sigma}^{t})+\mu_{3} \mathbf{u}\end{array}\right)=\left(\begin{array}{c}0\\ \mathbf{f}\end{array}\right),\quad\forall\mathbf{x}\in\Omega\subset\mathbb{R}^{ 3}, \tag{77}\]

where \(\boldsymbol{\rho}=(\mu_{1},\mu_{2},\mu_{3})\in[100,1000]^{2}\times[1,1]=\mathcal{ P}\subset\mathbb{R}^{3}\) and \(\mathbf{f}=(0,-1,0)\). The system can be rewritten as FS as in (20). We define the boundaries

\[\Gamma_{D}=\partial\Omega\cap\{z=0\},\qquad\Gamma_{N}=\partial\Omega\setminus \Gamma_{D}. \tag{78}\]

Mixed boundary conditions are applied with the boundary operator (24): homogeneous Dirichlet boundary conditions are imposed on \(\Gamma_{D}\) and homogeneous Neumann boundary conditions on \(\Gamma_{N}\).

We consider two parametric spaces:

\[\boldsymbol{\rho} =(\mu_{1},\mu_{2})\in[100,1000]^{2}=\mathcal{P}_{1}\subset\mathbb{ R}^{2},\qquad\mbox{\bf(CLE1)} \tag{79a}\] \[\boldsymbol{\rho} =(\mu_{1},\mu_{2},f_{1},f_{2})\in[100,1000]^{2}\times[-2,2]^{2}= \mathcal{P}_{2}\subset\mathbb{R}^{4},\qquad\mbox{\bf(CLE2)} \tag{79b}\]

where in the second case, the source term \(\mathbf{f}\) is now piecewise constant:

\[\mathbf{f}=\begin{cases}f_{1}\cdot(0,-1,0),&z<5,\\ f_{2}\cdot(0,-1,0),&z\geq 5.\end{cases} \tag{80}\]

We show two sample solutions for \(\mu_{1}=\mu_{2}=1000\) in Figure 8 for **CLE1** and \(\mu_{1}=\mu_{2}=1000\), \(f_{1}=1\) and \(f_{2}=-1\) for **CLE2**, on the left and on the right, respectively. The partitioned and repartitioned subdomains are shown in Figure 9. For the first case **CLE1** we employ a mesh of \(24\) cells and \(7776\) dofs, for the second **CLE2** a mesh of \(60\) cells and \(19440\) dofs.

Test case **CLE1**. This test case presents no region for which the restricted solutions are more or less approximable with a constant field, as would be detected by the variance indicator: as shown in Figure 10, the local relative \(L^{2}\)-reconstruction error in the region with _low variance_, assigned by \(I_{\mathrm{var}}\), deteriorates from the value \(2\cdot 10^{-3}\) of the abscissae \(0\%\) and \(100\%\) to \(1\cdot 10^{-2}\) of the abscissae \(4\%\). Nonetheless, despite the parametric solutions are not approximable efficiently with a constant field, they are well represented by a one dimensional linear subspace in the region located by the cellwise Grassmannian dimension indicator \(I_{G}\), for \(P_{l}=12\%\). The associated _low local Grassmannian dimension_ region for \(P_{l}=12\%\) is shown in Figure 9 in blue.

Also in this test case, the employment of DD-ROMs is not advisable, since there are little gains in the local relative \(L^{2}\)-reconstruction error for the _low local Grassmannian dimensional_ region (values around \(3\cdot 10^{-3}\), in orange for the abscissa \(P_{l}=12\%\), in Figure 10). The choice of local reduced dimensions \(r_{\Omega_{1}}=2\) and \(r_{\Omega_{2}}=3\) does not affect greatly the errors shown in Figure 11. Also in this case, we evaluate \(n_{\mathrm{train}}=20\) training full-order solutions and \(n_{\mathrm{test}}=80\) test full-order solutions,

Figure 8: **CLE. Left:** solution of the compressible linear elasticity FS **CLE1** with parameter values \(\mu_{1}=\mu_{2}=1000\). The cylindrical shell displacement \(\mathbf{u}\), and with a different colorbar also the field \(\sigma\mathbf{e}_{\mathbf{z}}\), named sigma_3, are shown. At the extremity close to \(z=0\) homogeneous Dirichlet boundary conditions are imposed. **Right:** solution of the test case **CLE2** with discontinuous values of the source terms along the computational domain \(\{z<5\}\cap\Omega\) and \(\{z\geq 5\}\cap\Omega\): \(\mu_{1}=\mu_{2}=1000\), \(f_{1}=1\) and \(f_{2}=-1\).
corresponding to a uniform independent sampling from the parametric domain \(\mathcal{P}\subset\mathbb{R}^{3}\). Also for these studies, we have fixed the local dimensions to \(r_{\Omega_{i}}=3,\ i=1,\ldots,K\) for \(K=4\), \(r_{\Omega}=3\) for the whole computational domain and \(r_{\Omega_{1}}=2,\ r_{\Omega_{2}}=3\) for the repartitioned case with \(k=2\).

Test case **CLE1**. Similarly to the previous case, we evaluate \(n_{\rm train}=20\) training full-order solutions and \(n_{\rm test}=80\) test full-order solutions, corresponding to a uniform independent sampling from the parametric domain \(\mathcal{P}_{2}\subset\mathbb{R}^{4}\). This time, if we vary the parameters \(f_{1}\) and \(f_{2}\) inside different subdomains \(\{z\geq 5\}\cap\Omega\) and \(\{z<5\}\cap\Omega\), we obtain the results shown in Figure 12. It can be seen that repartitioning \(\Omega\) in \(k=2\) DD-ROM subdomains with the local Grassmannian indicator \(I_{G}\) and \(P_{l}=50\%\) does not produce more accurate DD-ROMs compared to the case of a single reduced solution manifold for the whole computational domain and for the DD-ROM with \(k=4\). In this case, we kept the local dimension of DD-ROM repartitioned case with \(k=2\) equal \(r_{\Omega_{1}}=r_{\Omega_{2}}=3\). For this simple test case, there is not an appreciable improvement for some test parameters in the accuracy for \(k=2\) instead of \(K=4\) or a classical global linear basis ROM. The reason is that even if the parameters \(f_{1}\) and \(f_{2}\) affect different subdomains of \(\Omega\), the solutions on the whole domain are still well correlated. Differently from the previous test case **MS2** from section 5.3.1, this is a typical case for which DD-ROMs are not effective, even if the parametrization affects independently two regions of the whole domain \(\Omega\).

In Table 2, we list the computational times and speedups for a simulation with the different methods. For an error analysis with respect to the size of the reduced space, we refer to Appendix C.

Figure 10: **CLE1.** Local relative \(L^{2}\)-reconstruction errors of the snapshots matrix for elasticity equations restricted to the two subdomains of the repartitioning performed with the indicator \(I_{\rm var}\) (in red and light-blue), Definition 3, and \(I_{\rm G}\) (in orange and blue), Definition 4. The relative \(L^{2}\)-reconstruction error attained on the whole domain is shown in black for the indicator \(I_{\rm var}\) and in brown for the indicator \(I_{\rm G}\). The local reduced dimensions used to evaluate the local reconstruction errors is \(r_{\Omega_{i}}=3,\ i=1,2\).

Figure 9: **CLE. Left:** computational subdomains partitioned in \(K=4\) subdomains by \(\mathtt{petsc4py}\) inside \(\mathtt{deal.II}\). **Center:** test case **CLE1** repartition of the computational subdomain \(k=2\) with the cellwise Grassmannian dimension indicator \(I_{\rm G}\), Definition 4: \(12\%\) of the cells belong to the _low local Grassmannian dimension_ part, represented in blue inside the torus, and the other \(88\%\) belong to the _high local Grassmannian dimension_ part, represented in red. **Right:** test case **CLE2** repartition of the computational subdomain \(k=2\) with the cellwise Grassmannian dimension indicator \(I_{\rm G}\) and \(P_{l}=50\%\).
\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|c|} \hline \multicolumn{2}{|c|}{FOM} & \multicolumn{3}{c|}{ROM} & \multicolumn{3}{c|}{ROM} & \multicolumn{3}{c|}{DD-ROM} \\ \hline \(N_{h}\) & time & \(r\) & time & speedup & \(r_{i}\) & time & speedup \\ \hline \hline
7776 & 411.510 [ms] & 3 & 80.444 [\(\mu\)s] & \(\sim 5115\) & [3, 3, 3,3] & 85.108 [\(\mu\)s] & \(\sim 4835\) \\ \hline
19440 & 2.080 [s] & 3 & 69.992 [\(\mu\)s] & \(\sim 29718\) & [3, 3, 3,3] & 94.258 [\(\mu\)s] & \(\sim 22067\) \\ \hline \end{tabular}
\end{table}
Table 2: **CLE.** Average computational times and speedups for ROM and DD-ROM approaches for Maxwell equations. The speedup is computed as the FOM computational time over the ROM one. The FOM runs in parallel with 4 cores, so “FOM time” refers to wallclock time. The first row correspond to test case **CLE1**, the second to test case **CLE2**.

Figure 11: **CLE1.** Errors and estimators for elasticity equations corresponding to the \(n_{\text{train}}=20\) uniformly sampled training snapshots corresponding to the abscissae \(0,5,10,\ldots,95\), and \(n_{\text{test}}=80\) uniformly sampled test snapshots, corresponding to the other abscissae. The reduced dimensions of the ROMs are \(\{r_{\Omega_{i}}\}_{i=1}^{K}=[3,3,3,3]\) for \(K=4\) partitions, \(r_{\Omega}=3\) for \(k=1\) partition, and \(\{r_{\Omega_{i}}\}_{i=1}^{k}=[2,3]\) for \(k=2\) partitions. For the case \(k=2\) we employed the cellwise local Grassmannian dimension indicator \(I_{G}\), Definition 4, with \(P_{l}=12\%\).
Figure 12: **CLE2. Errors and estimators corresponding to the \(n_{\text{test}=80}\) uniformly sampled test snapshots corresponding to the abscissae \(0,5,10,\ldots,95\), and \(n_{\text{train}=20}\) uniformly sampled training snapshots, corresponding to the other abscissae. The reduced dimensions of the ROMs are \(\{r_{\Omega_{i}}\}_{i=1}^{K}=[3,3,3,3]\) for \(K=4\) partitions, \(r_{\Omega}=3\) for \(k=1\) partition, and \(\{r_{\Omega_{i}}\}_{i=1}^{k}=[3,3]\) for \(k=2\) partitions. For the case \(k=2\) we employed the cellwise variance dimension indicator \(I_{\text{var}}\), Definition 3, with \(P_{l}=50\%\).**
#### 5.3.3 Scalar concentration advected by an incompressible flow (ADR)

We consider the parametric semi-linear advection diffusion reaction equation in \(d=2\) dimensions, with \(m=3\) equations, rewritten in mixed form:

\[\begin{cases}\kappa^{-1}\sigma+\nabla u=0,&\text{ in }\Omega,\\ \nabla\cdot\sigma+\mathbf{v}\cdot\nabla u+u=f,&\text{ in }\Omega,\\ \sigma\cdot\mathbf{n}=0,&\text{ on }\Gamma_{N}\cup\Gamma_{D,0},\\ u=\sum_{i=1}^{P}\mu_{i}\chi_{I_{i}},&\text{ on }\Gamma_{D},\end{cases} \tag{81}\]

where \(\kappa=0.05\) is fixed for this study,

\[\boldsymbol{\rho}=(\mu_{1},\ldots,\mu_{P})\in\mathcal{P}\subset\mathbb{R}^{P}, \qquad\mathcal{P}=\{\boldsymbol{\rho}\in\{0,1\}^{P}|\mu_{i}=1,\ \mu_{j}=0,\ \forall j\in\{0,\ldots,99\}\backslash\{i\}\}, \tag{82}\]

and \(\{\chi_{I_{i}}\}_{i=0}^{N_{\text{par}}}\) are the characteristic functions of the symmetric intervals \(I_{i}=0\times[-i0.01+1.5,i0.01+2.5]\), with \(N_{\text{par}}=99\). The domain is shown in Figure 13. The advection velocity \(\mathbf{v}\) is obtained from the following incompressible Navier-Stokes equation at \(t=2\)s:

\[\begin{cases}\partial_{t}\mathbf{v}+\mathbf{v}\cdot\nabla\mathbf{v}-\nu \Delta\mathbf{v}+\nabla p=\mathbf{0},&\text{ in }\Omega\\ \nabla\cdot\mathbf{v}=0,&\text{ in }\Omega\\ \mathbf{v}\times\mathbf{n}=0,\ p=0,&\text{ on }\Gamma_{N}\\ \mathbf{v}=0,&\text{ on }\Gamma_{D,0}\\ \mathbf{v}(t=0)=\mathbf{v}_{b},&\text{ on }\Gamma_{D}\end{cases} \tag{83}\]

with initial conditions on the boundary \(\Gamma_{D}\), \(\mathbf{v}_{b}=\mathbf{v}(x,y,t=0)=(6y(4.1-y)/4.1^{2},0)\in\mathbb{R}^{2}\) and \(\nu\in\mathbb{R}\) such that the Reynolds number is \(Re=100\). The implementation is the one of step-35 of the tutorials of the deal.II library [6].

Homogeneous Neumann boundary on \(\Gamma_{N}\cup\Gamma_{D,0}\) and Dirichlet non-homogeneous boundary conditions on \(\Gamma_{D}\) are applied with the boundary operator (30). A sample solution is shown in Figure 14 for \(\mu_{i}=0,\ i=0,\ldots,98\) and \(\mu_{99}=1\), \(\kappa=0.05\). We remark that, for the moment, we consider only fixed values of \(\kappa=0.05\). For a convergence of ROMs to vanishing viscosity solutions with graph neural networks, see Section 6.

Figure 14: **ADR. Left**: scalar concentration \(u\) of the advection diffusion reaction equations (81), with \(\mu_{i}=0,\ i=0,98\) and \(\mu_{99}=1\), \(\kappa=0.05\). **Right**: advection velocity employed for the FS (81), obtained as the velocity \(\mathbf{v}\) from the INS (83) at \(t=2\)s.

Figure 13: **ADR.** Computational domain of the advection diffusion reaction equation FS (81) and the incompressible Navier-Stokes equations (83). The boundary conditions specified for each system are reported in the text.
The FOM partitioned and DD-ROM repartitioned subdomains are shown in Figure 15. We choose the variance indicator to repartition the computational subdomain in two subset: \(21\%\) of the cells for the _low variance_ part and \(79\%\) for the _high variance_ part. With respect to the previous test cases, now it is evident the change in the order of magnitude of the local relative \(L^{2}\)-reconstruction error in Figure 16, especially for the cellwise variance indicator \(I_{\text{var}}\). We expect that lowering the local reduced dimension of the _low variance_ repartitioned region will not affect sensibly the accuracy.

We use for the monodomain approach \(r_{\Omega}=5\) reduced basis as well as \(r_{\Omega_{i}}=5\) for \(i=1,\ldots,K\) for the FOM partitioned subdomains. In the DD-ROM approach, we can use even \(r_{\Omega_{1}}=2\) and \(r_{\Omega_{2}}=5\) for the lower and higher variance subdomains, respectively, without affecting the error of the ROM solution, as we see in Figure 17. Indeed, the accuracy in terms of \(L^{2}\) and energy norms is essentially identical for all approaches, even with so little number of basis functions for the DD-ROM one.

Again, we evaluate \(n_{\text{train}}=20\) training full-order solutions and \(n_{\text{test}}=80\) test full-order solutions, corresponding to the parameter choices \(\mu_{i}=1\) and \(\mu_{\bar{i}}=0\), for \(i=0,\ldots,99\), with fixed viscosity \(\kappa=0.05\), where \(\bar{i}\) represents all the indices in \(\{0,\ldots,99\}\) except from \(i\). So, the training snapshots correspond to \(i=0,5,10,\ldots,95\). For these studies we have fixed the local dimensions to \(r_{\Omega_{i}}=5,\ i=1,\ldots,K\) for \(K=4\), \(r_{\Omega}=5\) for the whole computational domain and \(r_{\Omega_{1}}=2,\ r_{\Omega_{2}}=5\) for the repartitioned case with \(k=2\), as mentioned.

In Table 3, we list the computational times and speedups for a simulation with the different methods.

## 6 Graph Neural Networks approximating Vanishing Viscosity solutions

In this section, we want to highlight how the well-known concept of vanishing viscosity solutions can be related to FS. In hyperbolic problems, the uniqueness of the weak solution is not guaranteed, already for very simple problems, e.g. inviscid Burgers' equations. In order to filter out the physically relevant solution, the concept of vanishing viscosity solution has been introduced, _inter alia_[37], and, consequently, vanishing viscosity methods have been developed, e.g. [64, 26].

Figure 16: **ADR.** Local relative \(L^{2}\)-reconstruction errors of the snapshots matrix of the advection diffusion reaction equation restricted to the two subdomains of the repartitioning performed with the indicator \(I_{\text{var}}\) (in red and light-blue), Definition 3, and \(I_{\text{G}}\) (in orange and blue), Definition 4. The relative \(L^{2}\)-reconstruction error attained on the whole domain is shown in black for the indicator \(I_{\text{var}}\) and in brown for the indicator \(I_{\text{G}}\). The local reduced dimensions used to evaluate the local reconstruction errors is \(r_{\Omega_{i}}=3,\ i=1,2\).

Figure 15: **ADR.** Domain of advection diffusion reaction equation. **Left**: computational subdomains partitioned in \(K=4\) subdomains by petsc4py inside deal.II. **Right**: DD-ROM repartition of the computational subdomain \(k=2\) with the cellwise variance indicator \(I_{\text{var}}\), Definition 3: \(21\%\) of the cells belong to the _low variance_ part, represented in blue, and the other \(79\%\) belong to the _high variance_ part in red.
\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline \multicolumn{2}{|c|}{FOM} & \multicolumn{2}{c|}{ROM} & \multicolumn{2}{c|}{ROM} & \multicolumn{2}{c|}{DD-ROM} \\ \hline \(N_{h}\) & time & \(r\) & time & speedup & \(r_{i}\) & time & speedup \\ \hline \hline
131328 & 3.243 [s] & 5 & 79.112 [\(\mu\)s] & \(\sim 40992\) & [5, 5, 5, 5] & 59.912 [\(\mu\)s] & \(\sim 54129\) \\ \hline \end{tabular}
\end{table}
Table 3: **ADR.** Average computational times and speedups for ROM and DD-ROM approaches for Maxwell equations. The speedup is computed as the FOM computational time over the ROM one. The FOM runs in parallel with 4 cores, so “FOM time” refers to wallclock time.

Figure 17: **ADR.** Errors and estimators for advection diffusion reaction equation corresponding to the \(n_{\text{train}=20}\) uniformly sampled train snapshots corresponding to the abscissae \(0,5,10,\ldots,95\), and \(n_{\text{test}=80}\) uniformly sampled test snapshots, corresponding to the other abscissae. The reduced dimensions of the ROMs are \(\{r_{\Omega_{i}}\}_{i=1}^{K}=[5,5,5,5]\) for \(K=4\) partitions, \(r_{\Omega}=5\) for \(k=1\) partition, and \(\{r_{\Omega_{i}}\}_{i=1}^{k}=[2,5]\) for \(k=2\) partitions. For the case \(k=2\) we employed the cellwise variance indicator \(I_{G}\), Definition 3, with \(P_{l}=21\%\).
We will consider the topic of vanishing viscosity solutions from the different perspective of model order reduction. It is known that slow decaying Kolmogorov n-width solution manifolds result in ineffective linear reduced order models. The origin of this problem rests theoretically on the regularity of the parameter to solution map [18, 19], and with less generality on the nature of some PDEs (e.g. advection dominated PDEs, nonlinearities, complex dynamics), on the size of the parameter space, and on the smoothness of the parametric initial data or parametric boundary conditions [5], mainly. A possible way to obtain more approximable solution manifolds is through regularization or filtering [88, 85], e.g. adding artificial viscosity. Heuristically, the objective is to smoothen out the parametric solutions of the PDEs, for example removing sharp edges, local features, complex patterns, with the aim of designing more efficient ROMs for the filtered solution manifolds. Then, the linear ROMs will be applied to different levels of _regularization_, still remaining in the regime where they have good approximation properties. Finally, the original (vanishing viscosity) solutions will be recovered with a regression method from the succession of filtered linear ROMs. This is realized without the need to directly reduce with a linear reduced manifold the original solution manifold, thus avoiding the problem of its approximability with a linear subspace and the slow Kolmogorov n-width decay.

In our case, we consider regularization by viscosity levels: the vanishing viscosity solutions \(u_{\nu}\) with viscosity \(0\leq\nu\ll 1\), will be recovered as the limit \(\lim_{i\to\infty}u_{\nu_{i}}=u_{\nu}\) of a potentially infinite succession of viscosity levels \(\{\nu_{i}\}_{i=0}^{\infty},\;\nu_{0}>\nu_{1}>\cdots>0\), each associated to its efficient reduced order model. In practice, \(\{\nu_{i}\}_{i=0}^{\infty}\approx\{\nu_{i}\}_{i=0}^{q}\), where \(q\) is the number of additional viscosity ROMs. It is clear the connection with multi-fidelity and super-resolution methods [35, 56]. The rationale of the approach is supported by the proofs of convergence to vanishing viscosity solutions of hyperbolic PDEs under various hypotheses [65, 57, 25, 38].

The framework is general and can be applied in particular to FS. We will achieve this for the advection-diffusion-reaction problem changing the viscosity constant \(\mathbb{R}\ni\kappa>0\) in (81). While this choice is specific for the model we are considering, a more general approach could consist in adding a viscous dissipative term to the generic FS obtaining another FS:

\[\begin{cases}Au=f+\mathbf{\kappa\Delta u},&\text{in }\Omega\\ (\mathcal{D}-\mathcal{M})(u-g)=0,&\text{on }\partial\Omega\end{cases} \to\begin{cases}\begin{cases}\kappa^{-1}\sigma+\nabla u=0\\ \nabla\cdot\sigma+Au=f\\ (\mathcal{D}-\mathcal{M})(u-g)=0,&\text{on }\partial\Omega,\end{cases}\end{cases} \tag{84}\]

recalling that the additional degrees of freedom are needed only for the high viscosity ROMs and FOMs (to collect the snapshots) and not the full-order vanishing viscosity solutions. This is only an example of how the procedure could be applied to any FS. In fact, the methodology is not designed specifically for FS.

The overhead of the methodology is related to the evaluation of the snapshots, the assembling of each level of viscosity \(\{\nu_{i}\}_{i=0}^{q}\), and the computational costs of the regression method. We remark that the full matrices of the affine decomposition of each \(\{ROM_{\nu_{i}}\}_{i=0}^{q}\) are the same. This is the price necessary to tackle the realization of reduced order models of parametric PDEs affected by a slow Kolmogorov n-width decay with our approach.

With respect to standard techniques for nonlinear manifold approximation, the proposed one is more interpretable as a mathematical limit of a succession of solutions to the vanishing viscosity one. Moreover, it has a faster training stage relying on the efficiency of the \(\{ROM_{\nu_{i}}\}_{i=0}^{q}\). To the authors' knowledge, cheap analytical ways to obtain the vanishing viscosity solution from a finite succession of high viscosity ones are not available, so we will rely on data-driven regression methods.

### Graph neural networks augmented with differential operators

Generally, machine learning (ML) architectures are employed in surrogate modelling to approximate nonlinear solution manifolds, otherwise linear subspaces are always preferred. The literature is vast on the subject and there are many frameworks that develop surrogate models with ML architectures. They promise to define data-driven reduced order models that infer solutions for new unseen parameters provided that there are enough data to train such architectures. This depends crucially on the choice of the encoding and inductive biases employed to represent the involved datasets: the training computational time and the amount of training data can change drastically.

On this matter, convolutional autoencoders (CNN) are one of the most efficient architectures to approximate nonlinear solution manifolds [59] for data structured on Cartesian grids, mainly thanks to their shift-equivariance property. For fields on unstructured meshes the natural choice are Graph neural networks (GNNs). Since their employment, GNNs architectures from the ML community have been enriched with physical inductive biases and other tools from numerical analysis. We want to test one of the first implementations and modifications of GNNs [81]. We also want to remark that in the literature, there are still very few test cases of ROMs that employ GNNs with more than \(\geq 50000\) degrees of freedom. The difficulty arises when the training is performed on large meshes, thus the need for tailored approaches.

The majority of GNNs employed for surrogate modelling are included in autoencoders [33, 68] or are directly parametrized to infer the unseen solution with a forward evaluation. These architectures may become heavy, especially for non-academic test cases. One way to tackle the problem of parametric model order reduction of slow Kolmogorov n-width solution manifolds is to employ GNNs only to recover the high-fidelity solution in a multi-fidelity setting, through super-resolution. Since efficient
ROMs are employed to obtain the lower levels of fidelity (high viscosity solutions in our case), the solution manifold dimension reduction is performed only at those levels, avoiding the costly and heavy in memory training of autoencoders of GNNs.

We describe the implementation of augmented GNNs as in [81], with the difference that we need to train only a map from a collection of DD-ROMs solutions to the full-order vanishing viscosity solution, and not an autoencoder with pooling and unpooling layers to perform dimension reduction. The GNN we will employ is rather thin with respect to autoencoder GNNs used to perform dimension reduction. Its details are reported in Table 4.

We represent with

\[\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{W}),\quad\mathcal{V}\in\mathbb{ R}^{n_{\text{nodes}}\times f},\quad\mathcal{E}\in\mathbb{R}^{n_{\text{edges}}\times 2}, \quad\mathcal{W}\in\mathbb{R}^{n_{\text{attr}}\times d}, \tag{85}\]

a graph with node features \(\mathcal{V}\), edges \(\mathcal{E}\) and edge attributes \(\mathcal{W}\). The number \(f\) represents the nodal features dimension. We denote with \(\mathbf{e}_{ij}=(i,j)\in\mathbb{N}^{2}\) the edge between the nodes \(\mathbf{n}_{i},\mathbf{n}_{j}\in\mathbb{R}^{f}\): \(\mathbf{e}_{ij}\) corresponds to a row of \(\mathcal{E}\), and \(\mathbf{n}_{i},\mathbf{n}_{j}\) correspond to the \(i\)-th and \(j\)-th rows of \(\mathcal{V}\), for \(i,j=1,\ldots,n_{\text{nodes}}\). Similarly, \(\boldsymbol{\omega}_{ij}\) represents the edge attributes of edge \(\mathbf{e}_{ij}\). We have \(n_{\text{edges}}=n_{\text{attr}}\). For their efficiency, GNNs rely on a message passing scheme composed of propagation and aggregation steps. Supposing that the graph is sparsely connected their implementation is efficient.

When the graph is supported on a mesh, it is natural to consider the generalized support points of finite element spaces as nodes of the graph and the sparsity pattern of the linear system associated to the numerical model as the adjacency matrix of the graph. We employ only Lagrangian nodal basis of discontinuous finite element spaces, but the framework can be applied to more general finite element spaces. As edge attributes \(\boldsymbol{\omega}_{ij}\), we will employ the difference \(\boldsymbol{\omega}_{ij}=\mathbf{x}_{i}-\mathbf{x}_{j}\in\mathbb{R}^{d}\) between the corresponding spatial coordinates associated to the nodes \(\mathbf{n}_{i},\mathbf{n}_{j}\in\mathbb{R}^{f}\). The nodes adjacent to node \(\mathbf{n}_{i}\) are represented with the set \(\mathcal{N}_{\text{neigh}(i)}\) for all \(i=1,\ldots,n_{\text{nodes}}\).

We consider only the two following types of GNN layers: a continuous kernel-based convolutional operator \(l_{\text{N}\text{Nconv}}\)[36, 78] and the GraphSAGE operator \(l_{\text{SAGEconv}}\)[40],

\[\mathcal{V}_{\text{out}} =l_{N\text{N}\text{conv}}(\mathcal{V}_{\text{inp}},\mathcal{E}, \mathcal{W})=\mathcal{V}_{\text{inp}}W_{3}+\text{Avg}_{1}(\mathcal{V}_{\text{inp }},h(\mathcal{W}))+\mathbf{b}_{3},\quad h(\mathcal{W})=\text{ReLU}(\mathcal{W} W_{1}+\mathbf{b_{1}})W_{2}+\mathbf{b_{2}}, \tag{86}\] \[\mathcal{V}_{\text{out}} =l_{SAGE\text{conv}}(\mathcal{V}_{\text{inp}},\mathcal{E}, \mathcal{W})=\mathcal{V}_{\text{inp}}W_{6}+\text{Avg}_{2}(\text{ReLU}( \mathcal{V}_{\text{inp}}W_{4}+\mathbf{b_{4}}))W_{5}+\mathbf{b}_{5}, \tag{87}\]

with weight matrices dimensions,

\[W_{1}\in\mathbb{R}^{2\times l},\ W_{2}\in\mathbb{R}^{l\times(f _{\text{inp}}\times f_{\text{out}})},\ W_{3}\in\mathbb{R}^{f_{\text{inp}} \times f_{\text{out}}},\ W_{4}\in\mathbb{R}^{f_{\text{inp}}\times f_{\text{inp }}},\ W_{5},W_{6}\in\mathbb{R}^{f_{\text{inp}}\times f_{\text{out}}}, \tag{88}\] \[\mathbf{b}_{1}\in\mathbb{R}^{l},\ \mathbf{b}_{2}\in\mathbb{R}^{(f_{ \text{inp}}\times f_{\text{out}})},\ \mathbf{b}_{3},\mathbf{b}_{5}\in\mathbb{R}^{f_{\text{out}}},\ \mathbf{b}_{5}\in\mathbb{R}^{f_{\text{inp}}},\] (89) \[h(\mathcal{W})\in\mathbb{R}^{n_{\text{edges}}\times(f_{\text{inp }}\times f_{\text{out}})},\quad\mathcal{W}=\{W_{s}^{h}\}_{s=1}^{n_{\text{ enders}}},\ W_{s}^{h}\in\mathbb{R}^{f_{\text{inp}}\times f_{\text{out},s}},\ \forall s=1,\ldots,n_{\text{edges}}, \tag{90}\]

with the following average operators used as aggregation operators,

\[(\text{Avg}_{1}(\mathcal{V},\{W_{s}^{h}\}_{s=1}^{n_{\text{edges}}}))_{i}=\frac {1}{\mathcal{N}_{\text{neigh}(i)}}\sum_{s\in\mathcal{N}_{\text{neigh}(i)}}W_{s }^{h}\mathbf{n}_{s},\qquad(\text{Avg}_{2}(\mathcal{V}))_{i}=\frac{1}{\mathcal{N }_{\text{neigh}(i)}}\sum_{s\in\mathcal{N}_{\text{neigh}(i)}}\mathbf{n}_{s}, \qquad\forall i=1,\ldots,n_{\text{nodes}}, \tag{91}\]

where \(\mathcal{V}_{\text{inp}}\in\mathbb{R}^{n_{\text{nodes}}\times f_{\text{inp}}}, \mathcal{V}_{\text{out}}\in\mathbb{R}^{n_{\text{nodes}}\times f_{\text{out}}}\) are the input and output nodes with feature dimensions \(f_{\text{inp}},f_{\text{out}}\). We remark that, differently from graph neural networks with heterogeneous layers, i.e., with changing mesh structure between different layers, in this network the edges \(\mathcal{E}\) and edge attributes \(\mathcal{W}\) are kept fixed, only the node features change. The feed-forward neural network \(h:\mathbb{R}^{n_{\text{edge}}\times d}\rightarrow\mathbb{R}^{f_{\text{inp}}, \times f_{\text{out},s}}\) defines a weight matrix \(W_{s}^{h}\in\mathbb{R}^{f_{\text{inp}},\times f_{\text{out},s}}\) for each edge \(s=1,\ldots,n_{\text{edges}}\). The number \(l\) is the hidden layer dimension of \(h\).

The aggregation operators are defined from the edges \(\mathcal{E}\) that are related to the sparsity pattern of the linear system of the numerical model. So, the aggregation is performed on the stencils of the numerical scheme chosen for every layer of the GNN architecture in Table 4. Many variants are possible, in particular, we do not employ pooling and unpooling layers to move from different meshes: we always consider the same adapted mesh.

Since our GNNs work on the nodal features, a good strategy is to augment their dimensions as proposed in [81]. In fact, in the majority of applications of GNNs for physical models the input features dimensions is the dimension of the physical fields considered and it is usually very small. Considering FS, the fields' dimension is \(m\). To augment the input features, we will filter them with some differential operators discretized on the same mesh in which the GNN is supported. We consider the following differential operators

\[\Delta:V_{h}(\Omega) \to V_{h}(\Omega),\] (Laplace operator) (92) \[\mathbf{v}\cdot\nabla:V_{h}(\Omega) \to V_{h}(\Omega),\] (Advection operator) (93) \[\nabla_{x}:V_{h}(\Omega) \to V_{h}(\Omega),\] (Gradient x-component) (94) \[\nabla_{y}:V_{h}(\Omega) \to V_{h}(\Omega),\] (Gradient y-component) (95)

for a total of four possible feature augmentation operators, where, in our case, \(\mathbf{v}\) is the advection velocity from the incompressible Navier-Stokes equations (83). We employ the representation of the previous differential operators with respect to
the polynomial basis of Lagrangian shape functions, so they act on the vectors of nodal evaluations in \(\mathbb{R}^{N_{h}}\). As in [81], we consider three sets of possible augmentations:

\[\mathcal{O}_{1} =\{\mathbb{I}_{N_{h}},\Delta,\mathbf{v}\cdot\nabla,\nabla_{x},\nabla _{y}\}, \tag{96}\] \[\mathcal{O}_{2} =\{\mathbb{I}_{N_{h}},\nabla_{x},\nabla_{y}\}\] (97) \[\mathcal{O}_{3} =\{\mathbb{I}_{N_{h}}\} \tag{98}\]

where \(\mathbb{I}_{N_{h}}\) is the identity matrix in \(\mathbb{R}^{N_{h}}\), \(|\mathcal{O}_{1}|=5=n_{\text{aug}}\), \(|\mathcal{O}_{2}|=3=n_{\text{aug}}\) and \(|\mathcal{O}_{3}|=1=n_{\text{aug}}\). We will reconstruct only the scalar concentration \(u\) with the GNN, so, in our case, the field dimension is 1, which is the output dimension. The input dimension depends on the number of high viscosity DD-ROMs employed that we denote with \(q\). Given a single parametric instance \(\boldsymbol{\rho}\in\mathbb{R}^{P}\) the associated solutions of \(\{\text{D-ROM}_{\kappa_{1}}\}_{i=1}^{q}\) are \(\{\mathbf{u}^{\text{RB}}(\boldsymbol{\rho}_{i})\}_{i=1}^{q}\).

We divide the snapshots \(\{\mathbf{u}^{\text{RB}}(\boldsymbol{\rho}_{i})\}_{i=1}^{n_{\text{train}}+n_{ \text{test}}}\) in training \(\{\mathbf{u}^{\text{RB}}(\boldsymbol{\rho}_{i})\}_{i\in I_{\text{train}}}\) and test snapshots \(\{\mathbf{u}^{\text{RB}}(\boldsymbol{\rho}_{i})\}_{i\in I_{\text{test}}}\), with \(|I_{\text{train}}|=n_{\text{train}}\) and \(|I_{\text{test}}|=n_{\text{test}}\). We have decided to encode the reconstruction of the vanishing viscosity solution \(\mathbf{u}_{q+1}\) learning the difference \(\mathbf{u}_{q+1}(\boldsymbol{\rho})-\mathbf{u}_{q}^{\text{RB}}(\boldsymbol{ \rho})-\overline{\mathbf{u}_{q+1}}\text{train}\) with the mesh-supported-augmented GNN (MSA-GNN) described in Table 4:

\[\mathbf{u}_{q+1}(\boldsymbol{\rho})=ReLU\left(\mathbf{u}_{q}^{\text{RB}}( \boldsymbol{\rho})+\text{MSA-GNN}\left(\{O_{a}\{\mathbf{u}_{i}^{\text{RB}}( \boldsymbol{\rho})\}_{i=1}^{q},\{\mathbf{u}_{i}^{\text{RB}}(\boldsymbol{\rho} )-\mathbf{u}_{i-1}^{\text{RB}}(\boldsymbol{\rho})\}_{i=1}^{q-1}\}\right)_{a=2 }^{n_{\text{aug}}}\right)+\overline{\mathbf{u}_{q+1}}_{\text{train}}\right), \tag{99}\]

where \(\overline{\mathbf{u}_{q+1}}_{\text{train}}=\frac{1}{n_{\text{train}}}\sum_{i=1 }^{n_{\text{train}}}\mathbf{u}_{q+1}(\boldsymbol{\rho}_{i})\). Learning the difference instead of the solution itself helps in getting more informative features. The input dimension is therefore \(3n_{\text{aug}}=15\) for \(\mathcal{O}_{1}\) and \(3n_{\text{aug}}=9\) for \(\mathcal{O}_{2}\).

### Decomposable ROMs approximating vanishing viscosity (VV) solution through GNNs

In this section, we test the proposed multi-fidelity approach that reconstruct the lowest viscosity level with the GNN. We consider the FS (81), with three levels of viscosity, from highest to lowest: \(\kappa_{1}=0.05\), \(\kappa_{2}=0.01\) and \(\kappa_{3}=0.0005\). We want to build a surrogate model that efficiently predicts the parametric solutions of the FS (81) for unseen values of \(\boldsymbol{\rho}\in\mathbb{P}\) with fixed viscosity \(\kappa_{3}=0.0005\). These solutions will be referred to as _vanishing viscosity_ solutions. The other two viscosity levels are employed to build the D-ROM\({}_{\kappa_{1}}\) and D-ROM\({}_{\kappa_{2}}\) with viscosities \(\kappa_{1}=0.05\) and \(\kappa_{2}=0.01\), respectively. The parametrization affects the inflow boundary condition and is the same as the one described in section 5.3.3, see equation (82). We also employ the same number of training 20 and test 80 parameters.

The DD-ROMs provided for \(\kappa_{1}=0.05\) and \(\kappa_{2}=0.01\) can be efficiently designed with reduced dimensions \(\{r_{\Omega_{i}}\}_{i=1}^{K}[5,5,5,5]\). To further reduce the cost, we employ an even coarser mesh for ROM\({}_{\kappa_{1}}\) and ROM\({}_{\kappa_{2}}\) and a finer mesh for the vanishing viscosity solutions. The former is represented on the left of Figure 18, the latter on the right. The degrees of freedom related to the coarse mesh are 43776, while the ones on the fine one are 175104.

For the training of the GNN we use the open source software library PyTorch Geometric [32]. The employment of efficient samplers that partition the graphs on which the training set is supported is crucial to lower the otherwise heavy memory burden [40]. We preferred samplers that partition the mesh with METIS [54] as it is often employed in this context. We decided to train the GNN with early stopping at 50 epochs as our focus is also in the reduction of the training time of the NN architectures used for model order reduction. It corresponds on average to less than 60 minutes of training time. The batch size is 100 and we clustered the whole domain in 100000 subgraphs in order to fit the batches in our limited GPU memory. Each augmentation strategy and additional fidelity level, do not affect the whole training time as they only

\begin{table}
\begin{tabular}{l|c|c|c} \hline \hline Net & Weights \([f_{\text{inp}},f_{\text{out}}]\) & Aggregation & Activation \\ \hline \hline Input NNConv & \([3n_{\text{aug}},18]\) & Avg\({}_{1}\) & ReLU \\ \hline SAGEconv & [18, 21] & Avg\({}_{2}\) & ReLU \\ \hline SAGEconv & [21, 24] & Avg\({}_{2}\) & ReLU \\ \hline SAGEconv & [24, 27] & Avg\({}_{2}\) & ReLU \\ \hline SAGEconv & [27, 30] & Avg\({}_{2}\) & ReLU \\ \hline Output NNConv & [30, 1] & Avg\({}_{1}\) & - \\ \hline \hline NNConvFilters & First Layer \([2,l]\) & Activation & Second Layer \([l,f_{\text{inp}}f_{\text{out}}]\) \\ \hline \hline Input NNConv & [2, 12] & ReLU & [12, \(3n_{\text{aug}}\cdot 18\)] \\ \hline Output NNConv & [2, 8] & ReLU & [8, 30] \\ \hline \hline \end{tabular}
\end{table}
Table 4: Mesh supported augmented GNN
increase the dimension of the input features from a minimum of \(1\) (\(1\) fidelity, no augmentation) to a maximum of \(15\) (all augmentations \(\mathcal{O}_{3}\), \(2\) fidelities). As optimizer we use ADAM [55] stochastic optimizer. Every architecture is trained on a single GPU NVIDIA Quadro RTX 4000.

Figures 19, 20 and 21 show the results of the algorithm for parameters with index \(i\in\{0,50,99\}\). In particular, we show on the left columns the FOM simulations, in the center column the ROM simulations and the error in the right column. Moreover, in the different rows, we have different viscosity levels. The first three rows use the classical DD-ROM approach. We can immediately see that the vanishing viscosity \(\kappa=\kappa_{3}=0.0005\) level shows strong numerical oscillations along the whole solution, which are not present in the FOM method. This phenomenon is observable also for higher viscosity levels but it is less pronounced and concentrated on the left of the domain, where the discontinuity are imposed as boundary conditions (see error plots). Finally, in the last row, we show the results of the GNN approach, which uses the first two viscosity levels to predict the vanishing viscosity one. Contrary the DD-ROM, we do not observe many numerical oscillations in the reduced solutions and they are much more physically meaningful. Thinking about extending this approach for more complicated problems, as Euler's equations, one could guarantee the presence of the correct amount of shocks and the right location or maintaining the positivity of density and pressure close to discontinuities.

In Figure 22, we show a quantitative measure of the error of the reduced approaches presented in terms of relative \(L^{2}\) error. Overall, we can immediately see that the new GNN approach can always reach errors of the order of \(1-2\%\) for the

Figure 19: **VV.** Scalar concentration advected by incompressible flow for \(i=0\). Comparison of ROM approach at different viscosity levels \(\kappa\in\{0.05,0.01,0.0005\}\) and GNN for \(\kappa=0.0005\). FOMs on the left, reduced solution at the center and error on the right.

Figure 18: **VV.** Left part of the computational domain partitioned in \(4\) for distributed parallelism: coarse mesh (**left**), fine mesh (**right**). The solution with viscosity \(\kappa\in\{\mathbf{0.05,0.01}\}\) are evaluated on the coarse mesh with \(4868\) cells and \(43776\) dofs, those with \(\kappa=\mathbf{0.0005}\) on the finer with \(19456\) cells and \(175104\) dofs.
Figure 21: **VV.** Scalar concentration advected by incompressible flow for \(i=99\). Comparison of ROM approach at different viscosity levels \(\kappa\in\{0.05,0.01,0.0005\}\) and GNN for \(\kappa=0.0005\). FOMs on the left, reduced solution at the center and error on the right.

Figure 20: **VV.** Scalar concentration advected by incompressible flow for \(i=50\). Comparison of ROM approach at different viscosity levels \(\kappa\in\{0.05,0.01,0.0005\}\) and GNN for \(\kappa=0.0005\). FOMs on the left, reduced solution at the center and error on the right.
vanishing viscosity solutions, with few peaks in the extrapolatory regime of \(8\%\), while the classical DD-ROM on the vanishing viscosity solutions perform worse, with errors around \(6\)-\(10\%\). On the other hand, the DD-ROM for higher viscosity levels have lower errors around \(3\%\) for \(\kappa_{2}\) and \(0.5\%\) for \(\kappa_{1}\), hence, they are still reliably representing those solutions.

On the different GNN approaches, in Figure 22 at the top we compare the different augmentations \(\mathcal{O}_{1}\), \(\mathcal{O}_{2}\) and \(\mathcal{O}_{3}\) and how many levels of viscosity we keep into considerations to derive the vanishing viscosity solution. The usage of multiple fidelity levels (two viscosity levels) is a great improvement for all the augmentations proposed and it can make gain a factor of \(2\) in terms of accuracy. There are slight differences with the used augmentations and, in particular, we observe that the \(\mathcal{O}_{1}\) augmentation, with all operators, guarantee better performance, while there are no appreciable differences between \(\mathcal{O}_{2}\) and \(\mathcal{O}_{3}\). Clearly, one could come up with many other augmentation possibilities choosing more operators, but at a cost of increasing the dimensions of the GNN and the offline training costs. We believe that all the presented options already perform much better with respect to classical approaches and can already be used without further changes.

In Table 5, we compare the computational times necessary to compute the FOM solutions, the DD-ROM ones, the training time for the GNN and the online costs of the GNN. As mentioned before, we employ only one GPU NVIDIA Quadro

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline  & \multicolumn{2}{c|}{FOM} & \multicolumn{4}{c|}{DD-ROM} \\ \hline \(\kappa\) & \(N_{h}\) & time & \(r_{i}\) & time & speedup & mean \(L^{2}\) error \\ \hline \hline
0.05 & 43776 & 3.243 [s] & [5, 5, 5, 5] & 59.912 [\(\mu\)s] & 54129 & 0.00595 \\ \hline
0.01 & 43776 & 3.236 [s] & [5, 5, 5, 5] & 79.798 [\(\mu\)s] & 40552 & 0.0235 \\ \hline
0.0005 & 175104 & 9.668 [s] & [5, 5, 5, 5] & 95.844 [\(\mu\)s] & 100872 & 0.0796 \\ \hline \end{tabular} 
\begin{tabular}{|c|c|c|c|c|c|c|} \hline \(\kappa\) & GNN training time & Single forward GNN online time & Total online time & GNN speedup & mean \(L^{2}\) error \\ \hline \hline
0.0005 & \(\leq 60\) [min] & \multicolumn{4}{c|}{2.661 [s]} & \multicolumn{4}{c|}{17.166 [s]} & \(\sim 56\) & 0.0217 \\ \hline \end{tabular}
\end{table}
Table 5: **VV.** Computational costs for scalar advected by a incompressible flow problem with GNNs approximating vanishing viscosity solutions (VV). The speedup is computed as the FOMs computational time over the ROM one. The speedup of the GNN is with respect to the FOM with viscosity \(\nu=0.0005\). The FOM runs in parallel with \(K=4\) cores as the DD-ROMs, so “FOM time” and “DD-ROM time” refers to wallclock time. Regarding the GNN results, “Single forward GNN online time” refers to a single online evaluation while “Total online time” refers to the evaluation of the 100 training and test snapshots altogether with only two separate GNN forward evaluations with batches of 50 inputs each. The speedup is evaluated as “FOM time” over ”Total online time” divided by 100.

Figure 22: **VV.** Relative errors for the scalar conservation advected by incompressible flow problem. The parameters corresponding to the snapshots used for the GNNs and DD-ROMs training correspond to the abscissae \(0,5,10,\ldots,95\) the rest are test parameters. The dashed red background highlights the extrapolation range. **Top:** errors on train and test set with different GNN approaches given by the three augmentation \(\mathcal{O}_{1}\), \(\mathcal{O}_{2}\) and \(\mathcal{O}_{3}\) and by using either \(1\) viscosity level (\(1\) fidelity) or \(2\) (all fidelities) and errors for DD-ROM with the same viscosity level \(\nu=0.0005\). **Bottom:** errors for DD-ROM approaches at different viscosity levels. The reduced dimensions of the ROMs are \(\{r_{\Omega_{i}}\}_{i=1}^{K}=[5,5,5,5]\) with \(K=4\) partitions.
RTX 4000 with 8GB of memory. Typical GNNs applications that involve autoencoders to perform nonlinear dimension reduction are much heavier. The training time of the GNNs for the different choice of augmentation operators vary between 48 minutes and 60 minutes approximately. We believe that in the near future more optimized implementations will reduce the training costs of GNNs. The computational time of the evaluation of a single forward of the GNN is on average 2.661 seconds but vectorization ensures the evaluation of multiple online solutions altogether: with our limited memory budget we could predict all the 100 training and test snapshots with just 2 batches of 50 stacked inputs each. The "Total online time" computed as previously described is 17.166 seconds that is 171.66 milliseconds per online solution with a speedup of around 56 with respect to the 9.668 seconds for the FOM.

Although the speedup for the GNN simulations are not as remarkable as for the DD-ROM, we want to highlight that the accuracy of the GNN solutions are qualitatively much better than the DD-ROM for that viscosity level, and physically more meaningful. This aspect is a major advantage with respect to classical linear ROMs that is probably worth the loss of computational advantage. In perspective, when dealing with nonlinear and more expensive FOM for different equations, the GNN approach will not require any extra computational costs, while FOMs and ROMs model might need special treatments for the nonlinearity that would make their costs increase.

## 7 Conclusions

We argue that Friedrichs' systems represent a valuable framework to study and devise reduced order models of many parametric PDEs at the same time: among them the ones studied in this work and others, like mixed elliptic and hyperbolic problems, complex and time-dependent FS and also nonlinear PDEs whose linearization results in FS, e.g. the Euler equations. The advantages include the availability of _a posteriori_ error estimators and the easy to preserve mathematical properties of positivity and symmetry from the full-order formulations to the reduced-order ones. We underlined in section 4.2 how optimally stable reduced-order models can be obtained from the ultraweak formulation. A more efficient numerical solver for Friedrichs' systems is the hybridized discontinuous Galerkin method [17]. These are possible future directions of research.

Working with discontinuous Galerkin discretizations is not only crucial from the possibly mixed elliptic and hyperbolic nature of Friedrichs' systems, but also to design domain decomposable reduced-order models with a minimum effort: in fact, penalties at the subdomains interfaces are inherited directly from the full-order models. We demonstrated with numerical experiments the limits and the ranges of application of domain decomposable ROMs: generally, with respect to single domain ROMs, there are benefits only when the model under study is truly decomposable, that is when the parameters affect independently different subdomains and the respective solutions are poorly correlated for unseen parametric instances. The results we showed in our academic benchmarks were obtained with the aim to tackle more complex multi-physics models like fluid-structure interaction systems. A typical application of DD-ROMs for FS is represented by parametric PDEs with a mixed elliptic and hyperbolic nature and possibly solution manifolds more and less linearly approximable respectively. The repartitioning strategies we developed are suited to adapt the reduced local dimension of the linear approximants, especially when the parameters influence only a limited region like in test case ADR 5.3.3. The implementation of _ad hoc_ physics inspired indicators can be a future direction of research.

The Friedrichs' systems formulation itself does not solve the problems caused by a slow decaying Kolmogorov n-width. DD-ROMs can help in this regard, isolating regions with a slow Kolmogorov n-width for which nonlinear approximants can be employed and regions with a fast decaying Kolmogorov n-width for which classical linear projection-based ROMs provide efficient and reliable predictions. Related to this subject and motivated also by the heavy computational resources that graph neural networks require when employed for model order reduction, we introduced a new paradigm for surrogate modeling: the inference with GNNs of vanishing viscosity solutions from a succession of higher viscosity projection-based ROMs. The approach is, of course, general and can be applied to PDEs that are not FS. The crucial hypotheses underneath this approach is the approximability with linear spaces of the solution manifolds corresponding to higher viscosity levels. We showed that the additional computational costs are not too large in our test case in section 6. Possible directions of research include more complex problems and different regularization or filtering choices, other than additional viscous terms.
Transformation into dissipative system

In some cases, the term \(\mathcal{A}^{0}=0\) or property (2b) is not satisfied, but there is a way to recover the previous framework. We want to recover a dissipative [61] or accretive system [52]. For example the linearized Euler equations in entropy variables [79] have \(\mathcal{A}^{0}=0\).

The condition of uniform positive definiteness

\[\exists\mu_{0},\quad\mathcal{A}^{0}+(\mathcal{A}^{0})^{t}-\mathcal{X}\geq 2 \mu_{0}\mathbb{I}_{m}\text{ a.e. in }\Omega, \tag{100}\]

is still valid if there exist \(\boldsymbol{\xi}\in\mathbb{R}^{d},\ \|\boldsymbol{\xi}\|=1\) and \(\beta\in\mathbb{R},\ \beta>0\) such that after the transformation

\[v(x)=e^{-\beta\xi\cdot x}z(x), \tag{101}\]

the resulting system

\[\sum_{i}A^{i}\partial_{i}v(x)+\beta\sum_{i=1}^{d}\xi_{i}A^{i}v(x)=e^{-\beta( \xi\cdot x)}f, \tag{102}\]

satisfies, with the newly found \(A^{0}=\beta\sum_{i=1}^{d}\xi_{i}A^{i}\),

\[\exists\mu_{0},\quad A^{0}+(A^{0})^{t}-\mathcal{X}=2\beta\sum_{i=1}^{d}\xi_{i }A^{i}-\mathcal{X}\geq 2\mu_{0}\mathbb{I}_{m}\text{ a.e. in }\Omega. \tag{103}\]

In some cases, such \(\xi\) and \(\beta\) exist, for example if the symmetric matrix \(\sum_{i=1}^{d}\xi_{i}A^{i}\) has at least one positive eigenvalue for some \(\xi\) for almost every \(x\in\Omega\), then taking \(\beta\) sufficiently large is enough to satisfy the condition. It is also sufficient that \(\sum_{i=1}^{d}\xi_{i}(x)A^{i}(x)\) has at least a positive eigenvalue for almost every \(x\in\Omega\) where \(\xi=\xi(x)\), see [52, Example 28].

_Remark 2_.: A more general transformation is

\[v(x)=w(x)z(x), \tag{104}\]

so that the positive definiteness condition becomes

\[\exists\mu_{0},\quad\mathcal{A}^{0}+(\mathcal{A}^{0})^{t}-\mathcal{X}=2\sum_ {i=1}^{d}\partial_{i}(-\log w)\mathcal{A}^{i}-\mathcal{X}\geq 2\mu_{0} \mathbb{I}_{m}\text{ a.e. in }\Omega. \tag{105}\]

## Appendix B Constructive method to define boundary operators

We report a procedure to define a boundary operator \(M\in\mathcal{L}(V,V^{\prime})\) starting from some specified boundary conditions. We exploit Theorem 4.3, Lemma 4.4 and Corollary 4.1 from [31]. It can be seen that the most common Dirichlet, Neumann and Robin boundary conditions can be found for some FS [28, 29, 24], following this procedure.

**Lemma 1** (Theorem 4.3, Lemma 4.4 and Corollary 4.1 from [31]).: _Let us assume that \((V_{0},V_{0}^{*})\) satisfy (10) and that \(V_{0}+V_{0}^{*}\subset V\) is closed. We denote with \(P:V\to V_{0}\) and \(Q:V:\to V_{0}^{*}\) the projectors onto the subspaces \(V_{0}\subset V\) and \(V_{0}^{*}\subset V\) of the Hilbert space \(V\), respectively. Then, the boundary operator \(\mathcal{M}\in\mathcal{L}(V,V^{\prime})\) defined as_

\[\begin{split}\langle\mathcal{M}u,v\rangle_{V^{\prime},V}& =\langle\mathcal{D}Pu,Pv\rangle_{V^{\prime},V}-\langle\mathcal{ D}Qu,Qv\rangle_{V^{\prime},V}+\\ &\langle\mathcal{D}(P+Q-PQ)u,v\rangle_{V^{\prime},V}-\langle \mathcal{D}u,(P+Q-PQ)v\rangle_{V^{\prime},V}\end{split} \tag{106}\]

_is admissible and satisfies \(V_{0}=\ker(\mathcal{D}-\mathcal{M})\) and \(V_{0}^{*}=\ker(\mathcal{D}+\mathcal{M}^{*})\). In particular,_

1. _If_ \(V=V_{0}+V_{0}^{*}\)_, then_ \(\mathcal{M}\) _is self-adjoint and_ \[\langle\mathcal{M}u,v\rangle_{V^{\prime},V}=\langle\mathcal{D}Pu,Pv\rangle_{V^{ \prime},V}-\langle\mathcal{D}Qu,Qv\rangle_{V^{\prime},V}\,.\] (107)
2. _If_ \(V_{0}=V_{0}^{*}\)_, then_ \(\mathcal{M}\) _is skew-symmetric and_ \[\langle\mathcal{M}u,v\rangle_{V^{\prime},V}=\langle\mathcal{D}Pu,v\rangle_{V ^{\prime},V}-\langle\mathcal{D}Pv,u\rangle_{V^{\prime},V}\,.\] (108)

We remark that, for fixed \((V_{0},V_{0}^{*})\), admissible boundary operators \(\mathcal{M}\in\mathcal{L}(V,V^{\prime})\) that satisfy \(V_{0}=\ker(\mathcal{D}-\mathcal{M})\) and \(V_{0}^{*}=\ker(\mathcal{D}+\mathcal{M}^{*})\) are not unique. The boundary operator defined in Lemma 1 is just a possible explicit definition, in general.
As an exercise, we show how to find the definition of the operator \(\mathcal{M}\) for our linear compressible elasticity FS, from Section 2.1.2. We want to impose the boundary conditions \(\mathbf{u}_{|\Gamma_{D}}=0\) and \((\boldsymbol{\sigma}\cdot\mathbf{n})_{|\Gamma_{N}}=0\), so,

\[V_{0}=V_{0}^{*}=\{(\mathbf{u},\boldsymbol{\sigma})\in V\mid\mathbf{u}_{|\Gamma _{D}}=0,\quad(\boldsymbol{\sigma}\cdot\mathbf{n})_{|\Gamma_{N}}=0\}=H_{ \boldsymbol{\sigma},\Gamma_{N}}\times[H_{\Gamma_{D}}^{1}(\Omega)]^{d}, \tag{109}\]

since we defined \(V=H_{\boldsymbol{\sigma}}\times[H^{1}(\Omega)]^{d}\), with \(H_{\boldsymbol{\sigma}}=\{\boldsymbol{\sigma}\in[L^{2}(\Omega)]^{d\times d}\mid \nabla\cdot(\boldsymbol{\sigma}+\boldsymbol{\sigma}^{t})\in[L^{2}(\Omega)]^{d}\}\), the traces \(\gamma_{D}:[H^{1}(\Omega)]^{d}\to[H^{\frac{1}{2}}(\Gamma_{D})]^{d}\) and \(\gamma_{N}:H_{\boldsymbol{\sigma}}^{d}\to[H^{-\frac{1}{2}}(\Gamma_{N})]^{d}\) on \(\Gamma_{D}\) and \(\Gamma_{N}\) are well-defined. In particular,

\[H_{\boldsymbol{\sigma},\Gamma_{N}}=\{\boldsymbol{\sigma}\in H_{\boldsymbol{ \sigma}}\mid\gamma_{N}(\boldsymbol{\sigma})=(\boldsymbol{\sigma}\cdot \mathbf{n})_{|\Gamma_{N}}=0\},\qquad[H_{\Gamma_{D}}^{1}(\Omega)]^{d}=\{ \mathbf{u}\in[H^{1}(\Omega)]^{d}\mid\gamma_{D}(\mathbf{u})=\mathbf{u}_{| \Gamma_{D}}=0\}. \tag{110}\]

Moreover, \((V_{0},V_{0}^{*})\) satisfy the properties of cone formalism (10). Thus, we can use the definition (108) of Lemma 1:

\[\begin{split}\langle\mathcal{M}(\boldsymbol{\sigma},\mathbf{u}), (\boldsymbol{\tau},\mathbf{v})\rangle_{V^{\prime},V}=&\langle \mathcal{D}P(\boldsymbol{\sigma},\mathbf{u}),(\boldsymbol{\tau},\mathbf{v}) \rangle_{V^{\prime},V}-\langle\mathcal{D}P(\boldsymbol{\tau},\mathbf{v}),( \boldsymbol{\sigma},\mathbf{u})\rangle_{V^{\prime},V}\\ =&-\langle\frac{1}{2}(\boldsymbol{\sigma}+ \boldsymbol{\sigma}^{t})\cdot\mathbf{n},\mathbf{v}\rangle_{-\frac{1}{2},\frac {1}{2},\Gamma_{D}}+\langle\frac{1}{2}(\boldsymbol{\tau}+\boldsymbol{\tau}^{t} )\cdot\mathbf{n},\mathbf{u}\rangle_{-\frac{1}{2},\frac{1}{2},\Gamma_{D}}+\\ &\langle\frac{1}{2}(\boldsymbol{\sigma}+\boldsymbol{\sigma}^{t} )\cdot\mathbf{n},\mathbf{v}\rangle_{-\frac{1}{2},\frac{1}{2},\Gamma_{N}}- \langle\frac{1}{2}(\boldsymbol{\tau}+\boldsymbol{\tau}^{t})\cdot\mathbf{n}, \mathbf{u}\rangle_{-\frac{1}{2},\frac{1}{2},\Gamma_{N}},\end{split} \tag{111}\]

where \(P:V=H_{\boldsymbol{\sigma}}\times[H^{1}(\Omega)]^{d}\to V_{0}=H_{\boldsymbol{ \sigma},\Gamma_{N}}\times[H_{\Gamma_{D}}^{1}(\Omega)]^{d}\) is the projector into the subspace \(V_{0}\) of the Hilbert graph space \(V\) with scalar product:

\[((\boldsymbol{\sigma},\mathbf{u}),(\boldsymbol{\tau},\mathbf{v}))_{V}=( \mathbf{u},\mathbf{v})_{[L^{d}(\Omega)]^{d}}+(\boldsymbol{\sigma},\boldsymbol {\tau})_{[L^{2}(\Omega)]^{d\times d}}+(A(\boldsymbol{\sigma},\mathbf{u}),A( \boldsymbol{\tau},\mathbf{v})). \tag{112}\]

## Appendix C ROM convergence studies

In this section, we validate the DD-ROM implementation, checking the convergence towards the FOM solutions with respect to the dimension of the reduced space. Uniform local reduced dimensions are employed \(\{r_{\Omega_{i}}\}_{i=1}^{K}\) and \(\{r_{\Omega_{i}}\}_{i=1}^{k}\). For each convergence study 20 uniformly independent samples are used as training dataset and 50 uniformly independent samples as test dataset.

Figure 23: **MS1. The convergence of DD-ROMS with uniform local reduced dimensions \(\{r_{\Omega_{i}}\}_{i=1}^{K}\) and \(\{r_{\Omega_{i}}\}_{i=1}^{k}\) is assessed. The uniform value of the local reduced dimensions is reported in the abscissae. For this test case an improvement of the accuracy with respect to the single domain reduced basis is not observed.**
In Figure 23, we show the \(L^{2}\)-error, the \(R\)-error and the energy error decay and their respective error estimators for the Maxwell equations test case **MS1**, section 5.3.1, with constant parameters \(\mu\) and \(\sigma\) on the whole domain. We clearly see an exponential behavior in the error as we add basis functions. On the other hand, we do not observe strong differences between the ROM, DD-ROM with repartition and DD-ROM with deal.II subdomains, for this simple test case. Similar results can be observed in Figure 24, where the same analysis is applied for the compressible linear elasticity test **CLE1** from section 5.3.2.

From these results, it should be clear that the employment of local reduced basis is not always useful to increase the accuracy of the predictions. Nonetheless, it may be used to locally reduce the dimension of the linear approximants. Possible benefits include the adaptation of the computational resources (higher dimensional reduced basis are chosen only where it is necessary) and the possibility to speedup parametric studies and non-intrusive surrogate modelling thanks to the further reduced local dimensions [86, 87].

Typical cases where DD-ROMs are effective to increase the accuracy of the predictions are truly decomposable systems where the parameters affect independently different regions of the computational domain, as in test case **MS2** in section 5.3.1.

## Acknowledgements

This work was partially funded by European Union Funding for Research and Innovation -- Horizon 2020 Program -- in the framework of European Research Council Executive Agency: H2020 ERC CoG 2015 AROMA-CFD project 681447 "Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics" P.I. Professor Gianluigi Rozza. We also acknowledge the PRIN 2017 "Numerical Analysis for Full and Reduced Order Methods for the efficient and accurate solution of complex systems governed by Partial Differential Equations" (NA-FROM-PDEs). Davide Torlo has been funded by a SISSA Mathematical fellowship within Italian Excellence Departments initiative by Ministry of University and Research.
