A local hybrid surrogate‐based finite element tearing interconnecting dual‐primal method for nonsmooth random partial differential equations

A domain decomposition approach for high‐dimensional random partial differential equations exploiting the localization of random parameters is presented. To obtain high efficiency, surrogate models in multielement representations in the parameter space are constructed locally when possible. The method makes use of a stochastic Galerkin finite element tearing interconnecting dual‐primal formulation of the underlying problem with localized representations of involved input random fields. Each local parameter space associated to a subdomain is explored by a subdivision into regions where either the parametric surrogate accuracy can be trusted or where instead one has to resort to Monte Carlo. A heuristic adaptive algorithm carries out a problem‐dependent hp‐refinement in a stochastic multielement sense, anisotropically enlarging the trusted surrogate region as far as possible. This results in an efficient global parameter to solution sampling scheme making use of local parametric smoothness exploration for the surrogate construction. Adequately structured problems for this scheme occur naturally when uncertainties are defined on subdomains, for example, in a multiphysics setting, or when the Karhunen–Loève expansion of a random field can be localized. The efficiency of the proposed hybrid technique is assessed with numerical benchmark problems illustrating the identification of trusted (possibly higher order) surrogate regions and nontrusted sampling regions.


F I G U R E 1 Left: Realization of composite material
A from Example 1; Right: Realization with multiple separated inclusions of Example 2 to solution map with the parameters determining the randomness. 1,2 However, in many science and engineering applications deviating from the common linear benchmark setting with smooth Karhunen-Loève random fields, 3,4 higher smoothness cannot be assumed globally in the parameter domain and physical space. This severely limits the use of techniques that rely on sparsity of a nonlinear approximation.
As a prototypical application example, which is for instance of relevance in material science, we have in mind a matrix composite material with random nonperiodic inclusions as depicted in Figure 1 (right). This setting exhibits discontinuities in the parameter dependence, rendering it basically intractable to functional approximations with global basis functions in parameter space as commonly used in stochastic Galerkin finite element method, stochastic collocation, and other nonintrusive projection approaches. [5][6][7][8][9][10][11] We note that in many applications, random composite materials are often analyzed by some type of homogenization or upscaling techniques. [12][13][14] By contrast to this, the requirements of the method developed in this article are determined by an engineering application which aims at analyzing the maximum principal stresses in a glue layer. This quantity of interest, while being of a local nature, is used as a stability measure for the underlying material and motivates the resolution of material imperfections which otherwise would get lost with homogenization or upscaling techniques.
While Monte Carlo (MC) sampling and its modern variants, for example, multilevel (MLMC) 3,15,16 and quasi-Monte Carlo 17 are widely applicable and robust with respect to stochastic dimensions, smoothness can only be exploited to a limited extend, resulting in slow convergence in the number of samples. Nevertheless, for problems with low parameter regularity, a sampling approach might seem the only option to pursue. It is a central aim of this work to overcome this restriction as (computationally) beneficial as possible. Sampling approaches have been analyzed for random material coefficients spatially being piecewise continuous in Reference 15 or essentially bounded in Reference 18. In both references, the approximation of the material coefficients (e.g., due to the discretization process) was avoided, leaving out the effect of discontinuities not being resolved perfectly in the numerical framework, see Examples 1 and 2.
The first part of this article attempts to close this gap within the random framework inspired by Bonito et al. 19 introducing a so-called p * -condition. More precisely, given higher-order spatial integrability of (weak) derivates of the solution, we are able to allow for an approximation of nonsmooth material coefficients in weaker norms such as L Q (Ω; L q (D)) with q < ∞ in contrast to the case of q = ∞ as in Reference 18. The discussion of the new theory in the context of MLMC is open and is out of the scope of this article.
A notion not frequently used for this kind of problems is a localization of randomness, for example, by means of a domain decomposition (DD) method as employed here. In fact, the guiding principle is to make the high-dimensional problem more accessible by considering smaller physical domains, locally leading to a dependence on fewer relevant random variables.
To achieve this, we rely on the well-known theory of DD techniques that are powerful for solving large-scale partial differential equations (PDEs) based on parallel computing. 20,21 A subclass of DD methods, also known as iterative substructuring methods, 20 is based on nonoverlapping subdomains which allows for two important features. First, it turns out to be flexible and efficient in handling deterministic elliptic problems with large jumps in the coefficients, [22][23][24] including extensions with adaptive coarse space techniques 25,26 or with handling irregular domains. 27 Second, it gives a natural framework for decoupled problems useful in the localization of randomness. Nonoverlapping DD methods might be classified based on the treatment of the interface solution compatibility. In primal methods such as the basic Schur complement method 20 or the balancing domain decomposition, 28 the unique unknown interface solution is computed by enforcing the continuity constraint on the subdomain interfaces. Alternatively, in dual methods such constraints might be handled via Lagrange multipliers, for example, yielding the classical finite element tearing and interconnecting (FETI) 29 class. A combination of these approaches leads to the concept of hybrid methods such as the finite element tearing interconnecting dual-primal method (FETI-DP). 25,27,30,31 The application of iterative substructuring methods in the framework of random PDEs has been presented by several authors, mainly using a Karhunen-Loève expansion (KLE) to model the material law relying on polynomial chaos expansions (PCE) in the arising approximation tasks. In References 32,33 the FETI-DP framework was applied to a random PDE based on a global (KLE) and a global PCE ansatz. This approach however is limited by the global random dimension since no localization of randomness was applied. A first localization approach was proposed in Reference 34 in the setting of the classical FETI method for two subdomains based on a separate representation in canonical tensor format (CP) with iterative updating via an alternating Rayleigh-Ritz method. From an approximation point of view the CP format is not closed 35 and may introduce instabilities in the numerical scheme. No localization of the random input field was considered. To our knowledge the combination of a localization of a KLE of a Gaussian random field with a substructuring method was first introduced in Reference 36. There, the Schur complement method using several local PCEs to approximate the local smooth parametric Schur complement matrices was applied to accelerate a global to local sampling process. A higher-order local PCE was needed to avoid losing the positive definiteness of the parametric Schur complement matrix. As those ideas inspired some aspects in our work, we give more details and discussion in Example 4, Remark 7, and Section 4.2. During the development of this article, the work 37 was published, which is based on the Schur complement method for the smooth KLE case for colored and white noise using a reduced basis method.
With our developed method, we attempt to extend the localization approach to the FETI-DP method for random fields that are nonsmooth in both the parameter and the physical space, which is the case with matrix composites with uncertain shape or position of the inclusions. Up to the author's knowledge this is the first application of FETI-DP to the setting of local nonsmooth random fields. In the proposed method we rely on several local surrogates to accelerate a global to local sampling stage. Due to the lack of regularity, the construction of surrogates of local operators is not straightforward. In particular, a classical PCE is not suitable since it would grow too large to stay feasible in practice. To overcome this severe obstacle, we introduce the concept of trust and no-trust regions for a prescribed error tolerance 38 locally in which a highly efficient surrogate can be evaluated ("trusted") or where one has to fall back to standard pointwise sampling ("nontrusted"). Local surrogates can be generated in parallel and depend only on local random coordinates. In order to get the largest gain from the construction, trust regions are required to cover as large an area as possible of the parameter domain related to each subdomain. For this, we introduce a local generalized multielement discretization in the local parameter spaces for which an hp-adaptive refinement procedure based on error indicators is presented.
The proposed hybrid local surrogate-based stochastic FETI-DP method allows for fast global sampling on the basis of local smoothness exploitation by appropriate surrogates. We would like to point out that the presented method is a single level method in terms of the physical discretization and might be used as an accelerated sampling emulator for modern sampling schemes such as MLMC.
The article is structured as follows: The following section is devoted to the theoretic framework for nonsmooth random PDEs. In particular we provide a new perturbation and convergence theory for a linear random model problem. Section 3 examines the construction of surrogate models, which are used locally in the proposed method. Special attention is paid to generalized multielement PCEs and a partition of unity (PoU) interpolation in order to treat nonsmooth dependence, leading to the notion of "trust zones." The employed parametric DD method is derived in Section 4, resulting in the central parametric FETI-DP method. Eventually, Section 5 demonstrates the accuracy and adaptive refinement behavior of the hybrid sampling approach with several benchmark problems.

