Kathrin Smetana
email: ksmetana@stevens.edu

Tommaso Taddei
email: tommaso.taddei@inria.fr

Localized model reduction for nonlinear elliptic partial differential equations: localized training, partition of unity, and adaptive enrichment

Keywords: parameterized partial differential equations, model order reduction, domain decomposition

HAL is

Introduction 1.Component-based model reduction for parameterized PDEs

Numerical modeling and simulation is of paramount importance to predict the response, improve the design, and monitor the structural health of engineering systems. Several problems of interest involve repeated solutions of a partial differential equation (PDE) for many values of the model parameters or require real-time responses: these tasks are prohibitively expensive for standard (e.g., finite element) methods. Parametric model order reduction (pMOR, [START_REF] Haasdonk | Reduced basis methods for parametrized PDEs -a tutorial[END_REF][START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF]) aims to reduce the marginal cost associated with the solution to parameterized systems over a range of parameters. The goal of this paper is to develop a pMOR procedure for large-scale nonlinear elliptic PDEs with parameter-induced topology changes.

pMOR techniques may rely on an offline/online decomposition to reduce marginal costs. During the offline phase, we rely on several high-fidelity (HF) solves to generate a reduced-order model (ROM) for the solution field. During the online phase, given a new value of the parameter, we query the ROM to estimate the solution field and relevant quantities of interest. Monolithic pMOR methods rely on HF solves at the training stage, which might be unaffordable for very large-scale problems. Furthermore, pMOR methods rely on the assumption that the solution field is defined over a parameter-independent domain or over a family of diffeomorphic domains: they thus cannot deal with problems for which parametric variations induce topology changes.

To address these issues, several authors have proposed component-based pMOR procedures (cf. [START_REF] Huynh | A static condensation reduced basis element method: approximation and a posteriori error estimation[END_REF] and the review [START_REF] Buhr | Localized model reduction for parameterized problems[END_REF]). During the offline stage, we define a library of archetype components and we build local reduced-order bases (ROBs) and local ROMs; then, during the online stage, we instantiate components to form the global system and we estimate the global solution by coupling local ROMs. CB-pMOR strategies consist of two distinct building blocks: (i) a rapid and reliable domain decomposition (DD) strategy for online global predictions, and (ii) a localized training strategy exclusively based on local solves for the construction of the local approximations.

CB-pMOR shares important features with multiscale methods [START_REF] Babuška | Optimal local approximation spaces for generalized finite element methods with application to multiscale problems[END_REF][START_REF] Efendiev | Multiscale finite element methods for porous media flows and their applications[END_REF][START_REF] Bris | An msfem type approach for perforated domains[END_REF][START_REF] Strouboulis | The design and analysis of the generalized finite element method[END_REF][START_REF] Owhadi | Bayesian numerical homogenization[END_REF][START_REF] Owhadi | Multigrid with rough coefficients and multiresolution operator decomposition from hierarchical information games[END_REF][START_REF] Owhadi | Localized bases for finite-dimensional homogenization approximations with nonseparated scales and high contrast[END_REF][START_REF] Owhadi | Polyharmonic homogenization, rough polyharmonic splines and sparse superlocalization[END_REF][START_REF] Målqvist | Localization of elliptic multiscale problems[END_REF][START_REF] Målqvist | Numerical homogenization by localized orthogonal decomposition[END_REF][START_REF] Smetana | Optimal local approximation spaces for component-based static condensation procedures[END_REF][START_REF] Ma | Novel design and analysis of generalized fe methods based on locally optimal spectral approximations[END_REF][START_REF] Grasedyck | The AL basis for the solution of elliptic problems in heterogeneous media[END_REF][START_REF] Verfürth | Numerical homogenization for nonlinear strongly monotone problems[END_REF][START_REF] Liu | Iterated numerical homogenization for multi-scale elliptic equations with monotone nonlinearity[END_REF][START_REF] Chen | Solving and learning nonlinear pdes with gaussian processes[END_REF][START_REF] Chen | Randomized sampling for basis function construction in generalized finite element methods[END_REF][START_REF] Calo | Randomized oversampling for generalized multiscale finite element methods[END_REF][START_REF] Chen | Manifold learning and nonlinear homogenization[END_REF][START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF][START_REF] Hughes | Multiscale phenomena: Green's functions, the dirichlet-to-neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods[END_REF]. Similarly to CB-pMOR, multiscale methods rely on local solves to build suitable approximation spaces that are tailored to the problem of interest. The emphasis in CB-pMOR is to devise and then exploit a library of inter-operable archetype components and associated ROMs that can be used for a broad range of potentially parameter dependent problems in a specific domain of interest.

Domain decomposition strategies within CB-pMOR

Since the seminal work by Maday and Rønquist [START_REF] Maday | A reduced-basis element method[END_REF] -that proposed a non-overlapping non-conforming reduced basis element method based on Mortar DD -several authors have combined DD methods with model reduction methods to devise effective CB-ROMs. As discussed in detail in the review [START_REF] Buhr | Localized model reduction for parameterized problems[END_REF], we can distinguish between conforming non-overlapping approaches [START_REF] Huynh | A static condensation reduced basis element method: approximation and a posteriori error estimation[END_REF][START_REF] Eftang | Port reduction in parametrized component static condensation: approximation and a posteriori error estimation[END_REF][START_REF] Smetana | Optimal local approximation spaces for component-based static condensation procedures[END_REF], non-conforming non-overlapping approaches based on Lagrange multipliers [START_REF] Hoang | Domain-decomposition least-squares petrov-galerkin (DD-LSPG) nonlinear model reduction[END_REF][START_REF] Iapichino | A reduced basis hybrid method for the coupling of parametrized domains represented by fluidic networks[END_REF][START_REF] Maday | A reduced-basis element method[END_REF][START_REF] Pegolotti | Model order reduction of flow based on a modular geometrical approximation of blood vessels[END_REF], non-conforming non-overlapping approaches based on discontinuous Galerkin (DG) coupling [START_REF] Albrecht | The localized reduced basis multiscale method[END_REF][START_REF] Antonietti | A discontinuous Galerkin reduced basis element method for elliptic problems[END_REF][START_REF] Ohlberger | Error control for the localized reduced basis multi-scale method with adaptive on-line enrichment[END_REF][START_REF] Riffaud | The DGDD method for reduced-order modeling of conservation laws[END_REF], and overlapping methods [START_REF] Bergmann | A zonal Galerkin-free POD model for incompressible flows[END_REF][START_REF] Buhr | Randomized local model order reduction[END_REF]. The vast majority of existing contributions (with few recent exceptions [START_REF] Baiges | A domain decomposition strategy for reduced order models. Application to the incompressible Navier-Stokes equations[END_REF][START_REF] Bergmann | A zonal Galerkin-free POD model for incompressible flows[END_REF][START_REF] Hoang | Domain-decomposition least-squares petrov-galerkin (DD-LSPG) nonlinear model reduction[END_REF][START_REF] Pegolotti | Model order reduction of flow based on a modular geometrical approximation of blood vessels[END_REF]) is restricted to linear PDEs.

In this work, we rely on the partition-of-unity method (PUM) to devise global approximation spaces from local reduced spaces, and on Galerkin projection to compute the global state estimate. PUM was proposed by Babuška and Melenk in [START_REF] Babuška | The partition of unity method[END_REF][START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF] and further developed and analyzed in the framework of generalized finite element methods for multiscale problems (cf. [START_REF] Babuška | Multiscale-spectral gfem and optimal oversampling[END_REF]); PUM was also considered in the pMOR literature for linear elliptic and parabolic problems [START_REF] Buhr | Randomized local model order reduction[END_REF][START_REF] Schleuß | Optimal local approximation spaces for parabolic problems[END_REF]. In the CB-pMOR framework, PUM offers a general (i.e., independent of the underlying PDE) framework with strong theoretical guarantees.

Localized training based on oversampling and randomization

Given the domain Ω associated with a given archetype component, oversampling methods consist in (i) defining a patch U ⊃ Ω and a suitable local PDE problem in U , (ii) solving the local PDE for several choices of the boundary conditions on ∂ U and then restricting the solution to Ω, and finally (iii) exploiting the results to build a local approximation space for the solution in Ω. Randomized methods rely on independent identically distributed (iid) samples of the boundary conditions on (a subset of) ∂ U : they thus require the introduction of a probability density function for the functions defined on ∂ U .

