Error estimation and model adaptation for a stochastic-deterministic coupling method based on the Arlequin framework

SUMMARY: The paper deals with the issue of accuracy for multiscale methods applied to solve stochastic problems. It more precisely focuses on the control of a coupling, performed using the Arlequin framework, between a deterministic continuum model and a stochastic continuum one. By using residual‐type estimates and adjoint‐based techniques, a strategy for goal‐oriented error estimation is presented for this coupling and contributions of various error sources (modeling, space discretization, and Monte Carlo approximation) are assessed. Furthermore, an adaptive strategy is proposed to enhance the quality of outputs of interest obtained by the coupled stochastic‐deterministic model. Performance of the proposed approach is illustrated on 1D and 2D numerical experiments. Copyright © 2013 John Wiley & Sons, Ltd.

multiscale methods; stochastic mechanics; Arlequin method; goal‐oriented error estimation; a posteriori error estimation

Numerical simulation has become an essential tool in design processes, providing for more flexibility and reduced production costs. Recent advances in this field include the development and consideration of multiscale and stochastic models. Multiscale models, on the one hand, aim at introducing appropriate models at each of the scales considered in a general problem and coupling them together. Stochastic models, on the other hand, aim at explicitly taking into account data uncertainties and modeling variabilities, which traditionally limit the predictability of deterministic numerical models. Such stochastic multiscale models, and the control of their accuracy, form the subject of the present paper.

We assume that there exists, for the study of a structure, a hierarchy of models defined at different scales. Although coarse models might yield effective results for the average response of the structure under distributed loading, these models may become ineffective when a localized defect arises (crack, hole, and local weakness) or when one is interested in local quantities. Finer‐scale models are then required, but usually involve higher computational costs, which renders their use unrealistic for industrial structures. The general idea of multiscale modeling is then to get around this dilemma by considering a coarse model wherever possible, while refining it locally around the zone of interest. Various multiscale strategies have been developed: (i) enrichment methods [1] , [2] in the FEM framework, in which appropriate fine‐scale basis functions are used around local defects; (ii) splitting methods [3] , [4] , [5] , [6] in which the solution field is split additively into a coarse field and a finer‐scale field, each being solved at its own scale, with fictitious boundary conditions at the lower scale; and (iii) superposition methods [7] , [8] , [9] , in which both coarse and fine‐scale models coexist in the local subdomain of interest and are solved concurrently by enforcing the matching of their solutions in an appropriate weak sense on an overlapping region. Extensions of these multiscale strategies to stochastic models have been proposed in the context of the heterogeneous multiscale method [10] , the multiscale finite element method [11] , and of the Arlequin method [12] , [13] . The latter is a superposition method that has been applied to the coupling of various models in the past decade (see [7] , [14] , [15] , [16] , [9] , [17] , [18] , [19] , [20] , among others) and on which we will concentrate in this paper [21] , [13] .

The purpose of stochastic multiscale modeling is to provide reliable results for engineering problems, with explicit account of data uncertainties and modeling variabilities at a reasonable computational cost. However, to reach that goal, the numerical errors associated to their solution must be evaluated and controlled. For more than 30 years, and particularly in the context of the FEM, a posteriori error estimators have been proposed (see [22] , [23] , [24] for reviews). These techniques allow to complement FE solutions with an estimation of their accuracy. After initial works on global error estimators, based on global energy norms of the solution, goal‐oriented error estimators have been developed [25] , [26] , [27] to assess the accuracy of specific outputs of interest (for example: stress average in a given zone or displacement of a point). Although most of the studies on error estimation are concerned with monoscale models, some specific developments have been proposed for coupled multi‐scale models. In particular, in the Arlequin framework, several papers have dealt with the evaluation and control of the accuracy of coupled atomistic and continuum models assessing both discretization and modeling errors, in the deterministic case [18] , [28] , [29] . A preliminary work has also been presented in the context of atomistic‐continuum stochastic models [12] but by using a different type of coupling from that considered in the present paper.

The main objective of the paper is to develop a framework for goal‐oriented error estimation when the Arlequin method is used to couple deterministic and stochastic models, as developed in [21] , [13] . To control the quality of approximate local quantities, we use residual‐type estimates and classical adjoint‐based techniques to get bounds on the local error. A second objective of the paper is to propose a method to split the different error sources, that is, modeling error, spatial discretization error, and stochastic error. The splitting is based on the introduction of intermediate reference models and on the construction of error estimates with respect to these new models, following [24] , [30] , [31] . The idea of evaluating and controlling modeling errors, independently of the discretization errors, was already discussed in [5] , [32] , [33] , [28] , [29] by comparing, in some adequate manner, the solutions obtained for the same problem evaluated at different scales. Extension to stochastic problems of such methodologies is, to the best of our knowledge, the first of its kind. This separation of sources enables to adapt the Arlequin model effectively by choosing, in a greedy manner, parameters of the numerical model that lead to the computation of the output of interest within a prescribed accuracy.

The paper is organized as follows. [NaN] introduces the reference (stochastic) model and its approximation in the context of the multiscale Arlequin method. In [NaN] , which deals with the first objective of the paper, the framework for goal‐oriented error estimation as applied to the stochastic‐deterministic Arlequin coupling is introduced. In [NaN] , we address the second objective of the paper. We propose a technique to split the error into different error sources based on the introduction of intermediate models, and we set up the associated adaptive modeling process. Finally, numerical examples are presented in [NaN] to discuss key numerical issues and illustrate the effectiveness of the error estimates for adaptive modeling purposes. Conclusions are drawn in [NaN] .

In this section, we describe the different models used in this paper. The reference monomodel is a continuum stochastic linear elliptic boundary value problem. This model is approximated using the Arlequin method to yield a surrogate problem, called the coupled model, in which a continuum deterministic model and a stochastic one interact with each other.

We consider here a stochastic linear elliptic Poisson problem representing a wide range of physical applications (heat equation, membrane, and Darcy equation). The problem is defined over an open‐bounded domain Ω with boundary ∂Ω, split into ΓD and ΓN such that ∂Ω=ΓD∪ΓN¯ (see Figure [NaN] ). Considering (Θ,F,P) a complete probability space with Θ a set of outcomes, F a σ‐algebra of events, and P:F→[0,1] a probability measure, the problem reads the following:

Find u such that

∇⋅K(x,θ)∇u(x,θ)+f(x)=0, almost everywhere (a.e.) in Ω, and almost surely (a.s.) ,

with the following:

u(x,θ)=0 a.e. on ΓD and a.s.,K(x,θ)∇u⋅n=g(x) a.e. on ΓN and a.s..

The stochastic behavior is driven by parameter K(x,θ) considered here as a stochastic field, P measurable on Θ ⊗ B(Ω) where B(Ω) is the Borel σ‐algebra generated by the open subsets of Ω. We assume that K(x,θ) is bounded and uniformly coercive [34] ; that is,

0

We denote W the functional space of admissible displacement fields as

W=v∈L2(Θ,H1(Ω)),v=0 on ΓD, a.s.

The weak formulation of the problem ‐ reads as follows:

Find u∈W such that

A(u,v)=L(v),∀v∈W,

where the internal virtual work is given by the bilinear form A:W×W→R and the external virtual work is given by the linear form L:W→R. They are defined, respectively, by the following:

A(u,v)=E∫ΩK(x,θ)∇u(x,θ)⋅∇v(x,θ)dΩ,

and

L(v)=∫Ωf(x)Ev(x,θ)dΩ+∫ΓNg(x)Ev(x,θ)dΓ,

where E[ ∙ ] is the mathematical expectation. It can be proved that the problem described previously admits a unique solution (see for example [34] ). We refer to this solution as the exact solution and denote it by uex.

Let us notice that problem  is a continuous stochastic problem. Therefore, two classical approximations processes have to be used: one for the spatial dimension and one for the stochastic dimension. These approximations could induce high computational cost, especially with regard to the stochastic dimension. To get an accurate approximation of statistical moments of uex by using a Monte Carlo approach for instance, a huge number of realizations [35] may be required, especially where critical zones have to be considered. Another method based on the polynomial chaos decomposition [36] , [34] can also be followed but the size of the matrix which has to be inverted can increase dramatically. In practice, because of the problem size, surrogate problems derived from this reference model are used. In the next paragraph, we present a family of surrogate problems, in the context of the Arlequin approach.

Here, we assume that we are interested in a local quantity of interest defined within a prescribed region Ωs ⊂ Ω. To reproduce accurately the effect of randomness on this quantity, we choose to keep the stochastic description in the local zone Ωs only. The parameter K is thus kept in this area as a stochastic field denoted by Ks(x,θ) while it is changed to a deterministic parameter denoted by Kd in the remainder of the structure. The latter is derived from K by means of stochastic homogenization [37] , [38] . In the Arlequin framework, this multi‐model representation is made possible by considering a deterministic model defined over Ωd ≡ Ω (called the substrate) on which we superpose a stochastic model defined over Ωs (called the patch). For simplicity of notation, we consider the case where the patch Ωs is fully embedded in Ωd as shown in Figure [NaN] . The two models are blended together in a coupling zone Ωc ⊂ Ωs (Figure [NaN] ) and the mechanical energy is distributed between the two models by using weight functions (αd(x),αs(x)) (see [7] , [9] , [21] , [13] for more details). The Arlequin problem then reads as follows:

Find (ud,us,λ) in Vd×Ws×Wc such that:

ad(ud,vd)+C(λ,vd)=ℓd(vd),∀vd∈Vd,As(us,vs)−C(λ,vs)=Ls(vs),∀vs∈Ws,C(μ,ud−us)=0,∀μ∈Wc,

where Vd={v∈H1(Ωd),v(x)=0 on ΓD} and Ws=L2(Θ,H1(Ωs)). The internal and external virtual works ad:Vd×Vd→R, As:Ws×Ws→R, ℓd:Vd→R and Ls:Ws→R are defined, respectively, as follows:

ad(u,v)=∫Ωdαd(x)Kd∇u(x)⋅∇v(x)dΩ,As(u,v)=E∫Ωsαs(x)Ks(x,θ)∇u(x,θ)⋅∇v(x,θ)dΩ,

and

ℓd(v)=∫Ωdαd(x)f(x)v(x)dΩ+∫ΓNg(x)v(x)dΓ,Ls(v)=E∫Ωsαs(x)f(x)v(x,θ)dΩ.

The mediator space Wc is built as a space of functions with a spatially varying expectation and a perfectly spatially correlated randomness: Wc={ψ+θcIΩc|ψ∈H1(Ωc),θc∈L2(Θ,R)} with θc a random variable and where IΩc is the characteristic function of Ωc, namely IΩc=1 in Ωc and IΩc=0 elsewhere. This mediator space has been chosen as it is the smallest one that ensures the stability of the formulation. Other possible choices included the restrictions on Ωc of either one of Vd or Ws. The former choice leads to an unstable mixed problem, and the latter leads to a waste or resources as the fine‐scale solution locks onto the coarse‐scale solution in the coupling area. Moreover, the specific structure of the chosen mediator space imposes implicitly that the (ensemble) average of the field us should be equal to the field ud, almost everywhere in Ωc, and that the variability of the space average quantity ∫ΩcEus−usdΩ should cancel (see [13] for more details). The coupling operator C:Wc×Wc→R is defined as follows:

C(u,v)=E∫Ωcκ0uv+κ1∇u⋅∇vdΩ,

where the parameters κ0 and κ1 are of the same units as a string rigidity divided by a length and a material rigidity, respectively (see [9] for details). The well‐posedness of problem  can be found in [13] .

The surrogate Arlequin problem provides two fields: ud defined on Vd, and us defined on Ws. We additionally introduce the continuous field uarlc∈W from ud and us given by the following:

uarlc=ud in Ωd∖Ωs,αdud+αsus in Ωs.

Obviously, the solution of the Arlequin model uarlc and the exact solution uex are different, even in the patch. An error is thus introduced by the consideration of the Arlequin model, even before any discretization process is introduced.