MODEL PROBLEM AND A PRIORI ESTIMATES
In this section, we introduce the random linear elliptic model problem used henceforth for the derivation of the proposed DD method. Of particular interest is the effect of approximation errors of the coefficient, which is of particular relevance for simulations with random data models. We examine conditions required for stability results of the approximate solution, especially in the case of nonsmooth coefficients. A prototypical application we have in mind is the numerical treatment of stochastic composite materials. This type of problem is much more involved than the frequently considered case of random PDEs with smooth data in UQs as, for example, presented in References 4,9,10,39,40, where the dependence on the countable (possibly infinite) parameter vector is analytic. Such a smooth dependence allows for the derivation of best n-term approximations with optimal (exponential) convergence rates and, while still computationally costly, the implementation of efficient numerical methods. We point out that this section is devoted to the theoretical question of approximation of nonsmooth random fields. This should be applicable in a wider sense than to the parametric DD we have in mind. In particular, the theory of MLMC could be enhanced with the introduced aspects. For the reader mainly interested in the algorithmic aspects of the proposed hybrid sampling approach based on DD, we suggest to only skim through this section to understand the examined model problem and the structure of the considered nonsmooth problems by the examples provided at the end. In Section 2.1, we introduce the model problem and collect assumptions that are needed for the pathwise and Bochner space perturbation theory for nonsmooth random fields examined in Section 2.2. The end of this section is devoted to examples that illustrate the challenges and motivate the derived a priori theory.

Random PDE model
As model problem we consider a random linear PDE. Given a probability space (Ω,  , P), a Lipschitz domain D ⊂ R d for d = 2, 3, with Dirchlet and Neumann boundary segments Then, the stationary diffusion model equation reads pointwise for ∈ Ω a.e.
Here A denotes the diffusion field, f and g are source fields and n is the outer normal on Γ 1 w.r.t. D. In the following, we state pathwise and global assumptions in the setting of elliptic (pathwise) second-order PDEs. In the context of integrability, we use lower case letters for spatial integrability and upper case letters for stochastic integrability, for example, L p (D) or L P (Ω). Then there exist random variables c, c ∶ Ω → R such that there hold pathwise uniform bounds a.e. in , (A3) Uniform bounds: 0 < C, C < ∞ such that ess inf c( ) ≥ C and ess sup c( ) ≤ C.
Note that Assumptions (PA1), (PA2), and (A2) are standard. 18 They are needed to ensure existence and uniqueness of a pathwise solution. Assumption (A1) is the key assumption for the proposed theory relying on the p * -condition. We emphasize that no uniform bound of the underlying random fields as in (A3) is required for the theory to hold. This somewhat generalizes the class of admissible random fields but discards, for example, lognormal random fields. Note that the latter yield higher regularity in terms of parametric dependence and thus are not in the focus of this work since standard approximation theory applies. In particular, a p * -condition is not necessary. Furthermore, it turns out that Assumption (A3)-while implying Assumption ((A2))-is useful for interpolation arguments employed later in terms of the approximation of the diffusion field in weaker norms.
We define the pathwise bilinear form a associated to the weak formulation of (1) for a given diffusion field B satisfying pathwise Assumptions (PA1) and (PA2) as and the pathwise linear form Then the pathwise weak formulation given B for ∈ Ω reads Seek Existence, uniqueness and stochastic integrability of (3) is standard and summarized in the following lemma. ) pathwise being the unique weak solution u( ) of (3) with B = A for P-almost all ∈ Ω. In particular u is measureable w.r.t. to the -algebra (f , g, A) generated by f , g, and A.

Error estimate for coefficient approximate solutions
LetÂ be some perturbation of A in (1), for example, implicitly introduced by a quadrature or triangulation scheme in the discretization process. Let u andû be (pointwise) weak solutions with respect to A andÂ. Then, we are interested in the distance of u toû in the abstract sense of where the type of convergence has to be determined.

