Localized model order reduction and domain decomposition methods for coupled heterogeneous systems

We propose a model order reduction technique to accurately approximate the behavior of multi‐component systems without any a‐priori knowledge of the coupled model. In the offline phase, we construct independent surrogate models by solving the local problems with parametrized interface boundary conditions, while we combine them using domain decomposition techniques during the online phase. We show the potential of the proposed approach in terms of accuracy, computational performance and robustness in a series of test cases, including nonlinear and multi‐physics problems.


INTRODUCTION
Mutual interactions in coupled, multi-component, parametrized systems can cause emergent behavior. 1 This defines phenomena that cannot be observed while looking at the individual components separately, as they are a result of the interactions of the constituents. Examples include synchronization, resonance, hydrodynamic instabilities and bifurcations arising from variations in the coupling terms or physical parameters. 2 Therefore, it is of interest to develop numerical methods able to capture such phenomena. Given the component-wise structure of the system, domain decomposition methods have proven to be effective in this direction. 3 The problem is decomposed into smaller domains that are solved independently through iterated exchanges of boundary data. The main advantage is that there is no need for a global solver for the specific coupled problem, which is useful in, for example, multi-physics frameworks or general problems coupled through interfaces.
Practical applications, including control problems, often require accurate and computationally fast multi-query simulations. This can be achieved through model order reduction techniques, among which a widely known paradigm is the reduced basis method. 4,5 One of its main requirements is the availability of solutions of the full system of interest at selected parameter values. However, this assumption might not always hold. On one hand, the coupled problem might not be a-priori available. Indeed, at the stage of the generation of the surrogate model, one may only have access to local models, without knowing how they will be assembled in future stages. On the other hand, even if the coupled model is available, its simulation can be computationally expensive. This is particularly true when the complexity of the problem is high, the number of components is large, or when the parameter space is high-dimensional and many samples are required to construct an accurate reduced basis. Additional complexities arise if the components are modeled using different simulation tools, as they can be difficult to integrate. 6 We propose a model reduction technique that relies on local surrogates, in the form of reduced models constructed in a component-wise, decoupled fashion. In the offline phase, the components are considered separately, and the snapshots are collected through a parametrization of the interface data. These are subsequently used to construct a local reduced basis. During the online phase, the reduced models are combined using a domain decomposition approach.
Localized model reduction of multi-component systems has been investigated in a number of works. 7 A related approach is the Reduced basis, Domain decomposition, Finite elements (RDF) method, 8 which constructs local surrogates using a parametrization of the interface data with Lagrange or Fourier basis functions and a greedy algorithm. This method shares many features with our approach, but the use of the full model in a few regions of the computational domain and the application to a single steady-state diffusion problem limit its generalization potential. Built upon this, a more recent work 9 proposes an algebraically non-overlapping decomposition strategy based on a Petrov-Galerkin projection and the enforcement of suitable compatibility constraints. Here, the authors solve the global problem using a Least-Square Petrov-Galerkin formulation equipped with different basis types and coupling constraints. Unlike the RDF method, it can easily handle nonlinear problems and constructs a fully reduced model. Our method goes in the same direction, as we aim to tackle complex multi-physics problems, but solves the coupled problem with iterative methods and is fully local both in the offline and the online phase.
Similar methods rely on Lagrange multipliers [10][11][12][13] to impose the coupling conditions. Although they can be easily constructed and efficient solvers can be designed, the introduction of additional variables can affect the computational performance of the method and might require modifications of the local solvers to accommodate this technique. In this regard, we mention a recent application in hemodynamics. 14 The idea is to split arteries into a small number of building blocks, construct local basis functions, and couple them through spectral Lagrange multipliers. The method handles time-dependent geometrically parametrized problems, and is well suited for biological applications, as blood vessels are naturally constructed as a combination of these building blocks. However, this might not hold for general cases, in which the components have significantly different structures. Moreover, although the bases are entirely local, the training is done by solving the global problem on a number of artificial configurations, which may not be feasible in more complex cases.
Similarly, static condensation methods [15][16][17] aim to solve the problem through an interface equation obtained by a splitting of internal and boundary degrees of freedom. Here, separate reduced bases are constructed to capture the internal dynamics and the interface behavior, respectively. Although they are not well suited for nonlinear problems, they have been proven extremely effective in large-scale structural dynamics problems. In this field, the design of local reduced models is an active area of research. 18,19 The main idea is to construct a local basis that captures both the static and the dynamic behavior, by studying both the system response in the presence of unitary boundary displacements and internal eigenvalue problems. Additionally, interface reduction methods 20 aim to further reduce the dimensionality at the component interfaces by using basis functions related to a secondary interface eigenvalue analysis. In this spirit, our approach relies on similar ideas to construct the parametrization of the interface data. However, these methods are based on model-driven approaches, and explicitly split internal and boundary degrees of freedom. Extending them to complex multi-physics problems, possibly characterized by nonlinearities and a nontrivial parameter dependence, seems rather prohibitive. Our method can easily handle this, and constructs a purely data-driven basis, which extracts the most relevant features of the dynamics in a problem-independent fashion, without an a-priori separation between its static and dynamic behavior or internal and boundary degrees of freedom.
Alternative approaches rely on randomization and oversampling. 21,22 The local problems are solved in a domain slightly larger than its physical one with random boundary conditions, and their solutions, once restricted to the physical domain, are used to construct local bases. Exponential convergence in the reduced dimension can be proven for linear problems, but its extension to complex nonlinear problems can be challenging. Moreover, the training procedure requires a modification of the computational domain, which might not always be possible.
A final approach, based on domain truncation, constructs a local reduced model for a combustion problem. 23 This method is once again related to the one presented here, and notably analyses a complex time-dependent problem. However, the artificial parametrization is based on a theoretical knowledge of the problem of interest, which prevents the method to generalize other classes of multi-physics problems. Moreover, the final problem is not fully reduced, as the surrogate model is constructed for a single subproblem only.
We stress that, unlike most of past work, our approach requires a minimal set of modifications in both the offline and the online phase, is targeted to time-dependent heterogeneous problems, relies on partitioned schemes and supports mixed interface conditions. The remainder of this paper is structured as follows. In Section 2, we present the coupled system of interest and an overview of domain decomposition methods. After a general discussion on model reduction, Section 3 describes our method, highlighting its steps in the offline and the online phase, an a-priori error estimate, and a discussion on the computational cost. Section 4 presents several numerical results, showing the potential of the proposed approach in terms of accuracy, computational efficiency and robustness. A few concluding remarks are found in Section 5.

