The role of interface boundary conditions and sampling strategies for Schwarz-based coupling of projection-based reduced order models

The notion of coupling subdomain-local PROMs with each other and/or with FOMs has been gaining attention in the PMOR literature in recent years. Existing coupling approaches can be classified as falling into two categories: (i) monolithic coupling methods, and (ii) iterative coupling methods. Methods from category (i) are typically more intrusive to implement but can possess more favorable convergence properties and runtimes than methods falling into the category (ii).

The vast majority of monolithic coupling methods are based on formulations that first define appropriate interface compatibility constraints, and then enforce these constraints using different techniques. Among the earliest approaches for coupling subdomain-local PROMs with each other and with FOMs using Lagrange multipliers is the work of Lucia et al. [53] (for POD-based PROMs) and Maday et al. [55] (for PROMs constructed via the reduced basis method). Both [53] and [55] impose compatibility using a mortar-type method that “glues” together non-overlapping subdomains using Lagrange multipliers approximated using low-order polynomials. Several subsequent works approach the coupling problem as a means to “stitch” together a variety of composable PROM “tiles.” In [81], continuity between tiles is enforced by duplicating the degrees of freedom on the interfaces and constraining their normal components to be equal. A Lagrange multiplier-free approach that is conceptually similar to [81] is the reduced basis method with domain decomposition and finite elements (RDF) [41]. In this non-overlapping approach, traditional finite element basis functions are defined on the boundaries between PROM domains so as to ensure continuity and provide a sort of finite element enrichment, with the goal of gluing together networks of repetitive blocks. Another recent work that accomplishes monolithic PROM-FOM coupling without Lagrange multipliers is [3]. Here, a discontinuous Galerkin (DG) formulation is assumed and coupling is achieved through a careful definition of numerical fluxes at discrete cell boundaries. In [37], Hoang et al. develop an algebraically non-overlapping method for coupling subdomain-local LSPG PROMs with each other. Compatibility conditions at subdomain boundaries can be imposed either strongly or weakly using Lagrange multipliers, and subdomain training can be performed independently. In a similar vein, Chung et al. [13] recently proposed their so-called “data-driven finite element method” (DD-FEM), which trains spatially-local reduced order models (ROMs) corresponding to small-scale unit components, tiles a computational domain with these components, and employs ideas from DG methods to glue the components together. In this method, the interface constraints are only enforced weakly through interface operators introduced to penalize constraint violations, thereby eliminating the need to introduce Lagrange multipliers. In [20], Diaz et al. extends the approach in [37] to nonlinear manifold autoencoder-based ROMs, which can be more efficient for problems in the slowly-decaying Kolmogorov $n$-width regime. In [12, 18], de Castro et al. develop a Lagrange multiplier-based explicit synchronous partitioned scheme for performing PROM-PROM and PROM-FOM coupling. One of the advantages of this method is that, when combined with explicit time-integration, the subdomain-local problems can be solved independently through a clever utilization of the Schur complement. Moreover, when utilizing independently-constructed PROM bases for the interfacial and interface degrees of freedom, it is shown that the approach leads to a provably nonsingular Schur complement independent of the underlying mesh size and/or reduced basis dimension [18].

The Schwarz alternating method falls into a category of iterative coupling methods. As mentioned above, these methods are usually less intrusive to implement in HPC codes, but may be less efficient than monolithic methods. Among the earliest Schwarz-based DD approaches for PROM-PROM and PROM-FOM coupling focuses on Galerkin-free PROMs. In [8], Buffoni et al. propose three methods to perform Galerkin-free PROM-FOM and PROM-PROM coupling: (i) a Schur iteration where the solution of the PROM is obtained by a projection step in the space spanned by the POD modes, (ii) a Dirichlet-Dirichlet iteration in the frame of a classical Schwarz method, and (iii) an approach obtained by minimizing the residual norm of the canonical approximation in the space spanned by the POD modes. Galerkin-free ROM-FOM and ROM-ROM couplings are also explored by Cinquegrana et al. [14] and Bergmann et al. [6]. The former approach considers overlapping DD in the context of a Schwarz-like iteration scheme, but, unlike our approach, requires matching meshes at the subdomain interfaces [14]. The latter approach, known as zonal Galerkin-free POD, defines an optimization problem which minimizes the difference between the POD reconstruction and its corresponding FOM solution in the overlapping region between a ROM and a FOM domain [6]. A true POD-greedy/Galerkin non-overlapping Schwarz method for the coupling of PROMs developed for the specific case of symmetric elliptic PDEs is presented by Maier et al. in [56]. In the recent work [42], Iollo et al. demonstrated that the classical overlapping Schwarz iteration can be interpreted as an optimization-based coupling, which can be solved directly in the method known as “one-shot Schwarz.” Optimization-based domain decomposition frameworks for coupling subdomain-local PROMs/FOMs have also been considered in several other recent works, e.g., [68] and [36].