Pathwise error estimates
Let u( ) be the solution of (3) with B = A and letû( ) be the solution of (3) with B =Â. Denote byĉ ∶ Ω → R the lower bound random variable with respect toÂ satisfying Assumptions (PA1) and (PA2). For simplicity, we shall assume that the right-hand side linear form can be evaluated exactly. Then it holds Now taking v = u −û and due to the assumptions on A, the left-hand side can be estimated from below bŷ Let us assume that ∇u( ) is in L p (D) for some p ∈ [2, ∞]. Then for q = 2p/(p − 2) ∈ [2, ∞], we obtain pathwise Combining (7) and (8) yields Based on this derivation, we can refine the above statement to a decomposition of the spatial domain D with relaxed integrability requirements of ∇u in each subdomain.
For p s = 2 for s = 1, … , N SD , the perturbation result recovers the standard L ∞ estimate, w.r.t. A. We are left with the question of integrability of ∇u. Here, we state a pathwise p*-condition in the spirit of Bonito et al. 19 adapted to our random framework.
Definition 1 (p * -condition). Let the integrability constant p * = p * (A, d, D) > 2 be independent of and let 2 ≤ p < p * . In addition, if |Γ 1 | > 0 there is a conormal derivative trace space (D) and a random variable C p ( ) = C p (d, A, D)( ) the following estimate holds true, Here, X p (Γ 1 ) ⊂ W We emphasize the most common setting with p = 2 and q = ∞. In this Hilbertian case, C 2 ( ) = c( ) −1 C(D, Γ 0 ) with c from Assumption (PA1) and a deterministic constant C(D, Γ 0 ) determined by Poincaré and trace inequality constants.
For p > 2, the verification of this condition is somewhat more involved but fortunately still holds true. In the case of a purely homogeneous Dirichlet boundary, that is, |Γ 1 | = 0 and for the plain Laplacian with A ≡ I, for any Lipschitz domain it follows by Jerison and Kenig 45 such that for all 2 ≤ p ≤ p * , The result pointwise a.e. in extends to A( ) ≠ I by a perturbation argument (see References 19,46) and can in particular be translated to our pathwise framework. Define with (p) ∶= (1 − 2∕p)∕(1 − 2∕p * ) ∈ (0, 1) monotonously increasing. Let random variables p * , C p be defined as depending on the random variable bounds c, c of A from (PA1). An important observation is that p * K (t) is monotonously decreasing in t. With this construction in mind, let u( ) be the pathwise solution of (3) with B = A satisfying Assumption (PA1). Then for 2 ≤ p < p * K ( ) pathwise, a.e. it holds Hence, in order to satisfy the p * -condition in Definition 1, the quotient c( )∕c( ) in (14) needs to be bounded away from zero to obtain ess inf p Note that the uniform boundedness Assumption (A3) implies (A1). The case of mixed boundary conditions and L p estimates of solution gradients can be derived from results of Gröger. 47 W 1, p (D) estimates for Robin boundary conditions can be obtained from Amrouche et al. 48 with boundary right-hand side in W We end the pathwise discussion of this section with an examination of the (p, q) relation. There are two important cases: . This kind of convergence is at hand in case of • higher regularity: A(⋅, ) ∈ C k, (D) yields an error of at least (h k ) using an appropriate quadrature scheme.
• exact meshing: A(⋅, ) is only piecewise  k, , k ∈ N 0 , ∈ (0, 1) with (a priori known, for example, Reference 15) discontinuities exactly being resolved, for example, through adapted meshes with possible curved cells 2. q < ∞, thus p > 2: This case is important if A(⋅, ) lacks spatial regularity or cannot be resolved in a discretization step by someÂ such that the error in the L ∞ norm stays (1), as in Example 1. In this case, the q = ∞, p = 2 estimate might become meaningless. However, the p * -condition still ensures convergence of the perturbed solution based on weaker approximation requirements on the coefficient. At this point, one may ask for the approximation in the weakest norm possible, that is, q = 2 such that We note that this indeed is possible by an interpolation argument, assuming that A andÂ satisfy the uniform boundedness (PA1) and (A3). Then, with a constant C( , q) = C(q, ||A(⋅, ) −Â(⋅, )|| L ∞ ). As a consequence, in order to obtain a pathwise good approximationû( ) of u( ), it is sufficient to control the L 2 approximationÂ(⋅, ) of A(⋅, ). However, we note that by the employed interpolation this estimate introduces a reduced convergence order by a factor 1/q. Hence, ideally one can strive for an L q approximation to avoid a degeneration of the convergence order.
Proof. By (9) and the p * -condition, it holds pathwise that The result then follows by multiple applications of the Hölder inequality. Choose such that Q = T and denote by ′ its conjugate exponent. Then, skipping the pathwise dependence, this yields .
Further iterative application of Hölder estimates yield the desired result. ▪ The a priori perturbation estimation gives a qualitative statement of some closeness of the approximation to the true solution if the approximation error of the involved coefficient A is controlled. We point out that the important case of coefficients that are numerically approximated in L Q (Ω, L ∞ ) is included. However, note that in this case the need of a p * -condition can be relaxed.
While the assumptions in Corollary 1 are rather mild, the important case of random fields A representing composite random materials with random inclusions (to be examined in the following example) might nevertheless be excluded. If this type of materials is modeled pathwise by piecewise Hölder continuous or even smoother data then it fits in the setting of Theorem 2. We shall illustrate this line of thought with a small example.

Example 1 (Inclusion with random radius). Let
with the Euclidean centered ball of radius ( ) denoted by B ( ) (0) and an indicator function , see Figure 1. In noncurved mesh-based discretization schemes using standard Gaussian quadrature rules, the spherical random interfaces and thus the jump of the coefficient A(x, ) ∶= (x, )I cannot be approximated pathwise in L ∞ (D). In particular, for any such finite quadrature scheme, by the Gibbs phenomenon this error is (1) with respect to the L ∞ norm.

Example 2 (Nonperiodic random matrix composite material). Consider a partition of
For the parameter reference domain Ξ ∶= [−1, 1] d+1 and r = s, consider a bijective mapping r ∶ Ξ r →  r for some 0 < ≪ R s with Here r maps the parameter y r to a radius r and a point p r within D s modeling variability in the size and the position of a circular inclusion having a positive distance from the boundary D s . Furthermore, we distinguish between the index s for physical structures and r for the parametric enumeration as in Section 4.1. Then with y = (y r ) r and the mapping yields a parametric representation of a composite material with inclusions defined on nonoverlapping subdomains having variable positions and sizes, see Figure 1.
In a numeric simulation, we may approximate this setting, for example, by the following approaches •̂M C is a pathwise projection of (⋅, y) to piecewise constants on an underlying mesh, •̂P C is a cluster approximation, that is, a piecewise constant approximation with respect to the parameter y, see Section 4.2.4 •̂s urr implicitly approximates via local surrogates in parametric DD methods, see Sections 4.2 and 5.2.
All of the above approximation schemes are not converging in L Q (Ξ; L ∞ (D)).
Remark 1. The convergence assumption in L Q (Ω; L q (D)) and the requirements of the p * condition can be further relaxed in the fully discretized setting. In fact, it only has to hold In the discrete case, by Strang's lemma, one needs an estimate of Hence, for a fixed discretization level , we have ∇w ⋅ ∇v ∈ H and as a consequence it is sufficient to have thatÂ converges weakly to A in the H topology to control the error in (23). In particular, it suffices thatÂ → A strongly in H instead of stronger convergence in L 2 (Ω; L ∞ (D)) or L Q (Ω; L q (D), when starting from the discrete setting at the cost of an unknown convergence rate.