The Arlequin solution presented previously is an approximation of the reference monomodel solution. Its estimation through numerical techniques introduces two additional levels of approximation: a space approximation linked to the introduction of a mesh and corresponding finite element bases and one related to the stochastic dimension through the use of the Monte Carlo technique. We assume that we can construct from (uda,usa) a field uarl∈W. We then define the error as e=uex−uarl.

Let us recall here that the Arlequin model is used to estimate a local quantity of interest inside the patch Ωs. This quantity is a function of the monomodel solution, denoted q(u). We therefore seek to quantify the error η=q(uex)−q(uarl) by using goal‐oriented techniques.

In the case of deterministic Arlequin models, preliminary works dealing with goal‐oriented error estimation were conducted for 1D tests [18] . With the introduction of an adjoint problem and tools for modeling error control as defined in [33] , the error on a given quantity of interest was assessed. These works were extended to 2D and 3D applications in [39] , [28] , [29] for adaptive modeling of polymeric materials.

We seek to quantify the error on a quantity of interest q(u):W→R obtained with the Arlequin method. A quantity of interest could be, for example, the space average over a given subdomain of the mathematical expectation of u or of a given component of  ∇ u.

The idea is to express the output of interest in the global form (see [25] , [40] , [33] ):

q(u)=E∫ΩfqudΩ+∫ΓNgqudΓ+∫Ωpq⋅∇udΩ,

where quantities fq,gq and pq are called extractors.

We then define the so‐called adjoint problem, which reads in the case of a linear quantity of interest:

Find p∈W such that:

A(v,p)=q(v),∀v∈W.

In the remainder, we denote the exact solution of this problem pex. As shown later, the error η can then be estimated as a function of uarl and pex.

Defining the residual function R:W×W→R associated to  by the following:

R(u,v)=L(v)−A(u,v),

the use of the adjoint problem provides a tool to estimate the error η[27] and one can easily show that

η=R(uarl,pex).

Indeed, by using the definitions of η and pex, we have the following:

η=q(uex)−q(uarl)=q(uex−uarl) since q is linear=A(uex−uarl,pex) thanks to Equation (6) with v=uex−uarl=L(pex)−A(uarl,pex) since A is bilinear and uex verifies (4)

that proves .

In the case of a nonlinear quantity of interest, q(uex)−q(uarl) would usually be linearized around uarl by writing [41] :

q(uex)=q(e+uarl)=q(uarl)+q′(uarl,e)+o(e),?q(uex)−q(uarl)≈q′(uarl,e)

where o is the little o of the Landau notation and q′(u,v)=limθ→0q(u+θv)−q(u)θ is the tangent operator of q at u. The extractors fq,gq and pq are then defined by the following:

q′(uarl,e)=E∫Ωfq(uarl)edΩ+∫ΓNgq(uarl)edΓ+∫Ωpq(uarl)⋅∇edΩ,

In practice, as the adjoint problem is defined on the same space as the reference one, its solution pex is as intractable as the solution of the reference problem. It then has to be approximated by a surrogate problem. By using again the Arlequin framework, the approximate adjoint problem reads as follows:

Find (pd,ps,pλ) in Ṽd×W̃s×W̃c such that

ad(wd,pd)+C(wd,λ)=qd(wd),∀wd∈Ṽd,As(ws,ps)−C(ws,λ)=qs(ws),∀ws∈W̃s,C(pd−ps,μ)=0,∀μ∈W̃c,

where Ṽd, W̃s, and W̃c can be different from Vd, Ws, and Wc, respectively (see [28] , [29] in the case of a deterministic coupling between continuum and discrete models). In the case studied here, under the continuum form, we have Ṽd≡Vd, and W̃s may be defined over a domain Ω̃s that is larger than Ωs and included in Ω (Ωs⊂Ω̃s⊂Ω). Moreover, in their discretized form, the spaces Ṽd, W̃s and W̃c are chosen richer than Vd, Ws, and Wc, respectively. The number of elements that define the finite element spaces may thus be larger, and the description of the stochastic dimension may be finer also (by using more Monte Carlo trials for instance).

By using the FEM and Monte Carlo techniques for instance, an approximate solution can be computed for the coupled adjoint problem . We denote this solution (pda,psa,pλa) and additionally introduce parl similarly to uarl, assuming parl∈W. The error on q is then approximated by the following:

η≈R(uarl,parl),

Note here that, because we approximate R(uarl,pex) by R(uarl,parl), the approximation of parl has to be more accurate than the approximation of uex by uarl. This will be extensively discussed in the examples of [NaN] .

To estimate the error η, we need to evaluate R(uarl,parl)=A(uarl,parl)−L(parl). In practice, the fields uarl and parl are approximated by the discrete fields uarla and parla given by the following:

uarla=uda in Ωd∖Ωs,usa in Ωs.

Through simple Monte Carlo sampling and by using realizations of parla generated by the solution of the adjoint Arlequin system , L(parla) can be easily computed. Nevertheless, the Monte Carlo estimation of A(uarla,parla) involves realizations of product K∇uarla⋅∇parla. Moreover, realizations of K are not necessarily the same for the two systems. Realizations of pairs (K,uarla) on the one hand, and (K,parla) on the other hand, are obtained, as described in the previous sections by solving the Arlequin systems  and , respectively. In general, we thus cannot compute directly A(uarla,parla).

We present in the following an approach to generate, from realizations obtained by solving , new realizations of usa that correspond to the realizations of K being used for the adjoint Arlequin system , while following the same first‐order and second‐order statistics of the solution usa of the primal Arlequin system . The field uarla is then reconstructed using Equation  and we can thus compute the term A(uarla,parla).

To simplify the presentation, we first assume that both K(x,θ) and usa(x,θ) are Gaussian random fields. An extension to more general first‐order marginal distributions will be presented at the end of this section. Also, we consider the random vectors Ξ, U, and P that correspond to the random fields K(x,θ) and usa(x,θ) after space discretization. The objective is to generate new realizations Û of U so that we can use Ξ and Û to compute the residual. These new realizations are such that

E[Û]=E[U]:=U̲,E[(Û−U̲)(Û−U̲)T]=E[(U−U̲)(U−U̲)T]:=CovU,E[(Û−U̲)(Ξ−Ξ̲)T]=E[(U−U̲)(Ξ−Ξ̲)T]:=CovUΞ,

