# A non-overlapping optimization-based domain decomposition approach to component-based model reduction of incompressible flows

Tommaso Taddei, Xuejun Xu, Lei Zhang

► To cite this version:

Tommaso Taddei, Xuejun Xu, Lei Zhang. A non-overlapping optimization-based domain decomposition approach to component-based model reduction of incompressible flows. Journal of Computational Physics, 2024, 509, pp.113038. 10.1016/j.jcp.2024.113038. hal-04382123

HAL Id: hal-04382123

https://inria.hal.science/hal-04382123v1

Submitted on 6 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A non-overlapping optimization-based domain decomposition approach to component-based model reduction of incompressible flows

Tommaso Taddei¹, Xuejun Xu²,³, Lei Zhang².

¹ Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France
Inria Bordeaux Sud-Ouest, Team MEMPHIS, 33400 Talence, France, tommaso.taddei@inria.fr

² School of Mathematical Sciences, Tongji University, Shanghai 200092, China, 22210@tongji.edu.cn

³ Institute of Computational Mathematics, AMSS, Chinese Academy of Sciences, Beijing 100190, China,
xxj@lsec.cc.ac.cn

## Abstract

We present a component-based model order reduction procedure to efficiently and accurately solve parameterized incompressible flows governed by the Navier-Stokes equations. Our approach leverages a non-overlapping optimization-based domain decomposition technique to determine the control variable that minimizes jumps across the interfaces between sub-domains. To solve the resulting constrained optimization problem, we propose both Gauss-Newton and sequential quadratic programming methods, which effectively transform the constrained problem into an unconstrained formulation. Furthermore, we integrate model order reduction techniques into the optimization framework, to speed up computations. In particular, we incorporate localized training and adaptive enrichment to reduce the burden associated with the training of the local reduced-order models. Numerical results are presented to demonstrate the validity and effectiveness of the overall methodology.

**Keywords:** component-based model order reduction; optimization-based domain decomposition; non-overlapping methods; Navier-Stokes equations.

## 1 Introduction

Parameterized model order reduction (pMOR) techniques [1, 2, 3, 4] have gained widespread popularity in science and engineering to reduce the computational cost in scenarios that involve repetitive computational tasks, such as many-query and real-time applications. Given the parameter domain *P* and a parameterized partial differential equation (PDE) of interest, pMOR strategies rely on an offline/online computational decomposition: in the offline stage, which is computationally expensive and performed only once, a reduced basis (RB) approximation space is generated by exploiting several high-fidelity (HF) solutions (e.g., finite element, finite volume) to the parameterized PDE for properly chosen parameter values, and a reduced order model (ROM) is then devised; in the online stage, for any new parameter value, the ROM can be solved with computational cost independent of the HF discretization size *N*<sub>hf</sub>, to ensure significant computational savings. Efficient training algorithms, such as proper orthogonal decomposition (POD, [5, 6]) and the weak-Greedy algorithm [3] are available to construct the reduced order basis (ROB). Additionally, effective projection-based techniques [7, 8] can be employed to devise ROMs that are suitable for online calculations.

The combination of RB methods and domain decomposition (DD) methods offers further advantages [9, 10, 11]. First, localized pMOR techniques do not require global HF solutions over the whole domain: this feature has the potential to dramatically reduce the offline computational burden for large-scale systems. Second, localization simplifies the task of defining a parameterization of the problem and enables model reduction of systems with parameter-induced topology changes (cf. section 2.2). Third, the DD framework offers the flexibility to seamlessly integrate ROMs with full order models (FOMs, generated by the HF discretization) or to accommodate multi-physics applications based on independent software.

Various approaches have been proposed to combine RB methods and DD methods which differ in the way local ROMs are coupled at components' interfaces. In the reduced basis element (RBE) method [12, 13, 14], local ROMs are glued together using Lagrange multipliers. This method has been introduced in the context of the Laplace equation [12, 13] and subsequently applied to the Stokes equations [14]. A more recent application of the RBE method to the unsteady 3D Navier-Stokes equations can be found in [15], where a spectral Lagrange multiplier on the 2D interfaces is employed to couple local solutions. Another approach is the static condensation
RBE (scRBE) method [9, 16, 17], which ensures the component coupling through a static condensation procedure [18]. Additionally, approximation spaces for the interfaces (ports) between the components are also constructed [16, 17] to further reduce the computational complexity associated with the static condensation system. Another advantage of the scRBE method is the interchangeability of the components, which enables the study of different systems from a single library of parameterized archetype components. The RB-DD-finite-element (RDF) method [19] uses parametric boundary conditions in the local problems to define versatile local RB spaces for handling of networks composed of repetitive geometries characterized by different parameters. A detailed review of these methods can be found in [19].

Iterative techniques based on substructuring and the Schwarz alternating methods [20, 21] have been adapted to the pMOR framework [22, 10, 11, 23, 24]. In [22], both a non-overlapping Dirichlet–Neumann iterative scheme and a Schwarz method for overlapping sub-domains are proposed to ensure coupling between the FOM and the ROM. The coupling is achieved by ensuring the solution compatibility between the FOM solution trace and ROM solution trace at the interface. Specifically, only Galerkin-free ROMs are considered in the work of [22]. Galerkin-based ROMs are explored in the context of DD in [11], the authors develop a versatile coupling framework for both FOM-ROM coupling and ROM-ROM coupling, which can be applied to both overlapping and non-overlapping domains. Similarly, in [10] Galerkin-based ROMs are employed to speed up Dirichlet-Neumann DD iterations. A Dirichlet-Neumann DD-ROM is developed in [23] to handle non-conforming interfaces. Here, the Dirichlet and Neumann interface data are transferred using the INTERNODES method [25]. In [24], the authors present a DD-ROM technique which is designed for heterogeneous systems: in this approach, components are treated separately, and a parametrization of the interface data is used to generate HF snapshots.

Moreover, several authors have proposed to formulate the coupling problem as a minimization statement [26, 27]. In [26], the optimization problem is framed as the minimization of the difference between the ROM reconstruction and the corresponding FOM solution within the overlapping region between the ROM and the FOM domain. This approach adopts Galerkin-free ROMs and is applied to approximating incompressible flows, such as the interaction between an airfoil and a vortex, and the incompressible turbulent flow past a vehicle with varying geometry. The one-shot overlapping Schwarz method [27] consists in a constrained optimization statement that penalizes the jump at the interfaces of the components, while adhering to the approximate fulfillment of the PDE within each sub-domain. This approach has been validated for a steady nonlinear mechanics problem and also applied to an unsteady nonlinear mechanics problem with internal variables [28], in combination with overlapping partitions. The results of [27] showed that the minimization framework, which enables the application of effective optimization solvers for nonlinear least-square problems, ensures rapid convergence to the solution and is also robust with respect to the overlapping size.

In the present work, we aim to extend the method of [27] to incompressible flows in non-overlapping domains: our point of departure is the variational formulation proposed in [29] and further developed in [30, 31, 32]. As in [29], we formulate the DD problem as an optimal control problem where the control is given by the flux on the components' interfaces and the dependent variables are velocity and pressure in each subdomain; our formulation reads as a constrained minimization problem where the objective functional measures the jump in the dependent variables across the common boundaries between subdomains, while the constraints are the partial differential equations in each subdomain. We modify the formulation of [29] to incorporate an auxiliary control variable for the continuity equation which weakly ensures continuous finite-dimensional pressure across the interface; furthermore, we propose a specialized sequential quadratic programming (SQP) method to efficiently solve the optimization problem without resorting to Lagrange multipliers. We remark that non-overlapping techniques are of particular interest for heterogeneous DD [33] tasks that necessitate the combination of different discretization methods in each subdomain. Non-overlapping methods are also of interest for interface problems with high-contrast coefficients [34] and for fluid flows in repetitive networks [15, 19, 17] such as the vascular system.

We here consider two-dimensional steady-state simulations at moderate Reynolds number; however, our ultimate goal is to devise a flexible computational tool to simulate vascular flows in real, patient-specific geometries. We interpret complex networks as the union of a small number of parameterized components. In order to avoid expensive global solves at training stage, we propose a combined localized training and global enrichment strategy that exclusively uses local HF solves to create local approximations for the archetype components, thus avoiding the need for computationally demanding global HF solves during the training phase.

Our work is related to several previous contributions to component-based (CB) pMOR. First, the variational formulation is strongly related to the recent work by Prusak et al. [35]. The authors of [35] consider separate spaces for velocity and pressure and rely on pressure supremizer enrichment in combination with Galerkin projection to ensure stability of the local problems; furthermore, they resort to a Lagrangian multiplier and gradient-based methods as in [29] to solve the global optimization problem. Instead, we consider a single reduced space for velocity and pressure; we rely on both the Galerkin projection and a Petrov-Galerkin formulation for the local problems; and we rely on the Gauss-Newton and SQP methods for optimization without resorting to Lagrange multipliers. Finally, the authors of [35] do not discuss the problem of localized training, which is of paramount importance for the success of CB techniques. Second, we emphasize that several authors have previously developed CB-pMOR methods for incompressible flows in repetitive geometries [19, 17]; in particular,
the work by Pegolotti and coauthors [15] first considered a CB-pMOR for the unsteady incompressible Navier-Stokes equations in realistic three-dimensional geometries. Third, the localized training and global enrichment strategies are an extension of the method proposed in [36]: localized training strategies have been previously proposed in [16, 17, 24]; similarly, enrichment techniques have been considered in several efforts for linear elliptic PDEs (see, e.g., [37]).

This paper is organized as follows. In section 2, we introduce the optimization-based domain decomposition method and the model problem considered in this work. In section 3, we review the variational formulation introduced in [29]; we present our new formulation; and we discuss the solution method based on Gauss-Newton and sequential quadratic programming. Then in section 4 we discuss the integration of projection-based ROMs into the proposed optimization framework and the hybrid solver that combines both the FOM solver and the ROM solver. In sections 3 and 4 we illustrate the method for a simplified geometric configuration with two components. Section 5 is dedicated to the presentation of the localized training and the adaptive enrichment techniques. Finally, in section 6, we present numerical results that validate the effectiveness of our methodology.

## 2 Optimization-based domain decomposition method for the Navier-Stokes equations

In this work, we consider the incompressible Navier-Stokes equations:

$$ egin{cases} -\nu \Delta \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p = \mathbf{f} & \text{in } \Omega, \\ \nabla \cdot \mathbf{u} = 0 & \text{in } \Omega, \\ \mathbf{u}|_{\Gamma_{\text{dir}}} = \mathbf{u}_{\text{in}}, \quad \mathbf{u}|_{\Gamma_{\text{dir}}^{0}} = 0, \quad (\nu \nabla \mathbf{u} - p \mathbf{I}) \mathbf{n}|_{\Gamma_{\text{neu}}} = 0, \end{cases} \quad (1) $$

where $\nu > 0$ denotes the kinematic viscosity of the fluid, $\Omega$ is a bounded Lipschitz domain; the open sets $\Gamma_{\text{dir}}$, $\Gamma_{\text{dir}}^{0}$, $\Gamma_{\text{neu}}$ constitute a partition of $\partial\Omega$, which are associated to non-homogeneous Dirichlet boundary conditions, homogeneous Dirichlet boundary conditions and Neumann boundary conditions, respectively. We consider two-dimensional problems; the extension to the three-dimensional case and to unsteady problems is beyond the scope of this paper.

### 2.1 Optimization-based domain decomposition

<figure><img src="image_2.png" /><figcaption>Figure 1: The domain $\Omega$ and a partition into two non-overlapping sub-domains.</figcaption></figure>

For the purpose of clarity, we introduce the optimization-based domain decomposition method in the case of two sub-domains. Note that this approach can be readily extended to accommodate many sub-domains, as discussed in the subsequent sections. Consider a non-overlapping partition of $\Omega$ into two open sub-domains $\Omega_1$ and $\Omega_2$ such that $\bar{\Omega} = \bar{\Omega}_1 \cup \bar{\Omega}_2$, as illustrated in Figure 1. The interface that separates the two sub-domains is denoted by $\Gamma_0$ so that $\Gamma_0 = \bar{\Omega}_1 \cap \bar{\Omega}_2$. The vectors $\mathbf{n}_i$, $i = 1, 2$, are the unit outward normals of $\Omega_i$ on $\Gamma_0$ (we thus have $\mathbf{n}_1 = -\mathbf{n}_2$). We define the local Dirichlet and Neumann conditions for each component $\Omega_i$, $i = 1, 2$ as

$$ \Gamma_{i,\text{dir}} = \Gamma_{\text{dir}} \cap \partial\Omega_i, \quad \Gamma_{i,\text{dir}}^{0} = \Gamma_{\text{dir}}^{0} \cap \partial\Omega_i, \quad \Gamma_{i,\text{neu}} = \Gamma_{\text{neu}} \cap \partial\Omega_i, \quad (2a) $$

and the spaces

$$ \mathcal{X}_i := \{ (\mathbf{v}, q) \in [H^1(\Omega_i)]^2 \times L^2(\Omega_i) : \mathbf{v}|_{\Gamma_{i,\text{dir}}^0} = 0 \}, \quad \mathcal{X}_{i,0} := \{ (\mathbf{v}, q) \in \mathcal{X}_i : \mathbf{v}|_{\Gamma_{i,\text{dir}}} = 0 \}, \quad \mathcal{G} := [L^2(\Gamma_0)]^2. \quad (2b) $$

The local solution $(\mathbf{u}_i, p_i) \in \mathcal{X}_i$ is fully determined by the flux $\mathbf{g}$ at the interface $\Gamma_0$: as in [29], we thus refer to $\mathbf{g}$ as the control. Given the control $\mathbf{g} \in \mathcal{G}$, the velocity-pressure pair $(\mathbf{u}_i, p_i)$ satisfies $\mathbf{u}_i|_{\Gamma_{i,\text{dir}}} = \mathbf{u}_{\text{in}}|_{\Gamma_{i,\text{dir}}}$ and

$$ \mathcal{R}_i(\mathbf{u}_i, p_i, \mathbf{v}, q) + \mathcal{E}_i(\mathbf{g}, \mathbf{v}) = 0 \quad \forall (\mathbf{v}, q) \in \mathcal{X}_{i,0}, \quad (2c) $$

where

$$ \mathcal{R}_i(\mathbf{u}_i, p_i, \mathbf{v}, q) = \int_{\Omega_i} (\nu \nabla \mathbf{u}_i : \nabla \mathbf{v} - p_i(\nabla \cdot \mathbf{v}) - q(\nabla \cdot \mathbf{u}_i) - \mathbf{f}_i \cdot \mathbf{v}) dx, \quad \mathcal{E}_i(\mathbf{g}, \mathbf{v}) = (-1)^i \int_{\Gamma_0} \mathbf{g} \cdot \mathbf{v} dx, \quad (2d) $$
for i = 1, 2. Here, the orientation of the flux g is chosen to be the same as $n_1$, i.e., from $\Omega_1$ to $\Omega_2$; the choice of the orientation is completely arbitrary. Note that an arbitrary choice of the control g does not guarantee that the local solutions $(u_i, p_i)$ are solutions to (1); however, if $(u_1 - u_2)|_{\Gamma_0} = 0$, we find that the field $(u, p)$ such that $(u|_{\Omega_1}, p|_{\Omega_1})$ and $(u|_{\Omega_2}, p|_{\Omega_2})$ satisfy (2c) is a weak solution to the global problem (1). The optimal control g should hence guarantee velocity equality at the interface $\Gamma_0$.

Gutzburger and coauthors [29, 32] proposed the following optimization-based domain-decomposition formulation to compute the desired control and the local solutions:

$$ \min_{\substack{(u_1, p_1) \in \mathcal{X}_1; \\ (u_2, p_2) \in \mathcal{X}_2; \\ g \in \mathcal{G}}} \frac{1}{2} \int_{\Gamma_0} |\mathbf{u}_1 - \mathbf{u}_2|^2 dx + \frac{\delta}{2} \int_{\Gamma_0} |\mathbf{g}|^2 dx \quad \text{s.t.} \quad \left\{ \begin{array}{l} \mathcal{R}_i(\mathbf{u}_i, p_i, \mathbf{v}, q) + \mathcal{E}_i(\mathbf{g}, \mathbf{v}) = 0 \quad \forall (\mathbf{v}, q) \in \mathcal{X}_{i,0} \\ \mathbf{u}_i|_{\Gamma_{i,\text{dir}}} = \mathbf{u}_{\text{in}}|_{\Gamma_{i,\text{dir}}}, \end{array} \right. \quad i=1,2. \qquad (3) $$

The second term in the objective function of (3) is a regularizer that is designed to penalize controls of excessive size; the positive constant $\delta$ is chosen to control the relative importance of the two terms in the objective. The proofs of the well-posedness of the optimization formulation, as well as the convergence of the optimal solution to the solution to (1) as the regularization parameter $\delta$ approaches 0, can be found in [32].

## 2.2 Model problem

As in [15], we assume that the geometry of interest can be effectively approximated through instantiations of the elements of a library of prototype components; the instantiated components are obtained by low-rank geometric transformations of the prototype components. As in [15], we consider a library with two prototype components: “junction” and “channel”; the two prototype components are depicted in Figure 2, where a number is assigned to each component edge. These edge numbers indicate boundary face groups that are associated with the ports and the different types of boundary conditions. Specifically, for the junction, edge numbers {1, 4, 7} denote the ports and edge numbers {2, 3, 5, 6, 8, 9} indicate homogeneous Dirichlet boundaries; while for the channel, edge numbers {1, 2} represent the ports and edge numbers {3, 4} correspond to homogeneous Dirichlet boundaries.

A system can then be constructed by instantiating the two prototype components as follows:

$$ \bar{\Omega} = \bigcup_{i=1}^{N_{dd}} \bar{\Omega}_i, \quad \text{where} \quad \Omega_i = \Phi^{L_i}(\tilde{\Omega}^{L_i}, \mu_i), \quad i = 1, \dots, N_{dd}, $$

where $L_i \in \{1, 2\}$ denotes the label of the i-th component of the system, $\tilde{\Omega}^1$, $\tilde{\Omega}^2$ represent the two prototype components, $\Phi^{L_i}$ encompasses geometric transformations such as rotation, translation and non-rigid deformation that are applied to the prototype component to obtain the corresponding instantiated component that appears in the target system. The deformation of the i-th component is governed by the geometric parameter $\mu_i$; the vector $\mu_i$ includes a scaling factor $\gamma$, the angle $\theta$ and a shift $\mathbf{x}_{\text{shift}}$ that characterize the linear map that ensures the exact fitting of consecutive elements at ports. For the junction component, the vector $\mu_i$ also includes the angle $\alpha$, which represents the angle between the main vessel and the branch vessel, as shown in Figure 2(a); for the channel, the vector $\mu_i$ includes the constant $h_c$, which is used in the parameterization of the bottom boundary of the channel as $y = -h_c(4t(1-t))^{\alpha_c}$, with $t \in [0, 1]$ and $\alpha_c = 4$.

We prescribe a parabolic (Poiseuille) profile at the left boundary $u_{\text{in}}$ and we prescribe homogeneous Neumann conditions at the other boundary ports. In conclusion, the complete system configuration is uniquely prescribed by (i) the component labels $\{L_i\}_{i=1}^{N_{dd}}$ and the geometric parameters $\mu = \text{vec}(\mu_1, \dots, \mu_{N_{dd}})$, and (ii) the Reynolds number $Re$ at the inlet. We define the Reynolds number as $Re = \frac{H u_0}{\nu}$, where $H = 1$ denotes the diameter of the vessel at the inlet, $u_0$ represents the centerline velocity imposed at the inlet, and $\nu$ is the kinematic viscosity. In the numerical implementation, we set $\nu = \frac{1}{Re_{\text{ref}}}$ in all the components of the network, and we consider the parametric inflow condition $u_0(Re) = \frac{Re}{Re_{\text{ref}}}$.

Figure 3 illustrates two examples of target system, which consist of 3 and 4 components, respectively: the red numbers indicate the indices of the components, while the blue numbers indicate the internal ports. Note that the two systems are not isomorphic to each other: parameter variations hence induce topology changes that prevent the application of standard monolithic pMOR techniques.

**Remark 1.** We here observe that each component includes mixed Dirichlet-Neumann boundary conditions: the presence of Neumann conditions prevents the problem of pressure indeterminacy (up to an additive constant), and the existence of Dirichlet conditions eliminates the need for any additional compatibility condition [30] concerning the control variable g.
<figure><img src="image_3.png" /><figcaption>Figure 2: archetyp e components.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: two example of target systems.</figcaption></figure>

**Remark 2.** We observe that the boundary face group 1 for the two archetyp e components either corresponds to an internal interface or to the inlet Dirichlet condition (for the first component of the network). In order to handle this scenario, we can either modify the objective function to include the extra-term $\int_{\Gamma_{\text{dir}}} |\mathbf{u} - \mathbf{u}_{\text{in}}|^2 dx$ or to distinguish between inflow and internal channel and junction components. The latter option leads to a library with ($N_c = 4$) archetyp e components. We here opt for the second strategy.

3 High-fidelity discretization

3.1 Finite element spaces

We proceed to discretize the optimization statement (3). Towards this end, we introduce the HF spaces $V_i^{\text{hf}} \subset [H_{0,\Gamma_{i,\text{dir}}}^1 (\Omega_i)]^2$, $Q_i^{\text{hf}} \subset L^2(\Omega_i)$. We further define the tensor product spaces $\mathcal{X}_i^{\text{hf}} = V_i^{\text{hf}} \times Q_i^{\text{hf}}$ and the lifted space $\mathcal{X}_{i,0}^{\text{hf}} = V_{i,0}^{\text{hf}} \times Q_i^{\text{hf}}$ with $V_{i,0}^{\text{hf}} = \{\mathbf{v} \in V_i^{\text{hf}} : \mathbf{v}|_{\Gamma_{i,\text{dir}}} = 0\}$ for $i = 1, 2$. We denote by $\{\varphi_{i,j}\}_{j=1}^{N_i^u}$ a basis of $V_i^{\text{hf}}$ and by $\{\psi_{i,j}\}_{j=1}^{N_i^p}$ a basis of $Q_i^{\text{hf}}$; we use notation $\bullet$ to indicate the FE vector associated with the FE field $\bullet$. We further define the trace spaces $\Lambda_i^{\text{hf}} = \{\tau_{\Gamma_0} \mathbf{v} : \mathbf{v} \in V_i^{\text{hf}}\}$ and $\Xi_i^{\text{hf}} = \{\tau_{\Gamma_0} q : q \in Q_i^{\text{hf}}\}$, where $\tau_{\Gamma_0} \bullet := \bullet|_{\Gamma_0}$ indicates the trace of the field $\bullet$ on $\Gamma_0$. We here consider conforming meshes such that nodes at the interface shared by the two sub-domains coincide, that is $\Lambda_1^{\text{hf}} = \Lambda_2^{\text{hf}} = \Lambda^{\text{hf}}$ and $\Xi_1^{\text{hf}} = \Xi_2^{\text{hf}} = \Xi^{\text{hf}}$; this assumption is exploited in the technical result of B; nevertheless, the formulation can be trivially extended to non-conforming grids. We further define the global spaces $\mathcal{X}_i^{\text{hf}} = V_i^{\text{hf}} \times Q_i^{\text{hf}}$ and $\mathcal{X}_0^{\text{hf}} = V_0^{\text{hf}} \times Q_i^{\text{hf}}$ with $V_0^{\text{hf}} = \{\mathbf{v} \in V_i^{\text{hf}} : \mathbf{v}|_{\Gamma_{\text{dir}}} = 0\}$.

In this work, we adopt a stabilized FE formulation that incorporates the Streamline Upwind/Petrov-Galerkin (SUPG) [38, 39] and the Pressure-Stabilized Petrov-Galerkin (PSPG) [40] stabilizations. The PSPG technique allows the use of the same polynomial degree for both pressure and velocity discretizations; the SUPG technique enhances robustness for high Reynolds numbers. The detailed description of these stabilization formulas is given in A. In conclusion, we consider the following local problems, which are the counterpart of (2c):

$$
\left\{
\begin{array}{ll}
\mathcal{R}_i^{\text{hf}}(\mathbf{u}_i, p_i, \mathbf{v}, q) + \mathcal{E}_i(\mathbf{g}, \mathbf{v}) = 0 & \forall (\mathbf{v}, q) \in \mathcal{X}_{i,0}^{\text{hf}} \\
\mathbf{u}_i|_{\Gamma_{i,\text{dir}}} = \mathbf{\Phi}_{i,\mathbf{u}_{\text{in}}} & i = 1, 2.
\end{array}
\right.
\qquad (4a)
$$

where $\mathbf{\Phi}_{i,\mathbf{u}_{\text{in}}} \in V_i^{\text{hf}}$ is the interpolant of the nodal values of $\mathbf{u}_{\text{in}}$ on $\Gamma_{i,\text{dir}}$ [20, p. 174]. In view of the discussion
below, we rewrite the HF residual as

$$ \mathcal{R}_i^{\text{hf}}(\mathbf{u}_i, p_i, \mathbf{v}, q) = \mathcal{R}_{i,u}^{\text{hf}}(\mathbf{u}_i, p_i, \mathbf{v}) + \mathcal{R}_{i,p}^{\text{hf}}(\mathbf{u}_i, p_i, q); \quad (4b) $$

the first term corresponds to the residual of the momentum equation (1)1, while the second term corresponds to the residual of the continuity equation (1)2.

## 3.2 Variational formulation

Exploiting the previous notation, we can introduce the HF counterpart of the optimization formulation (3):

$$ \min_{ \substack{ (\mathbf{u}_1, p_1) \in \mathcal{X}_1^{\text{hf}}; \\ (\mathbf{u}_2, p_2) \in \mathcal{X}_2^{\text{hf}}; \\ \mathbf{g} \in \Lambda^{\text{hf}} } } \frac{1}{2} \int_{\Gamma_0} |\mathbf{u}_1 - \mathbf{u}_2|^2 dx + \frac{\delta}{2} \int_{\Gamma_0} |\mathbf{g}|^2 dx \quad \text{s.t.} \quad \left\{ \begin{array}{l} \mathcal{R}_i^{\text{hf}}(\mathbf{u}_i, p_i, \mathbf{v}, q) + \mathcal{E}_i(\mathbf{g}, \mathbf{v}) = 0 \quad \forall (\mathbf{v}, q) \in \mathcal{X}_{i,0}^{\text{hf}} \\ \mathbf{u}_i|_{\Gamma_{i,\text{dir}}} = \mathbf{\Phi}_{i,\text{u}_{\text{in}}}, \quad i = 1, 2. \end{array} \right. \quad (5) $$

This formulation coincides with the statement considered in [32] and also [35] — with the minor difference that we here rely on a stabilized FE formulation for the local problems. In the remainder of this section, we discuss an alternative HF formulation that will be used to define the reduced-order model.

Formulation (5) does not ensure the continuity of pressure across the internal interfaces: we prove this result rigorously in B; here, we provide a sketch of the proof that justifies our new DD statement. If we denote by $(\mathbf{u}^{\text{hf}}, p^{\text{hf}}) \in \mathcal{X}^{\text{hf}}$ the solution to the global problem such that $\mathcal{R}^{\text{hf}}(\mathbf{u}^{\text{hf}}, p^{\text{hf}}, \mathbf{v}, q) = 0$ and we neglect for simplicity the stabilization term, we obtain

$$ \mathcal{R}_p^{\text{hf}}(\mathbf{u}^{\text{hf}}, q) = \int_{\Omega_1} \nabla \cdot \mathbf{u}^{\text{hf}} q \, dx + \int_{\Omega_2} \nabla \cdot \mathbf{u}^{\text{hf}} q \, dx = 0 \quad \forall q \in Q^{\text{hf}}. $$

Since $Q^{\text{hf}}$ is a space of continuous functions, it is in general false that $\mathcal{R}_{i,p}^{\text{hf}}(\mathbf{u}^{\text{hf}}|_{\Omega_i}, q) = \int_{\Omega_i} \nabla \cdot \mathbf{u}^{\text{hf}} q \, dx = 0$ for all $q \in Q_i^{\text{hf}}$, $i = 1, 2$; nevertheless, it is possible to show that there exists $h^* \in \Xi^{\text{hf}}$ such that

$$ \mathcal{R}_{i,p}^{\text{hf}}(\mathbf{u}^{\text{hf}}|_{\Omega_i}, q) + (-1)^i \int_{\Gamma_0} h^* q \, dx = 0 \quad \forall q \in Q_i^{\text{hf}}. $$

Similarly, there exists $\mathbf{g}^* \in \Lambda^{\text{hf}}$ such that

$$ \mathcal{R}_{i,u}^{\text{hf}}(\mathbf{u}^{\text{hf}}|_{\Omega_i}, p^{\text{hf}}|_{\Omega_i}, \mathbf{v}) + (-1)^i \int_{\Gamma_0} \mathbf{g}^* \cdot \mathbf{v} \, dx = 0 \quad \forall \mathbf{v} \in V_i^{\text{hf}}, \quad i = 1, 2. $$

We conclude that the tuple $(\mathbf{u}^{\text{hf}}|_{\Omega_1}, p^{\text{hf}}|_{\Omega_1}, \mathbf{u}^{\text{hf}}|_{\Omega_2}, p^{\text{hf}}|_{\Omega_2}, \mathbf{g}^*, h^*)$ is a solution to the minimization problem

$$ \begin{array}{ll} \min_{ \substack{ (\mathbf{u}_1, p_1) \in \mathcal{X}_1^{\text{hf}}; \\ (\mathbf{u}_2, p_2) \in \mathcal{X}_2^{\text{hf}}; \\ \mathbf{g} \in \Lambda^{\text{hf}}, h \in \Xi^{\text{hf}} } } & \frac{1}{2} \int_{\Gamma_0} |\mathbf{u}_1 - \mathbf{u}_2|^2 dx + \frac{1}{2} \int_{\Gamma_0} (p_1 - p_2)^2 dx \\ \text{s.t.} & \left\{ \begin{array}{l} \mathcal{R}_i^{\text{hf}}(\mathbf{u}_i, p_i, \mathbf{v}, q) + \mathcal{E}_i(\mathbf{g}, \mathbf{v}) + (-1)^i \int_{\Gamma_0} hq \, dx = 0 \quad \forall (\mathbf{v}, q) \in \mathcal{X}_{i,0}^{\text{hf}}, \quad i = 1, 2. \\ \mathbf{u}_i|_{\Gamma_{i,\text{dir}}} = \mathbf{\Phi}_{i,\text{u}_{\text{in}}}, \end{array} \right. \end{array} $$

This discussion suggests to consider a modified formulation that explicitly penalizes the jump of the pressure field. We introduce the state $\mathbf{w}_i := \text{vec}(\mathbf{u}_i, p_i)$, $i = 1, 2$ and the control $\mathbf{s} := \text{vec}(\mathbf{g}, h)$; we introduce the control space $\mathcal{S}^{\text{hf}} = \Lambda^{\text{hf}} \times \Xi^{\text{hf}}$ equipped with the norm

$$ \|\mathbf{s} = \text{vec}(\mathbf{g}, h)\|^2 = \int_{\Gamma_0} |\nabla_{\Gamma_0} \mathbf{g}|^2 + |\mathbf{g}|^2 + h^2 \, dx, \quad (6a) $$

where $\nabla_{\Gamma_0} \mathbf{g}$ denotes the gradient of $\mathbf{g}$ in the tangential direction; we use notation $\mathbf{w}_i(1:2)$ to indicate the first two components of the vector-valued function $\mathbf{w}_i$. Then, we introduce the variational formulation:

$$ \min_{ \substack{ \mathbf{w}_1 \in \mathcal{X}_1^{\text{hf}}; \\ \mathbf{w}_2 \in \mathcal{X}_2^{\text{hf}}; \\ \mathbf{s} \in \mathcal{S}^{\text{hf}} } } \mathcal{F}_{\delta}(\mathbf{w}_1, \mathbf{w}_2, \mathbf{s}) \quad \text{s.t.} \quad \left\{ \begin{array}{l} \mathcal{R}_i^{\text{hf}}(\mathbf{w}_i, \mathbf{z}) + \mathcal{E}_i^{\text{hf}}(\mathbf{s}, \mathbf{z}) = 0 \quad \forall \mathbf{z} \in \mathcal{X}_{i,0}^{\text{hf}}, \\ \mathbf{w}_i(1:2)|_{\Gamma_{i,\text{dir}}} = \mathbf{\Phi}_{i,\text{u}_{\text{in}}}, \end{array} \right. \quad i = 1, 2; \quad (6b) $$

where

$$ \mathcal{F}_{\delta}(\mathbf{w}_1, \mathbf{w}_2, \mathbf{s}) := \frac{1}{2} \|\mathbf{w}_1 - \mathbf{w}_2\|_{L^2(\Gamma_0)}^2 + \frac{1}{2} \delta \|\mathbf{s}\|^2, \quad (6c) $$
and $\mathcal{E}_i^{\text{hf}}(\mathbf{s}, \mathbf{v}, q) = \mathcal{E}_i(\mathbf{g}, \mathbf{v}) + (-1)^i \int_{\Gamma_0} hq \, dx$. Note that we replaced the $L^2$ norm for the control $\mathbf{g}$ with the $H^1$ norm: as discussed in section 6.1, we empirically observe that the use of the $H^1$ norm significantly reduces the oscillations in the profile of $\mathbf{g}$.

Some comments are in order. First, the addition of the pressure jump and of the control $h$ ensures that the optimal pressure is continuous in the limit $\delta \to 0$. Note that at the continuous level the test space $Q = L^2(\Omega)$ is discontinuous; therefore, the control $h$ is unnecessary. Similarly, if we rely on a P0 discretization for the pressure field [41], the pressure jump is also unnecessary. Second, since velocity and pressure have different units and might also have very different magnitudes, it might be necessary to rescale the objective function to avoid stability issues (see, e.g., [42]). In our numerical experiments, we solve the equations in non-dimensional form, and we do not include any scaling factor.

## 3.3 Solution methods for (6)

As in [27] and also [32], we resort to a gradient-based optimization method to find local minima of (6). In more detail, we consider the Gauss-Newton method (GNM) and sequential quadratic programming (SQP) [43]. As discussed below, both methods rely on static condensation to devise a reduced system for the control $\mathbf{s}$.

### 3.3.1 Gauss-Newton method

We define the local solution map $\mathcal{H}_i : \mathcal{S}^{\text{hf}} \to \mathcal{X}_i^{\text{hf}}$ such that $\mathcal{H}_i(\mathbf{s})(1:2)|_{\Gamma_{i,\text{dir}}} = \mathbf{\Phi}_{i,\text{u}_{\text{in}}}$ and

$$ \mathcal{R}_i^{\text{hf}}(\mathcal{H}_i(\mathbf{s}), \mathbf{z}) + \mathcal{E}_i^{\text{hf}}(\mathbf{s}, \mathbf{z}) = 0 \quad \forall \mathbf{z} \in \mathcal{X}_{i,0}^{\text{hf}}, \quad i = 1, 2. \quad (7) $$

Then, we rewrite (6) as an unconstrained optimization problem:

$$ \min_{\mathbf{s} \in \mathcal{S}^{\text{hf}}} \mathcal{F}_{\delta}^{\text{gn}}(\mathbf{s}) = \mathcal{F}_{\delta}(\mathcal{H}_1(\mathbf{s}), \mathcal{H}_2(\mathbf{s}), \mathbf{s}). \quad (8) $$

If we define the space $\mathcal{X}^{\text{hf}} = \Lambda^{\text{hf}} \times \mathcal{S}^{\text{hf}}$ equipped with the norm $\|\mathbf{r} = \text{vec}(\mathbf{w}, \mathbf{g}, h)\|_{\mathcal{X}^{\text{hf}}}^2 = \|\mathbf{w}\|_{L^2(\Gamma_0)}^2 + \|\mathbf{g}\|_{H^1(\Gamma_0)}^2 + \|h\|_{L^2(\Gamma_0)}^2$ and the operator $F_{\delta}: \mathcal{S}^{\text{hf}} \to \mathcal{X}^{\text{hf}}$ such that $F_{\delta}(\mathbf{s}) = \text{vec}(\tau_{\Gamma_0}(\mathcal{H}_1(\mathbf{s}) - \mathcal{H}_2(\mathbf{s})), \sqrt{\delta}\mathbf{s})$, we can rewrite (8) as a nonlinear least-square problem, that is

$$ \min_{\mathbf{s} \in \mathcal{S}^{\text{hf}}} \mathcal{F}_{\delta}^{\text{gn}}(\mathbf{s}) = \frac{1}{2} \|F_{\delta}(\mathbf{s})\|_{\mathcal{X}^{\text{hf}}}^2. \quad (9a) $$

The unconstrained problem (9a) can be solved efficiently using GNM: given the initial condition $\mathbf{s}^{it=0}$, we repeatedly solve, for $it = 0, 1, \dots$,

$$ \mathbf{s}^{it+1} = \arg \min_{\mathbf{s} \in \mathcal{S}^{\text{hf}}} \frac{1}{2} \| F_{\delta}(\mathbf{s}^{it}) + \frac{\partial F_{\delta}(\mathbf{s}^{it})}{\partial \mathbf{s}} (\mathbf{s} - \mathbf{s}^{it}) \|_{\mathcal{X}^{\text{hf}}}^2. \quad (9b) $$

with the termination condition

$$ \frac{\|\mathbf{s}^{it+1} - \mathbf{s}^{it}\|}{\|\mathbf{s}^{it}\|} \le tol, \quad (9c) $$

where $tol > 0$ is a predefined tolerance.

We observe that GNM requires the explicit calculation of $F_{\delta}$ and the gradient of $F_{\delta}$ with respect to the control at $\mathbf{s}^{it}$: the former involves the solution to the local problems (7) for all components, while the latter is given by

$$ \frac{\partial F_{\delta}(\mathbf{s}^{it})}{\partial \mathbf{s}} = \left[ \tau_{\Gamma_0} \left( \frac{\partial \mathcal{H}_1(\mathbf{s}^{it})}{\partial \mathbf{s}} - \frac{\partial \mathcal{H}_2(\mathbf{s}^{it})}{\partial \mathbf{s}} \right) \right] \quad \text{with} \quad \frac{\partial \mathcal{H}_i(\mathbf{s})}{\partial \mathbf{s}} = - \left( \frac{\partial \mathcal{R}_i^{\text{hf}}(\mathcal{H}_i(\mathbf{s}))}{\partial \mathbf{w}_i} \right)^{-1} \mathcal{E}_i^{\text{hf}}, \quad (10) $$

and $id$ is the identity map. We notice that the evaluation of $\frac{\partial F_{\delta}(\mathbf{s}^{it})}{\partial \mathbf{s}}$ involves the solution to $N^s$ linear systems where $N^s$ is the cardinality of the space $\mathcal{S}^{\text{hf}}$; it is hence computationally feasible only if the dimension of the control is moderate: this observation highlights the importance of port reduction [16] for optimization-based methods. Conversely, we remark that the computation of $\mathcal{H}_1(\mathbf{s}^{it})$, $\mathcal{H}_2(\mathbf{s}^{it})$ and their derivatives is embarrassingly parallel with respect to the number of components: as discussed in [27], GNM enables effective parallelization of the solution procedure if compared to standard multiplicative Schwartz iterative methods, provided that the computational cost is dominated by the solution to the local problems (7). Finally, we remark that the least-square problem in (9b) can be solved by explicitly assembling the normal equations; alternatively, we might employ the QR factorization [7]. We omit the details.
### 3.3.2 Sequential quadratic programming (SQP)

The SQP method solves a sequence of optimization subproblems, each of which optimizes a quadratic model of the objective subject to a linearization of the constraints. Since the objective (cf. (6c)) is quadratic, we hence find the iterative method

$$ (\mathbf{w}_1^{it+1}, \mathbf{w}_2^{it+1}, \mathbf{s}^{it+1}) = \arg \min_{\substack{\mathbf{w}_1 \in \mathcal{X}_1^{hf}; \\ \mathbf{w}_2 \in \mathcal{X}_2^{hf}; \\ \mathbf{s} \in \mathcal{S}^{hf}}} \mathcal{F}_{\delta}(\mathbf{w}_1, \mathbf{w}_2, \mathbf{s}) \quad \text{s.t.} \quad \left\{ \begin{array}{l} \mathcal{R}_{i,it}^{hf}(\mathbf{z}) + \mathcal{J}_{i,it}^{hf}(\mathbf{w}_i - \mathbf{w}_i^{it}, \mathbf{z}) + \mathcal{E}_i^{hf}(\mathbf{s}, \mathbf{z}) = 0 \\ \mathbf{w}_i(1:2)|_{\Gamma_{i,\text{dir}}} = \mathbf{\Phi}_{i,\mathbf{u}_{in}}, \quad \forall \mathbf{z} \in \mathcal{X}_{i,0}^{hf}, \quad i=1,2; \end{array} \right. \tag{11a} $$

where the linear forms {$\mathcal{R}_{i,it}^{hf}$}$_i$ and the bilinear forms {$\mathcal{J}_{i,it}^{hf}$}$_i$ are given by

$$ \mathcal{R}_{i,it}^{hf}(\mathbf{z}) = \mathcal{R}_i^{hf}(\mathbf{w}_i^{it}, \mathbf{z}), \quad \mathcal{J}_{i,it}^{hf}(\mathbf{w}, \mathbf{z}) = \frac{\partial \mathcal{R}_i^{hf}}{\partial \mathbf{w}_i} [\mathbf{w}_i^{it}](\mathbf{w}, \mathbf{z}), \quad \forall \mathbf{w} \in \mathcal{X}_i^{hf}, \mathbf{z} \in \mathcal{X}_{i,0}^{hf}, \quad it = 0, 1, \dots \tag{11b} $$

In the numerical experiments, we consider the same termination condition (9c) used for GNM.

The optimization problem (11a) is quadratic with linear constraints. The solution to (11a) hence satisfies

$$ \left\{ \begin{array}{l} \mathbf{s}^{it+1} = \arg \min_{\mathbf{s} \in \mathcal{S}^{hf}} \| \tilde{F}_{\delta}^{it} + \tilde{\mathcal{J}}_{\delta}^{it} (\mathbf{s} - \mathbf{s}^{it}) \|_{\mathcal{X}^{hf}}^2; \\ \mathbf{w}_i^{it+1} = \mathbf{w}_i^{it} - (\mathcal{J}_{i,it}^{hf})^{-1} (\mathcal{R}_{i,it}^{hf} + \mathcal{E}_i^{hf} \mathbf{s}^{it+1}), \quad i=1,2; \end{array} \right. \tag{12a} $$

where

$$ \tilde{F}_{\delta}^{it} = \begin{bmatrix} \tau_{\Gamma_0}(\mathbf{w}_1^{it} - \mathbf{w}_2^{it}) \\ \sqrt{\delta} \mathbf{s}^{it} \end{bmatrix}, \quad \tilde{\mathcal{J}}_{\delta}^{it} = \begin{bmatrix} \tau_{\Gamma_0}((\mathcal{J}_{1,it}^{hf})^{-1} \mathcal{E}_1^{hf} - (\mathcal{J}_{2,it}^{hf})^{-1} \mathcal{E}_2^{hf}) \\ \sqrt{\delta} \mathbf{id} \end{bmatrix}. \tag{12b} $$

In our implementation, we rely on (12) to solve (11a).

As for GNM, we obtain a least-square problem for the control by applying static condensation: while in the previous section we first derived the unconstrained statement (cf. (8)) and then we applied the optimization method, here we first optimize using SQP and then we apply static condensation at each iteration of the optimization algorithm.

Since the underlying PDE model is nonlinear, GNM requires to perform Newton subiterations to solve the local problems (7) (see also the definition of $F_{\delta}(\mathbf{s}^{it})$ in (9b)); conversely, SQP does not involve subiterations. The cost per iteration of SQP is hence significantly inferior to the cost of GNM. We empirically observe that the SQP approach mitigates the potential convergence issues of the sub-iterations for the local problems, particularly at the very early stages of the optimization loop.

We observe that (9b) and (12a)$_1$ are formally equivalent, while (10) and (12b) share the same structure. We conclude that the SQP and GNM approaches can be implemented using the same data structures and can be parallelized in the same way. We omit the details.

**Remark 3.** For high-Reynolds number flows, it is important to enhance the robustness of our approach by resorting to pseudo transient continuation (PTC) [44]. PTC introduces an additional pseudo-temporal integration with adaptive time step, that is performed until convergence to a steady-state solution. If we resort to the backward Euler scheme for the discretization of the time derivative, at each PTC step we solve the relaxed problem:

$$ \min_{\substack{\mathbf{w}_1 \in \mathcal{X}_1^{hf}; \\ \mathbf{w}_2 \in \mathcal{X}_2^{hf}; \\ \mathbf{s} \in \mathcal{S}^{hf}}} \mathcal{F}_{\delta}(\mathbf{w}_1, \mathbf{w}_2, \mathbf{s}) \quad \text{s.t.} \quad \left\{ \begin{array}{l} \displaystyle \frac{1}{\Delta t_k} \int_{\Omega_i} (\mathbf{w}_i(1:2) - \mathbf{w}_i^k(1:2)) \cdot \mathbf{v} \, dx + \mathcal{R}_i^{hf}(\mathbf{w}_i, \mathbf{z}) + \mathcal{E}_i^{hf}(\mathbf{s}, \mathbf{z}) = 0 \\ \displaystyle \mathbf{w}_i(1:2)|_{\Gamma_{i,\text{dir}}} = \mathbf{\Phi}_{i,\mathbf{u}_{in}}, \quad \forall \mathbf{z} = (\mathbf{v}, q) \in \mathcal{X}_{i,0}, \quad i=1,2. \end{array} \right. \tag{13} $$

where the index k refers to the temporal loop and $\Delta t_k$ is chosen adaptively based on the residual of the steady-state equations. We refer to [44] and to the references therein for further details. Note that (13) is formally equivalent to (6c): it can hence be solved using the same procedure outlined above. As discussed in A, the time derivative should also be included in the SUPG and PSPG stabilization terms.

## 4 Projection-based reduced order formulation

We rely on the formulation (6b) to define the CB-ROM. Towards this end, first, we identify a low-rank approximation of the control $\mathbf{s}^{hf}$ and the local states $\mathbf{w}_1^{hf}, \mathbf{w}_2^{hf}$; second, we devise local ROMs for the approximation of the solution maps (7); third, we devise specialized GNM and SQP methods for the formulation (6b) based on
approximate solution maps. We conclude the section by discussing the implementation of hybrid formulations that combine full-order and reduced-order local solution maps. We remark that in order to further enhance online performance we should also reduce the online costs associated with the computation of the $||\cdot||_{\chi^{hf}}$ norm in (9b) and (12a)$_1$ (cf. [27]): we do not address this issue in the present work.

## 4.1 Construction of the local basis

We denote by $\{\mu^{(k)} = \text{vec}(\mu_1^{(k)}, \mu_2^{(k)})\}_{k=1}^{n_{\text{train}}}$ a set of $n_{\text{train}}$ global configurations; we further denote by $\{\mathbf{w}_{i,k}^{\text{hf}} : i = 1, 2, k = 1, \dots, n_{\text{train}}\}$ and $\{\mathbf{s}_k^{\text{hf}} : k = 1, \dots, n_{\text{train}}\}$ the corresponding HF state and control estimates based on (6b). We resort to POD to devise a low-dimensional approximation space for the local solution manifolds and for the control

$$ \begin{cases} [\mathcal{Z}_i = \text{span}\{\zeta_{i,j}\}_{j=1}^n] = \text{POD} (\{\mathbf{w}_{i,k}^{\text{hf}} - \mathbf{\Psi}_{i,\mathbf{u}_{\text{in}}}^{(k)}\}_{k=1}^{n_{\text{train}}} \|\cdot\|_{\chi_i, n}); \\ [\mathcal{W} = \text{span}\{\eta_j\}_{j=1}^m] = \text{POD} (\{\mathbf{s}_k^{\text{hf}}\}_{k=1}^{n_{\text{train}}} \|\cdot\|_{\mathcal{H}}, m). \end{cases} \quad (14) $$

Here, the function POD($\mathcal{D}$, $||\cdot||$, $n$) returns the POD space of dimension $n$ associated with the snapshot dataset $\mathcal{D}$ and the norm $||\cdot||$ using the method of snapshots [45]. To ease the presentation, the integers $n$ and $m$ are here chosen a priori: in practice, we should choose $n, m$ based on the energy criterion. The fields $\mathbf{\Psi}_{1,\mathbf{u}_{\text{in}}}^{(k)}, \mathbf{\Psi}_{2,\mathbf{u}_{\text{in}}}^{(k)}$ satisfy the boundary conditions in (6b); we refer to section 5 for the explicit expression; this implies that the local space $\mathcal{Z}_i$ is contained in $\mathcal{X}_{i,0}$, for $i = 1, 2$. In the remainder, we further use notation $\mathcal{Z}_i^{\text{dir}} = \{\mathbf{\Psi}_{i,\mathbf{u}_{\text{in}}} (\mu) + \zeta_i : \zeta_i \in \mathcal{Z}_i\}$ to identify the affine approximation spaces that incorporate Dirichlet boundary conditions. Furthermore, given $\mathbf{w}_i \in \mathcal{Z}_i^{\text{dir}}$ and $\mathbf{s} \in \mathcal{W}$, we define the generalized coordinates $\alpha_1, \alpha_2 \in \mathbb{R}^n$ and $\beta \in \mathbb{R}^m$ such that

$$ \mathbf{w}_i(\alpha_i; \mu) = \mathbf{\Psi}_{i,\mathbf{u}_{\text{in}}}(\mu) + \sum_{j=1}^{n} (\alpha_j)_j \zeta_{i,j}, \quad i = 1, 2; \quad \mathbf{s}(\beta) = \sum_{j=1}^{m} (\beta)_j \eta_j. \quad (15) $$

## 4.2 Construction of the local reduced-order models

We rely on (Petrov-)Galerkin projection to define the local ROMs.

**Galerkin projection.** We consider the local solution maps $\hat{\mathbf{H}}_i^g : \mathcal{W} \to \mathcal{Z}_i^{\text{dir}}$ such that

$$ \mathcal{R}_i^{\text{hf}}(\hat{\mathbf{H}}_i^g(\mathbf{s}), \mathbf{z}) + \mathcal{E}_i^{\text{hf}}(\mathbf{s}, \mathbf{z}) = 0 \quad \forall \mathbf{z} \in \mathcal{Z}_i, \quad i = 1, 2. \quad (16) $$

It is useful to rewrite (16) in fully-algebraic form. Towards this end, we define the discrete residuals $\hat{\mathbf{R}}_i^g : \mathbb{R}^n \to \mathbb{R}^n$ and $\hat{\mathbf{E}}_i^g \in \mathbb{R}^{n \times n}$ such that

$$ (\hat{\mathbf{R}}_i^g(\alpha))_j = \mathcal{R}_i^{\text{hf}}(\mathbf{w}_i(\alpha_i), \zeta_{i,j}), \quad (\hat{\mathbf{E}}_i^g)_{j,k} = \mathcal{E}_i^{\text{hf}}(\eta_k, \zeta_{i,j}), \quad i = 1, 2, j = 1, \dots, n, k = 1, \dots, m; \quad (17a) $$

and the local algebraic solution maps $\hat{\mathbf{H}}_i^g : \mathbb{R}^m \to \mathbb{R}^n$ such that

$$ \hat{\mathbf{R}}_i^g(\hat{\mathbf{H}}_i^g(\beta)) + \hat{\mathbf{E}}_i^g\beta = 0, \quad i = 1, 2. \quad (17b) $$

**Least-square Petrov-Galerkin (LSPG, [7]) projection.** Given the reduced space $\mathcal{Y}_i \subset \mathcal{X}_{i,0}$, we introduce the local solution maps $\hat{\mathbf{H}}_i^{\text{pg}} : \mathcal{W} \to \mathcal{Z}_i^{\text{dir}}$ such that

$$ \hat{\mathcal{H}}_i^{\text{pg}}(\mathbf{s}) = \arg \min_{\zeta \in \mathcal{Z}_i^{\text{dir}}} \sup_{\mathbf{z} \in \mathcal{Y}_i} \frac{\mathcal{R}_i^{\text{hf}}(\zeta, \mathbf{z}) + \mathcal{E}_i^{\text{hf}}(\mathbf{s}, \mathbf{z})}{\|\mathbf{z}\|_{\chi_i}}. \quad (18) $$

For $\mathcal{Y}_i = \mathcal{X}_{i,0}$, (18) is referred to as minimum residual projection. In view of the derivation of the algebraic counterpart of (18), we denote by $\{\mathbf{v}_{i,k}\}_{k=1}^{j_{\text{es}}}$ an orthonormal basis of $\mathcal{Y}_i$; then, we define the algebraic residuals

$$ (\hat{\mathbf{R}}_i^{\text{pg}}(\alpha_i))_\ell = \mathcal{R}_i^{\text{hf}}(\mathbf{w}_i(\alpha_i), \mathbf{v}_{i,\ell}), \quad (\hat{\mathbf{E}}_i^{\text{pg}})_{\ell,k} = \mathcal{E}_i^{\text{hf}}(\eta_k, \mathbf{v}_{i,\ell}), \quad (19a) $$

with $i = 1, 2, \ell = 1, \dots, j_{\text{es}}, k = 1, \dots, m$; and the local algebraic solution maps $\hat{\mathbf{H}}_i^{\text{pg}} : \mathbb{R}^m \to \mathbb{R}^n$ such that

$$ \hat{\mathbf{H}}_i^{\text{pg}}(\beta) = \arg \min_{\alpha \in \mathbb{R}^n} |\hat{\mathbf{R}}_i^{\text{pg}}(\alpha) + \hat{\mathbf{E}}_i^g\beta|, \quad i = 1, 2. \quad (19b) $$

We observe that (19b) reads as a nonlinear least-square problem that can be solved efficiently using GNM; the combination of LSPG ROMs within the DD formulation (6b) is challenging: we address this issue in the next section.
The ROM (18) depends on the choice of the test space $\mathcal{Y}_i$. Following [46, 47], we propose to construct the test space $\mathcal{Y}_i$ using POD. Given the snapshots $\{\mathbf{w}_{i,k}^{\text{hf}}\}_k$ and the ROB $\{\zeta_{i,j}\}_{j=1}^n$, we compute the Riesz elements $\psi_{i,j,k} \in \mathcal{X}_{i,0}^{\text{hf}}$ such that

$$ (\psi_{i,j,k}, \mathbf{z})_{\mathcal{X}_i} = \frac{\partial \mathcal{R}_i^{\text{hf}}}{\partial \mathbf{w}_i} [\mathbf{w}_{i,k}^{\text{hf}}] (\zeta_{i,j}, \mathbf{z}), \quad \forall \mathbf{z} \in \mathcal{X}_{i,0}^{\text{hf}}, \qquad (20a) $$

for $i=1, 2, j=1, \dots, n, k=1, \dots, n_{\text{train}}$. Then, we apply POD to find the low-dimensional bases $\mathcal{Y}_1$ and $\mathcal{Y}_2$,

$$ [\mathcal{Y}_i = \text{span}\{\mathbf{v}_{i,j}\}_{j=1}^{j_{es}}] = \text{POD}(\{\psi_{i,j,k} : j=1, \dots, n, k=1, \dots, n_{\text{train}}, \|\cdot\|_{\mathcal{X}_i, j_{es}}\}, i=1, 2. \quad (20b) $$

As in [46, 47], we choose $j_{es} = 2n$; we refer to [46, Appendix C] for a rigorous justification of the choice of the test space for linear inf-sup stable problems.

**Remark 4.** The solution to (16) and (18) is expensive due to the need to evaluate the HF residual and its Jacobian at each iteration. To reduce the computational burden, several authors have proposed to resort to hyper-reduction strategies [48] to speed up assembly costs at prediction stage. We refer to the recent review [49] for a detailed presentation of the subject. Since the local problems (16) and (18) fit in the framework of monolithic pMOR, standard hyper-reduction techniques can be employed. We refer to a future work for the development and the assessment of hyper-reduction techniques for the DD formulation of this work.

## 4.3 Global formulation

We first introduce the algebraic counterpart of the objective (6c). We denote by $\{(\mathbf{x}_q^{\Gamma}, \omega_q^{\Gamma})\}_{q=1}^{N_{\Gamma}}$ the FE quadrature rule of $\int_{\Gamma_0}[\cdot] dx$ and we define the matrices $\mathbf{A}_1, \mathbf{A}_2 \in \mathbb{R}^{3N_{\Gamma} \times n}$ and the vector $\mathbf{b} \in \mathbb{R}^{3N_{\Gamma}}$ such that

$$ (\mathbf{A}_i)_{q+(\ell-1)N_{\Gamma}, j} = \sqrt{\omega_q^{\Gamma}} (\zeta_{i,j}(\mathbf{x}_q^{\Gamma}))_{\ell}, \quad (\mathbf{b})_{q+(\ell-1)N_{\Gamma}} = \sqrt{\omega_q^{\Gamma}} (\Psi_{1,\mathbf{u}_{\text{in}}}(\mathbf{x}_q^{\Gamma}) - \Psi_{2,\mathbf{u}_{\text{in}}}(\mathbf{x}_q^{\Gamma}))_{\ell} \quad (21a) $$

with $q=1, \dots, N_{\Gamma}, \ell=1, 2, 3, j=1, \dots, n$; then, we rewrite the objective function as

$$ \mathcal{F}_{\delta}(\alpha_1, \alpha_2, \beta) = \mathcal{F}_{\delta}(\mathbf{w}_1(\alpha_1), \mathbf{w}_2(\alpha_2), \mathbf{s}(\beta)) = \frac{1}{2}|\mathbf{A}_1\alpha_1 - \mathbf{A}_2\alpha_2|^2 + \frac{\delta}{2}|\beta|^2. \quad (21b) $$

For the Galerkin local ROMs, the DD ROM can be obtained by simply projecting (6c) onto the reduced spaces, that is

$$ \begin{array}{ll} \min_{\mathbf{w}_1 \in \mathcal{Z}_i^{\text{dir}}; \mathbf{w}_2 \in \mathcal{Z}_i^{\text{dir}}; \mathbf{s} \in \mathcal{W}} & \mathcal{F}_{\delta}(\mathbf{w}_1, \mathbf{w}_2, \mathbf{s}) \quad \text{s.t.} \quad \mathcal{R}_i^{\text{hf}}(\mathbf{w}_i, \mathbf{z}) + \mathcal{E}_i^{\text{hf}}(\mathbf{s}, \mathbf{z}) = 0 \quad \forall \mathbf{z} \in \mathcal{Z}_i, \quad i=1,2. \quad (22a) \\ & \mathbf{w}_1 \in \mathcal{Z}_i^{\text{dir}}; \mathbf{w}_2 \in \mathcal{Z}_i^{\text{dir}}; \mathbf{s} \in \mathcal{W} \end{array} $$

Note that non-homogeneous Dirichlet conditions are encoded in the choice of the ansatz. Exploiting the previous notation, we obtain the algebraic counterpart of (22a).

$$ \min_{\substack{\alpha_1, \alpha_2 \in \mathbb{R}^n; \\ \beta \in \mathbb{R}^m}} \mathcal{F}_{\delta}(\alpha_1, \alpha_2, \beta) \quad \text{s.t.} \quad \hat{\mathbf{R}}_i^g(\alpha_i) + \hat{\mathbf{E}}_i^g \beta = 0, \quad i=1,2. \quad (22b) $$

Problem (22b) can be solved using either GNM or SQP; as for the HF model, the methods require the computation of the derivatives of the local solution maps (17b), which satisfy

$$ \frac{\partial \hat{\mathcal{H}}_i^g}{\partial \beta}(\beta) = - \left( \frac{\partial \hat{\mathbf{R}}_i^g}{\partial \alpha_i} \left[ \hat{\mathcal{H}}_i^g(\beta) \right] \right)^{-1} \hat{\mathbf{E}}_i^g. \quad (22c) $$

Note that (22c) can be computed using standard FE routines that are readily available for the full-order model.

The combination of (6b) with the LSPG ROM (19b) is more involved since the resulting component-based ROM cannot be interpreted as the projection of (6b) onto suitable low-dimensional spaces. We here rely on an approximate SQP procedure. At each iteration, given the triplet $(\alpha_1^{it}, \alpha_2^{it}, \beta^{it})$, we compute

$$ \hat{\mathbf{R}}_i^{\text{pg},it} = \hat{\mathbf{R}}_i^{\text{pg}}(\alpha_i^{it}) \in \mathbb{R}^{j_{es}}, \quad \hat{\mathbf{J}}_i^{\text{pg},it} = \frac{\partial \hat{\mathbf{R}}_i^{\text{pg}}}{\partial \alpha_i}(\alpha_i^{it}) \in \mathbb{R}^{j_{es} \times n}; \quad (23a) $$

then, we solve the minimization problem

$$ \min_{\substack{\alpha_1, \alpha_2 \in \mathbb{R}^n; \\ \beta \in \mathbb{R}^m}} \mathcal{F}_{\delta}(\alpha_1, \alpha_2, \beta) \quad \text{s.t.} \quad (\hat{\mathbf{J}}_i^{\text{pg},it})^{\top} (\hat{\mathbf{R}}_i^{\text{pg},it} + \hat{\mathbf{J}}_i^{\text{pg},it} (\alpha_i - \alpha_i^{it}) + \hat{\mathbf{E}}_i^{\text{pg}} \beta) = 0, \quad i=1,2. \quad (23b) $$
We observe that for $n = j_{es}$ the constraints imply that $\hat{\mathbf{R}}_{i}^{\text{pg},it} + \hat{\mathbf{J}}_{i}^{\text{pg},it} (\alpha_{i} - \alpha_{i}^{it}) + \hat{\mathbf{E}}_{i}^{\text{pg}} \beta = 0$ for $i = 1, 2$. We hence recover the standard SQP procedure.

A thorough convergence analysis of the SQP procedure (23) is beyond the scope of the present work. Here, we observe that if $\alpha_{i}^{it} \to \alpha_{i}^{*}$ for $i = 1, 2$ and $\beta^{it} \to \beta^{*}$, the constraints in (23b) reduce to

$$ \left( \frac{\partial \hat{\mathbf{R}}_{i}^{\text{pg}}}{\partial \alpha_{i}} (\alpha_{i}^{*}) \right)^{\top} \left( \hat{\mathbf{R}}_{i}^{\text{pg}} (\alpha_{i}^{*}) + \hat{\mathbf{E}}_{i}^{\text{pg}} \beta^{*} \right) = 0, \quad i = 1, 2. $$

Given $i \in \{1, 2\}$, the latter implies that $\alpha_{i}^{*}$ is a stationary point of the function $\alpha_{i} \mapsto |\hat{\mathbf{R}}_{i}^{\text{pg}}(\alpha_{i}) + \hat{\mathbf{E}}_{i}^{\text{pg}}\beta^{*}|^{2}$; provided that (19b) admits a unique solution, we hence find that $\alpha_{i}^{*} = \hat{\mathbf{H}}_{i}^{\text{pg}}(\beta^{*})$.

## 4.4 Enrichment of the trial space

In (14), we construct the state and control spaces independently. We might hence obtain that the matrices $\frac{\partial \hat{\mathbf{H}}_{i}^{\text{pg}}}{\partial \beta}(\beta)$ and $\frac{\partial \hat{\mathbf{H}}_{i}^{\text{pg}}}{\partial \alpha_{i}}(\beta)$ are rank-deficient: as empirically shown in the numerical examples, rank deficiency of the sensitivity matrices leads to instabilities of the ROM and to poor approximations of the control $\mathbf{s}$. To address this issue, we propose to enrich the trial spaces $\mathcal{Z}_1, \mathcal{Z}_2$ with the perturbed snapshots $\{\tilde{\mathbf{w}}_{i,j,k}\}_{i,j,k}$

$$ \mathcal{R}_{i}^{\text{hf}}(\mathbf{w}_{i,k}^{\text{hf}}, \mathbf{z}) + \frac{\partial \mathcal{R}_{i}^{\text{hf}}}{\partial \mathbf{w}_{i}} [\mathbf{w}_{i,k}^{\text{hf}}] (\tilde{\mathbf{w}}_{i,j,k} + \mathbf{\Psi}_{i,\mathbf{u}_{\text{in}}} - \mathbf{w}_{i,k}^{\text{hf}}, \mathbf{z}) + \mathcal{E}_{i}^{\text{hf}}(\boldsymbol{\eta}_{j}, \mathbf{z}) = 0 \quad \forall \mathbf{z} \in \mathcal{X}_{i,0}^{\text{hf}}. \quad (24) $$

In more detail, given the snapshots $\{\mathbf{w}_{i,k}^{\text{hf}} : i = 1, 2, k = 1, \dots, n_{\text{train}}\}$ and the reduced spaces $\mathcal{Z}_1, \mathcal{Z}_2, \mathcal{W}$, we compute the perturbations $\{\tilde{\mathbf{w}}_{i,j,k}\}_{j,k} \subset \mathcal{X}_{i,0}^{\text{hf}}$ for $i = 1, 2$, and then we update the reduced spaces $\mathcal{Z}_1$ and $\mathcal{Z}_2$ as follows:

$$ \mathcal{Z}_{i}^{\text{new}} = \mathcal{Z}_{i} \oplus \mathcal{Z}_{i}', \quad \text{with} \quad [\mathcal{Z}_{i}'] = \text{POD} \left( \{\Pi_{\mathcal{Z}_{i}^{\perp}} \tilde{\mathbf{w}}_{i,j,k} : j = 1, \dots, m, k = 1, \dots, n_{\text{train}}\}, \|\cdot\|_{\mathcal{X}_{i}}(n'), \right) \quad (25) $$

where $\Pi_{\mathcal{Z}_{i}^{\perp}} \bullet$ denotes the projection of $\bullet$ onto the orthogonal complement of the space $\mathcal{Z}_{i}$ and $n'$ is a given integer.

Some comments are in order. The hierarchical construction of the state approximation space (25) has been proposed in a similar context in [50]. The integer $n'$ should be sufficiently large to ensure stability of the DD formulation; we further comment on the selection of $n'$ in the numerical experiments. Finally, in C, we provide a formal justification of the enrichment strategy for a linear problem.

## 4.5 Hybrid solver

In the introduction, we anticipated the importance of developing a DD formulation that enables the seamless coupling of local, independently generated models. We here illustrate how to combine the HF model introduced in section 3 with the local ROM introduced in section 4. To provide a concrete reference, we assume that the HF model (7) is solved in $\Omega_1$ and that the LSPG ROM (18) is solved in $\Omega_2$.

We set $N_1^w = N_1^u + N_1^p$ and we define the basis (cf. section 3.1)

$$ \{\boldsymbol{\xi}_{1,j}\}_{j=1}^{N_1^w} = \{\text{vec}(\varphi_{1,1}, 0), \dots, \text{vec}(\varphi_{1,N^u}, 0), \text{vec}(0, 0, \psi_{1,1}), \dots, \text{vec}(0, 0, \psi_{1,N_1^p})\}. $$

We introduce the vector-valued representation of the lifted state field $\hat{\mathbf{w}}_1 = \mathbf{w}_1 - \mathbf{\Psi}_{1,\mathbf{u}_{\text{in}}} \in \mathbb{R}^{N_1^w}$. Then, we introduce the matrices (see (19a) and (21a)) $\mathbf{A}_1^{\text{hf}} \in \mathbb{R}^{3N_1^\Gamma \times N_1^w}$ and $\mathbf{E}_1^{\text{hf}} \in \mathbb{R}^{N_1^w \times m}$ such that

$$ (\mathbf{A}_1^{\text{hf}})_{q+(\ell-1)N_1^\Gamma, j} = \sqrt{\omega_q^\Gamma} (\boldsymbol{\xi}_{1,j}(\mathbf{x}_q^\Gamma))_\ell, \quad (\hat{\mathbf{E}}_1^{\text{hf}})_{j,k} = \mathcal{E}_i^{\text{hf}}(\boldsymbol{\eta}_k, \boldsymbol{\xi}_{1,j}). $$

Then, we can state the SQP method for the hybrid coupled problem:

$$ \min_{\substack{\hat{\mathbf{w}}_1 \in \mathbb{R}^{N_1^w} \\ \boldsymbol{\alpha}_2 \in \mathbb{R}^n \\ \boldsymbol{\beta} \in \mathbb{R}^m}} \frac{1}{2} |\mathbf{A}_1^{\text{hf}} \hat{\mathbf{w}}_1 - \mathbf{A}_2 \boldsymbol{\alpha}_2 + \mathbf{b}|^2 + \frac{\delta}{2} |\boldsymbol{\beta}|^2 \quad \text{s.t.} \quad \left\{ \begin{array}{l} \mathbf{R}_1^{\text{hf},it} + \mathbf{J}_1^{\text{hf},it} (\hat{\mathbf{w}}_1 - \hat{\mathbf{w}}_1^{it}) + \hat{\mathbf{E}}_1^{\text{hf}} \boldsymbol{\beta} = 0; \\ (\hat{\mathbf{J}}_2^{\text{pg},it})^\top (\hat{\mathbf{R}}_2^{\text{pg},it} + \hat{\mathbf{J}}_2^{\text{pg},it} (\boldsymbol{\alpha}_2 - \boldsymbol{\alpha}_2^{it}) + \mathbf{E}_2^{\text{pg}} \boldsymbol{\beta}) = 0. \end{array} \right. \quad (26a) $$

where

$$ (\mathbf{R}_1^{\text{hf},it})_j = \mathcal{R}_1^{\text{hf}}(\mathbf{w}_1^{it}, \boldsymbol{\xi}_{1,j}), \quad (\mathbf{J}_1^{\text{hf},it})_{j,k} = \frac{\partial \mathcal{R}_1^{\text{hf}}}{\partial \mathbf{w}_i} [\mathbf{w}_1^{it}] (\boldsymbol{\xi}_{1,k}, \boldsymbol{\xi}_{1,j}), \quad j, k = 1, \dots, N_1^w, \quad (26b) $$

with $\mathbf{w}_1^{it} = \hat{\mathbf{w}}_1^{it} + \mathbf{\Psi}_{1,\mathbf{u}_{\text{in}}}$.
Problem (26a) can be solved using the static condensation procedure outlined in (12a) and (12b). Note that for (26a) the least-square problem (12a)$_1$ is of size $m$: the computational cost is hence independent of $N_1^w$. On the other hand, the cost to assemble the least-square problem in (12a)$_1$ is dominated by the cost of computing $(\mathbf{J}_1^{\text{hf},it})^{-1}\mathbf{E}_1^{\text{hf}}$, which requires the solution to $m$ linear systems of size $N_1^w$. We emphasize that the local models in (26a) only communicate through the vehicle of the control $\mathbf{s}$ (or equivalently through the generalized coordinates $\boldsymbol{\beta}$) and the matrices $\mathbf{A}_1^{\text{hf}}, \mathbf{A}_2$ in the objective function: the implementation of the local models is hence agnostic to the discretization that is employed in the neighboring subdomain.

# 5 Localized training and adaptive enrichment

In section 4 we devised the CB-ROM based on the DD formulation (6b). The major limitation of the approach is the need for global HF solves to generate the reduced spaces (cf. (14)). In this section, we propose a general strategy to adaptively construct the reduced space for state and control, for the model problem of section 2.2. First, in section 5.1, we present the general multi-component DD formulation and relevant quantities that are employed in the adaptive procedure. Then, in sections 5.2 and 5.3, we present the localized training strategies for the control $\mathbf{s}$ and for the local states. Finally, in section 5.4 we present the adaptive enrichment strategy that allows the correction of the local approximations based on global reduced-order solves.

## 5.1 Multi-component formulation

Given the archetypal components $\{\tilde{\Omega}^k\}_{k=1}^{N_c}$ and the reference port $\tilde{\Gamma}$, we introduce the instantiated system $\Omega \subset \mathbb{R}^2$ such that $\tilde{\Omega} = \bigcup_{i=1}^{N_{\text{dd}}} \tilde{\Omega}_i$ with $\Omega_i = \Phi^{L_i}(\tilde{\Omega}^{L_i}, \mu_i)$ for $i = 1, \dots, N_{\text{dd}}$ and ports $\{\Gamma_j\}_{j=1}^{N_f}$ such that $\Gamma_j = \Psi_j(\tilde{\Gamma})$, where $\mu_1, \dots, \mu_{N_{\text{dd}}}$ are geometric parameters associated with the elemental mapping and $\Psi_1, \dots, \Psi_{N_{\text{dd}}}$ are the mappings associated with the ports; we further introduce the union of all ports $\Gamma := \bigcup_{j=1}^{N_f} \Gamma_j$. For $i = 1, \dots, N_{\text{dd}}$, we denote by $\mathcal{I}_i^\Gamma \subset \{1, \dots, N_f\}$ the set of the indices of the ports that belong to $\partial\Omega_i$. We further denote by $\mathbf{n}_j^+$ the positive normal to the port $\Gamma_j$. We denote by $\tilde{\Psi}_{k, \mathbf{u}_{\text{in}}}$ the HF solution to the Navier-Stokes equations in $\tilde{\Omega}_k$ with inflow condition $u_0(\text{Re}_{\text{ref}})$ for some $\text{Re}_{\text{ref}} > 0$ and Neumann boundary conditions on the remaining ports; then, we introduce the parametric field $\tilde{\Psi}_{k, \mathbf{u}_{\text{in}}}(\text{Re}) = \frac{\text{Re}}{\text{Re}_{\text{ref}}} \tilde{\Psi}_{k, \mathbf{u}_{\text{in}}}$.

We introduce the FE spaces $\tilde{\mathcal{X}}_k^{\text{hf}}$ and $\tilde{\mathcal{X}}_{k,0}^{\text{hf}}$ associated with the domain $\tilde{\Omega}^k$ (cf. section 3.1) for $k = 1, \dots, N_c$; furthermore, we introduce the reduced spaces $\tilde{Z}_k \subset \tilde{\mathcal{X}}_{k,0}^{\text{hf}}$ and the affine spaces $\tilde{Z}_k^{\text{dir}}(\text{Re}) := \tilde{\Psi}_{k, \mathbf{u}_{\text{in}}}(\text{Re}) + \tilde{Z}_k$ —to shorten notation, we omit the dependence of $\tilde{Z}_k^{\text{dir}}$ on the Reynolds number. The choice $\tilde{Z}_k = \tilde{\mathcal{X}}_{k,0}^{\text{hf}}$ corresponds to considering the HF discretization in all components of type $k$. Then, we define the global discontinuous approximation space over $\Omega$

$$ \mathcal{X}^{\text{dd}} := \{ \mathbf{w} \in [L^2(\Omega)]^3 : \mathbf{w}|_{\Omega_i} \circ \Phi^{L_i}(\cdot, \mu_i) \in \tilde{Z}_{L_i}^{\text{dir}}, i = 1, \dots, N_{\text{dd}} \}. \quad (27) $$

We denote by $[\mathbf{w}] \in [L^2(\Gamma)]^3$ the jump of the field $\mathbf{w}$ on the interfaces of the partition

$$ [\mathbf{w}](\mathbf{x}) = \mathbf{w}^{+}(\mathbf{x}) - \mathbf{w}^{-}(\mathbf{x}) \quad \forall \mathbf{x} \in \Gamma_j, \quad \mathbf{w}^{\pm}(\mathbf{x}) := \lim_{\epsilon \to 0^+} \mathbf{w}(\mathbf{x} \mp \epsilon \mathbf{n}_j^{+}(\mathbf{x})), \quad j = 1, \dots, N_f. \quad (28) $$

Given the port reduced space $\tilde{\mathcal{W}} \subset [L^2(\tilde{\Gamma})]^3$, we also introduce the global port space over $\Gamma$

$$ \mathcal{W}^{\text{dd}} := \{ \mathbf{s} \in [L^2(\Gamma)]^3 : \mathbf{s}|_{\Gamma_j} \circ \Psi_j \in \tilde{\mathcal{W}}, j = 1, \dots, N_f \}. \quad (29) $$

We handle geometry deformations using the discretize-then-map approach (cf. [47]). Given the FE field $\mathbf{w} \in \mathcal{X}_i^{\text{hf}}$, we denote by $\tilde{\mathbf{w}} \in \tilde{\mathcal{X}}_{L_i}^{\text{hf}}$ the corresponding field in the reference configuration; the two fields share the same FE vector. We introduce norms in the reference components

$$ \| \mathbf{w} = \text{vec}(\mathbf{u}, p) \|_{\tilde{\mathcal{X}}_k}^2 = \int_{\tilde{\Omega}^k} \nabla \mathbf{u} : \nabla \mathbf{u} + |\mathbf{u}|^2 + p^2 \, dx, \quad \| \mathbf{s} = \text{vec}(\mathbf{g}, h) \|_{\tilde{\Gamma}}^2 = \int_{\tilde{\Gamma}} |\nabla_{\tilde{\Gamma}} \mathbf{g}|^2 + |\mathbf{g}|^2 + h^2 \, dx, \quad (30) $$

for $k = 1, \dots, N_c$. Then, we define the corresponding norms for the instantiated components that are obtained by applying the prescribed deformation

$$ \| \mathbf{w} \|_{\mathcal{X}_i} := \| \mathbf{w} |_{\Omega_i} \circ \Phi^{L_i}(\cdot, \mu_i) \|_{\tilde{\mathcal{X}}_{L_i}} , \quad i = 1, \dots, N_{\text{dd}}, \quad \| \mathbf{s} \|^2 = \sum_{j=1}^{N_f} \| \mathbf{s} |_{\Gamma_j} \circ \Psi_j \|_{\tilde{\Gamma}}^2. \quad (31) $$

Note that the algebraic norms associated with (31) are independent of the geometric parameters that enter in the mappings $\{\Phi^{L_i}(\cdot, \mu_i)\}_i$: there exist indeed $N_c$ matrices $\mathbf{X}_1, \dots, \mathbf{X}_{N_c}$ such that $\|\mathbf{w}\|_{\mathcal{X}_i} = \sqrt{\mathbf{w}^\top \mathbf{X}_{L_i} \mathbf{w}}$ for
i = 1, ..., N<sub>dd</sub>. This observation simplifies the implementation of the dual residual norm used in the adaptive strategy (cf. (35)). Similarly, the variational forms associated with the PDE problem are defined for each archetype component and then mapped to obtain the variational forms for each instantiated component. We define the forms $\tilde{\mathcal{R}}_k^{\text{hf}} : \tilde{\mathcal{X}}_k^{\text{hf}} \times \tilde{\mathcal{X}}_{k,0}^{\text{hf}} \times \mathcal{P}_k \times \mathbb{R}_+ \to \mathbb{R}$ such that

$$ \tilde{\mathcal{R}}_{L_i}^{\text{hf}}(\tilde{\mathbf{w}}, \tilde{\mathbf{z}}; \mu_i, \text{Re}) = \mathcal{R}_i^{\text{hf}}(\mathbf{w}, \mathbf{z}), \quad \forall \mathbf{w} \in \mathcal{X}_i^{\text{hf}}, \quad \mathbf{z} \in \mathcal{X}_{i,0}^{\text{hf}}. \qquad (32) $$

We further define the boundary form

$$ \tilde{\mathcal{E}}_{L_i, \ell}^{\text{hf}}(\tilde{\eta}, \tilde{z}, \mu_i) = \int_{\Gamma_{j_i, \ell}} |\|\mathbf{w}\||^2 dx + \frac{\delta}{2} \|\mathbf{s}\|^2 \quad \text{s.t.} \quad \mathcal{R}_i^{\text{hf}}(\mathbf{w}, \mathbf{z}) + \mathcal{E}_i^{\text{hf}}(\mathbf{s}, \mathbf{z}) = 0 \quad \forall \mathbf{z} \in \mathcal{Z}_i, \quad i = 1, \dots, N_{\text{dd}}; \qquad (33) $$

where $j_{i,\ell} \in \{1, \dots, N_f\}$ is the index (in the global numbering) of the $\ell$-th port of the $i$-th component of the system.

We have now the elements to present the DD Galerkin formulation:

$$ \min_{\mathbf{w} \in \mathcal{X}^{\text{dd}}, \mathbf{s} \in \mathcal{W}^{\text{dd}}} \frac{1}{2} \int_{\Gamma} |\|\mathbf{w}\||^2 dx + \frac{\delta}{2} \|\mathbf{s}\|^2 \quad \text{s.t.} \quad \mathcal{R}_i^{\text{hf}}(\mathbf{w}, \mathbf{z}) + \mathcal{E}_i^{\text{hf}}(\mathbf{s}, \mathbf{z}) = 0 \quad \forall \mathbf{z} \in \mathcal{Z}_i, \quad i = 1, \dots, N_{\text{dd}}; \qquad (34a) $$

where $\mathcal{Z}_i = \{\zeta \in [H^1(\Omega_i)]^3 : \zeta \circ \Phi^{L_i}(\cdot, \mu_i) \in \tilde{\mathcal{Z}}_{L_i}\}$ and

$$ \mathcal{E}_i^{\text{hf}}(\mathbf{s}, \mathbf{z}) = \sum_{j \in \mathcal{J}_i} \int_{\Gamma_j} \mathbf{s} \cdot \mathbf{z} \, dx, \qquad (34b) $$

for $i = 1, \dots, N_{\text{dd}}$. Formulation (34a) can be adapted to cope with Petrov-Galerkin ROMs using the strategy outlined in section 4.5: we omit the details.

Given the estimate $(\mathbf{w}^*, \mathbf{s}^*)$ of the solution to (34a), we devise two error indicators to assess its accuracy; the indicators are employed in section 5.4 to drive the enrichment strategy. First, we define the local errors

$$ e_i := \sup_{\mathbf{z} \in \mathcal{X}_{i,0}^{\text{hf}}} \frac{\mathcal{R}_i^{\text{hf}}(\mathbf{w}^*, \mathbf{z}) + \mathcal{E}_i^{\text{hf}}(\mathbf{s}^*, \mathbf{z})}{\|\mathbf{z}\|_{\mathcal{X}_i}}, \quad i = 1, \dots, N_{\text{dd}}. \qquad (35) $$

The quantity $e_i$ measures the performance of the $i$-th ROM to approximate the solution to the Navier-Stokes equations for the control $\mathbf{s}^*$. We further introduce the jump errors:

$$ e_j^{\text{jump}} := \sqrt{\int_{\Gamma_j} |\|\mathbf{w}\||^2 dx}, \quad j = 1, \dots, N_f. \qquad (36) $$

The indicator (36) controls the jump of the state estimate at the interfaces: the value of $e_j^{\text{jump}}$ can thus be interpreted as the measure of the ability of the control to nullify the jump at the $j$-th interface of the domain.

**Remark 5.** In order to enhance the compressibility of the local state and control manifolds, following [51], in the numerical experiments, we consider the approximation spaces

$$ \begin{align} \mathcal{X}^{\text{dd}} &:= \left\{ \mathbf{w} \in L^2(\Omega; \mathbb{R}^3) : \mathbf{A}(\theta_i)\mathbf{w}|_{\Omega_i} \circ \Phi^{L_i}(\cdot, \mu_i) \in \tilde{\mathcal{Z}}_{L_i}^{\text{dir}}, \ i = 1, \dots, N_{\text{dd}} \right\}; \tag{37a} \\ \mathcal{W}^{\text{dd}} &:= \left\{ \mathbf{s} \in L^2(\Gamma; \mathbb{R}^3) : \mathbf{A}(\omega_j)\mathbf{s}|_{\Gamma_j} \circ \Psi_j \in \tilde{\mathcal{W}}, \ j = 1, \dots, N_f \right\}; \nonumber \end{align} $$

where $\theta_i$ (resp., $\omega_i$) is the angle between the inlet port of the $i$-th deformed component $\Omega_i$ (resp., the $j$-th port $\Gamma_j$) and the $x_1$ axis, and

$$ \mathbf{A}(\theta) = \begin{bmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \end{bmatrix}. \qquad (37b) $$

We remark that several authors have considered more sophisticated (Piola) transformations to improve the compressibility of solution manifolds in internal flows, (e.g. [14]): in this respect, our choice is a compromise between accuracy and simplicity of implementation.

## 5.2 Pairwise training for the control variables

Following [16, 17], we pursue a pairwise-training approach to generate the port space $\tilde{\mathcal{W}}$. We perform HF simulations for systems of two components that represent all possible connections (channel-channel, channel-junction, junction-junction, junction-channel) based on random Dirichlet boundary conditions at the inflow,
random Neumann conditions at the outflow, and a random selection of the Reynolds number and the geometric parameters in prescribed parameter ranges (cf. Figure 4). The HF data for the ports are retrieved and stored, and finally the port space $\tilde{W}$ is constructed using POD. Recalling (37), the HF data are rotated using (37b) before applying the compression technique.

Similarly to [16, 52], we consider the inlet velocity

$$ \mathbf{u}_{\text{in}}(y) = -\frac{\text{Re}}{\text{Re}_{\text{ref}}} \left( u_0(y) + \delta_u \sum_{k=1}^{R} \frac{c_k}{k^2} P_k(-1+2y) \right) \mathbf{n}, \quad (38) $$

where $\text{Re} \sim \text{Uniform}(\text{Re}_{\text{min}}, \text{Re}_{\text{max}})$, $\{P_k\}_k$ are zero-flowrate weighted polynomials (cf. [17, section 3.1.1])

$$ P_k(y) = \begin{cases} (1 - y^2)y, & \text{if } k = 1, \\ (1 - y^2)(5y^2 - 1), & \text{if } k = 2, \\ (1 - y^2)\mathcal{L}_k(y), & \text{if } 3 \le k \le R, \end{cases} $$

and $\{\mathcal{L}_k\}_k$ are the Legendre polynomials. The coefficients of the expansion are sampled from a standard Gaussian distribution, $c_1, \dots, c_R \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$, $\mathbf{n}$ denotes the outward normal to $\Omega$ on the inlet boundary, $y \in (0, 1)$ is the curvilinear coordinate, $u_0(y) = 4y(1 - y)$ is the Poiseuille velocity profile, the coefficient $\delta_u$ is selected a posteriori to ensure that the inflow is positive for all $y \in (0, 1)$. Similarly, we prescribe the outward flux as

$$ \mathbf{g}_{\text{out}}(y) = \left( g_0 + \delta_g \sum_{k=1}^{R} c_k^{\text{out}} \mathcal{L}_k(-1+2y) \right) \mathbf{n}, \quad c_1^{\text{out}}, \dots, c_R^{\text{out}} \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1), \quad (39) $$

where $g_0 \sim \text{Uniform}(g_{0\text{min}}, g_{0\text{max}})$, and we choose the coefficient $\delta_g$ to prevent reverse flow.

<figure><img src="image_5.png" /><figcaption>Figure 4: channel-junction connection for the training of s.</figcaption></figure>

## 5.3 Localized training for the state variables

After having built the reduced space for the control, we repeatedly solve (34a) for several random configurations and several parameter values to acquire datasets of simulations for each phenotype component. Thanks to port reduction, the computational cost of the global problem is significantly reduced if compared with the full HF model; nevertheless, we choose to consider systems with a moderate number of components (up to four) to further reduce offline costs. The HF data for components of the same type are mapped in the reference configurations, rotated through (37b), and are then used to build the local reduced spaces $\tilde{Z}_1, \dots, \tilde{Z}_{N_c}$.

We observe that the training strategy is not fully local since it requires to assemble systems with up to four components. In our experience, the practical implementation of a fully localized training strategy for incompressible flows is extremely challenging due to the need to ensure that the fluid flows from left to right and that the prescribed Neumann conditions lead to physical velocities. The choice of considering global training based on a reduced control space for systems of moderate dimension represents a trade-off between offline efficiency and accuracy. The adaptive strategy presented in the next section provides a systematic way to improve the quality of the local reduced spaces.

## 5.4 Adaptive enrichment

In Algorithm 1, we present the full adaptive strategy for the construction of the reduced spaces. The procedure extends the method introduced in [36]; to clarify the presentation, we postpone two steps of the algorithm to sections 5.4.1 and 5.4.2.
<caption>**Algorithm 1** Adaptive enrichment procedure.</caption>

1: Generate the reduced space $\tilde{W}$ for the control through pairwise training (cf. section 5.2).
2: Generate the local spaces $\{\tilde{Z}_k\}_{k=1}^{N_c}$ for the state through global training (cf. section 5.3).
3: Enrich the reduced spaces $\{\tilde{Z}_k\}_{k=1}^{N_c}$ based on the port space $\tilde{W}$ (cf. section 5.4.2).
4: Sample $n_{\text{train}}^{\text{glo}}$ global configurations, $\mathcal{P}_{\text{train}} := \{\mu^j\}_{j=1}^{N_c}$.
5: (if LSPG projection is employed) Build the empirical test space (cf. section 4.2)
6: **for** $i = 1, \dots, \text{maxit}$ **do**
7: &nbsp;&nbsp;&nbsp;&nbsp;Initialize the datasets $\mathcal{D}_{(1)} = \dots = \mathcal{D}_{(N_c)} = \emptyset$ and $\mathcal{D}_s = \emptyset$.
8: &nbsp;&nbsp;&nbsp;&nbsp;**for** $\mu \in \mathcal{P}_{\text{train}}$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute the reduced solution using the CB-ROM solver (cf. (34a)).
10: &nbsp;&nbsp;&nbsp;&nbsp;Compute local residuals $\{e_i\}_{i=1}^{N_{\text{dd}}}$ (cf. (35)) and the jumps $\{e_j^{\text{port}}\}_{j=1}^{N_f}$ (cf. (36))
11: &nbsp;&nbsp;&nbsp;&nbsp;Mark the $m_w$ instantiated components with the largest residuals of each type $\{\Gamma_{\text{mark}}^{\mu,(k)}\}_{k=1}^{N_c}$.
12: &nbsp;&nbsp;&nbsp;&nbsp;Mark the $m_s$ instantiated ports with the largest port jumps of each type $\Gamma_{\text{mark}}^{\mu,\text{p}}$.
13: &nbsp;&nbsp;&nbsp;&nbsp;Update the datasets $\mathcal{D}_{(1)}, \dots, \mathcal{D}_{(N_c)}$ and $\mathcal{D}_s$ (cf. section 5.4.1)
14: **end for**
15: Update the port POD space $\tilde{W} = \tilde{W} \oplus \text{POD}(\{\Pi_{\tilde{W}^{\perp}} \tilde{s} : \tilde{s} \in \mathcal{D}_s\}, \|\cdot\|_{\bar{\Gamma}}, n^{\text{glo}})$.
16: Update the reduced spaces $\tilde{Z}_k = \tilde{Z}_k \oplus \text{POD}(\{\Pi_{\tilde{Z}_k^{\perp}} \tilde{w} : \tilde{w} \in \mathcal{D}_{(k)}\}, \|\cdot\|_{\tilde{X}_k^{\text{glo}}}, n^{\text{glo}})$, $k = 1, \dots, N_c$.
17: (Optional) Enrich the reduced spaces $\{\tilde{Z}_k\}_{k=1}^{N_c}$ based on the port space (cf. section 5.4.2).
18: (if LSPG projection is employed) Update the empirical test space (cf. section 4.2)
19: **end for**

As in [36], we add $m_w$ (resp., $m_s$) snapshots to the state (resp., control) datasets for each element of $\mu \in \mathcal{P}_{\text{train}}$, instead of selecting the marked elements after having computed the local indicators for all configurations: this choice avoids the storage of all reduced global solutions and ultimately simplifies the implementation. In our experience, the enrichment of the state spaces is only needed for localized training (Line 3 of the Algorithm) but not after each update of the control space $\tilde{W}$ (Line 17 of the Algorithm): a possible explanation is that the enrichment step inherently couples the construction of the two spaces. Further numerical investigations are necessary to investigate this aspect.

Algorithm 1 depends on several user-defined parameters. The localized training of the control space depends on (i) the sampling distributions for the Dirichlet inflow boundary condition (38) and for the Neumann outflow condition (39); (ii) the number $n_{\text{loc}}^{\text{s}}$ of samples; and (iii) the number $m_0$ of retained POD modes. The localized training for the state variables depends on (i) the number $N_{\text{dd}}$ components of the networks considered; (ii) the number $n_{\text{loc}}^{\text{w}}$ of samples; and (iii) the number $n_0$ of retained POD modes for each archetypic component. The enrichment strategy depends on (i) the number $n'$ of added modes (cf. section 5.4.2). The adaptive loop depends on (i) the number **maxit** of outlet loop iterations; (ii) the number $n_{\text{train}}^{\text{glo}}$ of global configurations; (iii) the numbers $m_w$ and $m_s$ of marked components and ports; (iv) the number $n_{\text{train}}^{\text{glo}}$, $m_{\text{train}}^{\text{glo}}$ of modes added at each iteration for state and control variables. We envision that the selection of several parameters can be automated: to provide a concrete reference, the parameters $n_{\text{train}}^{\text{glo}}$, $m_{\text{train}}^{\text{glo}}$ can be updated based on a energy/projection criterion. Nevertheless, further investigations are necessary to provide actionable guidelines to select all the parameters.

### 5.4.1 Computation of the local solutions

Given the sampled port $\Gamma_j$, we solve the HF model with flux boundary conditions given by the control $\mathbf{s}^*$ on the remaining port (cf. Figure 5) in the domain $\Omega^* = \Omega_j^+ \cup \Omega_j^-$ where $\Omega_j^+, \Omega_j^-$ are the elements of the network that share $\Gamma_j$. Given the sampled component $\Omega_i$, we consider two separate strategies: (i) we solve the global hybrid model in which we replace the local ROM with the local HF model in the sampled component, or (ii) we solve the HF model in the sampled component with boundary conditions prescribed by the control estimate $\mathbf{s}^*$. The first option is significantly less computationally expensive; however, we experienced some convergence issues for very inaccurate controls $\mathbf{s}^*$. For this reason, in the numerical experiments, we rely on global hybrid solves for the first iteration of the algorithm and to fully local solves for the subsequent iterations.
<figure><img src="image_6.png" /><figcaption>Figure 5: computation of local solution. Port update: if $\Gamma_1$ is the sampled port, we solve the HF model in the components $\Omega_1 \cup \Omega_2$ with Neumann boundary conditions on the port $\Gamma_2$ given by the predicted control $\mathbf{s}^*$ (cf. Line 9, Algorithm 1). State update: if $\Omega_2$ is the sampled component, we either solve the global problem using the HF discretization in $\Omega_2$ and the ROM discretization in $\Omega_3$ (option 1), or we solve the HF model in $\Omega_2$ with Neumann boundary conditions on the ports $\Gamma_1$ and $\Gamma_2$ given by $\mathbf{s}^*$ (option 2).</figcaption></figure>

### 5.4.2 Enrichment of the state spaces

It suffices to generalize the procedure of section 4.4. We denote by $\{\tilde{\mathbf{w}}_\ell^k\}_{\ell=1}^{n_{\text{train}}^k}$ a dataset of snapshots associated with the $k$-th prototype component and the local parameters $\{\mu_\ell^k\}_{\ell=1}^{n_{\text{train}}^k}$ and $\{\text{Re}_\ell\}_{\ell=1}^{n_{\text{train}}^k}$. The dataset $\{\tilde{\mathbf{w}}_\ell^k\}_{\ell=1}^{n_{\text{train}}^k}$ is extracted by the global simulations performed in the internal loop (cf. Lines 8–14) of Algorithm 1 or from the simulations performed to generate the initial local space $\tilde{Z}_k$ (cf. Line 2). We denote by $\{\eta'_j\}_{j=1}^m$ the newly-added modes of the port space; we further recall the definitions of the local residuals and (32) boundary forms (33). Then, we define $\tilde{\mathbf{w}}_{\ell,j,q}^k$ such that (compare with (24))

$$ \tilde{\mathcal{R}}_k^{\text{hf}}(\tilde{\mathbf{w}}_\ell^k, \mathbf{z}; \mu_\ell^k, \text{Re}_\ell) + \frac{\partial \tilde{\mathcal{R}}_k^{\text{hf}}}{\partial \tilde{\mathbf{w}}_k} [\tilde{\mathbf{w}}_\ell^k, \mu_\ell^k, \text{Re}_\ell] (\tilde{\mathbf{w}}_{\ell,j,q}^k + \tilde{\mathbf{\Psi}}_{k,\text{in}}(\text{Re}_\ell) - \tilde{\mathbf{w}}_\ell^k, \mathbf{z}) + \tilde{\mathcal{E}}_{k,q}^{\text{hf}}(\eta'_j, \mathbf{z}) = 0 \quad \forall \mathbf{z} \in \tilde{\mathcal{X}}_{k,0}^{\text{hf}}, $$

for $\ell = 1, \dots, n_{\text{train}}^k$, $j = 1, \dots, m$, and $q = 1, \dots, N_{\text{port}}^k$ ($N_{\text{port}}^k = 2$ for the channel component, and $N_{\text{port}}^k = 3$ for the junction component). After having computed the snapshots $\{\tilde{\mathbf{w}}_{\ell,j,q}^k\}_{\ell,j,q}$, we update the reduced space $\tilde{Z}_k$ with $n'$ modes using POD (cf. (25)).

# 6 Numerical results

We present numerical results of the proposed method for the parameterized incompressible flow of section 2.2. The parameters are the Reynolds number and the geometric parameters $\alpha$ and $h_c$ introduced for each instantiated component. We consider a P2 FE discretization with 1281 degrees of freedom for the channel, and 3329 degrees of freedom for the junction. The regularization constant $\delta$ is set equal to $10^{-8}$.

## 6.1 HF solver

We present the HF results for the Reynolds number $Re = 100$ and the geometric configuration shown in Figure 3(b). In Figure 6(a)-(b)-(c), we show the solution to the global HF problem (i.e., without domain decomposition) for the x-direction velocity, y-direction velocity, and the pressure, respectively. Figures 6(d)-(e)-(f) illustrate the difference between the solution to the global problem and the solution to the (multi-component generalization of the) DD formulation (6b). Our new formulation exhibits high accuracy, with a pointwise error of the order of $10^{-6}$ for the three variables. Here, we employ the SQP method introduced in section 3.3.2; GNM (cf. section 3.3.1) does not converge for this value of the Reynolds number. For the solution to the DD problem, the global prediction at the interfaces is obtained by averaging the solution in the two neighboring sub-domains.
<figure><img src="image_7.png" /><figcaption>Figure 6: HF formulation. (a)-(b)-(c) Behavior of the solution to the monolithic FE problem. (d)-(e)-(f) difference between the monolithic FE solution and the DD solution based on (6b).</figcaption></figure>

In Figure 7, we present the comparison between the monolithic FE solution and the solution to the DD formulation (5). The results of Figure 7 show much larger pointwise errors for both velocity and pressure — the error for the pressure is $O(10^{-2})$ as opposed to $O(10^{-6})$. This result justifies the addition of the control $h$ for the continuity equation.

In Figure 8, we present the variable jump across the interfaces for the new formulation (6b) and the standard formulation (5). For (5), the jump of the velocity field is modest, but it is significant ($O(10^{-1})$) for the pressure. In contrast, for (6b), the jump of both velocity and pressure is extremely modest. These results further corroborate the introduction of the control $h$ for the continuity equation.

<figure><img src="image_8.png" /><figcaption>Figure 7: HF formulation. Difference between the monolithic FE solution and the DD solution based on (5).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: HF formulation. (a)-(b)-(c) interface jump of the solution to (6b). (d)-(e)-(f) interface jump of the solution to (5).</figcaption></figure>

Figure 9 investigates the effect of the choice of the penalization norm for the control. In more detail, we compare the behavior of the horizontal control $g_x$ for the first port $\Gamma_1$ in Figure 3(b) for both $L^2$ regularization and $H^1$ regularization. We observe that the use of the $H^1$ regularization dramatically reduces the spurious oscillations in the proximity of the boundaries of the domain. We further observe that, since $n = \text{vec}(1, 0)$, the control $g_x$ should equal the viscous flux $-p + \nu \frac{\partial u_x}{\partial x}$; provided that $p \gg |\nu \frac{\partial u_x}{\partial x}|$, we hence find that $g_x \approx -p$.

<figure><img src="image_10.png" /><figcaption>Figure 9: Comparison of the $H^1$ norm and the $L^2$ norm for the regularization term.</figcaption></figure>

## 6.2 MOR procedure for networks of moderate size

We now evaluate the performance of the ROM introduced in section 4 for the system configuration shown in Figure 3. Since the total number of degrees of freedom is relatively modest, we can afford to solve the multi-component generalization of (6b) with HF local models and HF control. This enables a rigorous assessment of the results. For the test cases presented in this section and in section 6.3, we choose the dimension of the original ROB (i.e., without ROB enrichment) for the state $n$ to be equal to the dimension of the ROB for the control $m$.
### 6.2.1 Performance for a fixed geometry

We freeze the value of the geometric parameters and we let the Reynolds number vary in the domain $P = [50, 150]$. We train the local ROMs based on $n_{train} = 60$ snapshots with equi-spaced parameters in $P$, and we assess the performance of the resulting CB-ROM based on $n_{test} = 10$ randomly-selected out-of-sample parameters. We measure performance of the ROMs in terms of the average out-of-sample relative prediction error for the four components:

$$E_{\text{avg}, i} := \frac{1}{n_{\text{test}}} \sum_{\mu \in \mathcal{P}_{\text{test}}} \frac{\|\mathbf{w}_i^{\text{hf}}(\mu) - \hat{\mathbf{w}}_i(\mu)\|_{\mathcal{X}_i}}{\|\mathbf{w}_i^{\text{hf}}(\mu)\|_{\mathcal{X}_i}}, \quad i = 1, \dots, N_{\text{dd}} = 4, \quad (40)$$

and the three ports:

$$E_{\text{avg}, j}^{\text{port}} := \frac{1}{n_{\text{test}}} \sum_{\mu \in \mathcal{P}_{\text{test}}} \frac{\|\mathbf{s}_j^{\text{hf}}(\mu) - \hat{\mathbf{s}}_j(\mu)\|_{\Gamma_i}}{\|\mathbf{s}_j^{\text{hf}}\|_{\Gamma_i}}, \quad j = 1, \dots, N_f = 3. \quad (41)$$

Figure 10 shows the prediction error $E_{\text{avg}, i}$ for the state $\mathbf{w}$ associated with three different local ROMs, Galerkin, Petrov-Galerkin, and minimum residual, for the four components of the network; Figure 11 shows the prediction error $E_{\text{avg}, j}^{\text{port}}$ for the control on the three ports for the same choice of the local ROM. In this test, we do not perform the enrichment of the state spaces described in section 4.4. The Galerkin method exhibits stability issues, while both the minimal residual and the Petrov-Galerkin methods perform equally well in terms of accuracy with a relative error of the order of $10^{-5}$ for $n = m = 20$. The prediction of the control variables is far less accurate: for $n = m = 20$, the relative error is $O(10^{-1})$ for port 3 and $O(10^{-2})$ for the other two ports. Nevertheless, we envision that the results can still be considered satisfactory, as illustrated by the profiles of the control $\mathbf{g}$ at port 3 shown in Figure 12, where $\xi$ is the local coordinate along the port 3.

<figure><img src="image_11.png" /><figcaption>Figure 10: performance for a fixed geometry. Behavior of the error (40) for the subdomains (no enrichment).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: performance for a fixed geometry. Behavior of the error (41) for the three ports (no enrichment).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: performance for a fixed geometry. Profile of the two components of the control g for one representative parameter value along the port 3 (no enrichment).</figcaption></figure>

In Figure 13, we illustrate the performance of the ROM when we employ the enrichment strategy discussed in section 4.4. To facilitate comparison, we include dashed lines representing the results obtained without employing ROB enrichment, which corresponds to the data presented in Figure 11. Here, the number of additional modes $n'$ (cf. section 4.4) is chosen to be equal to the dimension of the ROB of the control, $m$. The ROB enrichment strategy significantly reduces the prediction error for the control; the state prediction achieved with ROB enrichment is comparable to the case without ROB enrichment and is not provided below. We further remark that the enrichment does not contribute to increase the number of SQP iterations: to provide a concrete reference, for $m = 10$, SQP converges in six iterations for both cases.

<figure><img src="image_14.png" /><figcaption>Figure 13: performance for a fixed geometry. Behavior of the error (41) for the three ports (with enrichment).</figcaption></figure>
### 6.2.2 Performance for a parametric geometry

We incorporate the geometric parameters described in section 2.2, along with the Reynolds number. For each junction component in the network, we set $\alpha \in [\frac{\pi}{8}, \frac{\pi}{6}]$; for each channel component, we set $h_c \in [0.1, 0.3]$; finally, we consider $Re \in [50, 150]$ with $Re_{ref} = 100$. We train the ROMs based on $n_{train} = 120$ snapshots and assess performance based on $n_{test} = 10$ randomly-selected out-of-sample parameters. As for the previous test, we analyze the prediction error for both **w** and **s** associated with the different ROMs. Figure 14 illustrates the prediction error $E_{avg,i}$ for the four components, while Figure 15 shows the prediction error $E_{avg,k}^{port}$ for the three ports. Interestingly, the Galerkin method is as effective as the minimal residual and the Petrov-Galerkin methods. All three ROMs yield a state prediction relative error of approximately $O(10^{-4})$ for $n = 20$; on the other hand, the control prediction error is roughly $O(10^{-1})$ for the third port, and $O(10^{-2})$ for the other two ports, for $n = 20$. In Figure 16, we perform a comparison of ROM errors associated to the three ports, with and without the ROB enrichment strategy outlined in section 4.4. The dashed lines represent the results obtained in the absence of ROB enrichment, which correspond to the data shown in Figure 15. As for the previous test, the ROB enrichment strategy significantly improves the accuracy of the control prediction. Here, the number of additional modes $n'$ (cf. section 4.4) is chosen to be twice as large as the dimension of the ROB for the ports $m$.

<figure><img src="image_15.png" /><figcaption>Figure 14: performance for a parametric geometry. Behavior of the error (40) for the subdomains (no enrichment).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: performance for a parametric geometry. Behavior of the error (41) for the three ports (no enrichment).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16: performance for a parametric geometry. Behavior of the error (41) for the three ports (with enrichment).</figcaption></figure>

## 6.3 Localized training and adaptive enrichment

In the previous test cases, a distinct reduced space is employed for each instantiated component: the same configuration is used for both training and assessment. This approach is computationally demanding when dealing with systems that comprise a large number of components; it is also unfeasible in the presence of topology changes. To address this issue, we apply the localized training and adaptive enrichment algorithms developed in section 5.

### 6.3.1 Application to networks with four components

We apply the localized training strategy of sections 5.2 and 5.3, for the same test set of section 6.2.2. In order to build the reduced space for the control, we consider 60 randomly selected boundary conditions for each connection described in section 5.2; on the other hand, we generate the reduced space for the state using 20 randomly-sampled networks with four components and the reduced space for the control.

Figure 17 presents the prediction error $E_{\text{avg},i}$ for the four components, while Figure 18 shows the prediction error $E_{\text{avg},k}^{\text{port}}$ for the three ports; we do not rely on the enrichment of the state space (cf. section 4.4). The results are comparable to those obtained in section 6.2.2 with slight deterioration in accuracy. Figure 19 displays the ROM errors for the three ports using ROB enrichment ($n' = 2m$), as represented by the solid line. The results exhibit significant improvement when compared to those obtained without the use of the state space enrichment, as illustrated by the dashed lines, which correspond to the data shown in Figure 18.
<figure><img src="image_18.png" /></figure>

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

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

<figure><img src="image_21.png" /><figcaption>Figure 17: localized training for networks with four components. State prediction error for the four sub-domains.</figcaption></figure>

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

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

<figure><img src="image_24.png" /><figcaption>Figure 18: localized training for networks with four components. Control prediction error for the three ports (without enrichment).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 19: localized training for networks with four components. Control prediction error for the three ports (with enrichment).</figcaption></figure>

### 6.3.2 Application to networks with ten components

We apply the full training procedure described in Algorithm 1 to $n_{test} = 10$ randomly selected configurations with ten components. As for the previous test case, we consider independent geometric variations for each instantiated component and we consider $Re \in [50, 150]$. We only present results for local Galerkin ROMs: the results obtained using minimum residual projection are comparable and are hence omitted.

Figure 20 shows the local relative error for the state and for the control, over the test set for the CB-ROM based on localized training: we use the same dataset considered in section 6.3.1 with port-based enrichment (cf. section 4.4). We observe that the error is roughly 10% for both state and control and does not decrease as we increase the number of modes.

Figure 21 shows the results for the full application of Algorithm 1. We initialize the algorithm with a ROB of size $m_0 = 10$ for the control using localized training; we apply the strategy of section 5.3, together with port-based enrichment, to find reduced spaces for the state of size $n_0 = 10 + 10$ for each component. Then, we apply adaptive enrichment: we consider $n_{train}^{glo} = 50$ global randomly-selected configurations with ten components; we mark $m_s = 1$ port, and $m_w = 1 - 3$ components of each type (specifically, we mark 1 component with the largest error of each type, along with the 2 adjacent components of the marked port). Then, we augment the bases for state and control with $n^{glo} = m^{glo} = 10$ modes. We do not apply the port-based enrichment strategy after each iteration of the adaptive loop (cf. Line 17). If Figure 21, iteration *it* corresponds to local ROBs of size $m = 10(it + 1)$ and $n = 10(it + 2)$.

We observe that the enrichment strategy clearly enhances the performance of the CB-ROM. This result empirically demonstrates the importance of adaptive enrichment when dealing with nonlinear PDEs.

<figure><img src="image_26.png" /><figcaption>Figure 20: application to networks with ten components. Boxplots of the out-of-sample error for reduced spaces of several sizes obtained using localized training (without adaptive enrichment).</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 21: application to networks with ten components. Boxplots of the out-of-sample error for several iterations of Algorithm 1.</figcaption></figure>

# 7 Conclusions

We developed and numerically validated a component-based model order reduction procedure for incompressible flows governed by the Navier-Stokes equations. Our point of departure is the optimization-based formulation of [29]: we included an additional control variable $h$ for the continuity equation that weakly enforces the continuity of pressure at interfaces; furthermore, we modified the regularization term to damp spurious oscillations of the control. We relied on sequential quadratic programming to solve the nonlinear optimization problem: at each iteration of the procedure, we relied on static condensation of the local degrees of freedom to enable trivial parallelism of the local solves and avoid the introduction of Lagrange multipliers. We relied on projection-based (Galerkin and Petrov-Galerkin) ROMs to speed up the solution to the local subproblems and we exploited port reduction to reduce the cost of the global problem. Finally, we adapted the localized training and adaptive enrichment strategy of [36] to build the local approximation spaces without the need for expensive global HF solves.

We illustrated the many pieces of our methodology for a parametric steady Navier-Stokes problem at moderate ($O(10^2)$) Reynolds number. The new DD formulation enables much tighter control of the discrepancy between the FE monolithic solver and the DD solution. LSPG projection is superior to Galerkin projection in the absence of geometric variability; interestingly, Galerkin and LSPG projection show comparable performance for all the test cases that involve varying geometries. The port-based enrichment of the state space (cf. section 4.4) is key to adequately approximate the control variables. The localized training strategy discussed in this paper leads to poor reconstructions of the state; adaptive enrichment driven by local error indicators is hence necessary to achieve accurate reconstructions.

In the future, we plan to extend our method to a broader class of problems including multi-physics (fluid-structure interaction) and unsteady problems, and to more challenging (higher-Reynolds, three-dimensional) test cases. Towards this end, it is of paramount importance to devise effective hyper-reduction techniques to speed up local solves and also the assembly of the objective function. We also plan to combine first-principles models with data-fitted models to enhance the flexibility of the method.

# Acknowledgements

The work of Lei Zhang is supported by the Fundamental Research Funds for the Central Universities of Tongji University.

# A Stabilized FE formulation

For completeness, we review the stabilized finite element formulation employed in the numerical results; we refer to [38, 39] for a thorough review of stabilized FE methods for incompressible flows. We denote by $\{D_k^i\}_k$ the elements of the mesh of $\Omega_i$; we further denote by $h_{k,i}$ the size of the $k$-th element of the mesh, and by $r$ the degree of the polynomials.
We consider the residual:

$$ \mathcal{R}_i^{\text{hf}}(\mathbf{u}_i, p_i, \mathbf{v}, q) = \mathcal{R}_i(\mathbf{u}_i, p_i, \mathbf{v}, q) + \mathcal{R}_i^{\text{supg}}(\mathbf{u}_i, p_i, \mathbf{v}) + \mathcal{R}_i^{\text{pspg}}(\mathbf{u}_i, p_i, \mathbf{v}) + \mathcal{R}_i^{\text{lsic}}(\mathbf{u}_i, \mathbf{v}), \quad \forall (\mathbf{v}, q) \in \mathcal{X}_{i,0}. \quad (42a) $$

The form $\mathcal{R}_i$ corresponds to the local residual introduced in (2d), while the other three terms are designed to improve the stability of the discrete problem. The form $\mathcal{R}_i^{\text{supg}}$ corresponds to the Streamline upwind Petrov-Galerkin (SUPG, [53]) stabilization, which is designed to handle advection-dominated flows,

$$ \mathcal{R}_i^{\text{supg}}(\mathbf{u}, p, \mathbf{v}) = \sum_k \int_{\mathcal{D}_k^i} ((\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p - \nu \Delta \mathbf{u} - \mathbf{f}) (\tau_{\text{supg}}\mathbf{u} \cdot \nabla \mathbf{v}) dx; \quad (42b) $$

the form $\mathcal{R}_i^{\text{pspg}}$ is the Pressure-Stabilized Petrov-Galerkin (PSPG) term [40] that is added to the mass conservation equation to eliminate spurious modes in the pressure solution when considering the same polynomial order for pressure and velocity,

$$ \mathcal{R}_i^{\text{pspg}}(\mathbf{u}, p, q) - \sum_k \int_{\mathcal{D}_k^i} \tau_{\text{pspg}} ((\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p - \nu \Delta \mathbf{u} - \mathbf{f}) \cdot \nabla q \, dx; \quad (42c) $$

finally, $\mathcal{R}_i^{\text{lsic}}$ is the least-squares incompressibility constraint (LSIC) stabilization term that is added to the momentum equation to improve accuracy and conditioning of the discrete problem [54, 55, 56],

$$ \mathcal{R}_i^{\text{lsic}}(\mathbf{u}, \mathbf{v}) = \sum_k \int_{\mathcal{D}_k^i} (\nabla \cdot \mathbf{u}) \tau_{\text{lsic}} \nabla \cdot \mathbf{v} \, dx. \quad (42d) $$

In the numerical experiments, following [57], we select the parameters $\tau_{\text{supg}}$, $\tau_{\text{pspg}}$, and $\tau_{\text{lsic}}$ as $\tau_{\text{supg}} = \tau_{\text{pspg}} = \alpha_{\text{supg}} \left[ \left( \frac{2|\mathbf{u}_i|}{h_{k,i}} \right)^2 + 9 \left( \frac{4\nu}{h_{k,i}^2} \right)^2 \right]^{-\frac{1}{2}}$, $\tau_{\text{lsic}} = \frac{h_{k,i}^2}{r^2 \tau_{\text{supg}}}$, where $0 \le \alpha_{\text{supg}} \le 1$ is a constant that enables the adjustment of $\tau_{\text{supg}}$ for higher-order elements. In the PTC formulation (cf. (13)), we modify the coefficients $\tau_{\text{supg}}$ and $\tau_{\text{pspg}}$ to account for the time step $\tau_{\text{supg}} = \tau_{\text{pspg}} = \alpha_{\text{supg}} \left[ \left( \frac{2}{\Delta t} \right)^2 + \left( \frac{2|\mathbf{u}_i|}{h_{k,i}} \right)^2 + 9 \left( \frac{4\nu}{h_{k,i}^2} \right)^2 \right]^{-\frac{1}{2}}$.

## B Justification of the pressure jump in the minimization formulation

We consider the configuration depicted in Figure 1 and we assume that the meshes of $\Omega_1$ and $\Omega_2$ are conforming on $\Gamma_0$. We denote by $\{\phi_i\}_{i=1}^{N_w}$ the Lagrangian basis associated with the global space $\mathcal{X}^{\text{hf}}$; we denote by $\mathcal{I}_1, \mathcal{I}_2$ the degrees of freedom associated with the domains $\Omega_1$ and $\Omega_2$, respectively. We further denote by $\mathcal{I}_0 = \mathcal{I}_1 \cap \mathcal{I}_2$ the nodes on the interface $\Gamma_0$; we introduce the local and global Dirichlet nodes $\mathcal{I}_{1,\text{dir}}, \mathcal{I}_{2,\text{dir}} \subset \{1, \dots, N_w\}$ and $\mathcal{I}_{\text{dir}} = \mathcal{I}_{1,\text{dir}} \cap \mathcal{I}_{2,\text{dir}}$. By construction, $\mathcal{I}_{\text{dir}} \cap \mathcal{I}_0 = \emptyset$ (cf. Figure 1). Finally, we recall the definition of the global problem

$$ \mathbf{w}^{\text{hf}}(1:2)|_{\Gamma_{\text{dir}}} = \mathbf{\Phi}_{\mathbf{u}_{\text{in}}}, \quad \mathcal{R}^{\text{hf}}(\mathbf{w}^{\text{hf}}, \mathbf{z}) = 0 \quad \forall \mathbf{z} \in \mathcal{X}_0^{\text{hf}}, \quad (43) $$

and the two local problems

$$ \mathbf{w}_i^{\text{hf}}(1:2)|_{\Gamma_{i,\text{dir}}} = \mathbf{\Phi}_{i,\mathbf{u}_{\text{in}}}, \quad \mathcal{R}_i^{\text{hf}}(\mathbf{w}_i^{\text{hf}}, \mathbf{z}) + (-1)^i \int_{\Gamma_0} \mathbf{s} \cdot \mathbf{z} \, dx = 0 \quad \forall \mathbf{z} \in \mathcal{X}_{i,0}^{\text{hf}}, \quad i=1,2; \quad (44) $$

which depend on the control $\mathbf{s}$.

Since the meshes are conforming, it is possible to verify that

$$ \mathcal{X}_i^{\text{hf}} = \operatorname{span}\{\phi_j|_{\Omega_i} : j \in \mathcal{I}_i\}, \quad i=1,2. $$

Furthermore, the global residual can be expressed as<fnref n="1" />

$$ \mathcal{R}^{\text{hf}}(\mathbf{w}^{\text{hf}}, \mathbf{z}) = \mathcal{R}_1^{\text{hf}}(\mathbf{w}^{\text{hf}}|_{\Omega_1}, \mathbf{z}|_{\Omega_1}) + \mathcal{R}_2^{\text{hf}}(\mathbf{w}^{\text{hf}}|_{\Omega_2}, \mathbf{z}|_{\Omega_2}) \quad \forall \mathbf{z} \in \mathcal{X}_0^{\text{hf}}. \quad (45) $$

Identity (45) implies that

$$ \mathcal{R}_i^{\text{hf}}(\mathbf{w}^{\text{hf}}, \phi_j) = 0 \quad \forall j \in \mathcal{I}_i \setminus \mathcal{I}_0, \quad i=1,2; \quad (46) $$

therefore, since the bilinear form $a(\mathbf{w}, \mathbf{z}) = \int_{\Gamma_0} \mathbf{w} \cdot \mathbf{z} \, dx$ is coercive in $\mathcal{Y} := \operatorname{span}\{\phi_j : j \in \mathcal{I}_0\}$, there exists a unique $\mathbf{s}^* = \operatorname{vec}(\mathbf{g}^*, \mathbf{h}^*) \in \mathcal{Y}$ such that

$$ \mathcal{R}_i^{\text{hf}}(\mathbf{w}^{\text{hf}}, \mathbf{z}) + (-1)^i \int_{\Gamma_0} \mathbf{s}^* \cdot \mathbf{z} \, dx = 0 \quad \forall \mathbf{z} \in \mathcal{X}_{i,0}^{\text{hf}}, \quad (47) $$

for $i=1,2$.

Exploiting the previous discussion, we can prove the following result.

<footnote n="1">¹The proof of (45) exploits the expressions of the residuals (42) and (2d). We omit the details. We further emphasize that at the right hand side of (45) we should use notation $\mathcal{R}_i^{\text{hf}}(\mathbf{w}^{\text{hf}}|_{\Omega_i}, \mathbf{z}|_{\Omega_i})$ for $i=1,2$.</footnote>
**Lemma 1.** Let $\mathbf{w}^{\text{hf}}$ be a solution to (43). The following hold.

1. The triplet $(\mathbf{w}^{\text{hf}}|_{\Omega_1}, \mathbf{w}^{\text{hf}}|_{\Omega_2}, \mathbf{s}^*)$ where $\mathbf{s}^*$ satisfies (47) is a global minimum of (6c) for $\delta = 0$.

2. Any global minimum of (6c) for $\delta = 0$ solves (43); in particular, if the solution $\mathbf{w}^{\text{hf}}$ to (43) is unique, the optimization problem (6b) admits a unique solution for $\delta = 0$.

*Proof.* Equation (47) implies that the triplet $(\mathbf{w}^{\text{hf}}|_{\Omega_1}, \mathbf{w}^{\text{hf}}|_{\Omega_2}, \mathbf{s}^*)$ satisfies the constraints of (6b) (cf. (44)); since $\mathbf{w}^{\text{hf}}$ is continuous, the objective function of (6b) (cf. (6c)) is equal to zero for $\delta = 0$. Since the function (6c) is non-negative, we conclude that $(\mathbf{w}^{\text{hf}}|_{\Omega_1}, \mathbf{w}^{\text{hf}}|_{\Omega_2}, \mathbf{s}^*)$ is a global minimum of (6b).

Exploiting the first part of the proof, we find that any global minimum $(\mathbf{w}_1, \mathbf{w}_2, \mathbf{s})$ of (6b) satisfies $\mathcal{F}_{\delta=0}(\mathbf{w}_1, \mathbf{w}_2, \mathbf{s})$
$\mathcal{F}_{\delta=0}(\mathbf{w}^{\text{hf}}|_{\Omega_1}, \mathbf{w}^{\text{hf}}|_{\Omega_2}, \mathbf{s}^*) = 0$. This implies that the function $\mathbf{w}: \Omega \to \mathbb{R}^3$ such that $\mathbf{w}|_{\Omega_1} = \mathbf{w}_1$ and $\mathbf{w}|_{\Omega_2} = \mathbf{w}_2$
is continuous and belongs to $\mathcal{X}^{\text{hf}}$. Recalling (44), we have that $\mathbf{w}$ satisfies $\mathbf{w}(1:2)|_{\Gamma_{\text{dir}}} = \mathbf{\Phi}_{\mathbf{u}_{\text{in}}}$. Furthermore,
since $\mathbf{z}|_{\Omega_i} \in \mathcal{X}_{i,0}$ for any $\mathbf{z} \in \mathcal{X}^{\text{hf}}$, we have that

$$ \mathcal{R}^{\text{hf}}(\mathbf{w}, \mathbf{z}) \stackrel{(45)}{=} \mathcal{R}_1^{\text{hf}}(\mathbf{w}_1, \mathbf{z}|_{\Omega_1}) + \mathcal{R}_2^{\text{hf}}(\mathbf{w}_2, \mathbf{z}|_{\Omega_2}) \stackrel{(44)}{=} - \int_{\Gamma_0} \mathbf{s} \cdot \mathbf{z} \, dx + \int_{\Gamma_0} \mathbf{s} \cdot \mathbf{z} \, dx = 0, $$

which is (43). We conclude that $\mathbf{w}$ solves (43). If the solution to (43) is unique, exploiting the previous argument, any solution $(\mathbf{w}_1, \mathbf{w}_2, \mathbf{s})$ should satisfy $\mathbf{w}_1 = \mathbf{w}^{\text{hf}}|_{\Omega_1}$ and $\mathbf{w}_2 = \mathbf{w}^{\text{hf}}|_{\Omega_2}$. Furthermore, since the solution to (47) is unique, we also find $\mathbf{s} = \mathbf{s}^*$. In conclusion, (6b) has a unique global minimum. $\square$

Lemma 1 illustrates the connection between the monolithic problem and the solution to the optimization problem (6b); the well-posedness analysis in [32] shows that in the continuous limit (i.e., $N_{\mathbf{w}} \to \infty$) $h^* = 0$; nevertheless, in general $h^* \neq 0$ for finite-dimensional discretizations. To illustrate this fact, we consider the solution to the Stokes problem (see Figure 1 for the definitions of the boundary subdomains)

$$ \left\{ \begin{array}{ll} -\Delta \mathbf{u} + \nabla p = \text{vec}(x_1, \cos(x_2^2)) & \text{in } \Omega = (0, 1)^2, \\ \nabla \cdot \mathbf{u} = 0 & \text{in } \Omega, \\ \mathbf{u}|_{\Gamma_{\text{dir}}} = \text{vec}((1-x_2)x_2, 0), \quad \mathbf{u}|_{\Gamma_{\text{dir}^0}} = 0, \quad (\nabla \mathbf{u} - p\mathbf{I})\mathbf{n}|_{\Gamma_{\text{neu}}} = 0, \end{array} \right. $$

based on a P3-P2 Taylor-Hood discretization for three meshes of increasing size. Figure 22(a) shows the final mesh used for computations whereas the blue dots indicate the interface $\Gamma_0$; Figure 22(b) shows the behavior of $h^*$ for three meshes with $4^2, 8^2, 16^2$ global elements: as expected, as we increase the size of the mesh, the magnitude of $h^*$ decreases.

<figure><img src="image_28.png" /><figcaption>Figure 22: justification of the pressure jump; Stokes model problem.</figcaption></figure>

## C Justification of the enrichment strategy

We consider the algebraic problem:

$$ \min_{\mathbf{w} \in \mathbb{R}^N, \mathbf{s} \in \mathbb{R}^M} |\mathbf{C}\mathbf{w} - \mathbf{b}| \quad \text{s.t.} \quad \mathbf{A}\mathbf{w} + \mathbf{B}\mathbf{s} = \mathbf{f}, \qquad (48) $$
with $\mathbf{A} \in \mathbb{R}^{N \times N}$, $\mathbf{B} \in \mathbb{R}^{N \times M}$, $\mathbf{C} \in \mathbb{R}^{M \times N}$, $\mathbf{b} \in \mathbb{R}^{M}$, $\mathbf{f} \in \mathbb{R}^{N}$ and $N > M$. If $\mathbf{A}$ is full rank, any solution $(\mathbf{w}^{*}, \mathbf{s}^{*})$ to (48) satisfies $\mathbf{w}^{*} = \mathbf{A}^{-1}(\mathbf{f} - \mathbf{B}\mathbf{s}^{*})$ and $\mathbf{s}^{*} = \arg \min_{\mathbf{s} \in \mathbb{R}^{M}} |\mathbf{D}\mathbf{s} - \mathbf{c}|$, with $\mathbf{D} = \mathbf{C}\mathbf{A}^{-1}\mathbf{B}$ and $\mathbf{c} = \mathbf{b} - \mathbf{A}^{-1}\mathbf{f}$. Therefore, provided that $\mathbf{C}$ is full rank, (48) is well-posed if and only if $\mathbf{A}^{-1}\mathbf{B}$ is full rank.

Let $\mathbf{Z} = [\zeta_1, \dots, \zeta_n] \in \mathbb{R}^{N \times n}$, $\mathbf{W} = [\eta_1, \dots, \eta_m] \in \mathbb{R}^{M \times m}$ and $\mathbf{Y} \in \mathbb{R}^{N \times n}$ be orthogonal matrices with $n < N$ and $m < M$; exploiting these definitions, we define the projected problem

$$ \min_{\alpha \in \mathbb{R}^n, \beta \in \mathbb{R}^m} |\bar{\mathbf{C}}\alpha - \mathbf{b}| \quad \text{s.t.} \quad \bar{\mathbf{A}}\alpha + \bar{\mathbf{B}}\beta = \bar{\mathbf{f}}, \qquad (49) $$

with $\bar{\mathbf{C}} = \mathbf{C}\mathbf{Z}$, $\bar{\mathbf{A}} = \mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}$, $\bar{\mathbf{B}} = \mathbf{Y}^{\top}\mathbf{B}\mathbf{W}$ and $\bar{\mathbf{f}} = \mathbf{Y}^{\top}\mathbf{f}$. It is straightforward to prove the following result: here, $\text{col}(\mathbf{X})$ denotes the linear space spanned by the columns of the matrix $\mathbf{X}$, while $\text{orth}(\mathbf{X})$ is the orthogonal matrix that is obtained by orthogonalizing the columns of $\mathbf{X}$.

**Lemma 2.** Let $\text{col}(\mathbf{A}^{-1}\mathbf{B}\mathbf{W}) \subset \text{col}(\mathbf{Z})$ and let $\mathbf{Y} = \text{orth}(\mathbf{A}\mathbf{Z})$. Then, $\bar{\mathbf{A}}^{-1}\bar{\mathbf{B}}$ is full rank, and (49) is well-posed.

*Proof.* We first prove that $\bar{\mathbf{A}} \in \mathbb{R}^{n \times n}$ is invertible. By contradiction, there exists $\alpha \in \mathbb{R}^n$ such that $\bar{\mathbf{A}}\alpha = 0$. Since $\mathbf{Y} = \text{orth}(\mathbf{A}\mathbf{Z})$, there exists $\beta \in \mathbb{R}^n$ such that $\mathbf{A}\mathbf{Z}\alpha = \mathbf{Y}\beta$. We hence find

$$ 0 = \beta^{\top}\bar{\mathbf{A}}\alpha = (\mathbf{Y}\beta)^{\top}\mathbf{A}\mathbf{Z}\alpha = |\mathbf{A}\mathbf{Z}\alpha|^2. $$

The latter implies that $\mathbf{Z}\alpha$ is a non-trivial element of the kernel of $\mathbf{A}$: this is in contradiction with the hypothesis that $\mathbf{A}$ is invertible.

Exploiting the same argument, we prove that $\bar{\mathbf{B}}$ is full rank. By contradiction, there exists $\beta \in \mathbb{R}^m$ such that $\bar{\mathbf{B}}\beta = 0$. Since $\text{col}(\mathbf{Y}) = \text{col}(\mathbf{A}\mathbf{Z})$ and $\text{col}(\mathbf{A}^{-1}\mathbf{B}\mathbf{W}) \subset \text{col}(\mathbf{Z})$, there exist $\alpha, \alpha' \in \mathbb{R}^n$ such that $\mathbf{B}\mathbf{W}\beta = \mathbf{A}\mathbf{Z}\alpha' = \mathbf{Y}\alpha$. We hence find

$$ 0 = \alpha^{\top}\bar{\mathbf{B}}\beta = (\mathbf{Y}\alpha)^{\top}\mathbf{B}\mathbf{W}\beta = |\mathbf{B}\mathbf{W}\beta|^2. $$

The latter implies that $\mathbf{W}\beta$ is a non-trivial element of the kernel of $\mathbf{B}$: this is in contradiction with the hypothesis that $\mathbf{B}$ is full-rank. $\square$

Lemma 2 provides a rigorous justification of the enrichment strategy in section 4.4. The matrix $\tilde{\mathbf{W}} = -\mathbf{A}^{-1}\mathbf{B}\mathbf{W}$ corresponds to the derivative of the state $\mathbf{w}$ with respect to the control $\hat{\mathbf{s}} = \mathbf{W}\beta$; the columns $\tilde{\mathbf{w}}_1, \dots, \tilde{\mathbf{w}}_m$ of the matrix $\tilde{\mathbf{W}}$ satisfy

$$ \mathbf{A}\tilde{\mathbf{w}}_k + \mathbf{B}\eta_k = 0, \quad k = 1, \dots, m, $$

which corresponds to (24). Similarly, as discussed in [46, 47], the choice of the test space in Lemma 2 is consistent with (20a).

## References

[1] A. Quarteroni, A. Manzoni, F. Negri, Reduced basis methods for partial differential equations: an introduction, Springer, 2015.

[2] J. S. Hesthaven, G. Rozza, B. Stamm, Certified reduced basis methods for parametrized partial differential equations, Springer, 2016.

[3] G. Rozza, D. B. P. Huynh, A. T. Patera, Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations, Archives of Computational Methods in Engineering 15 (2008) 229–275. doi:10.1007/s11831-008-9019-9.

[4] P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, W. Schilders, L. M. Silveira, Model order reduction: volume 2, snapshots-based methods and algorithms, De Gruyter, 2021.

[5] G. Berkooz, P. Holmes, J. Lumley, The proper orthogonal decomposition in the analysis of turbulent flows, Annual review of fluid mechanics 25 (1) (1993) 539–575. doi:10.1146/annurev.fl.25.010193.002543.

[6] S. Volkwein, Model reduction using proper orthogonal decomposition, Lecture Notes, Institute of Mathematics and Scientific Computing, University of Graz. see math.uni-konstanz.de/numerik/personen/volkwein/teaching/POD-Vorlesung.pdf 1025 (2011).

[7] K. Carlberg, C. Bou-Mosleh, C. Farhat, Efficient non-linear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approximations: efficient non-linear model reduction, International Journal for Numerical Methods in Engineering 86 (2011) 155–181. doi:10.1002/nme.3050.
[8] K. Carlberg, M. Barone, H. Antil, Galerkin v. least-squares Petrov-Galerkin projection in nonlinear model reduction, Journal of Computational Physics 330 (2017) 693-734. doi:10.1016/j.jcp.2016.10.033.

[9] D. B. P. Huynh, D. J. Knezevic, A. T. Patera, A static condensation reduced basis element method: approximation and a posteriori error estimation, ESAIM: Mathematical Modelling and Numerical Analysis 47 (2013) 213-251. doi:10.1051/m2an/2012022.

[10] I. Maier, B. Haasdonk, A Dirichlet-Neumann reduced basis method for homogeneous domain decomposition problems, Applied Numerical Mathematics 78 (2014). doi:10.1016/j.apnum.2013.12.001.

[11] J. Barnett, I. Tezaur, A. Mota, The Schwarz alternating method for the seamless coupling of nonlinear reduced order models and full order models, arXiv:2210.12551 (2022).

[12] Y. Maday, E. M. Rønquist, A reduced-basis element method, Journal of Scientific Computing 17 (2002) 447-459. doi:10.1023/A:1015197908587.

[13] Y. Maday, E. M. Rønquist, The reduced basis element method: application to a thermal fin problem, SIAM Journal on Scientific Computing 26 (1) (2004) 240-258. doi:10.1137/S1064827502419932.

[14] A. E. Løvgren, Y. Maday, E. M. Rønquist, A reduced basis element method for the steady Stokes problem, ESAIM: Mathematical Modelling and Numerical Analysis 40 (3) (2006) 529-552. doi:10.1051/m2an:2006021.

[15] L. Pegolotti, M. R. Pfaller, A. L. Marsden, S. Deparis, Model order reduction of flow based on a modular geometrical approximation of blood vessels, Computer Methods in Applied Mechanics and Engineering 380 (2021) 113762. doi:10.1016/j.cma.2021.113762.

[16] J. L. Eftang, A. T. Patera, Port reduction in parametrized component static condensation: approximation and a posteriori error estimation, International Journal for Numerical Methods in Engineering 96 (5) (2013) 269-302. doi:10.1002/nme.4543.

[17] A. Benaceur, A. T. Patera, Port-reduced reduced-basis component method for steady state Navier-Stokes and passive scalar equations (2022). doi:https://doi.org/10.21203/rs.3.rs-1473869/v1.

[18] R. Craig, M. Bampton., Coupling of substructures for dynamic analyses, AIAA Journal 6 (7) (1968) 1313-1319. doi:10.2514/3.4741.

[19] L. Iapichino, A. Quarteroni, G. Rozza, Reduced basis method and domain decomposition for elliptic problems in networks and complex parametrized geometries, Computers & Mathematics with Applications 71 (1) (2016) 408-430. doi:10.1016/j.camwa.2015.12.001.

[20] A. Quarteroni, A. Valli, Domain decomposition methods for partial differential equations, Oxford University Press, 1999.

[21] A. Mota, I. Tezaur, C. Alleman, The Schwarz alternating method in solid mechanics, Computer Methods in Applied Mechanics and Engineering 319 (2017) 19-51. doi:10.1016/j.cma.2017.02.006.

[22] M. Buffoni, H. Telib, A. Iollo, Iterative methods for model reduction by domain decomposition, Computers & Fluids 38 (6) (2009) 1160-1167. doi:10.1016/j.compfluid.2008.11.008.

[23] E. Zappon, A. Manzoni, P. Gervasio, A. Quarteroni, A reduced order model for domain decompositions with non-conforming interfaces, arXiv:2206.09618v2 (2022).

[24] N. Discacciati, J. S. Hesthaven, Localized model order reduction and domain decomposition methods for coupled heterogeneous systems, International Journal for Numerical Methods in Engineering (2023). doi:https://doi.org/10.1002/nme.7295.

[25] S. Deparis, D. Fortia, P. Gervasio, A. Quarteroni, INTERNODES: an accurate interpolation-based method for coupling the Galerkin solutions of PDEs on subdomains featuring non-conforming interfaces, Computers and Fluids 141 (2016) 22-41. doi:10.1016/j.compfluid.2016.03.033.

[26] M. Bergmann, A. Ferrero, A. Iollo, E. Lombardi, A. Scardigli, H. Telib, A zonal Galerkin-free POD model for incompressible flows, Journal of Computational Physics 352 (2018) 301-325. doi:10.1016/j.jcp.2017.10.001.

[27] A. Iollo, G. Sambataro, T. Taddei, A one-shot overlapping Schwarz method for component-based model reduction: application to nonlinear elasticity, Computer Methods in Applied Mechanics and Engineering 404 (2023) 115786. doi:10.1016/j.cma.2022.115786.
[28] G. Sambataro, Component-based model order reduction procedures for large scale THM systems, Ph.D. thesis, University of Bordeaux (2022).

[29] M. D. Gunzburger, J. S. Peterson, H. Kwon, An optimization based domain decomposition method for partial differential equations, Computers & Mathematics with Applications 37 (10) (1999) 77-93. doi: 10.1016/S0898-1221(99)00127-3.

[30] M. D. Gunzburger, M. Heinkenschloss, H. K. Lee, Solution of elliptic partial differential equations by an optimization-based domain decomposition method, Applied Mathematics and Computation 113 (2000) 111-139. doi:10.1016/S0096-3003(99)00076-4.

[31] M. D. Gunzburger, J. Lee, A domain decomposition method for optimization problems for partial differential equations, Computers & Mathematics with Applications 40 (2000) 177-192. doi:10.1016/S0898-1221(00)00152-8.

[32] M. D. Gunzburger, H. K. Lee, An optimization-based domain decomposition method for the Navier-Stokes equations, SIAM Journal on Numerical Analysis 37 (5) (2000) 1455-1480. doi:10.1137/S0036142998332864.

[33] A. Quarteroni, F. Pasquarelli, A. Valli, Heterogeneous domain decomposition: principles, algorithms, applications, in: Fifth International Symposium on Domain Decomposition Methods for Partial Differential Equations, 1992, pp. 129-150.

[34] Y. Gorb, D. Kurzanova, Heterogeneous domain decomposition method for high contrast dense composites, Journal of Computational and Applied Mathematics 337 (2018) 135-149. doi:10.1016/j.cam.2018.01.008.

[35] I. Prusak, M. Nonino, D. Torlo, F. Ballarin, G. Rozza, An optimisation-based domain-decomposition reduced order model for the incompressible Navier-Stokes equations, Computers & Mathematics with Applications 151 (2023) 172-189. doi:10.1016/j.camwa.2023.09.039.

[36] K. Smetana, T. Taddei, Localized model reduction for nonlinear elliptic partial differential equations: localized training, partition of unity, and adaptive enrichment, SIAM Journal on Scientific Computing 45 (3) (2023). doi:10.1137/22M148402X.

[37] A. Buhr, C. Engwer, M. Ohlberger, S. Rave, ArbiLoMod, a simulation technique designed for arbitrary local modifications, SIAM Journal on Scientific Computing 39 (4) (2017) A1435-A1465.

[38] T. Tezduyar, Stabilized finite element formulations for incompressible flow computations, Advances in Applied Mechanics 28 (1991) 1-44. doi:10.1016/S0065-2156(08)70153-4.

[39] T. Tezduyar, Y. Osawa, Finite element stabilization parameters computed from element matrices and vectors, Computer Methods in Applied Mechanics and Engineering 190 (2000) 411-430. doi:10.1016/S0045-7825(00)00211-5.

[40] T. J. Hughes, L. P. Franca, M. Balestra, A new finite element formulation for computational fluid dynamics: v. circumventing the Babuška-Brezzi condition: a stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Computer Methods in Applied Mechanics and Engineering 59 (1) (1986) 85-99. doi:10.1016/0045-7825(86)90025-3.

[41] X. Na, X. Xu, An optimal Robin-Robin domain decomposition method for Stokes equations, Applied Numerical Mathematics 171 (2022) 426-441. doi:10.1016/j.apnum.2021.09.015.

[42] K. M. Washabaugh, M. J. Zahr, C. Farhat, On the use of discrete nonlinear reduced-order models for the prediction of steady-state flows past parametrically deformed complex geometries, in: 54th AIAA Aerospace Sciences Meeting, 2016, p. 1814.

[43] J. Nocedal, S. J. Wright, Numerical optimization, Springer, 2006.

[44] C. T. Kelley, D. E. Keyes, Convergence analysis of pseudo-transient continuation, SIAM Journal on Numerical Analysis 35 (2) (1998) 508-523. doi:10.1137/S0036142996304796.

[45] L. Sirovich, Turbulence and the dynamics of coherent structures. I. Coherent structures, Quarterly of applied mathematics 45 (3) (1987) 561-571.

[46] T. Taddei, L. Zhang, Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs, ESAIM: Mathematical Modelling and Numerical Analysis 55 (2021) 99-130. doi:10.1051/m2an/2020073.
[47] T. Taddei, L. Zhang, A discretize-then-map approach for the treatment of parameterized geometries in model order reduction, Computer Methods in Applied Mechanics and Engineering 384 (2021) 113956. doi:10.1016/j.cma.2021.113956.

[48] D. Ryckelynck, Hyper-reduction of mechanical models involving internal variables, International Journal for Numerical Methods in Engineering 77 (1) (2009) 75–89. doi:10.1002/nme.2406.

[49] C. Farhat, S. Grimberg, A. Manzoni, A. Quarteroni, 5 Computational bottlenecks for PROMs: pre-computation and hyperreduction, De Gruyter, Berlin, Boston, 2021, pp. 181–244. doi:doi:10.1515/9783110671490-005.

[50] B. Haasdonk, Chapter 2: reduced basis methods for parametrized PDEs—a tutorial introduction for stationary and instationary problems, in: Model Reduction and Approximation, 2017, pp. 65–136. doi:10.1137/1.9781611974829.ch2.

[51] D.-Q. Bui, P. Mollo, F. Nobile, T. Taddei, A component-based data assimilation strategy with applications to vascular flows, ESAIM: Proceedings and Surveys 73 (2023) 89–106.

[52] C. Hoang, Y. Choi, K. Carlberg, Domain-decomposition least-squares Petrov–Galerkin (DD-LSPG) non-linear model reduction, Computer Methods in Applied Mechanics and Engineering 384 (2021) 113997. doi:10.1016/j.cma.2021.113997.

[53] A. N. Brooks, T. J. R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering 32 (1) (1982) 199–259. doi:https://doi.org/10.1016/0045-7825(82) 90071-8.

[54] L. P. Franca, S. L. Frey, Stabilized finite element methods: ii. the incompressible Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering 99 (2-3) (1992) 209–233. doi:10.1016/0045-7825(92)90041-H.

[55] T. Gelhard, G. Lube, M. A. Olshanskii, J. H. Starcke, Stabilized finite element schemes with LBB-stable elements for incompressible flows, Journal of Computational and Applied Mathematics 177 (2) (2005) 243–267. doi:10.1016/j.cam.2004.09.017.

[56] M. Braack, E. Burman, V. John, G. Lube, Stabilized finite element methods for the generalized Oseen problem, Computer Methods in Applied Mechanics and Engineering 196 (4-6) (2007) 853–866. doi:10.1016/j.cma.2006.07.011.

[57] J. W. Peterson, A. D. Lindsay, F. Kong, Overview of the incompressible Navier–Stokes simulation capabilities in the MOOSE framework, Advances in Engineering Software 119 (2018) 68–92. doi: 10.1016/j.advengsoft.2018.02.004.