SURROGATE RESPONSE FOR NONSMOOTH RANDOM OPERATORS
The construction of surrogate models for problems with high-dimensional input (parameters) becomes essential when extensive sampling is computationally expensive in comparison to the construction of adequate functional approximations. These may even be viable for nonsmooth data since it can still be possible to exploit local smoothness, which then results in sufficiently accurate surrogates (in certain parameter subdomains).
In the literature there is a vast amount of surrogate types, including generalized polynomial chaos expansions (gPCE), for example, References 6,11, its multielement extension, for example, Reference 49, low-rank 35,50 and sparse grid techniques, for example, References 7,51 or neural networks.
A central motivation for the proposed approach is the treatment of parametric composite materials. With this in mind, our aim is to build (local) surrogates for (local) maps within a DD framework presented in Section 4. For this, we focus on two surrogate types based on a parameter space decomposition, namely, the multielement generalized polynomial chaos expansion and a PoU interpolation. The advantages of these approaches is the ability to explore areas in the parameter domain in which the underlying maps are more regular. In the case of low-rank structures, also hierarchical tensor representations might become useful, which we defer to future research.

Surrogates for matrix valued functions
In preparation of the framework of localized descriptions of randomness in a DD setup, we discuss different surrogate models for random matrices of the form for some indices r, s ∈ N. In our application, n and m depend on the number of interface degrees of freedom (d.o.f.s) or subsets of these. The choice of surrogate for the map (24) should be made dependent on the regularity of the involved map.
In particular, we have the following surrogates in mind that are based on (quasi-)best approximations or interpolations: 1. hPCE surrogates: A set of piecewise orthonormal polynomials w.r.t. an underlying measure is constructed, see Section 3.2. The coefficients in the hPCE series are computed via projection, or (sparse) quadrature schemes. 2. Hierarchical tensor surrogates: A number of samples of M s and a given underlying discrete tensor basis is provided.
Then, a tensor reconstruction as in the (nonintrusive) variational MC method 52 yields a L 2 -compressed low-rank representation. 3. PoU interpolation surrogates: Based on an adaptive mesh of Γ k , a discrete PoU basis with respect to the mesh is used. This can, for example, be obtained by a Lagrange basis with respect to mesh nodes. Each basis coefficient is computed by a single sample. More details are provided in Section 3.3. 4. Sparse grid interpolation surrogates are build by evaluating a realization of M s on each sparse grid point.
Remark 2. In the case of low regularity of the map (24) and no knowledge of a basis for an approximation scheme with (quasi)-optimal convergence in the sense of p -summability of coefficients for some p ≥ 1, 53 one may be restricted to lower local parametric dimensions to obtain an affordable and accurate surrogate.
While out of the scope of this article, a progressive learning of an approximate basis as, for example, via a neural network regression might be useful, given its training only involves a manageable number of samples.

Multielement gPCE
The following presentation is motivation by the approximation of nonsmooth functions, where standard functional approaches may lack efficiency due to the Gibbs phenomenon.
We consider a probability measure on Ξ ⊂ R M ′ and the space For the envisaged application, Ξ (or Ξ r ) is the image and will be the push-forward of P of an underlying parameter random vector. Let  = N and assume that there exists a family of orthonormal polynomials then such a family does indeed exist. 54 Example 3. The condition (25) holds for any bounded domain Ξ or for Ξ = R M ′ and any Gaussian measure including the nonindependent case. 55 In case that exhibits a product structure, that is, the independent case, the index set  can be reshaped into a multiindex tensor structure. The existence of a complete orthonormal polynomial basis then may be answered via the Hamburger moment problem. 54 Remark 3. There exist probability measures such that no dense polynomial subset exists. A classical example is the log-normal case. For M ′ > 1 and a nonseparable probability measure or a nontensorized domain Ξ, such a family of polynomials may exist but is not necessarily unique. 54 Let  Ξ ⊂ N be finite and let {Ξ k ∈ (Ξ), k ∈  Ξ } be a finite nonoverlapping partition of Ξ with 0 < (Ξ k ) ≤ 1. This gives rise to the decomposition

Lemma 3.
There is a family of orthonormal polynomials that spans a dense subset in Y k .
Proof. Fix v ∈ Y k and > 0. Let k denote the indicator function with respect to Ξ k , then ,k ∶= k for ∈  is a polynomial in Y k . Let Ψ k be an orthonormalized version of Π k ∶= { ,k , ∈ } with respect to the inner product in Y k . Since Ψ is dense in Y , there exists  ⊂ , | | < ∞ and a polynomial ∈ spanΨ with ||E k 0 v − || L 2 (Ξ,(Ξ), ) < . Here, Motivated by Lemma 3, we define an (weak) orthonormal polynomial set which spans a dense subset in Y k . We note that in the independend case the decomposition of Ξ has to respect the tensor structure to obtain a product structure of the polynomial chaos.

Lemma 4.
⨁ k∈ Ψ k spans a dense subset in Y with Ψ k .
As the partition of Ξ can be interpreted as a possibly nonregular meshing, we abbreviate the construction of multielement generalized polynomial chaos as gHPCE, motivated by hp-FEM in the standard Lebesgue spaces L 2 (D, d ).
For more details on the existence of dense generalized polynomial chaos we refer to References 6,11,54 and for the multielement generalized PCE to Reference 56.