PROBLEM FORMULATION
In the following, we introduce the mathematical formulation of the problems of interest. We consider multi-component systems of partial differential equations defined on different non-overlapping subdomains, so that the local problems correspond to the components of the system. To simplify the notation, we focus on the case of two time-dependent problems, denoted by  1 and  2 , defined on as many subdomains. We consider a time interval [0, T] and a spatial domain where Ω i is the domain on which  i is defined. The interface between the two, where the coupling conditions are enforced, is denoted by Γ ∶= Ω 1 ∩ Ω 2 . Thus, the coupled system can be written as where  i is a generic second order differential operator modeling the dynamics of problem  i . The initial conditions in each subdomain are denoted by u i,0 . The boundary conditions on the physical domain boundaries are encoded by the operators h i (u i ), whereas the interactions between the components are modeled through suitable functions f i and g i . Typically, the coupling functions are designed to enforce constraints on physical quantities among the components, such as continuity of the solutions and the normal fluxes. Instead of enforcing them separately as in more classical formulations, one can consider a linear combination of these constraints. 3 Among other advantages, this addresses concerns of well-posedness of the local problems and allows to design faster solution algorithms. 24 Denoting the fluxes and the unit normal vectors as i = i (u i ) and n i , respectively, we consider coupling conditions of the form where i ∈ R + are tunable parameters.
Remark 1. More sophisticated coupling conditions can be considered, at least when simulating the coupled problem. Examples include adding tangential derivatives 24 or time derivatives 25 in (2). However, (2) is a sufficiently general functional form to model the interactions in coupled systems and it is commonly used in practical applications. Therefore, in this work we employ this functional representation in both the offline and the online phase, as shown in Section 3. We conjecture that the reduced models are robust enough such that other boundary conditions can be employed in the online phase without re-running the offline phase entirely.
Remark 2. Higher order problems require more general coupling conditions, so that a careful definition of (2) is needed. An interesting example is the biharmonic operator, for which several options can be considered. 26 However, we believe that the method proposed in Section 3 can be extended to such cases while keeping the same spirit.
As our goal is to efficiently solve (1) promoting component independence, we simulate the problem using a non-overlapping Schwarz method. This allows us to decouple the local problems, for which surrogate models can be constructed. We consider the implicit Euler scheme with time step Δt and N steps = T∕Δt steps, and we denote with u n i the numerical solution of  i at time nΔt. The idea is to use fixed-point iterations, solving the local problems  i with educated guesses on the solutions of  j , j ≠ i, which are in turn updated until a suitable convergence criterion is met. Specifically, at each iteration of the temporal loop, we proceed as follows: 1. Take a guess on u n+1 i , for example, u n+1,0 i = u n i and set s = 0.

Solve
with boundary conditions h 1 ( u n+1,s+1 1 ) = 0 and 3. Solve with boundary conditions h 2 ( u n+1,s+1 2 ) = 0 and 4. Compute error, for example, 5. Repeat Steps 2, 3, 4 until convergence, for example, s+1 < tol. 6. Set u n+1 Convergence of the fixed point algorithm has been proven in a number of cases, including steady 3,27,28 and unsteady 25,29,30 problems. We remark that, to ensure and accelerate convergence, a relaxation step might be needed. 3 As it is found not to be necessary for our problems of interest, and more generally when Robin-type conditions (2) are used, we decided not to explicitly include it in the algorithm. We also note that in Step 3 the interface condition is based on the updated solution of  1 , so that the proposed algorithm can be interpreted as a block Gauss-Seidel method. Alternatively, one could consider its Jacobi counterpart, that relies on u n+1,s 1 instead. A Jacobi implementation is easier to implement and is embarrassingly parallel, but the convergence is slower. As our current implementation is purely serial, we prefer the Gauss-Seidel strategy. A detailed discussion on the two approaches can be found in other works. 6 Extending the presented formulation to more than two subdomains is rather straightforward. Denoting by  i the index set of the problems  j that share an interface Γ ij with  i , the coupling conditions of each  i can be written as for each j ∈  i . The non-overlapping Schwarz method can be generalized in a similar way. A Jacobi implementation can be directly inferred from the two-subdomain case, whereas in the Gauss-Seidel case one needs to take particular care of the ordering in which the elements are processed. This problem is strictly related to the generation of a directed acyclic graph and graph coloring. 6

MODEL ORDER REDUCTION
In this Section, we firstly outline the main features of the reduced basis method, that is, the technique used to construct the local surrogate models. 4,5 Secondly, we show the application to the problem of interest.