Some emerging technologies in PMOR coupling are applicable to both monolithic and iterative coupling formulations. The approaches summarized above assume a predefined domain decomposition (DD) and PROM/FOM subdomain assignment, which is consistent with the scope of this paper. It is worth noting, however, that PROM-PROM and PROM-FOM coupling methods with on-the-fly basis and/or DD adaptation are also possible and beginning to emerge in the PMOR literature. In [16, 15], Corigliano et al. develop an adaptive non-overlapping Lagrange multiplier-based coupling method for nonlinear elasto-plastic and multi-physics problems. Here, on-the-fly PROM adaptation and PROM/FOM switching is performed through a plastic check during the reduced analysis. Hybrid PROM–FOM coupling in the context of solid mechanics applications is considered also in [43, 44, 70]. In [43, 44], a local/global model reduction strategy for the simulation of quasi-brittle fracture is developed, whereas, in [70], an adaptive sub-structuring (or DD) approach is proposed that not only enables on-the-fly adaptation of the PROM basis, but also on-the-fly substructuring/DD changes. Finally, in [40], Huang et al. develop a component-based modeling framework that can flexibly integrate PROMs and FOMs for different components or domain decompositions, with an eye towards large-scale combustion problems. It is demonstrated that accuracy can be enhanced by incorporating basis adaptation ideas from [64, 39].

Researchers have additionally begun in recent years to integrate ideas from machine learning (ML) into PROM-PROM and PROM-FOM coupling formulations. In [1], Ahmed et al. present a hybrid PROM-FOM approach in which a long short-term memory (LSTM) network is introduced at the subdomain interfaces and utilized to perform the multi-model coupling. In order to improve efficiency, Bochev et al. [7] replace the Schur complement solve required in the explicit synchronous partitioned presented in [12, 18] with a dynamic mode decomposition surrogate trained from interface flux data. The second direction involves the application to non-intrusive PROMs. In a similar vein to [7], in [21], Discacciati et al. propose a non-intrusive PMOR approach, which approximates the boundary response maps arising from a non-overlapping DD method using a combination of dimensionality reduction techniques and interpolation/regression approaches (e.g., kernel interpolation methods and ANNs). In [25], Farcas et al. develop a methodology for the coupling of nonintrusive and subdomain-local operator inference (OpInf) models, which can both improve accuracy and reduce memory/computation requirements. A similar approach for coupling subdomain-local OpInf models with each other and with FOMs was recently introduced in [57]. In [49, 50], the Schwarz alternating method was extended to enable the DD-based coupling of physics-informed neural networks. The methods proposed in these works, known as D3M [49] and DeepDDM [50] respectively, inherit the benefits of DD-based ROM–ROM couplings, but are developed primarily for the purpose of improving the efficiency of the neural network training process and reducing the risk of over-fitting. A similar approach, but considering both weak and strong formulations of the Schwarz BCs, was explored by Snyder et al. in [77]. The Schwarz alternating method has also been used for online coupling of independently pre-trained subdomain-localized neural network-based models, e.g. in [80], which develops a transferable framework for solving boundary value problems (BVPs) via deep neural networks that can be trained once and used forever for various unseen domains and BCs.

While our manuscript builds on some of our earlier work in developing the Schwarz alternating method for FOM-FOM [58, 59] and PROM-PROM/PROM-FOM [4] coupling in solid mechanics, it includes several novel contributions.

Whereas our past work focused on solid mechanics and finite element-based discretizations, here we consider challenging two-dimensional (2D) nonlinear hyperbolic conservation law problems, including problems with strong shocks, discretized using a cell-centered finite volume method (CCFV). While Dirichlet boundary conditions can be imposed strongly in FEM-based discretizations and PROMs based on top of these discretizations [4, 34], all boundary conditions in CCFV methods are enforced weakly using ghost cells, auxiliary computational cells that extend the physical domain for the purpose of boundary condition imposition. It is well-known that the non-overlapping Schwarz alternating method is in general non-convergent if Dirichlet-Dirichlet transmission boundary conditions are specified at the Schwarz boundaries; to achieve convergence, one must use either alternating Dirichlet-Neumann [52] or Robin-Robin [84] transmission conditions. Implementing these conditions can be cumbersome in HPC codes, as the creation of carefully-constructed transmission operators is required [60]. We demonstrate that the use of ghost cells for boundary condition imposition, which has the effect of introducing overlap into an otherwise non-overlapping domain decomposition, can improve the flexibility and simplify the implementation of the non-overlapping version of the Schwarz alternating method by enabling the use of Dirichlet-Dirichlet coupling conditions at the Schwarz boundaries. This is a novel result, to the best of our knowledge.