where we have introduced notations U̲ and Ξ̲ for the expectation of U and Ξ respectively, and covariance matrices CovU and CovU Ξ . For latter considerations, we also introduce the auto‐covariance of Ξ:CovΞ=E[(Ξ−Ξ̲)(Ξ−Ξ̲)T] that is known because the stochastic model of K is given. All expectations and covariance matrices are estimated using standard statistical estimators and the original realizations of the solution of the primal Arlequin problem .

We introduce a unitary centered Gaussian random vector Θ with uncorrelated components ( E[θiθj] = δij) and independent from Ξ, and the vector

Û=U̲+CovUΞCovΞ−1(Ξ−Ξ̲)+CovU−CovUΞCovΞ−1CovUΞTΘ,

where the square root sign for a matrix C indicates that CCT=C. Simple algebra (remembering the independence of Θ and Ξ and E[Θ] = 0) indicates that this random vector indeed verifies the desired first‐order and second‐order statistics  sought.

As neither K(x,θ) nor usa(x,θ) are Gaussian random fields, an isoprobabilistic approach should be followed to transform the original fields into Gaussian fields. After construction of Gaussian vector Û, the inverse approach can be used to transform new realizations of U into the original first‐order marginal distribution.

With the procedure described in this section, it is possible to generate new realizations of uarla that correspond to both K and parla. Hence, it allows us to estimate A(uarla,parla) through Monte Carlo sampling, and then η by using Equation .

The error introduced by the use of the surrogate problem can stem from several sources. In fact, if we decompose the construction of the approximate solution into different steps, we can distinguish (Figure [NaN] ) the following:

the use of a surrogate problem described by the Arlequin method,

the spatial discretization of continuum models,

and the approximation of the statistical moments (by using the Monte Carlo technique for example).

Each step is driven by specific parameters. The Arlequin model is described by the size of the patch and that of the coupling zone. The spatial discretization is parametrized by the element size on each domain ( Ωd, Ωs, and Ωc). Finally, if we use the Monte Carlo technique to deal with the random dimension, the third step is driven by the number of Monte Carlo samples used to describe the material property. We propose here a methodology to assess the error induced by each step and separate η into corresponding error indicators.

To split different error sources, we have to consider intermediate problems. In the following, the term ‘continuous’ refers to the opposite of ‘discretized’ (by using FEM) or ‘approximate’ (by using FEM and Monte Carlo technique). We call the following:

Reference problem, the monomodel stochastic continuous problem (whose the solution is denoted by uex and associated residual function R was introduced in ). To estimate the residual R(uarl,parl), we follow the approach described in [NaN] bu using a fine mesh and a large number of Monte Carlo draws.

First intermediate problem, the Arlequin problem, defined by  and Figure [NaN] , coupling a deterministic continuous and a stochastic continuous models. We denote the solutions (ud,us,λ)∈Vdh×Wsh×Wch, the associated monomodel solution uarlc, and the associated residual function Rc. The associated residual function Rc((.,.,.);(.,.)) is then defined for all (vd,vs) as follows:Rc(ud,us,λ);(vd,vs)=ℓ(vd,vs)−a((ud,us,λ),(vd,vs)),where ℓ(vd,vs) = ℓd(vd) + Ls(vs) and a((ud,us,λ),(vd,vs,μ)) = ad(ud,vd) + As(us,vs) + C(λ,vd − vs). As previously, we follow the approach described previously to evaluate the residual. As this problem refers to a continuous Arlequin model, the mesh used has to be sufficiently fine, and the number of Monte Carlo draws sufficiently large.

Second intermediate problem, the Arlequin problem coupling a deterministic discretized and a stochastic spatially discretized models. We denote the solution (udh,ush,λh)∈Vdh×Wsh×Wch, the associated monomodel solution uarlh, and the associated residual function Rh. This problem is defined from the previous one by using a discretization technique along the space dimension, by means of the FEM, for instance. The associated residual function Rh((.,.,.);(.,.)) is also defined for all (vdh,vsh) as follows:Rhudh,ush,λh;vdh,vsh=ℓhvdh,vsh−ahudh,ush,λh,vdh,vsh,where ℓh and ah correspond to discretized forms of ℓ and a. Contrary to the previous model, this residual refers to a discretized model. Therefore, to evaluate this residual, we use the same mesh as the one used for the computable problem, a coarse one. Nevertheless, the stochastic space is not approximated yet and we thus have to use a large number of Monte Carlo draws.

Computable problem, the Arlequin problem coupling a deterministic and a stochastic models after discretization in space and solving by the Monte Carlo technique in the random dimension. The associated solution (uda,usa,λa) with monomodel solution uarl is the only one we compute in practice. We will see how this solution can be used to estimate different error sources.

We can notice that

η=q(uex)−q(uarl)=q(uex)−q(uarlc)+q(uarlc)−q(uarlh)+q(uarlh)−q(uarl)

The different errors on q are then defined as follows:

the stochastic error (mainly due to the approximation of the statistical moments):ηθ=q(uarlh)−q(uarl),

the discretization error (mainly due to the discretization along the spatial dimension):ηh=q(uarlc)−q(uarlh),

and the modeling error induced by the use of the Arlequin method (by use of a patch and by use of a deterministic material parameter instead of the stochastic one in Ωd):ηm=q(uex)−q(uarlc).

As the only data available are the solution of the primal problem (uda,usa,λa), the solution of the adjoint problem (pda,psa,pλa), and the model (residual function) associated to each intermediate problem, we can derive estimates of the previous errors by using :

for the stochastic errorηθ≃RhΠhuda,Πhusa,Πhλa;Πhpda,Πhpsa,

for the discretization errorηh=q(uarlc)−q(uarl)−q(uarlh)−q(uarl)≃RcΠcuda,Πcusa,Πcλa;Πcpda,Πcpsa−RhΠhuda,Πhusa,Πhλa;Πhpda,Πhpsa,

and for the modeling errorηm=q(uex)−q(uarl)−q(uarlc)−q(uarl)≃R(uarl,parl)−RcΠcuda,Πcusa,Πcλa;Πcpda,Πcpsa,