Surrogates based on PoU-interpolation
This section is devoted to the derivation of an adapted hybrid surrogate construction. The adaptivity is such that it explores areas of higher parameter regularity resulting in "trust zones," which allow for efficient surrogate models. In numerical applications, a local sample is checked if it is contained within a trusted area. In this case, the interpolation surrogate is evaluated. Otherwise, the method falls back to classical sampling. Algorithm 1 summarizes the proposed idea. Note that the sample generation with the DD approaches in the upcoming section requires the assembly and algebraic modifications (such as an inversion and a reordering) of the involved operators. This indicates the potential speed up when large trust zones can be recovered by the algorithm since in this case the evaluation of the surrogate basically evokes no computational costs. Let Ξ r ∶= img r and assume a nested set of discrete function spaces of level = 0, 1, … , L, defined by such that its elements form a PoU on Ξ r , that is, ∑N i k=1 s, k ≡ 1. To describe an important class of this spaces, we consider a family of cell partitions  r of level = 1, … , L of Ξ r with tree structure: for each cell in  r +1 there is a unique father cell in  r . Furthermore, assume that for each k = 1, … , N the function s, k has support in one cell in  r +1 only, that is, We associate with each such function in U r s, a global degree of freedom and a unique associated element in Ξ r denoted by global coordinate degree of freedom.
With these constructions we formulate an adaptive scheme to build up surrogates for the parametric matrix ensemble in Algorithm 1, used in our numerical experiments. The aim of the proposed algorithm is to balance the maximization the trust region while keeping the creation process as cheap as possible.
Remark 4. The hybrid structure, namely, letting the surrogate coincide with a sample on the nontrust region, allows to control the cost of creating the surrogate in L iterations. For large L, this leads to overrefinement and thus too many sampler calls.
Remark 5. Although the above representation is presented for matrices, the technique also extends to matrices given only implicitly: If we are interested in building a surrogate for ( )  → A( ( )) −1 , due to the PoU approach, for each global coordinate degree of freedom we store a LU factorization only. Then, the evaluation of the surrogate at a point p requires several forward-backward substitutions associated to all basis functions evaluating to nonzero at p. Note that here each LU decomposition may have its own sparsity pattern determined by its pivotization and scaling. • surrogate quality parameter tol > 0 • maximal level of refinement L Output: hybrid surrogate for r ( )  → M s ( r ( )). for all T ∈  r do 8: judge quality of current surrogate on T w.r.t tol and set T to be trusted / not trusted 9: update trust region / nontrust region with T 10: if T is not trusted then 11: mark T in ( r , U r s, ) for refinement, that is, add to markings 12: end if 13: end for 14: if markings = ∅ then break 15: else  r +1 , U r s, +1 ∶= refine( r , U r s, , markings) 16: end for return hybrid surrogate: defined in U r for input within trust region and else coincides with .

PARAMETRIC DD
This section is devoted to the introduction of substructuring methods with locally defined random input fields. Such local definitions can be obtained by a localization of random fields as in Reference 57 or is given naturally as in the case of matrix composites. We first present the abstract framework of local random fields. Second, based on a semidiscretization we translate the framework to the simple parametric Schur complement method similar to Reference 36, which is rather simple in terms of implementation but lacks weak scalability. Finally, we introduce our parametric FETI-DP method and give details of appearing local parameter dependend operators for both, the system operators and the preconditioners. As a special case of a zeroth-order multielement PCE, we discuss a cluster sampling based FETI-DP method.
The proposed method represents a special case of a global sampling routine, which can be embedded into advanced sample schemes such as MLMC. A global sample is then transformed into several local contributions (i.e., local samples). Based on the underlying substructuring method, local surrogates based on local decoupled samples are constructed in a first (offline) step. The second (online) stage is the sampling stage in which the locally constructed surrogates are applied. Recall that by the use of the trust zone based surrogates described in Section 3.3, local sample processes can be accelerated by avoiding assembling and algebraic modification. Our method is summarized in Algorithm 2.

Abstract random domain decomposed model
Consider a partition of D into mutual disjoint, nonempty connected Lipschitz subdomains D s ⊂ D, s = 1, … , N SD < ∞.
We are interested in random fields given in parameterized form and possibly localized with respect to (D s ) s . Such local representations of random fields may for instance occur in the following scenarios • Localization to D s via KLE, for example, References 58,59, in case of underlying Gaussian random fields, see Example 4.
• Local uncertainties, for example, the random composite material model from Example 2. Instead of an abstract random field A, we consider a parametrization by a given stochastic coordinate system represented by some ∶ Ω → R M ′ , M ′ ∈ N ∪ {∞}. Let ( ) be the sigma algebra generated by and (Ω, ( ), P) the considered probability space. Furthermore, assume a subenumeration of the M ′ -dimensional random vector into M ≤ N SD blocks of m r -dimensional (sub) random vectors r ∶ Ω → R m r with 1 ≤ m r ≤ M ′ and components r i ∶= ( r ) i , i = 1, … , m r , r = 1, … , M. We are then interested in a physical model described by a linear PDE with randomness modeled by where at most one input random subvector r acts on a subdomain D s as shown in Figure 2. If the latter is the case, we call the tuple (r, s) an active index, which we collect in The abstract equations reads for a.a. ∈ Ω as with ( ( ), x) = ( r ( ), x) for x ∈ D s for each (r, s) ∈  active from (30) and ( ( ), x) = (x) else. An analog structure holds for the boundary operator . Here, we assume existence and uniqueness of a solution u ∈ L 2 (Ω, ( ), P) ⊗  for some separable Hilbert space , for example,  = H 1 Γ D (D) in the case of (1). Note that the abstract model includes several important cases of random modeling: • Direct modeling with independent random components, that is, ( s ) s are mutually independend random vectors and  is defined by material parameters that depend locally on r in a possibly nonlinear manner.
• The setting r ≡ corresponds to global random contributions, which includes the nonlocalized (standard) KLE case.
• The local coordinates r are obtained by a localized KLE. In the case that the underlying random field is Gaussian, and r consist of mutually independend Gaussian random variable components, although now ( r ) r is not independent in general.
• The case of one localized uncertainty input M = 1 as in Reference 60.
Remark 6. In order to balance the workload for parallel computations in the DD schemes introduced by varying sizes of subdomains described below, one might decompose a large subdomain, for example, D s , with operator dependence only on r , (r, s) ∈  active into D s = D s ′ ∪ D s ′′ , thus updating  active with (r, s ′ ), (r, s ′′ ) instead of (r, s). This may require more involved system preconditioners since the material coefficient on the interface D s ′ ∩ D s ′′ might behave nontrivially. Assume that the global KLE truncated after M terms yields the desired accuracy. In the case that the local eigenvalues ( s j ) j have a (much) lower magnitude, the local KLE can be truncated after m s ≪ M terms, that is, a low-dimensional local representation of the random field. Furthermore, there exists a matrix T s ∈ R m s ,M such that ( s j )

Example 4 (Localized KLE, 58). Given the mean and covariance kernel, a random field A = I may be written via KLE as
In the case of Gaussian random fields, one gets the favorable property that ( s j ) are independent Gaussian random variables as well, enabling local dense polynomial chaos approximations. Note however that ( s ) s is not independent in general.
Remark 7. Contrary to Example 4, in the case that ( i ) M i=1 are independent non-Gaussian random variables, the distribution of ( s j ) m s j=1 can be arbitrarily complex since the involved linear map T s introduces a nonlinear mapping of distributions due to m s < M. Hence, one cannot expect ( s j ) m s j=1 to be independent and the existence of a dense polynomial chaos in L Q (Ω, ( s ), P) is not ensured in general, for example, if img s ⊂ R m s is unbounded. Note that the image might be of lower topological dimension, for example, a submanifold only.