Our past work involving Schwarz for coupling has focused on the classical “multiplicative” variant of the method [58, 59, 4]. A well-known disadvantage of multiplicative Schwarz stems from the fact that the method solves subdomain-local problems in sequence until convergence is reached, which limits computational acceleration. This may be remedied by a variant of the method known as “additive” Schwarz, where subdomain-local problems are solved concurrently and boundary data are communicated asynchronously. Our numerical results demonstrate that, with subdomain-local PROMs coupled using additive Schwarz, it is possible to achieve up to two orders of magnitude speedup over equivalent coupled FOM solutions and moderate speedups over analogous monolithic FOM solutions. Our numerical experiments additionally reveal that a PROM-PROM coupling of a given size is typically more accurate than a monolithic PROM of comparable size defined on the full problem domain. This result suggests that the proposed coupling methodology has the potential to improve PROM accuracy by enabling the spatial localization of these models via domain decomposition.

Lastly, while our previous work on Schwarz for PROM-PROM coupling [4] considered nonlinear PDEs and thus included hyper-reduction to maintain efficiency, attention was restricted to one-dimensional problems, for which each subdomain had at most two Schwarz boundary nodes. In moving to the multi-dimensional case, where it is possible to have hundreds or thousands of Schwarz boundary mesh elements, it is important to develop a hyper-reduction strategy for sampling the Schwarz boundary in a way that maintains both accuracy and efficiency. We demonstrate that the optimal boundary sampling rate for applying Schwarz to hyperbolic conservation law PROMs with hyper-reduction is problem-specific, and, for the problems considered here, varies between sampling every two boundary cells (for problems exhibiting traveling shocks) and every ten boundary cells (for problems with smooth solutions).

The remainder of this paper is organized as follows. In Section 2, we describe the Schwarz alternating method in the context of a general hyperbolic PDE, introducing several new ingredients relative to our past related work [58, 59, 4], namely some details regarding additive Schwarz and the application of the method to CCFV discretizations. In Section 3, we overview PMOR based on the POD/LSPG method with hyper-reduction, and present the details relevant to applying the Schwarz alternating method to PROMs constructed using this approach. We evaluate our method in the context of three 2D nonlinear hyperbolic conservation laws in Section 4: the shallow water equations, Burgers’ equation, and the Euler equations. Conclusions are offered in Section 5, which also includes a brief discussion of potential future work.

Consider the PDE (1) in $\Omega\times I_{n}$ for some specified time-interval $I_{n}=[t_{n},\ t_{n+1}]$, where $T\geq t_{n+1}>t_{n}\geq 0$ with initial condition $u(x,\ t_{n})=v(x)$ for a given function $v$. Assuming without loss of generality a two subdomain decomposition like those shown in Figure 1, the additive Schwarz algorithm concurrently solves the following subdomain-local problems

and

for Schwarz iteration $k=1,2,...$ subject to the aforementioned initial conditions $u_{1}(x,0)=v|_{\Omega_{1}}$ and $u_{2}(x,0)=v|_{\Omega_{2}}$. Here, $u_{i}$ for $i=1,2$ denotes the solution in subdomain $\Omega_{i}$, $\Gamma_{i}$ is the so-called Schwarz boundary (see Figure 1), and $\mathcal{B}_{ij}$ for $i,j=1,2$ is a linear operator acting along the Schwarz boundary $\Gamma_{i}$ on the solution $u_{j}$. It is common to set $u_{1}^{(0)}=v|_{\partial\Omega_{1}}$ on $\Gamma_{1}$ and $u_{2}^{(0)}=v|_{\partial\Omega_{2}}$ on $\Gamma_{2}$ when initializing the Schwarz iteration process, to ensure solution compatibility with the initial condition. The iterative process in (2) and (3) continues until a set of pre-determined criteria are met. In the present work, convergence criteria are based on the Euclidean norm of the solution differences between consecutive Schwarz iterations; that is, Schwarz is deemed converged when $\epsilon_{\text{abs}}^{(k)}<\delta_{\text{abs}}$ and $\epsilon_{\text{rel}}^{(k)}<\delta_{\text{rel}}$ for some pre-specified tolerances $\delta_{\text{abs}},\delta_{\text{rel}}>0$, where

and

for Schwarz iteration $k=1,2,...$.