where Π(c,h)uda,Π(c,h)usa,Π(c,h)λa and Π(c,h)pda,Π(c,h)psa are projections of (uda,usa,λa) and (pda,psa), respectively, in Vd(c,h),Ws(c,h),Wc(c,h).

For the spatial projection, linear interpolation is chosen, whereas the stochastic part is treated following the approach described in [NaN] . As the patch has the same size in all the intermediate problems, the projection of the primal fields is obvious. As the patch used for the adjoint problem is bigger (Ωs⊂Ω̃s), the stochastic adjoint is restricted to the patch used of the primal model. For instance, the fields Πhpda and Πhpsa are obtained from (pda,psa) by using the following:

Πhpda=pdaΠhpsa=psa|Ωs.

Let us remind that Rc(⋅,⋅) is computed using a very fine mesh, whereas Rh(⋅,⋅) is estimated using the same mesh as for the primal model.

The quantities ηθ, ηh, and ηm are only error indicators as the adjoint problem is approximated using the Arlequin solution. Let us notice that the different indicators can be estimated using only solutions uarl and parl (direct post processing), avoiding additional cost. This is in opposition with other possible approaches in which intermediate adjoint problems would be involved and solved.

The basic adaptation process, based on a greedy algorithm, is described in Figure [NaN] . The approach consists in the following steps:

A set of initial parameters is used to perform primal and associated adjoint problems that are both of Arlequin type.

The total error η associated with the quantity of interest is estimated.

If this error is higher than a given precision, we follow the technique described previously to split the error sources.

The dominating error source is identified (in modulus), and the corresponding parameter is refined.

We continue the process until the total error is below the required precision (in modulus) or when the parameter, which has to be refined, is already as fine as possible.

This choice of adaptation strategy is not unique; indeed, step 4 could be changed and we could identify two (or more) dominating errors and refine the corresponding parameters.

In this section, we present some numerical results. [NaN] and [NaN] discuss several numerical issues, whereas [NaN] and [NaN] present full applications (in 1D and 2D) that illustrate the theoretical results of the paper. [NaN] explores the convergence of the error estimate proposed with respect to several key parameters and compares that error estimate with the real error. [NaN] elaborates on the importance of adequately solving the adjoint problem for the error estimate to be precise. Finally, in [NaN] , [NaN] , and [NaN] , full 1D and 2D cases are treated to assess the methodology proposed in this paper, including the issue of splitting error sources to drive adaptation. All the numerical examples of this section have been computed using the software CArl (Code Arlequin) [42] .

We study here the influence of various parameters of the Arlequin strategy on a simple 1D problem (Figure [NaN] ). The problem consists in a bar of unit length under traction loading with prescribed Dirichlet conditions, as well as a unit bulk load. The material property is random and modeled as a uniform field with bounds 0.3013 and 2.3601 (arithmetic average 1.3307, harmonic mean 1 ∕ E [1 ∕ K(x,θ)] = 1, and standard deviation σK = 0.2), and exponential correlation with correlation length Lcor = 0.01.

It is approximated by an Arlequin model defined in Figure [NaN] where (αd,αs) are chosen linear in the coupling zone Ωc. The model is driven by several parameters: Ls and Lc for the definition of the patch, hd and hs for the size of the spatial discretization (for domains Ωd and Ωs, respectively) and MC for the number of samples used for the Monte Carlo description. The spatial discretization of Ωc is driven by the same size hs as for Ωs.

We study here two quantities of interest:

One related to the mean of the gradient in a zone near the middle of the bar:q∇(u)=10.05E∫0.450.5dudx(x)dx.The associated extractor, by using integration by parts, is defined by the following (Figure [NaN] ):q∇(v)=10.05E∫Ω{δ(x−0.5)−δ(x−0.45)}v(x)dΩ,

One related to the variance of the gradient mean near the middle of the bar (which is nonlinear with respect to u):qv(u)=Var10.05∫0.450.5dudx(x)dx.

Following the idea developed in , for the nonlinear quantity qv(u), we can write the following (Figure [NaN] ):

qv(u)−qv(uarl)≈20.052EΔeΔ(uarl)−EΔ(uarl)

where Δe = e(x = 0.5) − e(x = 0.45) and Δuarl=uarl(x=0.5)−uarl(x=0.45). The extractor is therefore defined by the following:

qv(v)=20.052E∫Ωδ(x−0.5)−δ(x−0.45)Δuarl−EΔuarlv(x)dΩ.

We now investigate the evolution of the relative error indicator:

ηr=R(uarl,pa)q(uarl),

where pa is the approximate solution of the adjoint problem (by using the FEM and Monte Carlo techniques for instance). This error indicator (drawn in dashed line and circles) will be compared for each of the two quantities of interest with the true relative error (drawn in solid line):

ηref=q(uref)−q(uarl)q(uarl).

where q(uref) is evaluated using a reference monomodel whose parameters are listed and computed using the FEM and the Monte Carlo technique. The large number of Monte Carlo draws used and the fine spatial discretization allow to neglect the error (the standard deviation of the true relative errors, taking 20 different sequences of 100,000 Monte Carlo draws is less than 10 − 3 for q ∇  and less than 5 × 10 − 3 for qv). In the following subsections, we study the evolution of the error related to q ∇ (u) (results in [NaN] , associated adjoint problem described in Figure [NaN] ), and qv(u) (results in [NaN] , associated adjoint problem described in Figure [NaN] ), with respect to two parameters: the number MC of Monte Carlo draws used for the primal model (left result in figures) and the half‐size of the stochastic patch Ls (right result in figures). As q(uarl) depends on the sequence of Monte Carlo draws, we show results for 20 different sequences of draws. The set of parameters used for the reference model, the primal model, and the adjoint model are defined in Table [NaN] and Table [NaN] . In these tables, when a monomodel is used, we give also the spatial discretization mesh size h.

Definition of model parameters used for the study of q  ∇  and qv with respect to MC. Twenty different sequences of Monte Carlo draws are used for the primal model.

Definition of model parameters used for the study of q  ∇  and qv with respect to Ls. Twenty different sequences of Monte Carlo draws are used for the primal model.