Oversampling methods exploit low-pass filtering properties of the differential operator to identify lowdimensional structures: we refer to [START_REF] Taddei | Model order reduction methods for data assimilation; state estimation and structural health monitoring[END_REF]Chapter 5] and [START_REF] Smetana | Optimal local approximation spaces for component-based static condensation procedures[END_REF]Remark 3.3] for two representative working examples. In detail, Caccioppoli type inequalities (see e.g. [START_REF] Giaquinta | Caccioppoli's inequality and Legendre-Hadamard condition[END_REF]) provide the theoretical foundations for the application of oversampling methods to a particular class of PDEs. Oversampling methods have been suggested and used extensively in the context of multiscale methods (see e.g. [START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF][START_REF] Babuška | Optimal local approximation spaces for generalized finite element methods with application to multiscale problems[END_REF][START_REF] Henning | Oversampling for the multiscale finite element method[END_REF][START_REF] Målqvist | Localization of elliptic multiscale problems[END_REF] and references therein) and then used as well in CB-pMOR [START_REF] Eftang | Port reduction in parametrized component static condensation: approximation and a posteriori error estimation[END_REF][START_REF] Smetana | Optimal local approximation spaces for component-based static condensation procedures[END_REF] for linear PDEs.

As suggested in [START_REF] Buhr | Randomized local model order reduction[END_REF], randomized oversampling methods for linear parameter-independent PDEs can be linked to randomized singular value decomposition (SVD) techniques developed and analyzed in randomized numerical linear algebra [START_REF] Halko | Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions[END_REF][START_REF] Martinsson | A randomized algorithm for the decomposition of matrices[END_REF][START_REF] Mahoney | Randomized algorithms for matrices and data[END_REF][START_REF] Sarlos | Improved Approximation Algorithms for Large Matrices via Random Projections[END_REF][START_REF] Rokhlin | A randomized algorithm for principal component analysis[END_REF][START_REF] Drineas | RandNLA: Randomized Numerical Linear Algebra[END_REF][START_REF] Mahoney | CUR matrix decompositions for improved data analysis[END_REF][START_REF] Woodruff | Sketching as a tool for numerical linear algebra[END_REF]: this link allows to extend methodological and theoretical contributions in randomized linear algebra to CB-pMOR. In particular, we can exploit concentration inequalities (see e.g. [START_REF] Boucheron | Concentration inequalities. A nonasymptotic theory of independence[END_REF]) to analyze the error of randomized techniques, and inform the choice of the sampling distribution. The influence of the choice of the sampling distribution for nonlinear PDEs remains an open question in CB-pMOR.

Contributions of the paper and outline

In this work, we propose a CB-pMOR procedure based on PUM for parametric nonlinear elliptic PDEs; we do not require the nonlinear operator to be monotone. The contributions of the paper are twofold. First, we propose a randomized data compression algorithm based on oversampling: the approach relies on random samples of local parameters and boundary conditions on the oversampling boundary. We propose a new sampler that controls the smoothness of the boundary condition, and we empirically demonstrate its effectiveness for a nonlinear diffusion problem. Second, we propose a basis enrichment algorithm that relies on global reduced solves to enrich the local reduced spaces. The algorithm relies on a local residual-based error indicator to identify boundary conditions for which the local ROM is inaccurate and a rigorous global a posteriori error bound as a termination criterion. We prove in-sample a priori exponential convergence of the enrichment algorithm for linear coercive problems; we further investigate performance for a nonlinear diffusion problem.

Our randomized algorithm reads as a randomized proper orthogonal decomposition [START_REF] Yu | A randomized proper orthogonal decomposition technique[END_REF] with respect to parameter and boundary conditions. On the other hand, the enrichment algorithm is closely related to the online enrichment strategy proposed in [START_REF] Ohlberger | Error control for the localized reduced basis multi-scale method with adaptive on-line enrichment[END_REF] for non-overlapping DG DD, and to the residual-based online enrichment algorithm considered in [START_REF] Buhr | Exponential convergence of online enrichment in localized reduced basis methods[END_REF] for linear problems. The major difference is that the enrichment is performed at training stage and aims to update the local approximation spaces associated with the archetype components, rather than during the online stage on the "instantiated components".

The outline of the paper is as follows. In section 2, we introduce the model problem considered throughout the paper to illustrate the main definitions and to numerically validate our proposal: the model problem involves a high-dimensional (O(10 2 )) parameterization and topology changes. In section 3, we discuss the DD strategy based on the PUM and we introduce local and global discrete approximation spaces; in section 4, we discuss the randomized localized data compression; in section 5, we present the enrichment strategy; and in section 6 we present thorough numerical investigations for the model problem. Section 7 concludes the paper.

Model problem: nonlinear diffusion

Given n dd ∈ N and H = 0.1, we define the domains

Ω i+(j-1)n dd = {[x 1 + H(i -1), x 2 + H(j -1)] : x 1 , x 2 ∈ (0, H)} , (1)
for i, j = 1, . . . , n dd , and Ω = N dd k=1 Ω k with N dd = n 2 dd . We further introduce P = [0.1, 0.2] × [START_REF] Grasedyck | The AL basis for the solution of elliptic problems in heterogeneous media[END_REF][START_REF] Bris | An msfem type approach for perforated domains[END_REF], the permeability coefficient κ : Ω × R × N dd i=1 P → R + such that κ Ωi = κ x; u, µ (1) , . . . , µ (N dd ) Ωi satisfies

κ Ωi = 36
µ (i) 2   u(1 -u) u 3 + 12 µ (i) 2 (1 -u) 3   2 + µ (i) 1 , i = 1, . . . , N dd ; (2a)
and the source term

f (x; i ) = 100 e -50 x-x c,i 2 2 1 Ω i (x). (2b)
Then, we introduce the PDE model in strong form: given µ = [µ (1) , . . . , µ (N dd ) , i ] ∈ P glo (n dd ) :=

N dd i=1 P × {1, . . . , N dd }, find u µ such that -∇ • (κ µ (u µ ) ∇u µ ) = f µ in Ω, u µ = 0 on ∂Ω. (2c)
In fig. 1(a)-(c)-(d), we show the domain Ω and selected snapshots for N dd = 100. We discretize the problem using a Q3 spectral element method based on a structured grid with 961 degrees of freedom in each subdomain Ω i . The PDE model (2c) has been previously considered in [START_REF] Smetana | Hierarchical model reduction of nonlinear partial differential equations based on the adaptive empirical projection method and reduced basis techniques[END_REF], and is inspired by the model for immiscible two-phase flows in porous media studied in [START_REF] Michel | A finite volume scheme for two-phase immiscible flow in porous media[END_REF]. While we focus on eq. (2c) in this paper to ease the exposition of ideas, we emphasize that the proposed methods can be readily applied to other nonlinear elliptic PDEs.

In the following, we devise a component-based reduced-order model (CB-ROM) for (2c): the ROM should handle arbitrary choices of N dd and of the parameters µ (1) , . . . , µ (N dd ) ∈ P and i . In view of the presentation of the methodology, we introduce the overlapping partition that will be used to define the partition of unity in section 3,

{ω i } N dd i=1 , ω i = x ∈ Ω : min y∈Ωi x -y ∞ < δ over , i = 1, . . . , N dd . (3)
Note that i ω i = Ω. As explained in section 3, for non-homogeneous Dirichlet or Neumann boundary conditions we shall require i ω i ⊃⊃ Ω. Here, δ over is the overlapping size; in the numerical experiments, we consider δ over = 0.1H. We further introduce the archetype components that we use to describe the system: the "corner" (co) component is associated with the corner elements of the partition {ω i } i ; the "edge" (ed) component is associated with the edge elements of {ω i } i ; the "internal" (int) component is associated with the internal elements of {ω i } i (see fig. 1(b)). We denote by Ω co , Ω int , Ω ed the master elements associated with the three archetype components. For edge and corner components, we denote by Γ ed dir , Γ co dir the Dirichlet boundaries; furthermore, we introduce the local HF spaces associated with the spectral element discretization

Y int ⊂ H 1 ( Ω int ), Y ed ⊂ H 1 0, Γ ed dir ( Ω ed ), Y co ⊂ H 1 0, Γ co dir ( Ω co ) (4)
and the corresponding (semi-) norms • • with • ∈ {co, ed, int} that are introduced in section 3. We denote by L ∈ {co, ed, int} N dd the set of labels that link the elements of {ω i } i to the corresponding component; we further denote by Φ i : Ω Li → ω i the mapping from the (appropriate) component to the i-th element of the partition. We remark that the mappings Φ i are simple translations for all internal components, while they are the composition of a rigid translation and a rotation for boundary (edge and corner) components to ensure that Φ i ( Γ Li dir ) ⊂ ∂Ω and thus to facilitate the imposition of Dirichlet conditions.

3 Partition of unity method for localized model reduction

Partition of unity

In [START_REF] Babuška | The partition of unity method[END_REF][START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF], Babuška and Melenk proposed the partition of unity method (PUM) to construct ansatz spaces with local properties. As discussed in [START_REF] Babuška | The partition of unity method[END_REF], the PUM is designed to include a priori knowledge about the PDE in the ansatz spaces, and it provides a systematic way to construct ansatz spaces of any desired regularity. In the framework of CB-pMOR, the PUM provides a systematic framework to construct global ansatz spaces and offers strong theoretical guarantess concerning approximation and robustness.

Given the overlapping cover of Ω, {ω i } N dd i=1 , we denote by M the minimum constant such that

∀ x ∈ Ω, card {i ∈ {1, . . . , N dd } : x ∈ ω i } ≤ M, (5a)
where card(A) denotes the cardinality of the discrete set A. Then, we define the partition of unity (PoU)

{φ i } N dd i=1 such that          supp (φ i ) ⊂ ω i , 0 ≤ φ i (x) ≤ 1, ∇φ i L ∞ (Ω) ≤ C i , N dd j=1 φ j (x) = 1, x ∈ Ω, i = 1, . . . , N dd . (5b)
We say that

{φ i } N dd i=1 is of degree m if {φ i } N dd i=1 ⊂ C m (Ω; R).
Then, we define the PUM spaces

X pum := N dd i=1 φ i ψ i : ψ i ∈ X i ⊂ H 1 0 (Ω), (6)
where

X i = {ζ • Φ -1 i : ζ ∈ Y Li }. Note that by construction φ i ζ • Φ -1 i ∈ H 1 0 (ω i
) and can thus be trivially extended to R d . Next, given the reduced spaces {Z • } • such that Z • ⊂ Y • , we define the global reduced space

Z gfem := N dd i=1 φ i ζ i • Φ -1 i : ζ i ∈ Z Li ⊂ X pum . (7)
Theorem 3.1 provides a rigorous upper bound for the approximation properties of the PUM space in Ω -the local approximation condition [START_REF] Bergmann | A zonal Galerkin-free POD model for incompressible flows[END_REF] provides the foundations for the localized data compression strategy proposed in section 4.

Theorem 3.1. ([5, Theorem 1]) Let u ∈ H 1 0 (Ω). Assume that there exist ζ 1 , . . ., ζ N dd such that ζ i • Φ i ∈ Z Li and u -ζ i L 2 (Ω∩ωi) ≤ i , ∇u -∇ζ i L 2 (Ω∩ωi) ≤ ∇,i , i = 1, . . . , N dd , (8)
for some positive constants { i } i and { ∇,i } i . Then, the function

u gfem = N dd i=1 φ i ζ i ∈ Z gfem satisfies                u -u gfem L 2 (Ω) ≤ √ M N dd i=1 2 i ; ∇u -∇u gfem L 2 (Ω) ≤ √ 2M N dd i=1 C 2 i 2 i + 2 ∇,i . (9)
In this work, we consider a piecewise tensorized bilinear PoU {φ i+(j-1)n dd (x) = φ 1d i (x 1 )φ 1d j (x 2 )} n dd i,j=1 where {φ 1d i } i is a PoU subordinate to the cover

{ω 1d i = ((i -1)H -δ over /2, iH + δ over /2)} n dd i=1 .
For this choice of the PoU, we have that i=1 ω i = Ω. Note also that the constant M in (5a) is equal to four. We observe that there exist φ int , φ co , φ ed such that

d dx φ 1d i L ∞ (Ω) =
φ i = φ Li • Φ -1 i : i = 1, . . . , N dd ,
for any choice of n dd ∈ N: this is due to the particular choice of the mappings {Φ i } i and of the archetype components. For arbitrary partitions {ω i } i and arbitrary mappings

{Φ i } i , given { φ • } • such that 0 ≤ φ • (x) ≤ 1 in R d , φ • (x) = 0 if x / ∈ Ω • , ∇ φ • L 2 (R 2 ) ≤ C • , we can show that the set {φ i } i such that φ i = 1 N dd j=1 φ Li • Φ -1 j φ Li • Φ -1 i , ∀ i = 1, . . . , N dd , is a PoU in Ω.

Discrete variational formulation and functional norms

We introduce the local semi-norms

w • = φ • w H 1 ( Ω • ) , • ∈ {int, co, ed}. (10)
Note that for this choice of the local norms, since the mappings

{Φ i } i are roto-translations, if {ζ • j } n j=1 are orthonormal bases with respect to • • , then {φ i ζ Li j } n j=1 is orthonormal in H 1 (ω i ), for i = 1, . . . , N dd . Given the spaces X i,0 := {φ i ζ • Φ -1 i : ζ ∈ Y Li } for i = 1, .
. . , N dd , we further introduce the inner products and induced norms

(w, v) 1,ωi = ωi ∇w • ∇v + wv dx, w 1,ωi = (w, w) 1,ωi , w, v ∈ X i,0 ; (11)
the global norm w 1,Ω = Ω ∇w 2 2 + w 2 dx; and the dual norms

f -1,ωi = sup v∈Xi,0 f (v) v 1,ωi , F -1,Ω = sup v∈Xpum f (v) v 1,Ω , i = 1, . . . , N dd , (12)
for f ∈ X i,0 and F ∈ X pum . Then, we introduce the HF problem: given µ ∈ P glo , find u µ ∈ X pum such that

G µ (u µ , v) = 0 ∀ v ∈ X pum , (13a)
where

G µ (w, v) := Ω η µ (x; w, v) dx with η µ (x; w, v) = κ µ (x; w)∇w • ∇v -f µ v (13b)
and G µ : X pum → X pum .

Residual assembly and algebraic formulation of the reduced-order model

We omit dependence of Ω and P glo on n dd to shorten notation. We consider the Galerkin ROM:

find u µ ∈ Z gfem : G µ ( u µ , v) = 0 ∀ v ∈ Z gfem . (14)
Given the local approximation spaces Z int , Z ed , Z co with bases1 {ζ • i } n i=1 , we define the basis of Z gfem {ζ i,j } i,j such that

ζ i,j = ζ Lj i • Φ -1 i φ j , i = 1, . . . , n, j = 1, . . . , N dd . (15a)
Given u ∈ Z gfem , we set N := nN dd and we denote by u ∈ R N the vector of coefficients such that

u = N dd j=1 n i=1 (u) i+(j-1)n ζ i,j . (15b)
Then, we introduce the discrete residual R : R

N × P glo → R N such that (R µ (u)) i+(j-1)n = G µ (u, ζ i,j ) (16a)
and the algebraic nonlinear problem associated to ( 14):

find u µ ∈ R N such that R µ ( u µ ) = 0. (16b)
In order to discuss the practical evaluation of the discrete residual R µ in (16a), we define

{ u i,µ = u µ φ i } i and Neigh i = {j : ω i ∩ ω j = ∅}. Then, we observe that G µ ( u µ , ζ i,j ) = ωi η µ   x; k∈Neigh i u k,µ , ζ i,j   dx = Ω L i η (i) µ   x; k∈Neigh i u k,µ • Φ i , ζ Lj i φ Lj   dx, (17a)
where

η (i) µ (x; w, v) = κ µ (Φ i (x); w)∇Φ -1 i ∇Φ -T i ∇w • ∇v -f µ v det(∇Φ i ), (17b)
with f µ = f µ • Φ i . Since {Φ i } i are roto-translations, (17b) reduces to:

η (i) µ (x; w, v) = κ µ (Φ i (x); w)∇w • ∇v -f µ v. (17c)
We observe that the Jacobian J µ (•) of the algebraic residual R µ (•) is sparse for large values of N dd . More precisely, exploiting [START_REF] Canuto | A posteriori error analysis of the reduced basis method for nonaffine parametrized nonlinear PDEs[END_REF], it is easy to verify that the number of non-zero elements of J µ (•) is bounded by

nnz (J µ (u)) ≤ N dd i=1 n 2 card(Neigh i ) = O n 2 N dd , ∀ u ∈ R N .
For the model problem considered in this work we have card (Neigh i ) ≤ 9 for i = 1, . . . , N dd .

Assembly of the residual in ( 17) is extremely expensive due to the need to integrate over all instantiated components { Ω Li }. To speed up computations, we should thus resort to hyper-reduction techniques [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF][START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF][START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Yano | An LP empirical quadrature procedure for reduced basis treatment of parametrized nonlinear PDEs[END_REF]. The choice of the hyper-reduction procedure strongly depends on the PDE model of interest, on the underlying high-fidelity numerical scheme, and on the geometrical parameterization: we refer to [START_REF] Taddei | A discretize-then-map approach for the treatment of parameterized geometries in model order reduction[END_REF] for a discussion on the treatment of geometry parameterizations. We further observe that evaluation of (17a) involves evaluation of u j,µ in the mapped quadrature points of the mesh Ω Li : this evaluation is extremely expensive for unstructured meshes and thus requires a specialized treatment. The development of specialized hyper-reduction techniques for CB-pMOR is part of ongoing research and is not addressed in the present work.

Data compression: randomized localized training

To highlight the main features of our methodology without unnecessary notation, we assume here that the system is described by a single archetype component. We denote by Z ⊂ Y the local approximation space. For the model problem in section 2, this would correspond to the case of Neumann or Robin boundary conditions on ∂Ω. In the supplementary materials, we discuss the extension to the case of multiple components and we provide further details for the particular test case considered.

The aim of this section is to devise an actionable procedure to build a local approximation space Z ⊂ Y such that min

ζ∈Z u µ ωi -ζ • Φ -1 i 1,ωi ≤ ε tol for i = 1, . . . , N dd , µ ∈ P glo (n dd ), (18)
where ε tol > 0 is a prescribed tolerance. Condition [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF] implies that the local space Z should approximate the manifold

M = u µ ωi • Φ i : i = 1, . . . , N dd , µ ∈ P glo (n dd ), n dd ∈ N ⊂ Y. ( 19
)
The computation of snapshots that belong to the manifold M requires to solve global problems and is thus unfeasible in our framework. Instead, in section 4.1, we propose to rely on oversampling to identify an actionable localized manifold M for which we can compute snapshots; then, in section 4.2, we propose a randomized training algorithm to construct local approximation spaces.

Oversampling

We fix i ∈ {1, . . . , N dd }, and we define the patch U ⊂ R 2 with input boundary Γ in ⊂ ∂ U . We extend the mapping Φ i to U and we define U i := Φ i ( U ) -for the considered model problem, the mappings {Φ i } i are linear maps that can be trivially extended to R 2 . As depicted in fig. 2, we consider U i = j∈Neigh i Ω j where Neigh i is introduced in section 3. We denote by u µ,i the restriction of the solution u µ to U i and we define u µ,i := u µ,i • Φ i . We observe that u µ,i solves the problem (cf. (17c)):

U κ µ ( u µ,i )∇ u µ,i • ∇v dx = U fµ v dx, ∀ v ∈ Y ovr i,0 , (20a)
with u µ,i | Γin = u µ,i •Φ i and Y ovr i,0 = {v •Φ i : v| Ui ∈ X pum , v| ∂U = 0} ⊂ H 1 0 ( U ).
We observe that u µ,i is a function of the subset of parameters {µ (j) } j∈Neigh i , and of the index i : we can then define the active set of parameters P co = N co dd i=1 P × {1, . . . , N co dd , 0} where N co dd = card(Neigh i ) and i = 0 means that the source term is outside the patch. The parameterization P co is associated to the archetype component of interest and is independent of the size of the system (i.e., the number of subdomains N dd ). Exploiting (20a), we define the transfer operator T : G × P co → Y such that T µ (g) = u| Ω where G ⊂ H 1/2 ( Γ in ), u satisfies (20a) with u| ∂ U \ Γin = 0 and u| Γin = g.

In the implementation, we replace Y ovr i with a standard discretization of H 1 0 ( U ). We define the (unknown) set G true ⊂ H 1/2 ( Γ in ) that contains all possible restrictions of the solution field to the input boundary for all instantiated components ω i , all parameters, and all choices of n dd ; clearly, we have M = {T µ (g) : g ∈ G true , µ ∈ P co }. If we introduce the "approximation" G of G true , we obtain the localized manifold: We observe that snapshots of M can be computed by solving local problems in the patch U for prescribed choices of the active parameters µ ∈ P co and the boundary conditions. The patch U should be significantly smaller than Ω to ensure rapid computations; at the same time, U should be large enough to ensure decay of highfrequency modes on Γ in .

M = {T µ (g) : g ∈ G, µ ∈ P co } . (20b)
The choice of the set of boundary conditions G is of paramount importance; clearly, G should be rich enough to ensure that sup w∈M dist(w, M) ≤ ε tol . Since the problem is nonlinear, generating a discrete representative approximation of the high-dimensional set G is also particularly challenging. In the next section, we directly prescribe a probability density function (pdf) p bc of the space of boundary conditions: the set G is thus defined as the support of the pdf p bc .

Algorithm 1 Randomized localized training

Inputs: ntrain size of training set, pµ, p bc pdfs.

Output: Z local approximation spaces.

1: Generate µ (k) iid ∼ p µ , g (k) iid ∼ p bc , k = 1, . . . , n train . 2: Compute u k = T µ (k) (g (k)
) for k = 1, . . . , n train .

3: Z = POD u k ntrain k=1 , (•, •), n .

Randomized training

We introduce notation

Z = POD {u (i) } ntrain i=1 , (•, •), n
to refer to the application of proper orthogonal decomposition (POD, [START_REF] Volkwein | Model reduction using proper orthogonal decomposition[END_REF]) to the snapshot set {u (i) } ntrain i=1 with inner product (•, •); here, n denotes the number of POD modes in the output space Z. We further denote by p µ and p bc the pdfs for parameter and boundary conditions.

Algorithm 1 illustrates the randomized training procedure. The algorithm reads as a randomized POD [START_REF] Yu | A randomized proper orthogonal decomposition technique[END_REF] with respect to parameter and boundary conditions: the inputs of the algorithm are the number of training points n train , the size of the sought reduced spaces n, and the pdfs {p µ , p bc } for the archetype component; the output is the reduced space Z.

It is well-known that POD is optimal in L 2 (p µ × p bc ) in the limit n train → ∞; however, since the pdfs p µ , p bc are chosen a priori, they might not be representative of the true distributions for the global systems. Provided that additional information on the class of global systems of interest is available, these observations motivate the enrichment strategy proposed in section 5.

Remark 4.1. Probabilistic a posteriori error estimation. Given n test additional simulations {u (i) } ntest i=1 ⊂ M and the space Z, we introduce the error indicator

E := 1 n test ntest i=1 u (i) -Π Z u (i) u (i) , (21)
which measures the average relative projection error on the test set {u (i) } i . Here, Π Z : Y → Z is the projection operator on Z. Provided that u (i) = T µ (i) (g (i) ) with µ (i) iid ∼ p µ and g (i) iid ∼ p bc , then (21) is an unbiased estimator of the expected relative projection error

E := E µ∼pµ,g∼p bc T µ (g) -Π Z T µ (g) T µ (g) . ( 22
)
Note that the error indicator provides a measure of the performance of Z for the particular choice of the sampling distribution.

Random boundary conditions

The oversampling domain U in fig. 2 contains N co dd subdomains (cf. fig. 2): in absence of prior information, we propose to set µ = µ (1) , . . . , µ (N co dd ) , i , µ (i) iid ∼ Uniform P ,

Pr (i = t) = psrc N co dd t = 1, . . . , N co dd 1 -p src t = 0 ( 23
)
where p src is the probability that a source term is present in the patch. If N dd is known a priori, we might set p s = N co dd N dd . In this work, however, we consider p src = 0.5. In view of the definition of p bc , we introduce the curvilinear coordinate s ∈ [0, 1] (cf. fig. 2); then, given N f ∈ N and α ∈ R + , we define the complex-valued random field g such that

g(s; c re , c im ) = N f -1 k=0 c re k+1 + ic im k+1 1 + (2πk) 2α e 2πksi , c re k , c im k iid ∼ N(0, 1). ( 24
)
Algorithm 2 Random sample generator of boundary conditions Inputs: N f , α (cf. ( 24)), ūmax ∈ (0, 1].

Output: g : [0, 1] → [0, 1) boundary condition.

1: Draw c re , c im ∈ R N f s.t. c re k , c im k iid ∼ N(0, 1). 2: Draw X 1 , X 2 iid ∼ Uniform(0, ūmax ), set a = min{X 1 , X 2 }, b = max{X 1 , X 2 }.
3: Set g (1) = Real g(•; c re , c im ) .

4: g = a +

b-a max g (1) -min g (1) g (1) -min g (1) .

Recalling that for any k, k = 0, . . . , N f -1 and α ∈ N, we have

1 0 e 2πksi e -2πk si ds = δ k,k , d α ds α e 2πksi = (2πki) α e 2πksi ; (2πki) α (-2πk i) α = (4π 2 kk ) α ; we find that g(•; c re , c im ) 2 H α (0,1) = g(•; c re , c im ) 2 L 2 (0,1) + g (α) (•; c re , c im ) 2 L 2 (0,1) = N f -1 k=0 c re k+1 + ic im k+1 c re k+1 -ic im k+1 = N f k=1 (c re k ) 2 + c im k 2 .
The latter implies that the random variable X := g(•; c re , c im ) 2 H α (0,1) is distributed as a χ 2 distribution with 2N f degrees of freedom: therefore, the parameter α in ( 24) controls (in a probabilistic sense) the Sobolev regularity of the datum g.

In order to choose p bc , exploiting a physical argument -the solution u µ to (2c) represents water saturation -we anticipate that u µ ∈ [0, ūmax ] for some ūmax < 1. Furthermore, we wish to devise samplers that reflect the Sobolev regularity of the datum g. For these reasons, we propose to consider the procedure in algorithm 2 to generate random samples of the boundary condition. We first generate a sample of the random field g in [START_REF] Efendiev | Multiscale finite element methods for porous media flows and their applications[END_REF] and we extract its real part (cf. Line 3). then, we simply rescale the datum to ensure that the image of g, Im[g], is contained in [0, ūmax ]. The strategy in algorithm 2 is not well-suited for sampling of boundary conditions for the corner and edge components due to the presence of Dirichlet boundaries; we postpone the description of the full sampling strategy to appendix B.

In the numerical experiments, we provide samples of the boundary conditions for various values of α ∈ R + and we investigate performance for the model problem considered. In particular, we discuss the impact of the choice of α. Note that the sampling strategy proposed in this section depends on several parameters -N f , α, ūmax in algorithm 2 and p src in [START_REF] Drineas | RandNLA: Randomized Numerical Linear Algebra[END_REF] -that might be difficult to tune. This observation justifies the use of few global reduced solves at training stage to improve performance of the CB-ROM.

Basis enrichment based on reduced global solves

In several contexts, it is possible to identify at the training stage a class of global configurations of interest. To provide a concrete reference for the model problem of section 2, we might be interested in solving the global PDE for (i) any choice of n dd ∈ {n dd,LB , . . . , n dd,UB } with n dd,LB , n dd,UB ∈ N, (ii) any µ (i) ∈ P, (iii) up to n src distinct sources. The aim of this section is to devise a localized training procedure with adaptive global enrichment that exploits prior knowledge about the global system to enrich the local spaces. In section 5.1, we present a residual-based error estimator that will be used to drive the enrichment strategy; in section 5.2, we present the training procedure; in section 5.3, we present an a priori convergence result for linear coercive problems. As in section 4, we assume that the system is described by a single archetype component to shorten notation.

Residual-based error estimation

Exploiting notation introduced in section 3.3, given i ∈ {1, . . . , N dd }, and u ∈ H 1 (Ω), we define the local Riesz elements ψ µ [u] ∈ X i,0 as

(ψ µ [u], v) 1,ωi = ωi η (i) µ (x; u, v) dx, ∀ v ∈ X i,0 , (25a)
and the dual residual

r (i) µ [u] := ψ µ [u] 1,ωi . (25b)
Next, lemma 5.1 provides an upper bound for the global dual residual in terms of the localized dual residuals {r (i) [•]} i . An analogous result for linear elliptic problems was proved in [START_REF] Buhr | ArbiLoMod, a simulation technique designed for arbitrary local modifications[END_REF]Proposition 5.1]. The proof of lemma 5.1 can be found in appendix A.

Lemma 5.1. Let {φ i } i be a PoU that satisfies [START_REF] Babuška | The partition of unity method[END_REF]. Then, given u ∈ X pum , we have

G µ (u, •) -1,Ω ≤ √ M max i=1,...,N dd C r i N dd i=1 r (i) µ [u] 2 =: R µ [u], (26)
with C r i := max{C i + C 2 i + 1, 2}. We will employ the local residuals (25b) to mark instantiated components of the partitions where the error is large; see section 5.2. Let us also note that as the inifinite-dimensional analogon of G µ as a map from H 1 0 (Ω) to H -1 (Ω) is not in C 1 , one cannot expect that the • -1,Ω -norm of the residual (see eq. ( 12) for the definitions) stays bounded if the mesh size goes to zero. As a remedy one may consider G µ as a mapping from W 1,p 0 (Ω) to W -1,p (Ω), p > 2; see [START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF][START_REF] Pousin | Consistency, stability, a priori and a posteriori errors for Petrov-Galerkin methods applied to nonlinear problems[END_REF][START_REF] Smetana | Hierarchical model reduction of nonlinear partial differential equations based on the adaptive empirical projection method and reduced basis techniques[END_REF] and appendix A. As this significantly complicates the calculations of the dual norms, we opt here for assuming that the dimension of the HF space is fixed and considering the • -1,Ω -norm. We may then define the error indicator

∆ µ = N dd,µ i=1 r i µ 2 . ( 27
)
For linear problems it is straightforward to derive a rigorous a posteriori bound based on R µ [•] (see e.g. [START_REF] Babuška | The partition of unity method[END_REF][START_REF] Buhr | Localized model reduction for parameterized problems[END_REF]).

Here, we combine lemma 5.1 with the Brezzi-Rappaz-Raviart (BRR) theory [START_REF] Brezzi | Finite-dimensional approximation of nonlinear problems. I. Branches of nonsingular solutions[END_REF][START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF] to derive a rigorous residualbased error bound for the global error; see in particular [START_REF] Canuto | A posteriori error analysis of the reduced basis method for nonaffine parametrized nonlinear PDEs[END_REF][START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: rigorous reduced-basis a posteriori error bounds[END_REF] for the application of the BRR theory in the context of model order reduction. To that end, we require that

0 < β 2,p := inf w∈Xpum |w|1,Ω =0 sup v∈Xpum |v|1,Ω =0 G µ ( u µ )w, v |w| 1,Ω |v| 1,Ω , (28)
and that there exist constants γ 2,p and L 2,p such that

G µ ( u µ )w, v ≤ γ 2,p |w| 1,Ω |v| 1,Ω , (29)
G µ ( u µ ) -G µ (w) ≤ L 2,p | u µ -w | W 1,p (Ω) (30)
for w ∈ B( u µ , R) ⊂ X pum and v ∈ X pum . Here, R is supposed to be sufficiently large and |w| 1,Ω := ∇w L 2 (Ω) and |w| W 1,p (Ω) := ∇w L p (Ω) . To obtain a proximity indicator [START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: rigorous reduced-basis a posteriori error bounds[END_REF][START_REF] Canuto | A posteriori error analysis of the reduced basis method for nonaffine parametrized nonlinear PDEs[END_REF], which is based on localized and easily computable residuals via the Riesz representation, we employ, as in [START_REF] Smetana | Hierarchical model reduction of nonlinear partial differential equations based on the adaptive empirical projection method and reduced basis techniques[END_REF], the finite dimensionality of X pum and define c h := sup v∈Xpum

|v| W 1,p (Ω) |v|1,Ω and τ µ,p := 2L 2,p c h β 2 2,p √ M max i=1,...,N dd C r i N dd i=1 r (i) µ [u] 2 . ( 31
)
The proximity indicator τ µ,p will be used to validate whether u µ is close enough to u µ within the adaptive algorithm 3. We obtain the following result, which is proved in appendix A.

Proposition 5.1 (Global a posteriori error bound). Let τ µ,p < 1 and (28), ( 29) and (30) be fulfilled. Then there exists a unique solution u µ ∈ B( u µ , β2,p L2,pc h ) ⊂ X pum of (13) and the error estimator

∆ µ,p := β 2,p L 2,p c h (1 -1 -τ µ,p ) (32)
satisfies | u µ -u µ | 1,Ω ≤ ∆ µ,p . (33)
Remark 5.1 (Discussion of result). It is well-known that for nonlinear PDEs the dual norm of the residual can only be used as an a posteriori error estimator if the approximation is already close to the high-fidelity solution (see e.g. [START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF][START_REF] Verfürth | A posteriori error estimates for nonlinear problems. Finite element discretizations of elliptic equations[END_REF][START_REF] Ortner | A posteriori existence in numerical computations[END_REF]). Relying solely on the dual norm of the residual can therefore be problematic as it may seem that the approximation error is acceptable even though that might not be the case. The proximity indicator τ µ,p eq. ( 31), which only relies on computable constants, can be used to assess, whether indeed the approximation u µ is close enough to u µ such that the error estimation eq. ( 33) is valid. While the proximity indicator τ µ,p eq. ( 31) and thus the a posteriori error estimator eq. (32) solely rely on the dual norms of local residuals that can be computed on the components and therefore do not require any global solutions, the constants L 2,p and β 2,p are global constants. We will discuss some strategies on how to estimate these constants in the next remark 5.2.

To the best of our knowledge even for linear elliptic PDEs there are no results in the conforming setting that solely rely on local constants (the a posteriori error estimators in [START_REF] Buhr | ArbiLoMod, a simulation technique designed for arbitrary local modifications[END_REF][START_REF] Smetana | A new certification framework for the port reduced static condensation reduced basis element method[END_REF] e.g. both contain the global coercivity constant). A fully localizable a posteriori error estimator for nonlinear non-monotone PDEs would therefore be at least a paper on its own and is thus beyond the scope of this paper.

Remark 5.2 (Estimation of constants). Regarding the estimation of the constant c h in the inverse inequality, we refer to classical results e.g. in [START_REF] Ern | Theory and Practice of Finite Elements[END_REF] noting that the global inverse inequality only requires the measure of Ω.

Estimating the constant L 2,p relies on estimates of the constant in the Poincaré inequality for L p , W ). The estimation of c E can be easily localized. An estimate of the constant in the Poincaré inequality involving the measure of Ω can be found in [29, (7.44)] for functions that are zero on ∂Ω. We hope that if the local reduced bases contain the constant function it is maybe possible to obtain localized and more precise estimates of the Poincaré constant. Finally, similar as in [START_REF] Smetana | A new certification framework for the port reduced static condensation reduced basis element method[END_REF], we propose to use a localized model order reduction approximation of β 2,p . In detail, we suggest using the following heuristic and hierarchical estimator

β app 2,p := inf w∈ Z gfem |w|1,Ω =0 sup v∈ Z gfem |v|1,Ω =0 G µ ( u µ )w, v |w| 1,Ω |v| 1,Ω ,
where Z gfem Z gfem ⊂ X pum . We conjecture that using a certain number of additional local basis functions per component might already yield an acceptable estimate of β 2,p .

Adaptive algorithm

We define the pdfs p µ , p bc for localized sampling and the pdf p glo µ that is used to generate global problems. In the numerical examples, we consider n dd ∼ Uniform({4, . . . , 12}), µ (i) iid ∼ Uniform( P) and we assume that exactly one source term is active in Ω (that is, n src = 1). Given the partition {ω i } i , we define the local solution operators

T (i) µ : X i → X i,0 s.t. G µ u + T (i) µ (u), v = 0 ∀ v ∈ X i,0 , (34)
for i = 1, . . . , N dd,µ . The particular choice of the operators {T

(i)
µ } i is motivated by the convergence analysis in section 5.3.

Algorithm 3 contains the data compression procedure. First, we initialize the local spaces using algorithm 1.

Then, we sample n glo train configurations P train = {µ (k) } n glo train k=1 with µ (k) iid ∼ p glo µ , and we proceed with the enrichment iterations. At the -th iteration, for each µ ∈ P train , we resort to the CB-ROM proposed in section 3 to estimate the solution u µ ; then, we compute the local residuals (25)

r i µ = r (i) µ [ u µ ], i = 1, . . . , N dd,µ ,
and we mark the m r % instantiated components with the largest residual, I µ mark ⊂ {1, . . . , N dd,µ } . Then, we solve [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF] to obtain

u i,µ = 1 φ i T (i) µ ( u µ | ωi ), ∀ i ∈ I µ mark ,
and we update the dataset of simulations associated with the marked elements D,

D = D ∪ {u i,µ • Φ i : i ∈ I µ mark }.
Note that u i,µ is not well-defined on ∂ω i (i.e., division of 0 by 0): however, since we are ultimately interested in the PUM space Z gfem [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF] and due to the choice of the local norm • (cf. [START_REF] Brezzi | Finite-dimensional approximation of nonlinear problems. I. Branches of nonsingular solutions[END_REF]), this issue does not affect our procedure. In view of the termination condition, we further compute the error estimator ∆ µ,p eq. ( 32) with approximate constants. At the end of the loop over the parameters, we update the local space using POD (cf. section 4.2)

Z = Z ∪ POD {w -Π Z w : w ∈ D}, (•, •), n glo ,
and we check if max µ∈Ptrain ∆ µ,p is below a user-defined tolerance. Compute u µ using the PUM-CB-ROM (cf. section 3).

6:

Compute local residuals [START_REF] Eftang | Port reduction in parametrized component static condensation: approximation and a posteriori error estimation[END_REF]

r i µ = r (i) µ [ u µ ] for i = 1, . . . , N dd,µ . 7:
Mark the m r % instantiated components with the largest residuals, {ω i } i∈I µ mark .

8:

Solve the local problems [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF] in

{ω i } i∈I µ mark , u i,µ = 1 φi T (i) µ ( u µ | ωi ).

9:

Augment the dataset

D = D ∪ {u i,µ • Φ i : i ∈ I µ mark }. 10:
Compute ∆ µ,p eq. ( 32) with approximate constants.

11:

end for

12:

Update the POD space Z = Z ∪ POD {w -Π Z w : w ∈ D}, (•, •), n glo .

13:

if max µ∈Ptrain ∆ µ,p < tol then end if 16: end for Several steps of the Algorithm are embarrassingly parallelizable: the loop over the configurations (cf. Lines 4 to 11), the computation of the residuals (cf. Line 6), the solution to the local problems (cf. Line 8). Note also that the solution to [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF] is performed over the domain ω i (or equivalently Ω ⊂ U ), and the Newton solver can be initialized with the null solution: for this reason, it is significantly cheaper than the solution to (20a). As discussed in the introduction, the enrichment algorithm is closely linked to the online enrichment strategy proposed in [START_REF] Ohlberger | Error control for the localized reduced basis multi-scale method with adaptive on-line enrichment[END_REF] and also to related approaches in the multiscale FE literature [START_REF] Chung | Fast online generalized multiscale finite element method using constraint energy minimization[END_REF].

A priori convergence analysis for coercive linear problems

We study the in-sample a priori convergence of the enrichment procedure in algorithm 3 : we consider the case of linear coercive problems, and we apply the simplified randomized procedure contained in algorithm 4; the proof follows the argument of [START_REF] Buhr | Exponential convergence of online enrichment in localized reduced basis methods[END_REF]Theorem 1]. We assume P train = {µ} and we omit dependency on µ; in A, we discuss the extension to multiple configurations. Then, we define the model problem:

find u ∈ X : G(u, v) = f (v) -a(u, v) = 0 ∀ v ∈ X, (35)
where

H 1 0 (Ω) ⊂ X ⊂ H 1 (Ω)
is a suitable Hilbert space on Ω. We also introduce the energy norm and the associated dual norm:

w a = a(w, w) ∀ w ∈ X, f X = sup v∈X f (v) v a ∀ f ∈ X . (36)
Given the partition {ω i } N dd i=1 , we further define the associated mappings {Φ i } N dd i=1 , the associated PoU {φ i } N dd i=1 , and the local spaces

X i = H 1 (ω i ) ∩ X and X i,0 = H 1 0 (ω i ).
Then, we define the local dual residual norms such that

r (i) [u] = sup v∈Xi,0 G(u, v) v a , i = 1, . . . , N dd . (37)
Algorithm 4 simplified randomized localized training with global enrichment

1: Initialize Z = Z 0 . 2:
Sample n glo train = 1 configurations µ ∼ p glo µ , P train := {µ} k 3: for = 0, . . . , maxit do 4:

Compute u using the PUM-CB-ROM (cf. section 3).

5:

Find k = arg max i=1,...,N dd r (i) [ u ].

6:

Solve the local problem: find u loc k ∈ X k,0 such that G( u + u loc k , v) = 0 for all v ∈ X k,0 . Finally, we denote by c pu the constant such that (see [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]):

G µ (u, •) X ≤ c pu N dd i=1 r (i) [u] 2 . ( 38
)
Proposition 5.2 shows that the reconstruction error decreases exponentially with respect to the iteration count for any choice of the initial reduced space.

Proposition 5.2. The sequence of PUM-CB-ROM solutions { u } =1,2,... satisfies u-u a ≤ 1 -1 N dd c 2 pu /2 u- u 0 a .
Next Lemma summarizes two standard results that will be used in the proof of proposition 5.2.

Lemma 5.2. Let Z gfem ⊂ X and let u ∈ Z gfem satisfy G( u, v) = 0 ∀ v ∈ Z gfem . Then, we have u -u a = inf ϕ∈Z gfem ϕ -u a ; (39a) u -u a = G(u, •) X . (39b)
Proof. (proposition 5.2). Exploiting [START_REF] Huynh | A static condensation reduced basis element method: approximation and a posteriori error estimation[END_REF] and then (39b), we find

r (k) [ u ] 2 ≥ 1 N dd N dd j=1 r (j) [ u ] 2 ≥ 1 N dd c 2 pu G( u , •) 2 X = 1 N dd c 2 pu u -u 2 a . (40)
By construction, u loc k in algorithm 4 belongs to

Z ( +1) gfem . As a result, if we consider ϕ = u + u loc k in (39a), we find u -u +1 2 a ≤ u -u -u loc k 2 a = u -u 2 a -2a u -u , u loc k + u k 2 a . Since a u -u , u loc k = G u + u loc k , u loc k =0 + u loc k 2 a = u loc k 2 a , and
u loc k a = sup v∈X k,0 a(u loc k , v) v a = sup v∈X k,0 G(u -u , v) v a = r (k) [ u ],
we obtain

u -u +1 2 a ≤ u -u 2 a -r (k) [ u ] 2 . (41)
Thesis follows by combining ( 40) and (41).

6 Numerical results

Performance of randomized training for a linear problem

We first provide numerical investigations for the linear advection-diffusion-reaction problem on the random field

-∇ • µ 1 κ∇u µ,g + [µ 2 , µ 3 ] T u µ,g + µ 4 u µ,g = 0 in U = (0, 0.3) 2 , u µ,g = g on ∂U =: Γ in ,
g • (x; c) := i∈I • dir f(c i , ūmax )φ fe,• i (x), c i iid ∼ N ūmax 2 , ū2 max 4 , f(c, u) = max{min{c, u}, 0}, (44)
where {I • dir } • denote the set of indices of the mesh on the patch input boundaries and {φ fe,• i } • are the Lagrangian bases associated with the high-fidelity discretization. We refer to the sampling procedure in algorithm 6 as smooth sampling; we refer to the sampling procedure based on (44) as Gaussian sampling.

We compute n test = 30 global solutions for n dd = 10 (N dd = 100) components; then, we define the test datasets {D • } •∈{co,ed,int} by extracting the solution in each element of V -card(D • ) = 1920 (resp., 120, 960) for the internal (resp., corner, edge) component. Finally, we introduce the localized error indicators

E • avg,rel (Z • ) = 1 card(D • ) w∈D • w -Π Z • w • w • , • ∈ {co, ed, int}, (45)
which are used to assess performance. Figure 6 shows random samples of the boundary conditions on Γ in for internal and edge components as provided by algorithm 6 for various choices of α and N f = 20 and ūmax = 0.5. As for the linear case, the value of α encodes the spatial smoothness of the samples. We further observe that algorithm 6 automatically enforces the proper condition at the extrema s = 0 and s = 1 -g(0) = g(1) = 0 for • ∈ {co, ed}, g(0) = g(1) for • = int.

Figure 7 shows the behavior of the relative errors [START_REF] Mahoney | CUR matrix decompositions for improved data analysis[END_REF] for the three components, for smooth sampling for three choices of α (N f = 20, ūmax = 0.5), and for Gaussian sampling [START_REF] Mahoney | Randomized algorithms for matrices and data[END_REF]. To provide a reference, we show also performance of the POD spaces based on the datasets {D • test } •∈{co,ed,int ("opt") generated using 30 additional global simulations with N dd = 100 components. We observe that smooth sampling outperforms Gaussian sampling for the boundary components: we believe that this is due to the presence of strong Dirichlet conditions on ∂ U • \ Γ • in . We further observe that results weakly depend on the choice of α. Figure 8 shows the performance of the CB-ROM based on PUM. In fig. 8(a), we show the average global L 2 and H 1 relative errors over the test set of n test = 30 global simulations and we also compare with the projection error. We here consider ūmax = 0.5, N f = 20, and α = 1. We observe that Galerkin projection is nearly optimal for all choices of n; we further observe exponential convergence with respect to n. In fig. 8(b), we compare the H 1 relative projection error for ūmax = 0.5, N f = 20, and α = 1 with the results obtained for ūmax = 0.75, N f = 20, and α = 4: we observe that results weakly depend on the choice of these two hyper-parameters.

Adaptive enrichment

We apply algorithm 3 with error indicator ∆ µ eq. ( 27) to the model problem of section 6.2. We consider n loc train = 30, n loc = 20, we set {p • µ } • as discussed in section 4.2 and we generate random samples of boundary conditions at input ports based on (i) algorithm 6 with N f = 20, ūmax = 0.5, α = 1 or (ii) on iid realizations of [START_REF] Mahoney | Randomized algorithms for matrices and data[END_REF]. We further consider n glo train = 50, n glo = 10, maxit= 3, and we generate global configurations using the strategy outlined in section 5.2. We assess performance based on n test = 20 out-of-sample randomly-chosen configurations.

Figure 9(a) and (b) show boxplots of the relative H 1 error after each iteration of the training algorithmiteration 0 corresponds to the performance of the CB-ROM without global enrichment. Iteration 0 corresponds to a reduced space of size n = 20; iterations it = 1, 2, 3 correspond to reduced spaces of size n = 20 + 10 • it. We observe that the enrichment iterations significantly improve performance of the CB-ROM and reduce the impact of the initial sampling distribution. (27) and the relative H 1 error on the test set for all iterations of the enrichment algorithm for smooth sampling; Figure 9(d) shows the effectivity of the error indicator η = ∆ µ /E rel for smooth sampling. We observe that the residual indicator is strongly correlated with the global error.

Conclusions and perspectives

We presented a CB-pMOR method for parameterized elliptic nonlinear PDEs. The approach relies on the definition of several archetype components and associated local ROBs and ROMs. CB-pMOR rely on two building blocks: (i) a localized training strategy for the construction of the local approximation spaces, and (ii) a DD strategy for online global predictions. In this paper, we proposed a localized data compression procedure based on oversampling and randomized sampling of boundary conditions of controlled smoothness, and we relied on the PUM to devise global approximation spaces and on Galerkin projection to estimate the global state. Finally, we proposed an adaptive enrichment procedure that exploits global CB-ROM solves to improve approximation properties of the local reduced spaces. Numerical results for a nonlinear diffusion problem show the impact of the sampling distribution on performance: given a class of nonlinear PDEs, it is thus necessary to devise an effective sampler that is informed by the problem of interest. The approach presented in this work (cf. algorithm 2) is simple to implement, and incorporates relevant features of the problem of interest -lower and upper bounds for the solution, Sobolev regularity, Dirichlet boundary conditions. However, it depends on several hyper-parameters that might be difficult to set a priori. In this respect, we numerically showed that the proposed enrichment strategy reduces the impact of the initial sampling distribution.

In the future, we wish to extend the approach in several directions. First, we wish to devise specialized hyper-reduction strategies for CB-pMOR methods based on PUM: hyper-reduction is key to reduce efficient online memory and computational costs. Second, we wish to develop rigorous a posteriori error estimators for nonlinear PDEs for online certification. Third, we wish to analyze performance of randomized algorithms for nonlinear operators: this analysis is key to provide mathematical foundations for randomized methods for nonlinear problems and also to inform the choice of the sampling distribution. A Proofs Lemma 5.1

Proof. Since the form G µ is linear with respect to the second argument, we find that

G µ (u, v) = N dd i=1 G µ (u, v i ) = N dd i=1 ωi η (i) µ (x; u, v i ) dx = N dd i=1 (ψ µ [u], v i ) 1,ωi ≤ N dd i=1 r (i) µ [u] v i 1,ωi , (46a)
where v i = vφ i . Note that in the second identity we used (25a), while in the last bound we applied Cauchy-Schwartz inequality. We observe that ωi

(vφ i ) 2 dx ≤ ωi v 2 dx and ωi ∇(vφ ) 2 2 dx = ωi v∇φ i + φ i ∇v 2 2 dx = ωi φ 2 i ∇v 2 2 + v ∇φ i 2 2 +2vφ i ∇v • ∇φ i dx ≤ ωi ∇v 2 2 + C 2 i v 2 + C i v 2 + C i ∇v 2 2 dx
The latter two bounds imply

v i 2 1,ωi ≤ max{C i + C 2 i + 1, 2} v 2 1,ωi . (46b)
Substituting (46b) in (46a), we obtain

G µ (u, v) ≤ N dd i=1 C r i r (i) µ [u] v 1,ωi .
Then, apply Cauchy-Schwartz inequality and the estimate in [5, Lemma 2], we obtain

G µ (u, v) ≤ max i=1,...,N dd C r i N dd i=1 r (i) µ [u] 2 N dd i=1 v 2 1,ωi ≤ √ M max i=1,...,N dd C r i N dd i=1 r (i) µ [u] 2 v H 1 (Ω) ,
which proves [START_REF] Ern | Theory and Practice of Finite Elements[END_REF].

A posteriori error estimation

As noted in section 5.1, the infinite-dimensional analogon of G µ as a map from H 1 0 (Ω) to H -1 (Ω) is not in C 1 ; this can be easily seen by determining the Fréchet derivative of G µ . One should thus consider G µ as a mapping from W 1,p 0 (Ω) to W -1,p (Ω), p > 2 which yields a C 1 -mapping. Note that the latter is crucial to apply the Brezzi-Rappaz-Raviart (BRR) theory [START_REF] Brezzi | Finite-dimensional approximation of nonlinear problems. I. Branches of nonsingular solutions[END_REF][START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF]. One may then derive a rigorous a posteriori error bound for the error in the | • | W 1,p (Ω) -norm (see e.g. [START_REF] Pousin | Consistency, stability, a priori and a posteriori errors for Petrov-Galerkin methods applied to nonlinear problems[END_REF][START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF][START_REF] Smetana | Hierarchical model reduction of nonlinear partial differential equations based on the adaptive empirical projection method and reduced basis techniques[END_REF]). Unfortunately, the error estimator relies on the W -1,p -norm, which is challenging to estimate even in the finite-dimensional setting. As in [START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF] we thus exploit that for all z ∈ B(u pum , R), G µ (z) : W 1,p (Ω) → W -1,p (Ω) can be continuously extended as an operator in L(H 1 (Ω), H -1 (Ω)). Furthermore, we require that eq. ( 28)-eq. ( 30) are satisfied. We can then prove proposition 5.1:

Proof. Thanks to lemma 5.1 the assumption τ µ,p < 1 implies that τ µ,p :=

2L2,pc h β 2 2,p G µ ( u µ , •) -1,Ω < 1.
The existence of a unique solution u µ ∈ B( u µ , β2,p L2,pc h ) of eq. ( 13) and

| u µ -u µ | 1,Ω ≤ β 2,p L 2,p c h (1 -1 -τ µ,p ) (47)
then follows using the standard arguments in the BRR theory (see [START_REF] Canuto | A posteriori error analysis of the reduced basis method for nonaffine parametrized nonlinear PDEs[END_REF][START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: rigorous reduced-basis a posteriori error bounds[END_REF][START_REF] Caloz | Numerical analysis for nonlinear and bifurcation problems[END_REF][START_REF] Pousin | Consistency, stability, a priori and a posteriori errors for Petrov-Galerkin methods applied to nonlinear problems[END_REF] and for this particular PDE [START_REF] Smetana | A dimensional reduction approach based on the application of Reduced Basis Methods in the context of Hierarchical Model Reduction[END_REF]). As the function t(x) := 1 -√ 1 -x is strictly increasing on (0, 1), applying lemma 5.1 to the right side of eq. ( 47) concludes the proof.

Extension of proposition 5.2 to multiple configurations

In Algorithm 5, we generalize algorithm 4 to the case of multiple configurations and we discuss the proof of a priori exponential convergence. Note that algorithm 5 can still be interpreted as a simplified version of algorithm 3 that is more amenable for the analysis. We assume here that N dd,µ ≤ N dd,max with probability one for some N dd,max . All constants introduced below depend on N dd,max and also on the size of the training set n glo train .

Algorithm 5 simplified randomized localized training with global enrichment Compute u ,µ using the PUM-CB-ROM (cf. section 3) for µ ∈ P train .

5:

Find (k, µ ) = arg max i=1,...,N dd,µ ,µ∈Ptrain r (i) µ [ u ,µ ].

6:

Solve the local problem: find u loc k ∈ X k,0 such that G µ ( u ,µ + u loc k , v) = 0 for all v ∈ X k,0 .

B Extension to multiple components

As discussed in the main body of the paper, in our numerical experiment, we consider three components. We should thus construct the local approximation spaces Z co , Z ed , Z int , with Z • ⊂ Y • , such that min

ζ∈Z L i u µ ωi -ζ • Φ -1 i
1,ωi ≤ ε tol for i = 1, . . . , N dd , µ ∈ P(n dd ), [START_REF] Martinsson | A randomized algorithm for the decomposition of matrices[END_REF] Condition [START_REF] Martinsson | A randomized algorithm for the decomposition of matrices[END_REF] implies that the local spaces Z co , Z ed , Z int should approximate the manifolds

             M int = u µ ωi • Φ i : µ ∈ P glo (n dd ), L i = int, n dd ∈ N , M ed = u µ ωi • Φ i : µ ∈ P glo (n dd ), L i = ed, n dd ∈ N , M co = u µ ωi • Φ i : µ ∈ P glo (n dd ), L i = co, n dd ∈ N . (49)
For the three components in our library, we consider the oversampling domains depicted in Figure 10. Note that U • comprises N

Random boundary conditions

In algorithm 6, we discuss the complete random sample generator of boundary conditions considered in the numerical experiments. Note that, if • ∈ {co, ed}, since by construction d k ds k g (1) (0) = d k ds k g (1) (1) for k ∈ N, we define g (2) (s) = g (1) (cs) with c = 0.72 (cf. Line 7); then, we enforce that g(s) = 0 for s ∈ {0, 1} and g ≥ 0 (cf. Line 8); finally, in Line 9, we ensure that Im[g] ⊂ [0, ūmax ].

Adaptive enrichment

Algorithm 7 is the generalization of algorithm 3 to the case of multiple archetype components. Note that we here enforce that the dimension of the local spaces is the same for all components to simplify the implementation of the CB-ROM; however, the Algorithm can be modified to allow local spaces of different size. b-a max g (1) -min g (1) g (1) -min g (1) . 6: else 7:

g (2) (s) = g (1) (0.7s).

8:

g (3) (s) = a + b-a max g (2) -min g (2) g (2) -min g (2) s(1 -s). 9: g = X3 max g (3) g (3) .

Figure 1 :

 1 Figure 1: nonlinear diffusion. (a) domain Ω and subdomains {Ω i } i for N dd = 100. (b) instantiated archetype components. (c)-(d) solution fields for N dd = 100 and two choices of the parameters.

Figure 2 :

 2 Figure 2: Nonlinear diffusion. Archetype component with corresponding oversampling domain U .

Algorithm 3 1 :

 31 randomized localized training with global enrichment Inputs (localized training): n loc train = number of solves, n loc = size of the POD spaces, pµ, p bc sampling pdfs. Inputs (enrichment): n glo train = number of global simulations per iteration, n glo = number of modes added at each iteration, maxit = maximum number of outer loop iterations, tol = tolerance for termination criterion, p glo µ = global configuration sampler, mr = percentage of marked components at each iteration. Outputs: Z local approximation space.Localized training Apply algorithm 1 to obtain the local space Z.

Enrichment 1 : 3 :

 13 Sample n glo train configurations µ (k) iid ∼ p glo µ , P train := {µ (k) } k 2: for = 1, . . . , maxit do Initialize the dataset D = ∅.

4 :

 4 for µ ∈ P train do 5:

7 :

 7 Define u = u loc k φ k and update the local space Z = Z ∪ span{u • Φ k }. 8: end for

Figure 4 :

 4 Figure 4: linear problem. Out-of-sample performance; comparison with deterministic training for N r = 100 choices of the random samples and for two fixed test sets. (a) -(b) smooth test set. (c) -(d) Gaussian test set.

Figure 5 :

 5 Figure 5: linear problem. Effectivity of the error indicator η = E(ntest=10) E(ntest=100) for 100 independent runs and for both Gaussian and smooth training.

= 3 Figure 6 :

 36 Figure 6: nonlinear problem. Samples of random boundary conditions for three choices of α (N f = 20, ūmax = 0.5). (a)-(b)-(c) internal component. (d)-(e)-(f) edge component.

  Figure 9(c) shows the correlation between the residual indicator

Figure 7 :

 7 Figure 7: nonlinear problem. Local approximation errors (45) for three choices of α (N f = 20, ūmax = 0.5), and for Gaussian sampling (44).

Figure 8 :

 8 Figure 8: nonlinear problem. Performance of PUM CB-ROM on n test = 30 global solutions for N dd = 100. (a) Galerkin error vs projection error. (b) projection error for two choices of the parameters in algorithm 6.

Figure 9 :

 9 Figure 9: nonlinear problem; adaptive enrichment. (a)-(b) boxplots of the relative H 1 error on the test set for smooth and Gaussian sampling of localized BCs. (b) correlation between ∆ µ and relative H 1 error (smooth sampling). (c) out-of-sample effectivity of the error indicator ∆ µ /E rel (27) (smooth sampling).

1 : 2 :

 12 Initialize Z = Z 0 . Sample n glo train configurations µ k iid ∼ p glo µ , P train := {µ k } k 3: for = 0, . . . , maxit do 4:

2 ≤

 2 update the local space Z = Z ∪ span{u • Φ k }. 8: end for Given µ ∈ P train , we denote by { µ j } j ⊂ {1, . . . , maxit} the set of indices corresponding to the iterations at which we select µ = µ at Line 5 of algorithm 5. Then, exploiting proposition 5.2, we findu µ j+1 ,µ -u µ aµ ≤ 1 -1 N dd,µ c 2 pu,µ j/2 u 0,µ -u µ aµ ≤ Ce -α n glo train jwhere C, α > 0 are constants that do not depend on the iteration count j.Given ∈ N, there exists µ ∈ P train such that #{ µ j : µ j ≤ } ≥ n glo train -after iterations of algorithm 5, the parameter µ is selected more than n glo train times. This implies that max i=1,...,N dd,µ ,µ∈Ptrain r (i) µ [ u ,µ ] ≤ C e -α and thus max µ∈Ptrain u ,µ -u µ aµ = max µ∈Ptrain G µ ( u ,µ , •) X µ ≤ max µ∈Ptrain c C N dd,µ max µ∈Ptrain c pu,µ e -α .

• dd = 9 (Figure 10 :

 910 Figure 10: nonlinear diffusion. Archetype components with corresponding oversampling domain.

1 :

 1 Inputs: N f , α (cf. (24)), ūmax ∈ (0, 1], • ∈ {co, ed, int}.Output: g : [0, 1] → [0, 1) boundary condition. Draw c re , c im ∈ R N f s.t. c re k , c im k iid∼ N(0, 1).

2 : 3 iid∼

 23 Draw X 1 , X 2 , X Uniform(0, ūmax ), set a = min{X 1 , X 2 }, b = max{X 1 , X 2 }.

3 :

 3 Set g(1) = Real g(•; c re , c im ) .

4: if • = int then 5 :

 5 g = a +

1 : 1 : 3 :

 113 10: end if Algorithm 7 randomized localized training with global enrichment Inputs (localized training): n loc train = number of solves, n loc = size of the POD spaces, {p • µ , p • bc }• sampling pdfs. Inputs (enrichment): n glo train = number of global simulations per iteration, n glo = number of modes added at each iteration, maxit = maximum number of outer loop iterations, tol = tolerance for termination criterion, p glo µ = global configuration sampler, mr = percentage of marked components at each iteration. Outputs: {Z • } •∈{co,ed,int} local approximation spaces. Localized training Apply algorithm 1 to obtain the local spaces {Z • } •∈{co,ed,int} . Enrichment Sample n glo train configurations µ (k) iid ∼ p glo µ , P train := {µ (k) } k 2: for = 1, . . . , maxit do Initialize the datasets D • = ∅ for • ∈ {co, ed, int}.

4 :

 4 for µ ∈ P train do 5:

  1,p and the Sobolev embedding inequality v C 0 (Ω) ≤ c E |v| W 1,p (Ω) (see e.g.[START_REF] Smetana | A dimensional reduction approach based on the application of Reduced Basis Methods in the context of Hierarchical Model Reduction[END_REF] Subsection 3.1.2]

Here, we choose n = dim(Z int ) = dim(Z ed ) = dim(Z co ). This choice simplifies notation and is also convenient for code vectorization. The extension to reduced spaces of arbitrary size is straightforward.

The choice c = 0.7 is not crucial for the methodology.

. The linear problem allows us to compare our randomized method with a previously-developed data compression algorithm. Note that the transfer operator T : (µ, g) → u µ,g | Ω is nonlinear due to the presence of parameters. We discretize the problem using the finite element method based on cubic (P3) polynomials, with N in = 360 degrees of freedom on the boundary Γ in .

We compare performance of our randomized algorithm with the approach in [START_REF] Taddei | A localization strategy for data assimilation; application to state estimation and parameter estimation[END_REF] (TE+POD): given the training set P train = {µ k } ntrain k=1 ⊂ P, we first solve n train independent transfer eigenproblems [START_REF] Babuška | Optimal local approximation spaces for generalized finite element methods with application to multiscale problems[END_REF] for each value of the parameter and then we use POD to combine the resulting spaces. We refer to [START_REF] Taddei | A localization strategy for data assimilation; application to state estimation and parameter estimation[END_REF] for further details and analysis, and we refer to [START_REF] Smetana | Optimal local approximation spaces for component-based static condensation procedures[END_REF] for a similar data compression algorithm. In the numerical experiments, we set n train = 100: this implies that TE+POD requires to solve n train • N in = 36000 PDEs. We envision that the total number of PDE solves can be reduced up to O(n • n train ) by resorting to Krylov methods to solve the transfer eigenproblem: we refer to the above-mentioned literature for further details.

We set p µ = Uniform(P) and we consider samples of the random field g = Real[ g(•; c re , c im )] (cf. ( 24)) with N f = 20. In fig. 3, we show random samples of g -g avg with g avg = 1 0 g(s) ds for various choices of α: we observe that, as α increases, the samples become increasingly smooth. Given the restriction of the finite element Lagrangian basis to the input boundary {φ fe i } i∈I dir , we further define the random field

which is used below for comparison. To assess performance, we compare the maximum relative projection error

for the two choices of p bc -"smooth" (with α = 1) and "Gaussian" (42) -and n test = 100. Figure 4 shows the results for smooth and Gaussian training and test sets. Here, we consider training sets of size n train = 50 in algorithm 2; furthermore, we compare errorbar plots based on 100 independent choices of the training set. We observe that our smooth sampling strategy is nearly as effective as TE+POD for n 40 for both smooth and Gaussian test sets. This result empirically demonstrates that randomized methods are extremely effective to identify dominant POD modes even for nonlinear transfer operators. We further observe that Gaussian sampling is clearly inferior when tested on smooth data, while it performs as accurately as smooth sampling on the Gaussian test set: we conjecture that this behavior is due to the low-pass filtering properties of the differential operator.

In fig. 5, we show the behavior of the error indicator E in remark 4.1: more precisely, we show boxplots of the approximate effectivity η = E(ntest=10) E(ntest=100) for 100 independent runs and for both Gaussian and smooth training. Note that, with very high-probability, η ∈ [0.5, 1.5]. Note also, however, that E strongly depends on the choice of the sampling distribution, which in practice is largely unknown.

Application to the nonlinear diffusion problem

We consider the application to the nonlinear diffusion problem introduced in section 6.2. We apply algorithm 2 with n train = 200; we set {p • µ } • as discussed in section 4.2 and we consider the smooth sampler described in algorithm 6 for N f = 20, various choices of α and ūmax -we recall that algorithm 6 is the generalization of algorithm 2 to the case of three components. We further compare performance with randomized training based Compute u µ using the PUM-CB-ROM (cf. section 3). Compute local residuals [START_REF] Eftang | Port reduction in parametrized component static condensation: approximation and a posteriori error estimation[END_REF]

Mark the m r % instantiated components of type • with the largest residual, {V i } i∈I µ mark,• .

9:

Solve the local problems [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF] in

µ ( u µ | ωi ). Compute ∆ µ [START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF]. if max µ∈Ptrain ∆ µ < tol then