Numerous past works have investigated the definitions of the operators $\mathcal{B}_{ij}$ to ensure well-posedness and convergence of each Schwarz subproblem and the Schwarz iteration process (the reader is referred to [28] for a nice summary). Table 1 summarizes several choices of boundary operators defining a set of convergent Schwarz transmission conditions on the boundaries $\Gamma_{i}$ for the specific case of a two subdomain domain decomposition. Whereas Dirichlet-Dirichlet transmission conditions on $\Gamma_{i}$ will give rise to a convergent Schwarz iteration process provided the overlap region is non-empty [51] ($\Omega_{1}\cap\Omega_{2}\neq\emptyset$; see Figure 1, non-overlapping domain decompositions ($\Omega_{1}\cap\Omega_{2}=\emptyset$; see Figure 1 typically require alternating Dirichlet-Neumann [84, 27, 17, 46] or Robin-Robin [52, 19, 54] transmission conditions on $\Gamma$. In Table 1, projection operators $P_{ij}$ from subdomain $\Omega_{j}$ onto boundary $\Gamma_{i}$ are introduced to enable application of the method to the case where $\Omega_{1}$ and $\Omega_{2}$ are discretized using non-conformal meshes. As described in more detail in [59, 60], the projection operator from one subdomain to another consists of the simple application of the existing finite element interpolation functions. That is, any time the value of a field is needed at a point at the boundary in one subdomain, a search is performed to determine the element containing that point in the other subdomain, and the field value is calculated by using the nodal values of that element and the corresponding interpolation functions.

While the PDEs considered herein are dynamic, our discussion thus far has focused primarily on the spatial aspects of Schwarz-based coupling. For long-time dynamic problems defined on a large time-interval $[0,\ T]$, it is highly inefficient to integrate each subdomain $\Omega_{i}$ to time $T$ before synchronizing boundary conditions at the Schwarz boundary. If each subdomain is advanced forward in time using the same time-step $\Delta t$, a natural mitigation to this difficulty is to converge the Schwarz iteration process on a time-step by time-step basis, as illustrated in Algorithm 1.

Note that we have specifically presented the additive Schwarz algorithm. This is in contrast to the original multiplicative Schwarz algorithm considered in our past works [58, 59, 60, 4], which solves each subdomain-local problem in sequence, rather than concurrently, thus incorporating the most up-to-date interface information. To modify Algorithm 1, line 8 would simply be altered to interpolate the solution $u^{(k)}_{1}$ onto $\Gamma_{2}$. Although it enables concurrent subdomain calculations, the additive Schwarz algorithm will in general require more iterations to converge than would the multiplicative algorithm [28]. However, it has been shown that, provided appropriate definition of the operators $\mathcal{B}_{ij}$ (see Table 1), the additive Schwarz iteration (2) and (3) will converge to the same solution of (1) as the multiplicative Schwarz iteration. All results shown in this paper utilize the additive Schwarz algorithm to enable concurrent calculations across all subdomains.

For spatial discretizations which define degrees of freedom at the endpoints, vertices, or edges of mesh elements, the application of non-overlapping Schwarz boundary projection operators (as outlined in Table 1) is rather straightforward. Often, this simply involves interpolation on a line segment for a 2D domain, or on a plane for a three-dimensional domain. Finite difference, finite element, and vertex-centered finite volume discretization schemes are readily adapted to non-overlapping Schwarz decompositions as a result. On the other hand, CCFV schemes define degrees of freedom at the centroid of the mesh volumes. For such schemes, boundary conditions can be applied by a “ghost cell” formulation, whereby fictitious cells are placed outside the computational domain, adjacent to cells whose faces compose the domain boundaries. These ghost cells need not have any geometry associated with them, though a centroid coordinate may be used for the sake of convenience. The system state in these ghost cells is specified to weakly enforce boundary conditions through the finite volume fluxes. For example, for a cell at an impermeable slip wall, the accompanying ghost cell state would be set to have the same density, pressure, and tangential velocity of the interior cell, but the opposite normal velocity. The fluxes at the boundary face would then be solved exactly the same as for any other cell face, and implicitly enforce zero transport through the wall. An overview of ghost cell boundary conditions for a variety of other physical boundaries can be found in many classical texts on finite volume methods, such as [48].

Applying the Schwarz alternating method to a domain decomposition with a CCFV discretization is largely indistinguishable from its application to any other discretization. Figure 2 shows a simplified graphical depiction of a uniform 1D CCFV overlapping decomposition. Solid lines indicate physical domains, while dashed lines indicate ghost cells. The Schwarz projection operator simply interpolates the interior state of one domain to the ghost cells of the neighboring subdomains. Solving the boundary face fluxes in the neighboring subdomain thus accounts for the effects of the solution as it is updated through the Schwarz algorithm.

The use of ghost cells in a finite volume discretization leads to a somewhat nonstandard treatment of non-overlapping Schwarz interfaces, as demonstrated visually in Figure 2. While the physical subdomains $\Omega_{1}$ and $\Omega_{2}$ do not overlap and only share a physical cell face at a single point (the midpoint between $x_{1,4}$ and $x_{2,1}$), the ghost cells create an effective overlap region having a width of two cells. This permits the application of a Dirichlet–Dirichlet Schwarz interface condition by assigning the state of the physical cells at the boundary to the corresponding ghost cell of the adjacent subdomain. The result is a more versatile method, as it is possible to obtain a convergent sequence of Schwarz iterations when utilizing a (more flexible) non-overlapping DD with Dirichlet transmission BCs at the Schwarz boundaries (row 1 of Table 1). This greatly simplifies the interface treatment by avoiding Dirichlet–Neumann or Robin–Robin conditions as in traditional non-overlapping Schwarz schemes.

While the scheme described above is not in truth a non-overlapping scheme, as there is an effective overlap region defined by the ghost cells, the resulting Schwarz variant possesses the same flexibility of a “true” non-overlapping method by enabling the independent discretization of neighboring subdomains and application to interface problems (see Remark 2). The fact that ghost cells are non-physical constructs hints at interesting possibilities for a simplified Schwarz decomposition for complex finite volume meshes. This treatment is not dissimilar from halo exchange regions utilized in parallel computing, where neighboring parallel partitions must communicate data at interfaces to propagate dynamics throughout the full domain. Unlike parallel partitioning, the Schwarz iterative method provides a simple and rigorous coupling mechanism for subdomains which may be governed by dissimilar physics, time integrators, or mesh topologies. The utility of this “non-overlapping” finite volume Dirichlet-Dirichlet interface will be demonstrated empirically in Section 4.

The convergence of the Schwarz alternating method has been studied by a number of authors for a variety of different PDEs, e.g., [51, 52, 84, 29, 58, 59]. Analysis of the method for the specific case of nonlinear conservation laws such as those considered here was examined by Gander and Rohdi [30], and is not repeated here.

While we are not aware of any publications discussing in general the convergence of the non-overlapping Schwarz alternating method with Dirichlet-Dirichlet transmission BCs for the specific case of CCFV discretizations, it is shown by Dolean et al. [22, 23] that the Dirichlet-Dirichlet form of the method can be convergent for certain linear hyperbolic PDEs regardless of the underlying spatial discretization. Attention is restricted to the Cauchy–Riemann [22] and Maxwell equations [23]. It is demonstrated that, when discretized in time, an initial boundary value problem (IBVP) having Dirichlet BCs for one of these PDEs is equivalent to an elliptic problem having Robin BCs. Since the Schwarz alternating method is provably convergent for elliptic problems with Robin-Robin transmission BCs [84, 28], the authors make the argument that the Dirichlet-Dirichlet Schwarz iteration will converge for the hyperbolic IBVPs from which the elliptic PDE was derived. We attempted to repeat the analysis in [22, 23] for the nonlinear hyperbolic conservation laws considered here, but were unable to show the same equivalence between the time-discretized versions of these PDEs and a corresponding elliptic PDE, due largely to the presence of nonlinear terms in the equations of interest. As discussed earlier in Section 2.2, we believe that the convergence of the non-overlapping Schwarz iteration with Dirichlet-Dirichlet transmission BCs observed here (see Section 4) is due to the small overlap introduced into the discretization via ghost cells used for BC enforcement (Figure 2). This convergence phenomenon can be explained using an intuitive argument. In a typical node-centered discretization, applying a Dirichlet-Dirichlet transmission BC on a non-overlapping boundary will not allow sufficient propagation of solution information from the interior of one subdomain onto the boundary of its neighboring subdomain. This is not the case, however, in discretizations where the primary degrees of freedom are defined at the cell centers and boundary data are communicated via ghost cells. Here, data from the interior of each subdomain (half an element in) are communicated through the Dirichlet BCs imposed on the Schwarz boundaries, and it is for this reason that the Dirichlet-Dirichlet Schwarz iteration is convergent for CCFV discretizations.

It is worth noting that the ghost cells used in CCFV methods introduce only a minimum amount of overlap (Figure 2). It is well-known [51, 28, 59, 78] that the convergence rate of the overlapping Schwarz alternating method is directly related to the size of the overlap region; faster convergence is expected when a larger overlap region is defined. We thus expect our Dirichlet-Dirichlet Schwarz method to require more iterations to reach convergence when used to couple non-overlapping subdomains discretized using CCFV discretizations (which are really overlapping domains having a minimal amount of overlap). Even with this inefficiency, we are often able to achieve speedups for our PROM-based couplings with respect to a corresponding monolithic FOM, as shown in Section 4. This is partially due to the fact that employing a non-overlapping DD eliminates redundant calculations for cells in the overlap region during the Schwarz iteration process. It may be possible to improve convergence of our non-overlapping Schwarz method through the development of optimized transmission BCs [28, 19, 54, 32, 61, 29, 35, 62], which can have a mix of boundary conditions, as discussed earlier in Remark 1. We do not attempt that here, but it may be an interesting extension to explore in future work.

To begin, PROMs propose an approximation of the solution $\mathbf{u}$ which is represented by a relatively small number $M\ll N$ of degrees of freedom. While nonlinear approximations have been proposed [47], we restrict our discussion to linear approximations of the form

where $\mathbf{\widetilde{u}}:[0,T]\rightarrow N$ is the approximate state, and $\mathbf{\overline{u}}\in\mathbb{R}^{N}$ is a centering (or reference) state. The trial basis $\mathbf{\Phi}:=[\boldsymbol{\phi}_{1},\ \ldots,\ \boldsymbol{\phi}_{M}]\in\mathbb{R}^{N\times M}$ is composed of $M$ trial basis vectors $\boldsymbol{\phi}_{i}$, while $\mathbf{\widehat{u}}:[0,T]\rightarrow M$ are the generalized coordinates (or modal coefficients) associated with the trial basis.

Inserting the approximate solution in (8) into the governing ODE (6) results in the formulation

which is still an $N$-dimensional ODE. Projection by an appropriate test basis $\mathbf{\Psi}\in\mathbb{R}^{N\times M}$ and rearranging terms arrives at the desired reduced-order model

where we define the term $\mathbf{\Pi}:=[\mathbf{\Psi}^{\top}\mathbf{\Phi}]^{-1}\mathbf{\Psi}^{\top}\in\mathbb{R}^{M\times N}$ for notational simplicity. Ignoring for the moment the fact that evaluating the non-linear function still scales with $N$, (10) is an $M$-dimensional ODE, and its integration in time should incur a much lower computational cost than that of (6). Naturally, the quality of the state approximation in (8) and the robustness of the PROM ODE in (10) are highly dependent on the choice of trial basis $\mathbf{\Phi}$ and test basis $\mathbf{\Psi}$. We address each in turn.

As noted in the literature [11], deliberately sampling mesh elements at boundaries promotes the incorporation of boundary condition effects on HPROM dynamics. While one may reasonably expect that certain boundary conditions permit zero sampling (e.g., supersonic outflows) or low sampling rates (e.g., far field boundaries), Schwarz interfaces which transmit complex interior dynamics often require careful sampling treatments. In the case of extremely sparse or non-existent sampling at Schwarz interfaces, empirical evidence (not presented here) suggests that the inability to accurately transmit interface information invariably leads to an unstable solution. This follows logically from the spatially-local nature of conservation laws, resulting from an inability to propagate waves across boundaries.

In this paper, we address this issue with the fairly naive solution of sampling all Schwarz interfaces at a fixed interval $N_{b}$. For example, if $N_{b}=5$, then every fifth cell at the boundary is sampled (separate from any additional random or greedy sampling). In the absence of a priori system knowledge which might inform more intelligent sampling (e.g., knowing where a shock forms), this is a simplistic yet effective sampling methodology. However, as will be shown in Section 4, the interface sampling interval $N_{b}$ must be balanced under a fixed sampling budget $N_{s}$ to ensure both robust interface transmission and interior solution accuracy. That is, decreasing the sampling interval necessarily draws sample points from the interior of each subdomain, often decreasing the accuracy of the unsteady solution. This behavior is illustrated in Figure 4.

The versatility of the Schwarz alternating method lies in its minimal intrusiveness, relying only on a means of projecting solution information from one subdomain onto any interfaces with neighboring subdomains. This flexibility extends to coupling PROMs and HPROMs with very little additional machinery. In the case of PROMs for nonlinear systems, the approximate full order solution $\mathbf{\widetilde{u}}$ must already be accessible at every Schwarz sub-iteration in order to compute the nonlinear terms $\mathbf{{N}}(\cdot)$. The application of the projection operators $P_{ij}$ (detailed previously in Table 1) is thus trivial, as they are unchanged from an equivalent FOM-FOM coupling scenario.

Coupling HPROMs requires slightly more care. Recall the concept of the sample mesh described in Section 3.1.3: any sampled degrees of freedom require access to adjacent mesh elements to compute the fully-discrete residual. Sample mesh elements located at a Schwarz interface thus require access to mesh elements in neighboring subdomains which may not “exist,” insofar as they are not included in the neighboring subdomain’s sample mesh set. Especially for large-scale applications with memory restrictions, any elements not explicitly included in the sample mesh will only be accessible in post-processing. Accommodating these additional neighboring degrees of freedom simply requires appending their indices to the sample set, updating the solution at those elements throughout the Schwarz algorithm, and communicating it to neighboring subdomains appropriately. These operations, and the modifications to the projection operators $P_{ij}$ acting on the sampling operator $\mathbf{S}$, largely only require some additional accounting steps. Otherwise, the Schwarz algorithm is unchanged.

As the Schwarz algorithm is fairly agnostic to the choice of model in each subdomain, arbitrary couplings between FOM, PROM, and HPROM subdomains are possible. However, the goal of reduced-order modeling is minimizing the cost of evaluating a given system. Any decomposed system which includes any FOM subdomains (or PROM subdomains for nonlinear systems) will therefore be extremely limited in its ability to achieve any computational speedups. This has been previously noted in [4], where FOM-PROM or FOM-HPROM coupling fails to achieve significant runtime reductions. In this work, we will restrict analysis to PROM-PROM coupling (as an intermediate step) and HPROM-HPROM coupling.

Before introducing hyper-reduction with the aim of achieving computational cost savings, it is important to investigate the characteristics of PROM coupling relative to monolithic PROM solutions. To begin, we examine the quality of the trial spaces induced by the monolithic and decomposed domains. The left column of Figure 7 shows, for each problem, the projection error of the solution onto both the monolithic and decomposed trial spaces, for a fixed number of trial basis modes $M$ in the monolithic domain and each subdomain. It is immediately clear that for a given basis size $M$, the decomposed trial bases more accurately represent the solution as a whole, generally halving the error relative to the monolithic projection error at a given parameter value. This is an unsurprising result, as the subdomain-local trial basis may be optimized to better represent spatially-local dynamics. The decrease in error is smaller at higher basis sizes and test parameter values, as shown in the left image of Figure 7 for the Euler equations with $M=200$. This indicates the limits of the trial basis in generalizing to unseen parametric data for which the spatial locality is unable to compensate.

Of primary concern in constructing Schwarz-based solvers is the amount of overlap between neighboring subdomains. Classical Schwarz wisdom dictates that under sufficient convergence criteria, the Schwarz solution should converge to the monolithic solution regardless of the overlap size. However, a larger overlap generally results in faster convergence with fewer Schwarz iterations, at the cost of redundant solutions for mesh elements in the overlapping regions. In the case of coupling PROMs, one should never expect the decomposed PROM solution to converge exactly to the monolithic PROM solution, as the decomposed trial space is not the same as the monolithic trial space. One should expect, however, that the overlap size has relatively little effect on the solution accuracy. This is confirmed in the right column of Figure 7, which shows, for each problem, the online coupled PROM error as the overlap size and number of trial basis modes are varied. Indeed, the overlap size has a minimal effect on the solution accuracy when varied from zero overlap to 20 overlapping cells. In general, the zero overlap solution achieves the lowest error for any given parameter value, though this distinction is marginal.

The effect of overlap size on Schwarz convergence speed is more significant, largely confirming the classical Schwarz wisdom, though this appears to be somewhat problem-dependent. The uppermost plots in Figure 8 display the average number of Schwarz iterations for decomposed PROMs of varying overlap sizes and fixed trial basis sizes. For a fairly smooth problem such as the SWE, the overlap size has little effect, achieving convergence in roughly three iterations for every overlap size and parameter value. For Burgers’ and the Euler equations, which are characterized by strong shocks propagating across Schwarz interfaces, Figures 8 and 8 reveal not only a generally slower convergence (approximately four iterations for $N_{o}>0$), but particularly a sizable slowing of convergence for $N_{o}=0$. On average, this accounts for roughly two additional Schwarz iterations per time step. The effect is even stronger for higher diffusion coefficients $D$ (Burgers’) and higher initial pressure $p_{4}$ in the upper right quadrant (Euler).

It is interesting to remark that the increase in the average number of Schwarz iterations does not directly translate to a more expensive simulation. Increasing the overlap size naturally increases the computational cost for a single Schwarz iteration as each subdomain contains commensurately more mesh elements. The lower plots in Figure 8 display normalized runtime measurements. For the SWE in Figure 8, the added cost of redundant cells in the overlap region outweighs any slim cost savings due to fewer Schwarz iterations, and the case where $N_{o}=0$ is both the most accurate and the fastest. For the Burgers’ and Euler equations, the case where $N_{o}=0$ is usually (though not always) the most costly. While for the Burgers’ equation at high $D$ this disparity is quite significant, for all other parameter values and for the Euler equations the costs are roughly comparable. On the whole, setting $N_{o}=0$ (the cell-centered finite volume Dirichlet-Dirichlet “non-overlapping” scenario) interestingly leads to an accurate coupled solution with minimal impact on computation expense, except in certain extremes. Thus, for the remainder of these results, we restrict analysis to the case for $N_{o}=0$.

Recall that PROMs for sufficiently complex non-linear systems rarely achieve runtime reductions relative to the full-order model due the cost of evaluating the nonlinear terms in the FOM dimension $N$. We now introduce hyper-reduction via collocation to achieve computational cost savings. Similar to the PROM results shown above, all coupled HPROM results are computed on a $2\times 2$ decomposed domain and utilize an HPROM model in every subdomain. For all HPROM domains (both monolithic and decomposed), the trial basis size is fixed at a sufficient level to at least guarantee low projection error. The SWE experiments are fixed at $M=80$, the Burgers’ equation at $M=150$, and the Euler equations at $M=200$.

In generating the sample meshes, the sample indices set $\mathcal{S}$ is first seeded with any Schwarz interface samples as dictated by the interface sampling interval $N_{b}$. The Burgers’ equation sample meshes are further seeded with the degrees of freedom selected by QDEIM [24], computed from a trial basis of size $M=200$. The remainder of the sample points under a fixed sample mesh size $N_{s}$ are selected randomly, with higher sampling rates simply appending additional sampled mesh elements to those selected at lower sampling rates, ensuring uniformity in all numerical experiments. A range of sample mesh sizes are investigated, ranging roughly logarithmically from 0.1% to 10% of the full mesh size, assuming the number of seed points does not already exceed this sample mesh size. In general, all HPROM solutions for the SWE are stable above $N_{s}=0.1\%\times N$, while solutions for the Burgers’ and Euler equations requires a sampling rate of more than $1.0\%\times N$. This is a natural consequence of the system dynamics. Whereas the SWE generate a fairly smooth solution which is easily represented by a sparse sample mesh, the Burgers’ and Euler equations require relatively dense sample meshes to accurately propagate shocks.

The body of hyper-reduction literature confirms the intuition that increasing the sample mesh size generally improves solution accuracy and robustness. For coupled HPROMs, however, the question of interface sampling is more nuanced. For all cases investigated here, failing to seed the sample mesh with interface mesh elements leads to an unstable solution, except at extremely high sampling random rates (10% $\times N$ for the SWE). We must then deliberately seed the sample mesh with elements at the Schwarz interface, in this case utilizing a simple fixed sampling interval $N_{b}$, as described in Section 3.2. As stated previously, decreasing the sampling interval naturally improves the ease of transmitting interface information between neighboring subdomains, but draws sample mesh elements away from the subdomain interior and degrades overall solution accuracy. This balancing act is illustrated in Figure 9, where an optimal interface sampling interval generates accurate and stable coupled HPROMs, while values of $N_{b}$ with are either too high or too low may generate unstable solutions. For the SWE, this optimum appears to be around $N_{b}=10$, while for the Burgers’ and Euler equations this is approximately $N_{b}=2$, or every other interface cell. The disparity between these optimums is again due to the smoothness of the SWE solution, in contrast to the shocks experienced by the Burgers’ and Euler solutions. The ability to accurately transmit a shock across the Schwarz interface reasonably requires a much higher sample mesh resolution.

Under sufficient sampling rates for a stable solution, it is useful to know whether the Schwarz coupling induces errors at interfaces. That is, as each subdomain has a unique trial space which does not guarantee perfect continuity across interfaces, are there inordinately high error levels at the interfaces? Or are the Schwarz convergence criteria sufficient to ensure accurate transmission between disparate trial spaces? For a fixed trial basis dimension $M$ and sample mesh size $N_{s}$, Figure 10 displays representative monolithic FOM, monolithic HPROM, and decomposed HPROM solutions at $t=T$. In a broad sense, the decomposed HPROM solution mirrors that of the monolithic HPROM closely: slight “ringing” effects are apparent in the Burgers’ solution, and more so in the Euler solution. The decomposed SWE solution roughly appears to approximate the FOM solution more closely than does the monolithic HPROM. Most notably, there are not any apparent discontinuities at the subdomain boundaries which would signal poor convergence at the interface.

To examine the error behavior more precisely, Figure 11 shows the time-average absolute error for the same monolithic and decomposed HPROM solutions. If the Schwarz coupling induces error at the interface, one might expect high error concentration along the subdomain boundaries. This does not appear to be the case, as not only are errors generally lower for the decomposed solution (again, the greater accuracy of subdomain-local trial bases), but error patterns are largely the same (dominated by features of the flow dynamics, such as shocks). This indicates that the Schwarz coupling algorithm has no difficulty in accurately transmitting information under sufficient interface sampling conditions.

Finally, we address the effect of coupling on the cost efficiency of HPROMs. Figure 12 plots the solution error as a function of either the computational speedup with respect to the all-FOM decomposed solution or the monolithic FOM solution. For the PROMs, each point represents a different test parameter value and trial basis size $M$. For the HPROMs, each point represents a different test parameter value and sample mesh size $N_{s}$. These cost measurements take into account the three additional cores used to parallelize the additive Schwarz algorithm.

Unsurprisingly, the PROMs are generally unable to achieve any computational cost savings compared to the monolithic PROM solution, and the decomposed PROM solution is invariably slower than the monolithic PROM solution due to the added cost of Schwarz iterations. The HPROM solutions, on the other hand, are all capable of achieving noticeable speedups relative to the decomposed FOM solution. The decomposed HPROM solutions broadly achieve lower error levels than the monolithic HPROM solutions, but do not reach the same maximum speedup levels. For the SWE, the monolithic HPROM is able to achieve speedups over 1,000 times faster than the decomposed FOM, while the decomposed HPROM is able to achieve roughly 100 times speedup at best. The same is true for the Burgers’ and Euler equations, where the monolithic and decomposed HPROM solutions achieve a maximum speedup of 100 times and 10 times, respectively, over the decomposed FOM. Note that especially for the shock-driven problems, the decomposed HPROMs are not able to achieve significant speedups over the monolithic FOM, as they require relatively dense sample meshes to capture the shock dynamics and low interface sampling intervals to transmit shocks across subdomain boundaries. That being said, these results indicate that the Schwarz alternating method can robustly couple HPROMs and achieve good cost savings relative to an equivalent coupled FOM solution. For systems which must be coupled, Schwarz coupling may offer appreciable computational speedups.