In Figure [NaN] , as the adjoint model is solved with the same refinement as for the reference model, the estimated error is superposed to the true relative error. As expected, the error decreases in both cases. On Figure [NaN] (right), when the half‐size of the stochastic patch evolves from 0.2 to 0.5, the error decreases from almost 2% down to less than 0.5%. On Figure [NaN] (left), with more than 500 Monte Carlo draws, for a half‐size of 0.3, the error is lower than 5%. Nevertheless, we remark that, when the number of Monte Carlo draws increases, the error does not tend to zero. This comes from the fact that the primal model is approximated using Monte Carlo draws but also by using the Arlequin method. That means that both Monte Carlo draws and half‐size of the stochastic patch (Ls = 0.3) contribute to the total error. Indeed, the asymptotic value of the error in the left picture corresponds to an error close to 1.4% for Ls = 0.3 in Figure [NaN] (right).

In Figure [NaN] , we remark that the true error and the estimated error decrease. Nevertheless, the error estimate is very optimistic (lower in modulus) for the evolution with respect to Ls and quite pessimistic (higher in modulus) for the evolution with respect to MC. Contrary to the cases involving q ∇ , here, the main difference is that the adjoint problem cannot be solved perfectly as its load is defined from the primal displacement solution. Moreover, the difference between the true error and the estimated one is also due to the linearization of the quantity of interest in . Again, a key point here is that the true error, in the picture on the left, seems to tend to a significant nonzero value. For the same reasons as those seen previously, it is because both Monte Carlo draws and half‐size of the stochastic patch contribute to the total error. We then propose in the next section to study the evolution of the error when we take an Arlequin model with the same half‐size of the stochastic patch as the reference.

We can remark that the parameter MC studied before drives the approximation of the stochastic space. If we follow the splitting description of the total process, we can observe that the stochastic approximation process occurs between two spatially discretized Arlequin models. That explains the fact that when we study the evolution of the error with respect to MC, we also evaluate the error induced by the use of the Arlequin method. To avoid that, we propose here to study the evolution of the error with respect to a reference model that is an Arlequin model with the same coupling and the same spatial mesh step as described in Table [NaN] .

Definition of models used for the study of q  ∇  or qv with respect to MC with Arlequin model as the reference. Twenty different sequences of Monte Carlo draws are used for the primal model.

Figure [NaN] shows the evolution of the true relative error ηref (solid line) and the error estimate ηr (dashed line and circles) with respect to the number of Monte Carlo draws used in the primal problem, for the evaluations of q ∇  (left) and qv (right). We remark in this case that we efficiently capture the error related to the parameter MC alone. This time, the true errors actually decrease down to a value close to zero. In fact, the last point proposed here corresponds to the Arlequin reference model used in the example. These curves correspond to a translation of the curves in Figures [NaN] and [NaN] (left) of a value corresponding to the error when a patch of half‐size Ls = 0.3 is used.

Definition of models used for the calibration of the adjoint problem (parameter Ls ) associated to q  ∇  ( u ). Twenty different sequences of Monte Carlo draws are used for the adjoint model.

For the quantities of interest studied here, the size of the stochastic patch or the number of Monte Carlo draws used for the primal problem have a strong influence on the error. When the adjoint model is solved with sufficient accuracy, the estimated error perfectly matches the true error in the case of the linear quantity of interest q ∇ . For the estimation of qv, the difference between the errors is because q(uex)−q(uarl) is estimated by q ′ (uarl,e) that corresponds to the first‐order term only (see ). Finally, [NaN] showed that the error for high MC and large Ls comes from the use of an Arlequin surrogate model for the solution of the adjoint problem.

In [NaN] , we propose to study the influence of the enrichment of the adjoint model, compared with the primal model one, when an approximate adjoint model is used to estimate the quantity of interest q ∇ .

As seen previously, in practice, the adjoint problem will also be approximated by the Arlequin method. Not to increase too much the computation time, parameters of the adjoint problem have to be set up correctly. In this part, we study the effectivity of the error estimates when adjoint model parameters vary. The parameters under study are as follows:

the half‐size of the stochastic patch (Lsa) used for the adjoint model (when the corresponding parameter Lsp varies in the primal model),

and the number of Monte Carlo draws (MCa) used for the adjoint model (when the corresponding parameter MCp varies in the primal model).

The aim is to get a feeling of how to enrich the adjoint model compared with the primal model to get an accurate estimation of the true error.

Table [NaN] shows the evolution of ηr and ηref with respect to the half‐size of the stochastic patch used for the primal model (Lsp). In this case, the increase of the size of the stochastic patch for the adjoint model does not seem to have a strong influence on the accuracy of the error estimate. We can remark a worsening of this accuracy when the size of the stochastic patch for the primal model increases. We can also notice that with a half‐size of patch of 0.5 for the adjoint model, the difference between the true error and the estimated error is low but still significant. That means that even with the largest patch, the adjoint problem is not solved as accurately as in the previous section. This is mainly because the coupling zone constrains the stochastic solution in a nonzero part of the problem.

Evolution of ηr [10  − 3 ] with respect to the half‐size of the stochastic patch of the primal model ( L sp ), for the evaluation of q  ∇  ( u ) with models defined by Table  . Standard deviations are also given between parenthesis.

For the adaptive examples described in [NaN] , we choose to set the patch size of the adjoint model with only one value (Lsa = 0.45) independently of the patch size of the primal model (Lsp).

Table [NaN] shows the evolution of ηr and ηref with respect to the number of Monte Carlo draws used for the primal model. We remark again that the estimation is only slightly modified when we take a large number of Monte Carlo draws for the adjoint model. In Table [NaN] , for MCp = MCa = 100,000, the relative estimated error 10.8 × 10 − 4 can also be read in Table [NaN] for Ls = 0.3 for the primal model and between Ls = 0.4 and Ls = 0.5 for the adjoint model. For the adaptive examples described in [NaN] , the number of Monte Carlo draws for the adjoint model (MCa) is set to 100,000.

Definition of models used for the calibration of the adjoint problem (parameter MC ) associated to q  ∇ . Twenty different sequences of Monte Carlo draws are used for the adjoint model.