General properties
Consider a parametrized dynamical system of the form where x ∈ R N is the state vector, ∈  ⊂ R P is a vector containing all parameters, and f ∈ R N is the function characterizing the system dynamics.
The reduced basis method relies on the assumption that the state of the system can be well approximated by a linear combination of appropriately chosen basis functions {u i } k i=1 for any parameter in . The number k of basis functions identifies the dimension of the reduced system, where ideally k ≪ N. Thus, we have where ∈ R k denotes the vector of the latent coordinates { i } k i=1 , and U ∈ R N×k is a matrix that contains the basis vectors as columns.
To construct an evolution equation for , we replace x in (9) with its approximationx to obtain where r denotes the residual due to the approximation (10). As the system (11) is overdetermined, a widely used closure technique is the Galerkin projection. Requiring the residual vector to be orthogonal to U and assuming U to be orthonormal, we obtain Other choices, on the line of Petrov-Galerkin projections, are possible, 31 but they are found not to be necessary to achieve an efficient reduction for the types of problems we consider. To construct the basis {u i } k i=1 different approaches can be considered, among which one of the most widely used is the Proper Orthogonal Decomposition (POD). Assuming that the solution x is available at appropriately chosen times and parameter instances { j } q j=1 , we define the snapshot matrix as the matrix having such solutions as columns, that is, The POD basis U minimizes the projection error of the snapshots onto the reduced space, that is, it is the solution to the optimization problem arg min where || ⋅ || F is the Frobenius norm and I k is the identity matrix of dimension k. The Eckart-Young theorem states that the solution of (14) is given by the first k left singular vectors of S. In particular, if is the singular value decomposition of S, then An offline-online splitting clearly emerges. The construction of the basis functions, although computationally expensive, is done in the offline phase. Conversely, the online phase is mainly constituted by the time integration of (12). Under parameter variations, only the new latent coordinates need to be found, whereas the basis remains unchanged. Therefore, a crucial aspect of the reduced basis method is the efficient integration of (12). Assume that the function f can be decomposed as where L is a parameter-independent, constant-in-time, linear operator, and f nonlin is a general nonlinear operator. Using (17), the right-hand-side of (12) reads By linearity, the operator U T LU ∈ R k×k can be precomputed in the offline phase, making the online computations independent of the full dimension N. Conversely, the evaluation of the nonlinear part requires the reconstruction of the solution, the function evaluation, and the projection onto the low-dimensional space, leading to a cost even larger than that of the full system (9). To overcome this bottleneck, several techniques have been developed. [32][33][34] Unless stated otherwise, in this work we rely on the Discrete Empirical Interpolation Method (DEIM), 35 whose main idea is to evaluate only a selected number of components of the nonlinearity. We seek to approximate the nonlinear function as a linear combination of suitable basis vectors {v i }̃k i=1 as where the vector c ∈ R̃k contains the coordinates of the approximation in the basis V. To compute this vector, we select k rows from the overdetermined system (19). Consider the matrix where Assuming that the matrix P T V is invertible, the vector c is found by solving the linear system The basis V is constructed by solving a problem similar to (14), replacing the solution snapshots with the evaluation of the nonlinear function f nonlin for each sample included in S. 36 The indices {p i }̃k i=1 correspond to the components of the largest magnitude of an appropriately defined residual vector. 35 After replacing f nonlin withf in (18), we obtain If f nonlin depends sparsely on x, the evaluation of P T f nonlin requires the reconstruction of a number (k) of components, in contrast to the N required by (18). As the operator U T V(P T V) −1 ∈ R k×k can be precomputed during the offline phase, the online computational cost becomes independent of N once again. Combining (12) and (22), the POD-DEIM approximation of (9) reads A generalization of (17) and (23) can be easily derived when the operator L is affine in the parameters or in the presence of vector-valued differential equations. In this second case, a separate basis is constructed for each variable (e.g., velocity and pressure in fluid problems), leading to a block-diagonal matrix U. We now adapt this general technique to our problem of interest (1), recalling that our main assumption is that the fully coupled model is not available when constructing the reduced model. In the following, we refer to (9) as the Full Order Model (FOM), while (12) and (23) constitute the Reduced Order Models (ROM).

Offline phase
The main objective of the offline phase is the generation of the basis functions. As the coupled problem (1) is not available, we consider the components separately and we build a basis for each of them. Specifically, each problem  i in (1) has the form where the subscripts i have been omitted for simplicity and g is a generic function of the form In (25), j are functions used to expand the boundary datum, N exp is their number, and j are the expansion coefficients. The idea is to collect snapshots of (24) from simulations with different parameter values. We have two classes of parameters: • Physical parameters and time, analogous to standard model reduction approaches. Samples can be collected using, for example, greedy algorithms 5 or quasi-random sampling. • The artificial parameters defining the boundary condition (24d) and (25). Samples are collected using the following method: -The Robin parameter is fixed a-priori. The only condition on its value is the well posedness of (24), which is typically ensured if > 0. -The functions j are also fixed. Inspired by similar works, 37 we choose them as the smallest eigenfunctions of the Laplace-Beltrami operator on the interface Γ, that is, with suitable boundary conditions. In order not to be restricted to a specific set of interface profiles, to increase the level of generality of the basis, and to have a larger insight on the local dynamics, we consider both homogeneous Dirichlet and homogeneous Neumann conditions. Moreover, we additionally found that including the zero function in (25) is beneficial for the reduced model, as it contributes to capture the dynamics of the completely decoupled problem. Thus, we select N exp = 2N basis = N Dir + N Neu + 1, where N basis = N Dir + 1 = N Neu is the number of the smallest eigenvalues to be kept. If the interface is one-dimensional with length L, this reduces to where is a curvilinear coordinate defined on the interface. -The coefficients used to generate the snapshots can vary. Assuming that we simulate the problem (24) N sampl This means that (24) is simulated with each basis function as a boundary condition, including the zero function. The choice of the scaling factor A is arbitrary. We only require A to scale as 1 (resp. ) in the limits of small (resp. large) , corresponding to Dirichlet and Neumann conditions. Thus, we consider A = max{1, }.
The main steps of the offline phase follow, that is, for each instance of the physical parameters, we simulate (24) using the parametrized Robin datum, and collect the solution snapshots. These are summarized as in Algorithm 1. The extension to vector-valued differential equations is rather straightforward, and can be done by parametrizing the boundary condition of each variable independently. Similarly, the extension to coupled problems with more than two subdomains can be easily derived. As  i can have multiple interfaces Γ ij and neighboring subdomains j ∈  i , we repeat the proposed procedure for each Γ ij , assigning homogeneous coupling conditions on the remaining Γ ij ,j ≠ j.
Remark 3. We choose to treat the physical and the artificial parameters in an independent, separate way. Although this can lead to a large offline computational cost, this is not prohibitive, especially in the case of a parallel implementation. Other choices can be considered, with the goal of jointly sampling their values and reducing the associated offline cost. An example can be found in Reference 8, where a greedy algorithm is used.
Remark 4. We stress that the choice (28) implies that the local problems are simulated using the eigenfunctions j individually with constant-in-time coefficients. Different choices can be considered, with the goal of improving the sampling of the boundary data (25). An option could be a random sampling of the coefficients ij , possibly combined with a suitable modal decay rate. More challenging would be the introduction of a time-dependence, although theoretically possible. In the absence of information on the coupled problem, an appropriate selection of the time interval and/or a basis expansion in the time domain can be difficult. Both these options would increase the offline cost, introduce additional hyperparameters, and impact on the reduction potential, and are found not to be necessary to achieve a satisfactory level of accuracy for our problems of interest.