Parametric hybrid DD based on semidiscretization
DD generally refers to the splitting of a PDE or an approximation thereof into coupled problems on smaller subdomains, forming a partition of the original domain, see Reference 20. We present the framework on a semidiscrete formulation with regard to some suitable discrete subspace of . Then, the structure of the decomposed random model (31) and (32) results in a random system of equations ( ( ))u( ( )) = F( ( )).
The DD approach for random PDEs is applied in the context of global and local KLEs in References 32 and 36 with a combination of model reduction techniques 61 for the Schur complement or the FETI-DP method. We note that the technique based on global random coordinates as in Reference 32 is based on a stochastic FE formulation without a sampling stage, which limits it to a moderate number of random coordinates due to the curse of dimensionality. The analysis and application of hierarchical tensor formats in this context is deferred to future research. To alleviate the issue of moderate dimensions, our aim is to rely on local random coordinates only, for example, based on canonical tensor representations as in Reference 34. In what follows, we present both techniques, the Schur complement and the FETI-DP method based on adaptively constructed local surrogates that enable accelerated sampling. The devised general hybrid approach is summarized in Algorithm 2.
Remark 8. If the local representation of ( r ( ), x) is given by means of a gPCE with independent ( r k ) then the involved local random operators may be represented in low-rank formats as in Reference 62, allowing for higher dimensional local input. This approach transfers to compressed structures of all involved parametric suboperators.
We note that in case that the quantity of interest QoI(u) is only local, not all decoupled local problems may need to be resolved in Algorithm 2.

Parametric Schur complement method
As a first approach to nonoverlapping DD methods, we consider the Schur complement method, see, for example [ 20, ch. 4-5] for details. Given the simple design of the Schur complement method, we present some details of the appearing local operators. This should serve the reader to translate the abstract framework of the preceding section to the FETI-DP method proposed in the following. Further details of the approach applied for random PDEs can be found in References 36,37. For simplicity, the presentation is based on the elliptic linear equation (1) with the random dependence setting introduced in the beginning of this section. Given a physical discretization space V h ⊂  ∶= H 1 Γ 0 (D) spanned by a nodal Lagrange basis on a matching triangle (tetrahedral for d = 3) mesh of D, we consider a semidiscretization with respect to the physical coordinate x ∈ D. Let us denote by Π the (primal) interface if surrogate creation cost < C $ then 3: .
Block Gaussian elimination leads to an equivalent system to (33) given bỹ ( ( ))u( ( )) =F( ( )), where ( ( )) ∶=  (30), it follows that these matrices are given by where each matrix with index s only depends on r ( ) if (r, s) ∈  active and A ΠL = A LΠ almost everywhere in Ω. In particular, for each active index (r, s), we have Algorithm 3. Realization as parametric Schur complement method • local system surrogates for ( )  → S ΠΠ ( ( )) • local surrogates for preconditioner of (parametric) S ΠΠ . Output: approximated realization of u k = u( k ( )) or of subdomain parts u s k = u k|D s . 1: Compute right-hand side realizationF k = [f L ,f Π,k ] =F( k ( )) of (34) via system surrogates. 2: Iteratively solve S ΠΠ,k u Π,k =f Π,k via PCG. The application of operator S ΠΠ,k and its preconditioner are based on evaluation of local surrogates. 3: Solve block diagonal system A LL,k u L,k = f L − A LΠ,k u Π,k in parallel.
In the deterministic setting, (34) is solved by first iteratively solving for u Π and second solving (in parallel) for interior contributions u s The full matrix S ΠΠ is never formed explicitly. Its action on a vector involves the application of the smaller matrices S s ΠΠ . In the deterministic case, the latter might be realized by a LU decomposition of A s LL and we build local surrogate modelsŜ s ΠΠ ≈ S s ΠΠ . Depending on the regularity of the maps and the input dimension of r , this might be realized by interpolation (e.g., sparse grid, Lagrange interpolation) or (quasi-)best approximation techniques (tensor reconstruction, hPCE). We note that the construction of the surrogate models may be expensive if the number of local d.o.f.s associated with Π or the random dimension of r grows or if there is a lack of regularity. In a practical implementation, this cost factor has to be compared with the cost of a full sampling approach with several backward and forward solves of the involved LU decompositions in the iterative process of the application of the Schur complement matrix. It is known that the condition number of a deterministic Schur complement matrix grows like (1∕Hh), where H is the diameter of the subdomains and h the maximal element size in the subdomains. 20 The proposed method is summarized in Algorithm 3.

Parametric FETI-DP
The FETI-DP method is known to successfuly balance the requirement of a minimal communication by a coarse space by keeping good convergence rates within the preconditioned conjugate gradient solution scheme. For the sake of simplicity, we stay on the algebraic matrix level subject to the discretization of the symmetric model problem (1). For interpretation in a Hilbertian framework, we refer to Reference 63 and for a formulation in intermediate approximation FE spaces to Reference 20.
After choosing the interior, dual and primal d.o.f.s indexed by the subindices I, Δ, and Π for a given physical discretization based on meshes on D s that are aligned with each other on the subdomain interfaces, the structure of the system (33) is given by Here where ( ( )) ∶= Here, it holds thatũ Π = with the random FETI-DP dual matrix F and right-hand side contributions such that Note that for better readability we omit the parametric dependence of matrices on the right-hand sides of the above equations. Examining the occurring matrices, their dimensions and their random dependence on local random vectors s ( ), we now discuss the construction of surrogate models. In order to avoid multiple LU forward and backward substitutions in the iteration of the system F =b , we aim for a lower total cost of computing and storing the surrogates for a given total number of samples. The parametric matrix F contains three contributions, which are discussed separately: • With the random dependence A s LL = A s LL ( r ( )), (r, s) ∈  active , recall that The inverse of A s LL can be seen as a 2 × 2 block matrix Due to the design of  s L , only information of B ΔΔ contributes. Consequently, we aim for a surrogate of the mappings • local-based surrogates for preconditioner of (parametric) F. Output: approximated realization of u k = u( k ( )) or subdomain parts u s k = u k|D s . 1: Compute r.h.s. sampleb k = [f L ,b Π,k ,b ,k ] =b( k ( )) (38) via surrogates (41)-(44). 2: Iteratively solve F k k =b ,k via a preconditioned CG method via application of F k and its preconditioner based on (40)- (42) and (50). 3: Solve Schur complement systemS ΠΠ,kũΠ,k =b Π,k +S Π ,k k via (42). 4: Solve block diagonal system A LL,k u L,k = f L −Ã LΠ,kũΠ,k −  T L k in parallel.
The surrogates of this mapping involve the largest coefficient matrices, depending only on the local dual d.o.f.s.
• We are now concerned with constructing surrogates to buildS Π . Note that with the notation in (39) it follows Thus, we aim for surrogates of the mapping whereas the coefficient matrix sizes only depend on local d.o.f.s associated to the coarse space and hence are rather small compared with the classic coarse space in the Schur complement method.
The right-hand side of the parametric system can also be computed via surrogates.
• The local surrogates employed for the evaluation ofb( ( )) are given by for all active indices (s, r).
• The computation ofb can be based on local surrogates for active indices (s, r) in addition to the surrogates forS Π ,S ΠΠ andb.
Note that all local surrogates can be computed in parallel and stored separately. The suggested approach is summarized in Algorithm 4.
The advantage of deterministic FETI-DP lies in its potential to lead to a weakly scalable algorithm that is enabled by a suitable choice of coarse space (associated with Π) and preconditioner P for the submatrix F. 20 In the context of DD techniques, the term scalability is associated with the iterative solution cost of the discrete system, which should not deteriorate when the number of subdomains grows. Let H denote the diameter of the subdomain and h the maximal diameter of the subdomain mesh cells. It can then be shown that where denotes the condition number, see, for example, Reference 20 and the references therein. The construction of the scaled lumped and Dirichlet preconditioner and its extension to the random case using surrogates is discussed in Section 4.2.3.