Evolution of ηr [10  − 3 ] with respect to the number of Monte Carlo draws used for the primal model ( MCp ), for the evaluation of q  ∇  with models defined by Table  . Standard deviations are also given between parenthesis.

As soon as the Arlequin model is used for the adjoint problem, the estimation of the error loses efficiency, specially concerning the number of Monte Carlo draws used for the adjoint model. Nevertheless, when the stochastic patch is long enough and when the number of Monte Carlo draws is sufficiently large, the difference between the reference error and the estimated error is quite small. In these cases, the estimated error tends to underestimate the true error.

In regard to the number of Monte Carlo draws used for the adjoint model, it seems that the estimated error is sufficiently accurate when this number is very large. For instance, if we take 100,000 Monte Carlo draws, the estimated error is very accurate as the patch size of the adjoint problem (Lsa = 0.45) still affects the approximation of the adjoint solution and therefore the quality of the estimation.

Finally, for model adaptation, we choose to approximate the adjoint problem by using the Arlequin method but with a large patch size (Ls = 0.45), and a large number of Monte Carlo trials (MC = 100,000). The other parameters are set to the following: hd = 0.01 and hs = 0.002. This choice is non‐unique and can be discussed and changed if needed (in the case where we have to refine the mesh size for example).

We study in this section the problem described in [NaN] . In particular, we present results from the adaptation process obtained when we investigate the two quantities of interest q ∇  and qv. The parameters are Ls (half‐size of the stochastic patch), MC (number of Monte Carlo draws used), the spatial mesh size hd of the deterministic model, whereas the spatial mesh size hs of the stochastic model is fixed at 0.002.

In Figure [NaN] and Table [NaN] , we show the results for the error associated with q ∇ . We observe that the main error sources are due to the use of a low number of Monte Carlo draws, then to the use of the Arlequin method. By increasing the corresponding parameters (Ls from 0.2 to 0.45 and MC from 5 to 2000), the relative estimate error ηr decreases from 8.7% to 6.6%. The corresponding true error decreases from 14% to 0.7%. More precisely, the algorithm stops after iteration 3 and the dominant error is due to the use of the Arlequin model. The last iteration is added manually to illustrate the possibility to decrease the total error by increasing the number of Monte Carlo draws (the second main source of the error) and therefore decreasing the stochastic error. Table [NaN] also shows the evolution of the true errors with respect to the patch size, the discretization and the Monte Carlo truncation. As we can see, the various error estimates (ηm, ηh, ηθ) give a good indication on the true errors, specially for ηθ and ηh. The bias on the estimation of the total error impacts only the estimation of the modeling error.

Evolution of parameters and relative errors (true and estimate) for the evaluation of q  ∇ .

In Figure [NaN] and Table [NaN] , we study the error associated to qv. In this case, the study shows that the main error sources are due to the following: (i) the use of a low number of Monte Carlo realizations; (ii) the size of the stochastic patch; and (iii) finally (for the last iteration) the size of the mesh of the deterministic model. By increasing the corresponding parameters alternately, the relative estimate error ηr decreases from more than 86% to less than 0.5%.

Evolution of parameters and relative error for the evaluation of qv.

This simple example shows the efficiency of the splitting error sources methodology. Indeed, the stochastic error and the discretization error correspond to the true errors as soon as the second iteration. The approximation of the adjoint problem by using an Arlequin model mainly affects the modeling error estimate. Finally, the use of this technique permits to decrease efficiently the error estimate by adapting the corresponding parameter for each step of the process.

We now consider a 2D sample Ω inscribed in the box [ − 3,3] × [ − 1,1] (Figure [NaN] ). The sample is submitted to a prescribed Dirichlet condition, with no bulk load. The boundary conditions are u(x = − 3,y) = 0,u(x = 3,y) = 1, and  ∇ u ⋅ n = 0 for the remaining edges, almost surely. The model is described by a random material property K(x,θ), modeled as a uniform field with bounds 0.3542 and 2.1938 (with geometric mean 1 ∕ E [1 ∕ K(x,θ)] = 1, and standard deviation σK = 0.2), and exponential correlation with correlation length Lcor = 0.05 in each direction. To estimate the residual numerically, the reference monomodel problem is spatially discretized by 33,792 elements, and we use MC = 5000 Monte Carlo draws to represent the stochastic field K.

The quantity of interest considered here is the component x of the gradient in a given zone Ωint. Considering the model given in Figure [NaN] , Ωint is located near the middle of the sample (where the loading is applied). This quantity of interest is defined by the following:

q2D(u)=E∫Ωint∇u⋅ixdΩ.

The reference problem is approximated by the Arlequin method with a centered patch. The corresponding adjoint problem is described in Figure [NaN] . It is loaded by pq = K(x,θ) ∇ u(x,θ) = ix with ix the unit vector of the x‐axis. It is spatially discretized by 8448 elements for the deterministic model and 6144 for the stochastic model in a patch of half‐size Ls = 1.8. Moreover, we use 5000 Monte Carlo draws for the adjoint model.

By using the adaptive strategy introduced in [NaN] , we investigate the absolute value of ηr :

ηr=R(uarl,pa)q2D(uarl).

Figure [NaN] and Table [NaN] show that the main error sources are due to the use of a weak number of Monte Carlo draws and due to the use of the Arlequin method. By increasing the corresponding parameters, the relative estimate error ηr decreases from 21% to 5%.

Evolution of parameters and relative errors for the evaluation of q 2D. Nd and Ns represent the number of elements of the deterministic discretized model and of the stochastic patch, respectively. Note that Ns increases only because Ls increases, the patch mesh is not finer.

As for the 1D adaptation example, the splitting of error sources technique permits to identify the dominating error source. The associated parameter can then be refined. For this 2D example, the Arlequin model with a relatively small number of Monte Carlo draws, and a medium‐sized patch, gives an estimation of the quantity of interest with less than 6% of error.

In dental restoration, dentists use bio‐compatible resin to replace the ill part of a tooth (damaged by caries for example). The adhesion of this resin with the tooth is a key issue of the treatment.