ALGORITHM 1. Offline phase
Require: Local problems Ensure: Local bases U i and local reduced operators for all system components i do for all instance of the physical parameter vector (greedy, random sampling, … ) do for all boundary data (25) and (28) do Simulate the problem and collect the snapshots end for end for Construct basis for the given component Compute reduced operators end for Remark 5. Effectively, the proposed method introduces a number of additional hyperparameters. The first is the number of interface basis functions N basis . In practical applications, this can be decided a-priori by fixing the number of frequencies that should be retained in the representation of the boundary data. The second is the Robin parameter . If available, its value is determined by the coupled problem. Otherwise, at least in the offline phase it can be chosen arbitrarily, as long as the local problem (24) is well-posed. As shown in Section 4, we numerically observed that its choice does not strongly affect our conclusions, as long as well posedness of the local problems and convergence of the fixed point iterations are guaranteed. Moreover, if these conditions are ensured, there is nothing that prevents one to choose → 0 and → ∞. Other hyperparameters include the number of samples of the boundary condition N sampl (see Remark 4) and the scaling factor A in (28).

Online phase
Once the local reduced models have been constructed, we combine them to approximate the solution of a given problem of interest. Mimicking the algorithm presented in Section 2, the general formulation of the domain-decomposition-based reduced order model reads as follows: 1. Take a guess on n+1 with boundary conditions h 1 ( 3. Solve 4. Compute error, for example, 5. Repeat Steps 2, 3, 4 until convergence, for example, s+1 < tol.
A few notes are in order. Firstly, enforcing the interface boundary conditions requires reconstructing a functional of the solution at the boundary nodes only. In the case of nonconforming meshes, this is followed by a projection/interpolation step. The reduced problem is then solved by a projection of the interface condition onto the reduced space. If the relevant operators are linear and can be accessed in the offline phase, their reduced counterparts can be precomputed, and the online phase does not require the aforementioned reconstruction step. Secondly, the norms required in Step 4 can be computed in practice without reconstructing the solution. The online phase is summarized in Algorithm 2. Note that, if the nonlinearities are treated appropriately, the computational cost of online phase becomes in principle independent of the full dimension, consistent with the discussion in Section 3.

ALGORITHM 2. Online phase
Require: Local reduced problems Ensure: Approximated solutions u i for all instance of the physical parameter vector do Simulate the coupled problem to find i Reconstruct the solution u i when required end for

Error analysis
We provide a simple but effective a-priori error analysis, identifying the different contributions to the total numerical error. It follows and generalizes the estimate derived for the RDF method. 8 Focusing again on the two-subdomain case, we denote with u = u 1 1 Ω 1 + u 2 1 Ω 2 the exact solution and with u ROM = u 1,ROM 1 Ω 1 + u 2,ROM 1 Ω 2 the reduced solution, that is, the solution of the reduced coupled problem. A similar notation is used for the other variables. Here, 1 Ω i denotes the characteristic function of Ω i . Focusing on the space-time L 2 -norm for simplicity, by the triangular inequality, we have where u h is the discrete solution and û is the discrete solution whose interface boundary conditions can be written as a linear combination of the selected Laplace-Beltrami eigenfunctions, that is, are of the form (25). The first term in (34) behaves as where h is a measure of the mesh size of the coupled problem. The convergence rate p can be derived using standard error estimates 4 and its value depends on the regularity of the problem, the selected discretization method, and the accuracy in computing the Robin data when enforcing the interface conditions. The second term in (34) is related to the approximation of a generic function using selected Laplace-Beltrami eigenfunctions. Splitting the error on the different subdomains and assuming that standard stability and trace estimates hold for each subproblem, 8 we have whereĝ i is the projection of g i using N exp functions as in (25). This projection error typically behaves 38,39 as The rate r depends on the regularity of the problem, its spatial dimension and the compatibility of the selected Laplace-Beltrami eigenfunctions with the boundary conditions at the physical boundaries. An exponential decay can appear, for example, when the problem is sufficiently regular and a boundary compatibility of the high-order derivatives is satisfied. 38,39 The third term in (34) is the reduction error and highly depends on the complexity of the local problems. This can be quantified with the Kolmogorov k-width, which represents the largest error that can be obtained by approximating the solution manifold with a k-dimensional linear space. In many cases, including diffusion equations, 40 it behaves as where k is the reduced dimension as in Section 3, at least assuming that the parameter space is appropriately sampled. For transport-dominated problems, polynomial rates appear instead. 41 Hence, the accuracy is mainly affected by the mesh size h, the number of boundary basis functions N basis and the reduced dimension k. The constants C 2, … ,6 can additionally depend on the hyperparameters defining the boundary conditions, including , A, N samples . However, as we assume these to be fixed, we have not explicitly included them in the analysis.

Computational cost
We conclude this Section by providing an estimate of the computational cost of our method, considering again Robin boundary conditions and two subdomains only. For simplicity, we focus on the linear case, so that no hyper-reduction is required. We assume that  i has size N i and that the meshes are conforming at the interface, which in turn has N Γ degrees of freedom. We denote by  the cost of each step of the algorithm.
Simulations of the coupled full order problem require, for each parameter sample, time step, and iteration of the domain decomposition algorithm: • Computing the Robin boundary condition for  1 , that is, assembling the right-hand-side of (4).  = (N Γ ).
• Solving  1 with a given Robin data, that is, computing the solution of (3) with the boundary condition computed above.  = (N 1 1 ). • Computing the Robin boundary condition for  2 , that is, assembling the right-hand-side of (6).  = (N Γ ).
• Solving  2 with a given Robin data, that is, computing the solution of (5) with with the boundary condition computed above.  = (N 2 2 ). • Computing the error (7).
The values of i and i ≤ i depend on the local underlying problems and discretization algorithms. At the reduced level, both the offline and the online phases contribute to the computational cost, although the former has to be done only once. The offline phase requires, for each system component, parameter sample, instance of the interface boundary datum, and time step: • Solving  i with a given Robin data, that is, computing the solution of (3) or (5) with (24d) as interface condition for a given value of the coefficients j in (25).
After all the snapshots of each component have been collected, the reduced bases have to be computed, as well as the linear reduced operators. This has a cost  = (N 2 i k i + N i k 2 i ). We assumed here that the Laplace eigenfunctions can be computed analytically. If that is not the case, an additional computational cost of  = (N 2 basis N Γ ) should be added for each system component. The online phase requires, for each parameter sample, time step, and iteration of the domain decomposition algorithm: • Computing the Robin boundary condition for  1 , that is, assembling the right-hand-side of (30), and project on basis . Similar to Section 3.3, if the boundary operators can be precomputed, the reconstruction step is not required and the computational cost is  = (k 1 k 2 ).
• Solving  1 with a given projected Robin data, that is, computing the solution of (29) with the projected boundary condition computed above.  = (k 1 1 ). • Computing the Robin boundary condition for  2 , that is, assembling the right-hand-side of (32), and project on basis U 2 .  = (N Γ + k 1 N Γ + k 2 N Γ ). As above, if the boundary operators can be precomputed, the reconstruction step is not required, and the computational cost is  = (k 1 k 2 ).
As in the full model, the values of i and i ≤ i depend on the local underlying problems and discretization algorithms.
To summarize, assuming that the boundary operators can be precomputed, the computational cost is: Here, the number of online parameter instances and time steps are denoted by N online and N steps,online , respectively, whereas the average number of domain decomposition iterations for the reduced order model is denoted by N iter,ROM,avg . The meaning of the remaining quantities can be directly inferred, except for the constant C, which acts as a normalization.
Indeed, the offline cost should be either amortized with the number of online evaluations of the reduced model, that is, C = N online N steps,online N iter,ROM,avg or simply ignored, that is, C → ∞. In this case, which is the one that is often considered in practice, the resulting speedup is: Although it is influenced by both the dimensionality reduction factor and the ratio in the number of domain decomposition iterations, we numerically observed that the former typically dominates. In the case of equal subdomains (N 1 = N 2 = N, 1 = 2 = ) and reduced dimensions (k 1 = k 2 = k, 1 = 2 = ), the expected speedup scales approximately as (N + N Γ )∕(k + k 2 ) ≈ (N )∕(k ).

NUMERICAL RESULTS
We now present a number of applications. For each coupled problem, we show that the reduced model efficiently recovers the full solution, supporting this by analyzing relevant quantities of interest, including the reduction error, the computational cost and the number of domain decomposition iterations. For each parameter value and time step, the errors between the full and the reduced solution are measured using the relative L 2 -norm in the appropriate spatial domain. This reads if a global error is considered, and if the local error of problem  i is considered instead. When necessary, (42) and (43) are averaged over time and and a suitable parameter test set. Unless stated otherwise, the tolerance for the Schwarz iterations is tol = 10 −6 .

Advection-diffusion-reaction
Our first system of interest is a coupled time-dependent advection-diffusion-reaction problem, in which the coefficients are allowed to vary between the components. The main goal of this test case is to validate our technique in a relatively simple framework, yet allowing the constituents to be characterized by different physical phenomena. We consider Ω = The governing equation for  i is The coefficients i , i , i control diffusion, advection and reaction, and are assumed to be constant in space and time. Homogeneous Neumann conditions are considered at the physical domain boundaries, whereas the initial condition is set as u 0 (x, y) = cos(2 x) cos(2 y) in the entire domain. We consider Robin coupling conditions at the interface, defined by for suitable parameters i . Firstly, we assume all relevant parameters to be fixed. We consider 1 = 1∕3, 1 = [3, 2] T , 1 = 1, 2 = 1, 2 = [0, 0] T , 2 = 0, while the Robin coefficients are set to 1 = 10, 2 = 10 in the offline and the online phase. The snapshots are generated from numerical simulations of the decoupled problems with N basis = 15 and a final time T = 0.05. The discretization is done using a second-order finite difference method in space with grid spacing Δx = Δy = 1∕40, and implicit Euler in time, with a time step obtained using a CFL condition with C = 1. The reduced dimension is set to k = 80 for each problem. The singular values of the local problems are reported in Figure 1, whereas the first four corresponding POD modes are shown in Figure 2. The exponential decay of the singular values indicates that both problems can be efficiently reduced, at least at a decoupled level. The most relevant modes capture the main features of the dynamics.
The first mode appears to be highly influenced by the initial condition, whereas the successive ones capture additional low-frequency characteristics. The solution of the coupled problem at the final time is shown in Figure 3 for both the full and the reduced model, together with the reduction error. This shows the high accuracy of the reduced model, and it qualitatively validates our method. Although the bases are generated independently of the coupling, the surrogates are able to capture the features of the dynamics induced by the interaction terms. Additionally, in Figure 4 we report relevant quantities of interest as a function of time. We observe that, for the selected parameter values, the relative error between the full and the reduced solutions is (10 −4 ). The number of iterations of the domain decomposition algorithm does not vary between the full and the reduced model, consistent with the high accuracy observed in Figure 3. This suggests that the computational speedup is uniquely determined by the dimensionality reduction factor, and it is in this case around 8. Although this is not extremely large, at this point we didn't aim to optimize the hyperparameters, as we simply wanted to validate the method. A lower reduced dimension could be used as well, without significantly compromising the accuracy. However, we should remark that the construction proposed in Section 3.2 is not specific to a given coupled problem, and the reduced basis needs to be robust with respect to different coupling functions. This implies that the speedup will be F I G U R E 1 Singular values of the advection-diffusion-reaction problem.

F I G U R E 3 Numerical solutions and reduction error of the advection-diffusion-reaction problem.
lower compared to more standard reduced basis approaches, 5 as we need to select a larger value of k to achieve a given error value. We mention that, consistent with a theoretical analysis, 30 the number of iterations depends on Δt in both the full and the reduced model. As our focus is on model reduction, we used the same timestep in the two cases. However, there is nothing specific about that, and it can be changed if required by a specific application. Now, we study the effects of varying hyperparameters on the reduction. Figure 5 reports the relative error, the computational cost and the number of iterations, averaged over time, for different values of the reduced dimension and boundary basis functions used to generate the reduced basis in the offline phase. Looking at both Figure 5A,D, we observe that an increase in the reduced dimension leads to a smaller reduction error, as expected. For a sufficiently large k, the reduction error saturates to a value related to the projection error of the Robin data onto the space spanned by the boundary basis functions. In turn, increasing N basis , we observe a trade-off. On one hand, as more snapshots are available, the dynamics of the full system can be better captured, and the accuracy of the reduced model improves. On the other, the heterogeneity in the snapshot matrix increases, leading to a larger variability in the reduced basis and an increase in the number of modes needed to accurately represent the solution. For a fixed number of latent variables, the former effect dominates when N basis is small, and adding more snapshots improves the reduced model. When too many basis functions are added, the reduced model starts to be affected by the increasing variability, resulting in a lower level of accuracy. For these reasons, in Figure 5A we can identify two regions. For large N basis and small k, the reduced basis error dominates, similar to more standard approaches. Instead, for small N basis and large k the boundary projection error is the dominant component, as the basis is sufficiently rich to capture the system dynamics but is unaware of the high-frequency components induced by higher boundary harmonics. The behavior of the computational cost ( Figure 5B,E) and number of iterations ( Figure 5C,F) can be interpreted in a more straightforward way. The number of boundary basis functions does not impact the computational cost per number of iterations, as N basis does not play a role in the online phase. As the number of domain decomposition iterations is not significantly altered when varying the hyperparameters, the total computational cost is only affected by k, up to minor variations.
Although the offline phase is based on a parametrization of Robin data with a fixed , the reduced basis is able to capture the dynamics of the system with other coupling functions, including, for example, Robin data with different coefficients. Figure 6 shows a comparison of different quantities of interest under variations in the Robin parameter used in the online phase. The reduction error is rather robust with respect to variations of , although minor differences are present. These are unavoidable, but the features of the dynamics are well captured even with a single offline instance of . Notably, the error decay as a function of k is not affected by changes in . The number of domain decomposition iterations, and consequently the computational cost, exhibit more significant changes. The former is related to the optimization of the interface conditions in domain decomposition methods, in which the value of is shown to impact the convergence rate. 25 The latter varies accordingly, as it is proportional to the number of iterations.
Additionally, in Figure 7 we report the trend observed when the Robin parameter used in the training phase is varied, while keeping the used in the testing phase constant. The reduction error is slightly affected, with lower accuracy for  To conclude, we mention that the reduced model can be used to speedup optimization problems. For instance, one could be interested in tuning the Robin parameters in order to minimize the number of domain decomposition iterations. Although constructing an appropriate control problem goes beyond the scope of this work, we argue that the surrogate model could be used in such a context. Figure 8 shows the iterations of the domain decomposition method under variations in the Robin parameters 1 and 2 for the two problems. The differences between the full and reduced model are mild. In turn, the associated computational cost, although not reported here, exhibits significant changes, similar to

Checkerboard diffusion
We now focus on another problem, a steady parametrized diffusion equation. We seek to show that the method can be efficiently applied to parametrized problems, proving its effectiveness also in high-dimensional parameter spaces. The computational domain is Ω = [0, 1] 2 , divided into two subdomains by a vertical interface Γ = We take the diffusion coefficients to be piecewise constant in the domain in a checkerboard fashion, that is, we assume We consider parametrized non-homogeneous Dirichlet boundary conditions at the left boundary of Ω 1 and at the right boundary of Ω 2 , that is, whereas homogeneous Dirichlet conditions are assumed on all the other physical boundaries. As in the previous example, at the interface we impose Robin boundary conditions of the form The parameters of interest are i 1 ∈ [0.1, 10] for i = 1, … , 6, i 2 ∈ [0.1, 10] for i = 1, 2, g 1 ∈ [1, 5], g 2 ∈ [1,5], so that the global parameter space has dimension 10. This splits into two local spaces of dimensions 7 and 3, respectively.
The reduced models are generated by sampling 50 parameter values using a Latin Hypercube Sampling in each subdomain and N basis = 15. The discretization is done using P 1 finite elements with grid spacing h ≈ √ 2∕60. The reduced dimension is set to k = 40 for each subdomain. The singular values of the local problems are reported in Figure 9, whereas the solution of the coupled problem and the corresponding reduction error for a single parameter instance not included in the training set are reported in Figure 10. As in the previous test case, the reduced models accurately approximate the high-fidelity solution.
We now study the reduction error, the computational cost and the number of iterations as a function of the reduced dimension. They are reported in Figure 11A-C for 20 and 11D-F 200 offline parameter samples per subdomain, respectively. In both cases, they have been computed by averaging the relative error (42), the cost and the number of iterations over 100 parameter values generated again using a Latin Hypercube Sampling, but different from the training samples. Increasing the reduced dimension improves the accuracy, until a saturation level, related to the projection of the boundary datum onto the set spanned by the Laplace-Beltrami eigenfunctions, is reached. In turn, increasing N basis improves the accuracy by reducing the saturation value, without impacting on the online computational cost.
Furthermore, Figure 11 allows us to compare the performance of the reduced bases generated using either the decoupled problems or the coupled problem. The second basis is generated from localized snapshots obtained from simulations of the coupled problem, keeping the other hyperparameters unchanged. In this way, the basis functions remain localized  to each component, but they are aware of the dynamics of the coupled problem. We observe that if the number of parameter samples is sufficiently large or a large reduction error is sufficient for the specific application, the basis generated with the coupled problem is more effective. This is tailored for that specific system, so it is expected to perform better than a general, system-independent, model. However, when the online parameters are chosen in such a way that the dynamics is very different from the one of the snapshots, using the decoupled systems allows one to achieve a low reduction error, thanks to the generality of the model. This is particularly evident when the number of offline parameter samples is small, as in Figure 11A. The sampling error in the coupled model is large, limiting the accuracy of the reduced model. This effect is mitigated by considering the systems independently, as the local parameter spaces can be sampled more effectively, F I G U R E 11 Comparison (number of basis functions and reduced dimension) of the parameter-averaged quantities of interest of the checkerboard diffusion problem. The terms 'Dec' and 'Cou' refer to the bases generated using the decoupled problems or the coupled problem, respectively. despite the need to parametrize the interface boundary conditions. As the dimension of the parameter space grows, this effect becomes more significant. Finally, we recall that in many applications the global coupled problem is not available when constructing the basis. Thus, although in some cases the basis generated from localized snapshots outperforms the proposed method, it might not be possible to construct it.

Multiple subdomains
To further increase the complexity of the problem, we consider a variation in the number of subdomains. The global spatial domain is Ω = [0, 1] 2 , which is split into M 2 square subdomains, obtained by dividing each direction into M intervals. Here, we focus on M = 2 and M = 3 only, as they already retain the properties that larger problems would feature. The problem we consider is a time-dependent diffusion equation, where the viscosity is assumed to be constant in each subdomain. Thus, problem  i is governed by The diffusion coefficients are chosen as shown in Figure 12, whereas f i = 1 in each subproblem. Moreover, homogeneous Neumann conditions are considered at the physical boundaries, the initial condition is u 0 (x, y) = cos(2 x) cos(2 y), and the final time is T = 0.2. At each interface Γ ij of problem  i , we impose Robin boundary conditions of the form where u i and u j denote the local and neighboring solution, respectively. We choose all the Robin parameters ij equal to 1. The discretization is done using P 1 finite elements with grid spacing h ≈ √ 2∕(20M) and a time step Δt = 0.01. This implies that the local number of degrees of freedom is independent of M, in the spirit of a weak scaling analysis. 6 F I G U R E 12 Viscosity distribution of the multiple subdomains problem. As a visual comparison of the full and reduced solution gives satisfactory results, in Figure 13 we directly report the relevant quantities of interest as a function of the reduced dimension for different N basis . The results are similar to the previous problems, to which we refer for a more detailed discussion. We note that increasing the number of subdomains and keeping the other parameter values unchanged leads to a slight increase in the error, especially for large N basis . This is expected, as the complexity of the local solutions increase due to the larger number of interfaces. An increase in the number of subdomains also leads to an increase in the number of the DD iterations. This is again consistent with our expectations, as more iterations are needed for information to propagate throughout the entire domain when it is decomposed into more components. To conclude, we mention that in our simulations we processed the elements using a two-color, checkerboard pattern. A different choice may lead to a different number of iterations and a different number of sequential step required in a parallel implementation, 6 both at the full and the reduced level.

FitzHugh-Nagumo
We now consider the two-dimensional FitzHugh-Nagumo model, as an example of a reaction-diffusion problem. Because of the nonlinear reaction terms, the solution develops Turing patterns. 42 The main goal is to extend the presented approach to both vector-valued differential equations and nonlinear problems. We consider Ω = [0, ] 2 , divided into two subdomains by a vertical interface Γ = The variables u i and v i represent the electric potential and a recovery variable, respectively, while the coefficients , , d, Γ govern different physical phenomena. They are chosen to have the same value in the two subdomains. Homogeneous Neumann conditions are considered at the physical domain boundaries, whereas the initial condition is set as where controls the amplitude. We consider Robin coupling conditions at the interface, defined by for suitable parameters i . Firstly, we focus again on a specific setting, by fixing the relevant parameters. Following a theoretical analysis, 42  The discretization is done using second-order finite differences in space with grid spacing Δx = Δy = 1∕40, and implicit Euler in time, with a time step obtained using a CFL condition with C = 40. The reduced dimension is set to k = 80 for each problem and variable. The nonlinear term is treated efficiently in the reduced model using DEIM, with a corresponding basis sizek = k. The first k singular values of the local problems are reported in Figure 14 for each variable. Again, they have a large reduction potential. Similarly, in Figure 15 we show the singular values associated to the nonlinear term and the indices selected by the DEIM as described in Section 3. The decay is again quite fast, and the effect of the nonlinearity appears to be stronger closer to the interface. This is not surprising, as in the offline phase we impose a large variability of the coupling conditions, which seems to affect the nonlinearity in a stronger way. The solution u of the coupled problem at the final time is reported in Figure 16 for both the full and the reduced model, together with the reduction error. Although the nonlinear term increases the complexity of the dynamics and the reduced model, the solution is well approximated with the reduced basis. Relevant indicators are shown in Figure 17 as a function of time. We observe that the reduced model accurately reproduces the solution (the relative error is around 1%), resulting in a speedup factor of 15. Again, the number of local iterations is practically independent of the dimension of the problem.
The effect of varying hyperparameters on the reduction is studied in Figure 18. We report the main performance indicators, averaged over time, for different values of the reduced dimension and boundary basis functions used in the offline phase. In terms of reduction error, the behavior as a function of the reduced dimension reflects our expectations, as the trend is exponentially decreasing. When the number of boundary basis functions is varied, we still experience the trade-off between higher accuracy and larger variability, at least when the number of latent variables is sufficiently large. The number of domain decomposition iterations has a larger variation compared to the previous test cases, and it is influenced by both the reduced dimension and the number of boundary basis functions. This reflects on the computational cost, which depends on the dimensionality reduction factor and the varying number of iterations.

Stokes-Darcy coupling
We now switch to a simple multi-physics problem, modeling the interaction between an incompressible flow governed by the unsteady Stokes equations and a porous medium described by the unsteady Darcy equation. 37 Although in a relatively simple context, this problem retains the main features of a multi-physics system, and is used to validate our technique with a more heterogeneous coupled system. We consider Ω = where (u, p) = ( u + u T ) − pI is the Cauchy stress tensor, = 0.04 the fluid viscosity and u and p are the fluid velocity and pressure. The model for  2 is consider together with a zero tangential velocity component for the fluid problem. We note that Neumann interface conditions are used to define the coupled problem. This is common practice for this problem, as they naturally arise in the weak form of the coupled problem and lead to a symmetric system. 43 By introducing suitable coefficients, they can be combined in a Robin form. As our main goal is to show that we can efficiently recover the solution of the coupled problem by constructing local reduced models, rather than optimizing the coupling conditions, we stick to this formulation for both the offline and online phase. Moreover, as shown for simpler problems, we can expect that a nonzero Robin parameter would not significantly alter the results.
Our numerical setting comprises a number of basis N basis = 10 and a final time T = 10. The spatial discretization is done using P 2 -P 1 , respectively P 2 , finite elements with a characteristic mesh size h ≈ √ 2∕40 for the Stokes, respectively Darcy, problem. An implicit Euler method with a timestep Δt = 0.5 is used for time integration. As this leads to a saddle point problem for the fluid component, particular care has to be given to the construction of the reduced model, as the vanilla model (23) may not be (inf-sup) stable. Specifically, integrating (9) requires solving a system of the form (59) With our choice of the finite element spaces for the approximation of u and p, (59) is well-posed. The reduced counterpart, constructed by computing the velocity and pressure bases U u and U p as described in Section 3, does not enjoy the same property. For this reason, we enrich the velocity basis U u with properly chosen (approximated) supremizer solutions. 44 The idea is to map each element of the pressure space to a suitable velocity element ensuring the inf-sup condition, and to include these in the velocity space. This leads to an enriched reduced space, making (60) stable. Consequently, we solve a problem of the form for each pressure snapshot p j , collect the solutions s(p j ), and apply the POD to obtain a number of supremizer modes that are added to the matrix U u . The matrix K in (61) can be chosen as the one arising from the discretization of the Laplace operator for u. We remark that the supremizer enrichment is not the only technique that ensures stability of the reduced system. Examples are the pressure Poisson formulation of the Navier-Stokes equations 45 and the artificial compressibility method, 46 which ensure stability at the full and the reduced level for more general velocity-pressure spaces without the need to enrich the velocity space. For this particular problem, the reduced dimension is set to k = 40 for each problem and variable, and the velocity space is enriched with 40 supremizer modes. The singular values of the local problems are reported in Figure 19 for each variable. The decay is quite fast, meaning that the intrinsic dimension of the local problems is quite small. The solution of the coupled problem at the final time is reported in Figure 20, together with the reduction error. All components are accurately reproduced at a reduced level. Relevant indicators are shown in Figure 21 as a function of time. For the selected hyperparameters, both the level of accuracy and the computational speedup reach high values. We note that, although the introduction of the supremizers negatively affects the computational cost, the speedup is still significant. The number of iterations is found to be independent of the dimension of the problem. Concerning variations in the hyperparameters of the reduced model, Figure 22 reports the relevant quantities of interest for different values of the reduced dimension and basis functions used to generate the snapshots in the offline phase. Most of the comments related to the previous test cases remain valid, even though the trend of the reduction error is less smooth as compared to Figure 5. This is possibly due to the larger heterogeneity and complexity of the problem, the addition of the supremizer modes, or a higher level of noise introduced by higher frequency components.

Fluid-structure interaction
Our final test case involves a more complex, heterogeneous multi-physics problem, selected as a prototype of general fluid-structure interaction problems. This test case is used to show the applicability of our reduction method in a real scenario, as the simplifying assumptions of the model are minimal. We consider Ω = [0, 1.2] × [0, 0.5], divided into two subdomains by a curved interface as depicted in Figure 23.
where f (u, p) = f ( u + u T ) − pI is the Cauchy stress tensor, f = 1.0 and f = 1.0 are the fluid density and viscosity, respectively, and u and p are the fluid velocity and pressure. Conversely, the governing equation for  2 is the linear elasticity equation is the elastic stress tensor. The physical parameters, s , E and represent the density of the solid, the Young's modulus and Poisson's ratio, respectively. For this test we have used 37 s = 1.0, E = 10 4 , = 0.48. Equation (63) corresponding to a Dirichlet-Neumann coupling. As in the Stokes-Darcy problem, they can be generalized to Robin conditions, but for the present problem this is not necessary. Thus, we maintain this formulation, parametrizing the Dirichlet and Neumann data respectively. We consider a number of basis N basis = 10 and a final time T = 0.5. The discretization of the fluid, respectively elastic, problem is done using P 2 -P 1 , respectively P 2 -P 2 , finite elements in space based on the mesh reported in Figure 23, and implicit Euler in time with Δt = 0.02. The reduced dimension is set to k = 80 for each problem and variable, and the velocity space is enriched with 80 supremizer modes.
The singular values of the local problems are reported in Figure 24 for each variable. The decay for the fluid problem, although exponential, is slower when compared to the previous test case. This is not surprising, as the current problem exhibits more complex features. In the elastic problem, the differences among the x-and y-component are minor, at least F I G U R E 25 Numerical solutions and reduction error (fluid) of the fluid-structure interaction problem.

F I G U R E 26
Numerical solutions and reduction error (structure) of the fluid-structure interaction problem.

F I G U R E 27
Quantities of interest of the fluid-structure interaction problem. Concerning variations in the hyperparameters of the reduced model, Figure 28 reports the relevant quantities of interest for different values of k and N basis . The results are consistent with the previous test cases, although the error decay is negatively impacted by the high degree of heterogeneity and complexity of the problem.

CONCLUSION
In this work, we propose a reduced basis method for coupled heterogeneous systems. Unlike most techniques, this approach requires no simulations of the coupled problem, as it relies on local solvers only. The modifications to the local models are minimal, making the approach straightforward and easily generalizable to other problems of interest. The offline phase introduces an artificial parametrization of the boundary conditions at the interface, which is used to collect snapshots of a given subsystem. The online phase couples the reduced models using a domain decomposition strategy inherited from the full order problem. We show that the solution of a given coupled problem can be efficiently approximated by the reduced model. The reduction error depends on both the number of boundary basis functions used in the offline phase and the reduced dimension, whose values could possibly be tuned to achieve a specified tolerance, at least for simple problems. The computational cost is mainly influenced by the dimensionality reduction factor, as the number of domain decomposition iterations is only mildly affected by the reduction. Finally, the reduced models are robust with respect to changes in the coupling functions and parameters. Possible future research directions include the extension to coupled systems at a larger scale, further increasing the complexity of the local models, the number of components, or the spatial dimension of the problem.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

ACKNOWLEDGMENT
Open access funding provided by Ecole Polytechnique Federale de Lausanne.