Preconditioner
In what follows we introduce surrogates for two classic FETI-DP preconditioners, namely, the lumped and the Dirichlet preconditioner P −1 lumped and P −1 Dir . These are given pointwise by Here, M s # = M s # ( r ( )) for all (r, s) ∈  active . The diagonal scaling weight mappings W s have a more involved parametric structure due to the coupling of neighbored random dependencies shown in (47).
We fix a subdomain D s and let D s ′ be any neighboring subdomain such that there exist dual d.o.f.s Δ s i and Δ s ′ j in the local spaces associated with local meshes on those subdomains that both account for the same global dual degree of freedom. We collect those D s ′ into D s Δ , noting that D s ∈ D s Δ per definition. Then, following Reference 23, the diagonal scaling W s is defined as With this construction, the diagonal scaling enables the preconditioner to take material heterogeneities into account. An important observation is that if the assembled local dual matrices A s ΔΔ and their neighborhood counterparts A s ′ ΔΔ depend on r and r ′ , respectively, then W s is a matrix valued function depending on r and r ′ . Hence, a direct surrogate construction may be not applicable due to the sum of involved local parameter dimensions. Since the underlying physical local meshes are fixed, we can work around this issue. For s ′ such that D s ′ ∈ D s Δ , we introduce diagonal matrices W s s ′ with Then, we can compactly write where W s s ′ = W s s ′ ( r ′ ) for all (r ′ , s ′ ) ∈  active . With this construction, we may generate surrogate of the following maps, which depend only on local parameters, and define the application of the preconditioner based on these local surrogates.
Remark 9. In summary, we presented a pointwise surrogate approach, which for each sample up to surrogate precision enables weak scalability based on the deterministic results. Alternatively, a fixed preconditioner such as the mean E[P −1 # ] can be considered as in Reference 36 for the Schur complement method.

Parametric FETI-DP for cluster sampling
In this subsection, we consider cluster sampling (as a special case) that corresponds to a piecewise constant approximation of the involved parametric assembled discretization in the spirit of zero-order gHPCE from Section 3.2. While the discretization is applicable to low local parametric dimensions only due to cost aspects, it turns out that the resulting surrogates have a very simple structure. In this case, the surrogates are not needed to be build explicitly and the method involving its preconditioner is a simple generalization of the classic deterministic FETI-DP in all its algorithmic details, like the use of precomputed LU decompositions of the involved inverse matrix applications. We shall describe the benefits illustrated with the random FETI-DP matrix F. Let (r, s) be active and consider the zero-order gHPCE for any involved local random matrix M s such that A key observation follows if M s j is invertible for all j ∈  . By the push-forward r = P # r and identification y r = r ( ) and Ξ r = ∪ j∈ Ξ r j , it holds pointwise As a consequence, for example, the inverse of the Schur complement structure S s ΠΠ ( r ( )) −1 has an easy form given as where only one summand contributes to a realization of r . With that structure, for the sth summand R s with indicator function Γ r j . Thus, for a global realization, the application of F relies on precomputed factorizations of A s LLj and storage of all involved coefficient matrices for s = 1, … , N SD , and j ∈ J. This indicates a limiting factor if |J| gets too large by the curse of dimensionality.
The piecewise constant approximation scheme in the stochastic space is suitable for domain-wise random fields  → a(x, ) = a(x, X( )) parameterized by a discrete random variable X with finite-dimensional image {x 1 , … , x |J| }. Such fields might be represented as with polynomials j with j (x i ) = ij and functions j in physical space. A special case of this is the checkerboard material with 0 (y) = y, 1 (y) = 1 − y, and j ≡ j.

NUMERICAL EXPERIMENTS
This section illustrates the numerical performance of the FETI-DP method with local hybrid surrogates using the PoU interpolation approach of Section 3.3. The choice of the surrogate is motivated by the fact that gHPC surrogates rely on F I G U R E 3 Adaptive surrogate construction based on macroelement anisotropic refinement procedure for s  → B s ΔΔ ( s ) associated to the inner square D s with tol = 10 −4 and 100% trusted zone after six iterations for the smooth problem in Section 5.1 integral computations (such as projections) in the parameter space. Assuming a nonsmooth underlying structure, these are not directly accessible via an exact numerical quadrature. The interpolation approach relies on point samples that are always available by the local sampler  s sample , which subsumes the involved operators, see Algorithm 2. We note that the choice of surrogate is independent of the underlying substructuring method. In particular, we do not discuss the Schur complement method, as it has been discussed in Reference 36 at least for the smooth KLE case based on gPCE. The implementation of the hybrid local surrogate FETI-DP method was carried out as part of the open source ALEA python library, 64 including the surrogate adaptivity, FE backend and all required algebraic modifications. We use the deterministic FE backend FEniCS 65 for assembling of local operators and UMFPACK 66 for the inversion of the involved A s LL from (39). The discrete random coordinate spaces are realized with a custom implementation of a hierarchical tree-based decomposition of M-dimensional cells such as tensorized hyperquadrilaterals. Each such element can be refined anisotropically by bisection separately in each coordinate direction. In addition, a cell is refined into macroelements (marked as gray in the figures below) while remaining a hyperquadrilateral if subareas are not necessary for further refinement in the marking process. This is realized by keeping the Lagrange basis of the preceding level to represent the solution on a macro element only while refined subdomains get enriched by further Lagrange basis functions of variable degree. This has the advantage of saving various d.o.f.s in the surrogate creation process, which determine the total amount of work in the offline phase.
Two experiments based on the model problem (1) are considered with homogeneous Dirichlet boundary conditions. In Section 5.1, we examine a smooth local parametric dependence with an additional challenge introduced by a nonlinear coupling of random and physical coordinates. A random cookie problem as described in Example 2 representing nonsmooth data is depicted in Section 5.2.

Smooth problem with nonlinear coupling
We consider D = [0, 1] 2 and a partition of D into N × N subsquares D s with N = 3, x = (x 1 , x 2 ). A smooth random field with a coupling of physical and stochastic coordinates is given by In the experiment we utilize a physical discretization with p = 1 FE on uniform 30 × 30 triangulations of D s , s = 1, … , N SD . The nonlinear coupling of random and physical coordinates introduces a layer of the involved matrix valued random maps that differs for all s = 1, … , 9, also depending on the size of values taken in x. We illustrate the different structures in Figure 6 by final meshes obtained by the adaptive scheme. We point out that = 2 , or on larger F I G U R E 5 Adaptive surrogate construction procedure for s  → S s PP ( s ) associated to an inner square D s with tol = 10 −4 and 100% trusted zone after five iterations for the smooth problem in Section 5.1 domains D and thus related larger range of x, the surrogate construction becomes more involved and finer layers have to be resolved.
In the numerical investigation, it is observed that given a tol > 0 for the local interface surrogates, the error between a sampled solution based on surrogates denoted byû k h and its analogue (MC) sample solution denoted by u k h is of higher order. The situation is shown in Figure 7 with a pointwise difference of order (10 −6 ) and ||û k h − u k h || H 1 0 (D) ∕||u k h || ∈ (6 × 10 −5 ).

5.2
Cookie problem, model problem for more general composite materials Remark 10. For the physical discretization, a fixed uniform (local) triangular mesh is used implicitly for each realization of the composite structure. In particular, given a realization of the composite, the mesh does not resolve the jump in the coefficient, which in a deterministic setting is critical from an accuracy point of view. In order to benefit from using such a fixed mesh, the FE basis functions can be enriched by an element-wise basis function with a jump in the gradient of the discrete basis function as described in Reference 67. With this, the presented technique exhibits a better approximation quality with respect to the solution u. On each square domain D s , the random field A = aI P-a.e. is described by (x, ) = (x, s ( )) = B( s ( )) (x) + 20(1 − B( s ( )) (x)), x ∈ D s .
Here, B s ( s ( )) denotes a random ball modeled by a random radius and a random x-position p s x 1 of a point p s = (p s x 1 , p s x 2 ) in D s s.t. ( , p s x ) = s ( s 1 ( ), s 2 ( )) with a map s chosen such that B( s ( )) is uniformly bounded away from D s by a given distance as in Example 2. By a push-forward we may identify y s = s ( ) and interpret y s  → (⋅, y s ) as an element in L 2 ([− 1, 1] 2 ; L ∞ (D)). We stress that ∉ ([−1, 1] 2 ; L ∞ (D)) but ∈ L ∞ ([−1, 1] 2 ; L ∞ (D)). The defined composite material coefficient lacks regularity in the random as well as the physical coordinates, satisfying the conditions of Section 2. Note that the full random dimension is M ′ = 2N H N W , for example, with M ′ = 200 for 20 × 5 squares as shown in Figure 8. Due to the same structure of the domain-wise nonperiodic material description for the homogeneous Dirichlet problem, it suffices to compute local surrogates for nine subdomains only (four associated to each corner and edge, one inner domain. The adaptive procedure is illustrated for the inner domain case in  We point out the refinement pattern toward the boundary of the parameter domain. The right interface corresponds to the case of maximal radius of the random ball, the top and bottom interface correspond to the maximal displacement of the ball. Due to the surrogate construction based on a semidiscretization, given sufficiently fine physical meshes, the intersection of the random ball and the support of the basic functions associated locally to the boundary shrinks. Thus, the dependence of the nonsmooth influence gets smaller. Furthermore, we observe a constant error for B ΔΔ in Figure 9 that accounts for the remaining support interaction and exhibits a notable very slow decrease. This interaction is further shown in Figure 12, when solving for a smaller tolerance tol = 10 −3 . In the numerical experiments, building the surrogate s  → A s LL ( s ) −1 as in Remark 5 leads to massive uniform refinements with an (arbitrarily) slow error decrease. However, the involved interface operators are observed to exhibit rather smooth subareas in the parametric space with remaining small nontrusted areas.
The underlying pathwise quadrature in the assembling of the system matrices is implemented via a piecewise constant representation of the matrix composite realization with values derived by areas of intersections of circles and triangles. Since the composite is of circular shape and a first-order FEM on triangles is used, this quadrature is exact up to machine precision without the need of adaptive quadrature schemes. Due to the design and meshing of the subdomain, surrogates for only nine subdomains (four touching the corner, four having intersection with left, right, top, and bottom side, and one interior) have to be created.

F I G U R E 11
Adaptive surrogate construction procedure for s  → S s PP ( s ) for an inner square D s with tol = 0.5 × 10 −2 . 100% trusted zone after five iterations for an interior subdomain of the nonsmooth problem in Section 5.2 F I G U R E 12 Adaptive procedure for s  → B s ΔΔ ( s ) associated to an inner domain D s with tol = 10 −3 . The maximal displacement areas and the maximal radius areas dominate the error and remain nontrusted. After six iterations a trusted area of 93.35% is obtained for an interior subdomain of the nonsmooth problem in Section 5.2

Discussion and conclusion
We developed a DD sampling approach based on local surrogates with the potential so significantly accelerate the sample generation by avoiding the assembly of operators and algebraic operations. The hybrid local surrogates are based on the exploration of "trust zones" and fall back to classic sampling in case that a local sample is untrusted, that is, a surrogate would not yield accurate approximations. For the surrogate construction, we presented an interpolation scheme. This turns out to be advantageous since local samples are accessible and no quadrature (as in projections for the gHPC) is needed. We point out that quadrature-based techniques might not be feasible at all for the considered nonsmooth random fields.
While the chosen FETI-DP framework needs more local surrogates than a Schur complement, it allows for weak scalability as a big advantage. In fact, the global stochastic dimension easily exceeds M > 100 for the considered application with random matrix composite random fields. This motivates the derivation of an efficient sampling algorithm which can be embedded into an advanced (multilevel) sampling scheme such as MLMC. An analysis can then be carried out by means of the derived perturbation theory relying on the p * , closing the theoretical gap of random field approximation in a weaker topology.