When caries reach the dentin, the adhesion is currently achieved by micromechanical seal of the resin with the demineralized dentin. The resin has to sufficiently infiltrate the dentin to obtain an effective adhesion. We study in this part the infiltration speed of a resin in the demineralized dentin. The studied problem is described in Figure [NaN] . The resin infiltrates first the tubule (represented by " holes" in the figure) and then a network of collagen fibers that is the main component of a demineralized dentin. The porosity of this medium is modeled by a stochastic field.

We assume that the flow follows the Darcy assumptions in a permanent regime. The pressure verifies the equation  where the random material property K(x,θ) is modeled as a uniform field with bounds [1.43 × 10 − 6,1.43 × 10 − 4]μm2(Pa ⋅ s) − 1 (with geometric mean 1 ∕ E [1 ∕ K(x,θ)] = 3.19 × 10 − 5μm2(Pa ⋅ s) − 1), and exponential correlation with correlation length Lcor = 1μm in each direction. Only one of the tubule is filled of resin, whereas the other tubule does not. This condition is modeled by a prescribed pressure of 1 bar on the bottom tubule (which corresponds to the pressure applied manually by the dentist), whereas the other edges are submitted to an over pressure of 0 bar. The problem is described using the Arlequin approach in which a square patch, centered on the sample, is used to model the tubules and the stochastic field, whereas the substrate does not contain tubules and heterogeneities. With these considerations, the smallest patch is a square of 18μm side length. We seek to adapt the different parameters of the method by considering as quantity of interest the average flow velocity along the y‐axis in a specific zone defined by Zint = {M(x,y) | (x,y) ∈ [ − 1 ∕ 3,1 ∕ 3] × [2 − 1 ∕ 3,2 + 1 ∕ 3]}:

qf(p)=E−1|Zint|∫ZintK(x,θ)∇p(x,θ)⋅iydΩ

The results (Figure [NaN] ) show that the flow is highly localized between the two tubules. Following the approach described on this paper, the quantity qf can be evaluated with an estimated accuracy of less than 0.2% for the absolute value of the total error ( 0.06% for the modeling error, − 0.18% for the stochastic error, and 0.001% for the discretization error) reached at the first increment.

The main contribution of the paper has been to propose a computational framework for addressing the control of errors coming from three different sources:

the modeling error, due to the use of the Arlequin method coupling a deterministic model with a stochastic one to approximate a full stochastic problem,

the discretization error, due to the discretization of the spatial dimension, by using the FEM for instance,

and the stochastic error, due to the use of the Monte Carlo technique to approximate statistical moments.

A goal‐oriented error estimation technique was introduced to quantify the capability of the Arlequin framework to evaluate some specific quantities of interest. Adaptivity was also considered. For that purpose, the corresponding error strategy was enriched to identify separately the different error sources. By using only the solution of the primal problem, the solution of the adjoint problem, and the definition of the residual associated to specific intermediate problems, we showed that it is possible to identify the major error source and refine the corresponding parameter to reduce the total error.

Future research will deal with the development of an optimal adaptive process where the error could be spatially analyzed, on one hand, and the coupling between two stochastic models, on the other hand.

This work was partially supported by the ANR project TYCHE (Advanced methods using stochastic modeling in high dimension for uncertainty modeling, quantification and propagation in computational mechanics of solids and fluids), with project number ANR‐2010‐BLAN‐0904, and by grants from DIGITEO and Région Ile‐de‐France, with project number 2009‐26D.

Graph: Reference stochastic monomodel with heterogeneous random coefficient K ( x , θ ).

Graph: Arlequin description of the problem with coupling between the two models in zone Ω c.

Graph: Approximation steps for the resolution of a 1D stochastic problem.

Graph: Greedy algorithm for the adaptation process.

Graph: 1D reference problem of a cantilever beam.

Graph: 1D surrogate model with the Arlequin method.

Graph: 1D adjoint monomodel associated to q ∇  ( u ).

Graph: 1D adjoint monomodel associated to qv ( u ) where Fv=20.052Δuarl−EΔuarl.

Graph: Evolution of relative errors for the evaluation of q ∇  with primal models defined by Table  and Table  : ηr and ηref (superposed) in gray markers for several sequences of draws with mean in solid line and standard deviations (error bars), with analytical rate of 1∕MC in dashed line with respect to the number of Monte Carlo draws [left] and with respect to the half‐size of the stochastic patch [right].

Graph: Evolution of relative errors for the evaluation of qv with primal models defined by Tables  and  : ηr in gray circles for several sequences of draws, with mean (solid line) and standard deviations (error bars) and of ηref in gray crosses for several sequences of draws with mean (dashed line) and standard deviations (error bars), with respect to the number of Monte Carlo draws [left] and with respect to the half‐size of the stochastic patch [right].

Graph: Evolution of relative errors for the evaluation of q ∇  [left] and qv [right] with models defined by Table  : ηr (gray circles) and ηref (gray crosses) for several sequences with mean in solid line and standard deviations (error bars), analytical rate of 1∕MC in dashed line.

Graph: Evolution of the relative estimates of the total error ηr (white), modeling error ηm (light gray), discretization error ηh (dark gray), and stochastic error ηθ (black) during adaptive process for the estimation of q ∇ .

Graph: Evolution of the relative estimates of the total error ηr (white), modeling error ηm (light gray), discretization error ηh (dark gray), and stochastic error ηθ (black) during adaptive process for the estimation of qv.

Graph: 2D case and associated mesh used for the estimation of the residual R.

Graph: 2D associated adjoint problem.

Graph: Evolution of relative estimates of the total error ηr (white), modeling error ηm (light gray), discretization error ηh (dark gray), and stochastic error ηθ (black) during adaptive process for the estimation of q2D.

Graph: Study of resin infiltration in demineralized dentin. Two dentinal tubules of diameter 3 μm are modeled.

Graph: Mean value of y ‐component of the flow velocity [ μm  ⋅  s − 1 ] in an Arlequin patch (size: 18 μm ) obtained from 100 Monte Carlo realizations and with a triangular mesh (whose the maximum diameter of the inscribed circle is 0.05 μm (corresponding to a total estimated error of − 0.2 % for the evaluation of qf ). The boundary of the substrate is visible in solid lines.