J P Rouse

P Kerfriden

M Hamadi

Keywords: Multiscale, Finite Element Methods, Error Estimation, Elasticity

come

A probabilistic hierarchical sub-modelling approach through a posteriori Bayesian state estimation of finite element error fields

A probabilistic hierarchical sub-modelling approach through a posteriori Bayesian state estimation of finite element error fields.

Introduction

Structural analysis, utilised to answer many key questions in the design process of complex engineering systems, often requires the abstraction of features at many different length scales. Without some level of length scale limit in model development practical analysis problems would quickly become intractable and solutions would not be furnished at a rate required by the iterative design process. Geometric features, material property definitions, and discretisation particulars (i.e. element formulations and refinement) will commonly be varied over the different length scales. Broadly speaking, large, coarse ("macro scale") models will be solved to address questions relating to "general" stiffness of a structure and to inform unknown driving boundary conditions for smaller scale, detailed ("micro scale") models that represent localisations where gradient quantities (likely stresses and strains) are of interest. Much has been written about the various approaches adopted for these multiscale problems (a brief summary is given later). All involve the transfer for information (stiffness and boundary conditions in most cases) between macro and micro models, either in a strongly coupled sense (where information flow is bidirectional) or in an weakly coupled sense (where information flow is unidirectional). The present work concerns itself mainly with weakly coupled multiscale problems (however extensions to strongly coupled problems are clear) and asks what are the ramifications of "error" in macro models? Bearing in mind the practical limitations associated with large scale analyses, how can error estimates in large scale models be generated quickly and used to their greatest effect, such that uncertainties can be propagated to local micro models? In the present work, we motivate, propose, and demonstrate a full error field recovery approach that utilises a small number of goal orientated error estimate (GOEE) observations in combination with a Gaussian Process state estimation framework. Importantly, a guided model order reduction technique, enabled by stochastic partial differential equation (PDE) representations of the prior, is proposed that alleviates computational bottlenecks in the Gaussian process steps, thereby expediting full error field estimation (in line with industry requirements). Full field error estimation, that is to say the estimation of a continuous error field over an entire domain rather than at several discrete locations, is now motivated through a simple aerospace example. Consider a large airframe assembly, such as the fuselage "barrel" section of civil aircraft. It is clear that the real assembly will consist of a multitude of components, from large skin, stringer, and stabiliser components that are vital in the approximation of stiffness, to small features such as holes, rivets, and fillet sections, which are prime failure initiation sites that must be given careful consideration if a safe design is to be produced. An approximation of global stiffness of the barrel section may well be required in order to anticipate the outcome of wing bend tests (among others). Clearly, small scale features (holes, rivets, and through thickness sections, for example) add little to global stiffness, justifying their omission from large scale models. Micro models, where small scale features (holes and composite material laminates, for example) are explicitly represented, may be used to answer concerns over the potential for premature failure. Let us consider a hole feature in a laminate material section by way of example. Discretisation is refined in this micro models to provide better approximation of terms related to gradients in the displacement field. The hole micro model will, of course, be much smaller than the macro barrel model, and it is highly unlikely that any loads or boundary conditions will be known for the hole model a priori, meaning that they will need to be extracted from the barrel model. One may ask, if the barrel model is coarse (and presumably relatively inaccurate in terms of local deformation fields), can the hole model which is driven by it be trusted? Confidence in micro model results could (conceptually) be gained by performing error estimation around the region of interest (RoI) in the barrel model and propagating the uncertainty in driving degrees of freedom (DoFs) to the hold model. This observation alone does not prompt full error field estimation, as it assumes the RoI in the barrel model is well defined. In reality, there will be many hole features in the barrel model at a multitude of locations and analysts may well be unsure as to which holes are critically important. Performing error analysis on all RoI boundaries would be both prohibitively time consuming and wasteful. Each GOEE observation demands its own simulation and, in many cases, observations that are close to one another will, conceivably, be very similar. The additional computation effort associated with this simplistic error field recovery process greatly limits its applicability and the method is not sensitive to the iterative design process. In practice, hole model positions and geometries may change as the design progresses (as refine cable routes are determined, for example). Given that the holes are not represented in the barrel model, it is justifiable to drive the new design hole models from the original barrel model. A brute force approach would be to repeat error estimation at the new RoI boundaries, however this is unsatisfying and impractical. In the present work, we resolve this challenge by considering potential RoI boundaries (to drive the placement of error estimate observations), however we recover a full error field, statistically with confidence intervals, such that any potential RoI boundary or micro model can be considered. Further discussion on the relationship between macro and micro models in aerospace design can be found in the authors' previous work [START_REF] Bonney | Goal oriented error estimation in multi-scale shell element finite element problems[END_REF]. A treatment of error estimation in mixed element (i.e. "shell" macro model to "continuum" micro model) multiscale problems is also presented in the author's previous work, however readers should note that this does not impact the novelty of the present work. Previous work was concerned with error estimation for a particular set of element formulations, whereas the present work takes a small set of these observations and considers how a stress analyst could use them effectively in practical multiscale problems.

The present work considers error estimation in finite element (FE) problems, as FE is by far the most commonly employed approach for engineering stress analysis. Approximation errors are the inaccuracies which are inherent to the discretisation methods that are required in order to approximate the solutions to mathematical models. This contrasts with modelling error, which is a measure of how well an abstract model approximates real physical phenomena. Estimates of approximation error may be termed a priori estimates or a posteriori estimates [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF], with the former utilising the problem definition and the discretisation to estimate error and the latter using the solution approximation itself to estimate error. A further distinction may also be made between the approximation of error bounds and the estimation of error itself [START_REF] Grtsch | A posteriori error estimation techniques in practical finite element analysis[END_REF]. The former can be guaranteed but may be inaccurate. Error bounds may be large, for example, although it is important to note that this is not always the case and bounds may be sufficiently narrow that they function as estimates. The latter is usually not guaranteed but the range of values associated with it will typically be narrower than for error bounds. Grätch and Blathe note that error estimates should possess several properties [START_REF] Grtsch | A posteriori error estimation techniques in practical finite element analysis[END_REF]. In addition to accuracy, error estimates should asymptotically tend to zero as the discretisation density increases, produce tight bounds for the error, be computationally inexpensive, be applicable in a wide range of potentially nonlinear applications, and inform mesh refinements such that approximate solution processes can be optimised. At present, no one estimator can satisfy all of these requirements however, for the applications considered in the present work, several relevant methods can be found in the literature (as will be discussed later). Multiscale modelling methods allow for the transfer of information across many length and time scales in structural analysis problems [START_REF] Marco | Multiscale modeling: A Bayesian perspective[END_REF][START_REF] Weinan | Principles of multiscale modeling[END_REF]. They are particularly effective in cases where there are multiple physical phenomena of interest (that span a range of length/time scales) and that cannot be practically rationalised in a single model. Multiscale methods can be viewed as a set of approaches which homogenise the heterogeneity observed in all real structures at some length scale. Multiscale modelling techniques have been applied to a wide variety of fields, including aerospace [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF][START_REF] Said | An iterative multiscale modelling approach for nonlinear analysis of 3D composites[END_REF][START_REF] Sturm | Multiscale modeling methods for analysis of failure modes in foldcore sandwich panels[END_REF], marine [START_REF] Lua | Multi-scale dynamic failure prediction tool for marine composite structures[END_REF], and civil [START_REF] Addessi | A multi-scale enriched model for the analysis of masonry panels[END_REF].

Many different multiscale methods can be found in the literature and several characterisation schemes have been proposed. Weinan bisects the field and categorises multiscale problems as either Type A, in which there are local defects/singularities that require a locally micro-scale model in an otherwise coarser global model, or Type B, which have micro-scale features throughout and require fine-scale modelling everywhere (with some form of computational homogenisation, for example) [START_REF] Weinan | Principles of multiscale modeling[END_REF]. Other authors, such as Geers in [START_REF] Marc | Homogenization methods and multiscale modeling: nonlinear problems[END_REF], have used similar categorisations but under different names. In [START_REF] Marc | Homogenization methods and multiscale modeling: nonlinear problems[END_REF], Type A methods are called "Hierarchical" methods while Type B methods are called "Concurrent" methods. There is inevitably some overlap between the techniques applicable to Type A and Type B problems, for example Kim and Swan used adaptive refinement of voxel meshes of representative volume elements within their numerical homogenisation approaches [START_REF] Hyung | Voxel-based meshing and unit-cell analysis of textile composites[END_REF]. Type A approaches include "classical" sub-modelling [START_REF] Gendre | A two-scale approximation of the schur complement and its use for non-intrusive coupling[END_REF], domain decomposition [START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF][START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF], and local mesh refinement including adaptive mesh refinement. The procedures used within the CleanSky2 project MARQUESS, wherein a pre-computed database of solutions enables a single sub-model to be applied to multiple instances of a recurring feature on a global FE model of a composite component, are particularly relevant to the present work [START_REF] Zou | A computationally efficient approach for analysing the onset of failure in aerospace composite structures[END_REF]. In the apporach, pre-computed sub-models are utilised in a bottom-up and top-down approach to identify local failure locations [START_REF] Zou | The identification of failure initiation hotspots in idealised composite material component models using a bottom-up database method[END_REF][START_REF] Zou | Comparison of interface modelling strategies for predicting delamination in composite l-angle sections under four-point bending. composite structures[END_REF][START_REF] Bonney | Bayesian reconstruction of goal orientated error fields in large aerospace finite element models[END_REF][START_REF] Bonney | Goal oriented error estimation in multi-scale shell element finite element problems[END_REF]. The methodology is particularly useful as the superposition principle used enables rapid approximation of micro model states. Displacements from the global model can be extracted and used to estimate local stress fields in the micro model by weighting pre-computed solutions that represent particular deformation modes. Type A sub-modelling approaches are the primary concern of the present work. As will be shown later, full error fields will be developed across a macro model that describe uncertainties in driving DoFs. These distributions will in turn be sampled to develop possible loading conditions for local micro models.

An example of a Type B problem is the multiscale modelling of complex fibre architectures, such as 3D woven textile composites. Many techniques can be identified for tackling this type of problem. Computational homogenisation, for example, attempts to determine equivalent properties via representative volume elements within a periodic displacement field. The FE2 approach, on the other hand, utilises a fine mesh discretization linked to the Gauss points of a coarser mesh. The multiscale FE method (MsFEM) uses a fine mesh substructure to replace each element in the global model [START_REF] Weinan | Principles of multiscale modeling[END_REF][START_REF] Paladim | Guaranteed error bounds in homogenisation: an optimum stochastic approach to preserve the numerical separation of scales[END_REF]. Examples from the recent literature include Liang et al.s use of voxel models (generated using the well-known TexGen software) to analyse woven textile composites with each element being assigned to a particular material component [START_REF] Liang | Multi-scale modeling of mechanical behavior of cured woven textile composites accounting for the influence of yarn angle variation[END_REF]. Shi et al. use a three-scale model involving representative volume elements at the micro and meso-scales to model the fracture of braided composites [START_REF] Shi | Multi-scale ageing mechanisms of 3D four directional and five directional braided composites impact fracture behaviors under thermo-oxidative environment[END_REF]. Liu et al. use a variant of a voxel FE mesh termed the inhomogeneous FE mesh to model a woven textile composite, with the material varying from integration point to integration point rather than with material boundaries being assumed to follow the mesh [START_REF] Liu | Multi-scale progressive failure simulation of 3D woven composites under uniaxial tension[END_REF].

The approximation of errors in elliptical partial differential equations (PDEs), such as those which govern linear elastic stress analysis, has received a great deal of attention [START_REF] Eriksson | Adaptive finite element methods for parabolic problems i: A linear model problem[END_REF][START_REF] Verfrth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF][START_REF] Bank | Some a posteriori error estimators for elliptic partial differential equations[END_REF][START_REF] Ladevéze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF]. In many cases, error measures are utilised to drive mesh refinement through, for example, polytree decomposition algorithms [START_REF] Verfrth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF][START_REF] Larsson | On two-scale adaptive fe analysis of microheterogeneous media with seamless scale-bridging[END_REF][START_REF] Zhang | Adaptive analysis using scaled boundary finite element method in 3d[END_REF][START_REF] Song | A novel error indicator and an adaptive refinement technique using the scaled boundary finite element method[END_REF][START_REF] Allendes | Stabilized finite element approximations for a generalized boussinesq problem: A posteriori error analysis[END_REF]. Error estimates have been developed for a wide range of solution approximation methods. In addition to "conventional" FE, error estimates have been derived for boundary element methods [START_REF] Zhang | Adaptive analysis using scaled boundary finite element method in 3d[END_REF], immersed surface methods [START_REF] Heltai | Error estimates in weighted sobolev norms for finite element immersed interface methods[END_REF], multigrid and composite FEA methods [START_REF] Larsson | On two-scale adaptive fe analysis of microheterogeneous media with seamless scale-bridging[END_REF][START_REF] Rech | Two-scale composite finite element method for dirichlet problems on complicated domains[END_REF][START_REF] Pramanick | Error estimates for two-scale composite finite element approximations of parabolic equations with measure data in time for convex and nonconvex polygonal domains[END_REF], extended FE (XFEM) [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] González-Estrada | Locally equilibrated stress recovery for goal oriented error estimation in the extended finite element method[END_REF], and stress singularity problems [START_REF] Cai | Error estimate of a finite element method using stress intensity factor[END_REF]. Oden and Prudhomme's [START_REF] Oden | Estimation of modeling error in computational mechanics[END_REF] and Larson and Runesson's [START_REF] Larsson | On two-scale adaptive fe analysis of microheterogeneous media with seamless scale-bridging[END_REF] contributions are particularly relevant to the present work due to its treatment of error estimation in multiscale problems. The interested reader may also refer to our previous work on submodelling based on Oden's work [START_REF] Rahimi | Cad model simplification error estimation for electrostatics problems[END_REF][START_REF] Paladim | Guaranteed error bounds in homogenisation: an optimum stochastic approach to preserve the numerical separation of scales[END_REF] and to [START_REF] Chamoin | Goal-oriented error estimation and adaptivity in msfem computations[END_REF] for a contribution related to the so-called multiscale FE method (MsFEM). Numerous reviews of a priori and a posteriori error estimators are available in the literature [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] Grtsch | A posteriori error estimation techniques in practical finite element analysis[END_REF]. In most cases, error estimators are categorised as energy norm based (including element residual, subdomain residual methods and the constitutive relation error) or recovery based estimators. The latter includes the well know Babuska and Rheinboldt estimator, the Kelly, Gage, Zienkiewicz and Babuska estimator, and the Zienkiewicz-Zhu patch recovery technique [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF]. In the present work GOEEs are utilised as, rather than providing a general indication of approximation error, they allow for the quantification of error in a specific quantity of interest (QoI) [START_REF] Oden | Estimation of modeling error in computational mechanics[END_REF][START_REF] Grtsch | A posteriori error estimation techniques in practical finite element analysis[END_REF]. GOEEs were proposed in the 1990's [START_REF] Lin | Enriched goal-oriented error estimation applied to fracture mechanics problems solved by xfem[END_REF] in the work of Prudhomme, Oden, and Ainsworth, [START_REF] Oden | Estimation of modeling error in computational mechanics[END_REF][START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] Ainsworth | A unified approach to a posteriori error estimation using element residual methods[END_REF], Ladevéze [START_REF] Ladevze | Local error estimators for finite element linear analysis[END_REF][START_REF] Chamoin | A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems[END_REF], Bathe [START_REF] Grtsch | A posteriori error estimation techniques in practical finite element analysis[END_REF], and Cirak [START_REF] Cirak | A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem[END_REF] to name a few of the most influential authors in this prolific area of research. In many cases, the GOEE approach has been used to represent the uncertainty in a physically relevant quantity that may be turned into an error on a failure criterion [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF][START_REF] Andrés | Mesh adaptivity driven by goal-oriented locally equilibrated superconvergent patch recovery[END_REF][START_REF] Cirak | A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem[END_REF]. In the multiscale setting, a posteriori error estimates have been applied to Type A problems in the work of Tirvaudey et al. [START_REF] Tirvaudey | A posteriori error estimation and adaptivity in non-intrusive couplings between concurrent models[END_REF], wherein a weighted residual based GOEE is applied in non-intrusive sub-modelling problems. A posteriori error estimates for MsFEM problems have been developed in the work of Chung and Chamoin [START_REF] Chung | Goal-oriented adaptivity for gmsfem[END_REF][START_REF] Chamoin | Goal-oriented error estimation and adaptivity in msfem computations[END_REF]. In the context of the present work, QoIs are driving DoFs extracted from macro models. If meaningful error estimates for those quantities can be derived, it is the straightforward to sample from the resulting distributions and propagate uncertainties to continuum element sub-models. The statistical recovery process alluded to here is the focus and novelty of the present work. The methods introduced here allows an analyst to compute corrected microscopic solutions by propagating the effect of driving DoFs corrections -the error estimatesfrom macroscopic to microscopic scales.

The present work concerns itself with an industry relevant problem. We seek to evaluate the influence of macroscopic mesh quality onto the solution of microscopic problem. The solution developed here is non-intrusive in terms of remeshing at the macroscale, as this is typically not feasible in real problems. Furthermore, the solution should allow for redefinition of microscale models and should minimise computational expense. The resolution proposed here uses Zhu-Zienkiewicz (ZZ) error estimates to calculate errors in macroscopic deformation fields at a limited number of locations. Gaussian processes are then trained using a stochastic differential equation interpretation of the Matérn prior [START_REF] Lassi Roininen | Whittle-Matérn priors for Bayesian statistical inversion with applications in electrical impedance tomography[END_REF] to interpolate between different observations of the error, i.e. the ZZ estimates are seen as partial measurements of the full error field in a state estimation framework [START_REF] Mohinder | Kalman Filtering: Theory and Practice with MATLAB[END_REF][START_REF] Chilés | Geostatistics: Modeling Spatial Uncertainty[END_REF]. Computational bottle necks in the generation of the Matérn prior are circumvented though a GOEE model order reduction method that makes use of hierarchical clustering (enabled through k-medoids type approaches). Novelty in the work is derived from the application of ZZ error estimates in the multiscale setting, the use of ZZ evaluated adjoint problems in the GOEE setting for assessing the quality of local deformations, the use of functional Gaussian processes in multiscale problems, and the model order reduction approaches to make the process tractable in the industrial setting.

1 Problem setting and Formulation

Linear elasticity

We begin with a general setting for linear problems we will consider in the present work. Let us consider the linear elasticity problem over domain Ω ∈ R d , with d ∈ {2, 3} (only two dimensional problems will be considered in the present work, however extensions to three dimensions are trivial and the authors' previous work considers "shell" element formulations [START_REF] Bonney | Goal oriented error estimation in multi-scale shell element finite element problems[END_REF]). We look for a displacement field u ∈ U = H 1 (Ω), the space of fields of R d with values in R d whose derivative are square-integrable in Ω. The solution is assumed to satisfy homogeneous Dirichlet boundary conditions u = 0 on ∂Ω u , which is of non-zero measure. The complementary part of boundary ∂Ω, over which Neumann boundary conditions are applied, is denoted ∂Ω t . The weak form of the problem of elasticity is as follows. We look for u ∈ U such that

∀v ∈ U, a(u, v) = l(v) (1)
In the previous equation, bilinear form a is defined by

a(u, v) := Ω ∇ s u : C : ∇ s v dΩ (2)
where ∇ s is the symmetrised gradient operator, C is the fourth-order Hooke tensor corresponding to isotropic elasticity, and linear form l is defined by

l(v) := Ω f • v dΩ + ∂Ωt T • v dΓ (3)
where f is a known source field, and T is a known field of prescribed boundary tractions.

Finite element approximation

We now introduce the usual "P1" Lagrange FE method to approximately solve the elasticity problem introduced above. Domain Ω is decomposed into a set T h of non-overlapping simplexes. We assume for simplicity that this decomposition is error-free. We now look for FE field

u h ∈ U h ⊂ U such that ∀v h ∈ U h , a(u h , v h ) = l(v h ) (4)
U h is the space of continuous, piecewise linear FE displacement fields that satisfy the homogeneous Dirichlet conditions. Mathematically

U h =    u h (x) ∈ U | ∃U ∈ R N , u h (x) = N i=1 d j=1 (φ i (x)e j ) U (i-1)d+j = N i=1 ψ i (x)U i    (5)
where the N = N/d shape functions φ i defined over Ω and with values in R are the standard piecewise linear "hat" functions associated with the vertices of tessellation T h , and e j is the j th canonical vector of R d . In the previous expression, only the hat functions that vanish over Dirichlet boundary ∂Ω d are considered so that the homogeneous boundary conditions are automatically enforced.

The FE model may be written in matrix form as follows

KU = F (6)
where K ij = a(ψ i , ψ j ) is an element of the stiffness matrix K and F i = l(ψ i ) is an element of the force vector F.

Finite element error. The FE error e

h = u -u h ∈ U satisfies the weak form ∀v ∈ U, a(e h , v) = r(v) := l(v) -a(u h , v) (7)
and verifies Galerkin orthogonality condition

∀v h ∈ U h , a(e h , v h ) = 0 (8)

Hierarchical local submodelling

We now consider a sub-modelling approach which draws on the global model solution to provide more comprehensive results in some RoI. For the present work, sub models (or micro models) are considered to possess sufficiently fine mesh densities such that, compared to the global model, they are error free. In practice, the GOEE approach developed here can easily be applied to micro models as well as macro models, however this is omitted in the present work as to limit tangents that do not add to reader comprehension. It should also be noted that, in practice, an analyst has far greater control over micro models that macro models. In the aerospace example, many engineers across multiple departments may be called upon to produce a global aircraft model. Remeshing in such a case would be unacceptable and the degree of discretisation may be limited by computational resources. Local feature models may well be generated by single analyst however, thereby providing capacity for refinemnt. An RoI Ω I ⊂ Ω is to be re-analyzed using a finer local model (e.g. nonlinear, time-dependent model, 3D extrusion of a global plate/shell model). To achieve this, the displacements computed using the coarse FE model are extracted from ∂Ω I , and are subsequently applied to the boundary of the submodel (see e.g. [START_REF] Oden | Estimation of modeling error in computational mechanics[END_REF][START_REF] Paladim | Guaranteed error bounds in homogenisation: An optimum stochastic approach to preserve the numerical separation of scales[END_REF]). In our examples, the model kinematics and FE meshes are compatible at the interface ∂Ω I . This is not a restriction of the approach but one of our current implementation of the proposed methodology (see Figure 2 for further details on this aspect). A continuous error field will be developed in the following sections, meaning Ω I can be located without restriction within Ω. Our goal is to estimate the error in the boundary field to be transferred from the coarse to the fine model, and propagate this error though the chosen sub-modelling technique.

Goal oriented error estimation 2.1 Gradient-recovery-based error estimation

We wish to estimate e h a = a(e h , e h ) (error in energy norm) without computing too expensive an estimate for exact error field e h . The ZZ recovery-based estimate of e a is the following [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineerng analysis[END_REF]

ν 2 zz = Ω ( -∇ s u h ) : C : ( -∇ s u h ) dΩ (9)
where is a continuous, smoothed recovered strain field that is obtained by post-processing the discontinuous, piecewise constant tensor field ∇ s u h . One may simply average the values of the strain tensor corresponding to all the elements connected to a node, for each node of the mesh, and use the FE shape functions to obtain a continuous strain field. More advanced techniques include the super-convergent patch recovery approach (SPR) [START_REF] Zienkiewicz | The superconvergent patch recovery (spr) and adaptive finite element refinement[END_REF], or equilibrium-informed recovery procedures [START_REF] González-Estrada | Locally equilibrated stress recovery for goal oriented error estimation in the extended finite element method[END_REF][START_REF] González-Estrada | On the role of enrichment and statical admissibility of recovered fields in a posteriori error estimation for enriched finite element methods[END_REF][START_REF] González-Estrada | Mesh adaptivity driven by goal-oriented locally equilibrated superconvergent patch recovery[END_REF].

Error estimates for engineering quantities of interest

The energy is not necessarily the only quantity whose accuracy is of interest to the FE practitioner. In the present context of multilevel modelling, we are interested in a limit-states quantities in hierarchically-defined solution enrichment (i.e. structural zooms), such as maximum stress measures. To address the limitation of a-norm error estimation, the classical GOEE methodology starts by formulating the adjoint (or dual) problem.

We begin by reformulating the original FE problem. We look to find z ∈ U such that

∀v ∈ U, a(v, z) = Q(v) (10)
The linear functional Q : U → R extracts one of the QoIs from an arbitrary solution field v ∈ U.

Many different forms of Q(u) can of course be chosen and the most appropriate choice will depend on which QoIs relate to the analyst's objectives. For instance,

Q(u) = 1 |Ω I | Ω I e 1 • u dΩ is the average of u in the first spatial dimension, in region Ω I ⊂ Ω. Q(u) = 1 |Ω I | Ω I Σ : C : ∇ s u
dΩ extracts a weighted average of the component of the stress field in Ω I . For the present work, we look to define errors in driving DoFs in a multiscale problem setting, therefore QoIs will be weighted averaged extractions of components of u. Now, we have that

∀v ∈ V, a(e, z) = Q(e) = Q(u) -Q(u h ) = r(z) (11)
The adjoint field z (the influence function in the work of Prudhomme and Oden [START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF]) is unknown and must be approximated by using the FE method. We will look for FE field

z h ∈ U h such that ∀v ∈ U h , a(v, z h ) = Q(v) (12)
We can therefore write

Q(e) = a(e, z h ) + a(e, z -z h ) = r(z h ) + a(e, e z ) (13)
The residual term vanishes due to Galerkin orthogonality. The second term contains both exact errors (in the direct and adjoint problems), and may be evaluated by a ZZ-type recovery-based error estimate, as follows

Q(u) -Q(u h ) ≈ Qh -Q(u h ) := Ω ( -∇ s u h ) : C : ( z -∇ s z h ) dΩ ( 14
)
where Qh is an estimate of the QoIs, recovered from the FE field, that approximates the QoI Q(u). If error analyses were limited to a single RoI with a small number of driving DoFs, the above treatment would be sufficient for error propagation. Pointwise QoIs, evaluated using a Dirac delta function for the development of Q(u), could be established along the boudnary of the RoI in the macro model and GOEE terms calcualted. In the present work, it is assumed that either no such restriction on the placement of ajoints in Ω exists, or that the number of ajoints required to sufficiently constrain all driving DoFs at all possible RoI boundaries is so high that implementing all of them would be prohibitively expensive. As a resolution to this challenge, we will now present a reconstruction of the macroscopic error field by training a Gaussian field using a sparse number of GOEE observations. The resulting probabilistic representation will yield a statistical microscopic error distribution wherever needed.

3 Full error field reconstruction through Bayesian state estimation 3.1 "Truth" reconstruction of finite element error fields using the goal oriented error estimation approaches

We now introduce a corrected model

ûh (x) = N i=1 ψ i (x) Ûi (15)
where

Û = U + ê (16)
In the previous expression, ê is a vector of FE DoF values corresponding to a FE representation of the exact error field e h = u -u h . This is to be understood in the sense that the corrected FE displacement should yield QoIs that are closer to that delivered by an infinitely refined FE model. That means that we aim to estimate the pollution error, whilst the interpolation error itself is assumed to be small enough to be ignored, or is to be estimated by other means (this will be discussed later on for the specific case of hierarchical local sub-modelling). Two relevant FE representations come to mind for correction field ûh (x)-

u h (x) = N i=1 ψ i (x)ê i • the FE interpolant I h e h of the exact error field e h in U h ,
• the L 2 (Ω)-projection of the exact error field e h onto U h .

It seems appropriate at this stage to require ê to be the vector of nodal values of the L 2 (Ω)projection of the exact error field onto the FE space. Indeed, any FE functional of the exact error, i.e Q(e) = Ω ξ h • e dΩ, where ξ h is a vector-valued FE field of U h , has the following property

Ω ξ h • e h dΩ = Ω ξ h • π h L2 e h dΩ (17)
where π h L2 is the L 2 (Ω)-projector and the above equality stems from standard operations on projectors. Hence, the above observation of the exact error is also an observation of its L 2 projection. As a consequence, it is possible to construct a ZZ-based approximation of the L 2 projection of the exact error onto the FE space by noting that:

∀i ∈ 1 N , Ω ψ i • e h dΩ = c i T Mê (18)
where:

M ij = Ω ψ i • ψ j dΩ (19)
and c i is the i th canonical vector of R d . Hence, knowing these N GOEE ZZ estimates, an inversion of the mass matrix M gives access to ê. A similar conclusion can be made if only part of the error field is to be reconstructed from GOEE. Of course, computing N GOEE ZZ estimates is out-of-the question in all but simplest of demonstrative models. We will instead proceed by Gaussian Process functional regression [START_REF] Nguyen | Gaussian functional regression for state prediction using linear PDE models and observations[END_REF][START_REF] Cuong Nguyen | Functional Regression for State Prediction Using Linear PDE Models and Observations[END_REF], to statistically estimate these N quantities based on the availability a selected few of them.

Bayesian reconstruction of finite element errors fields from sparse functional observations

Linear state estimation

Bayesian state estimation setting. We assume that we have calculated an array of M ZZ estimates of errors in scalar QoIs. These estimates are seen as noisy observations of the corresponding exact error in QoI, which reads as

∀i ∈ 1, M , d i = Q i T ê + i = ∆Q i + i ( 20
)
where

∆Q i ≈ (Q i (π h L2 u) -Q i (u h ))
is the true error in the i th QoI. It is vital to note that this is unavailable at this stage and will be given a probability distribution later on in this section, following Bayesian modelling of lack of knowledge. Above, i ∼ N (0, σ i ) is a white noise that allows us to represent the fact that the error estimate is not exact (we have, after all, used the ZZ estimate to this end). We will set σ i = |∆Q i | 2 , which qualitatively encodes the fact that "the typical error in the error estimate is of the order of the error estimate itself"1 .

We define the vector of errors in QoIs as

d = ∆Q + =     ∆Q 1 ∆Q 2 ... ∆Q M     + =      Q 1 T Q 2 T ... Q M T      ê + = Hê + ( 21
)
where ∼ N (0, Σ ) is a M -dimensional random vector containing the errors in ZZ-estimated errors in QoI, ∆Q is a M -dimensional vector containing the true errors in QoIs, and H is a M by N linear observation operator. In our examples, we will use FE shape functions as extractor of quantities of interest. In this case, Q i T = c I(i) T M, as shown in the previous section.

We further assume that the unknown error vector ê is zero-mean and Gaussian distributed, as follows ê ∼ N (0, Σ)

which, together with the distribution of , encodes all prior knowledge about epistemic sources of uncertainty in our state estimation problem.

Posterior distributions. We can now look for the posterior probability distribution of ê, given its M noisy and partial measurements. Following standard techniques in Gaussian processes inference, the distribution of the state and the noisy observations is jointly and consistently Gaussian, summarised by

d ê ∼ N 0 0 , HΣH T + Σ HΣ ΣH T Σ (23)
The posterior probability of ê is simply obtain by Gaussian conditioning 2 , i.e.

ê|d ∼ N (e , Σ ) (25)
where e = ΣH T (HΣH

T + Σ ) -1 d Σ = Σ -ΣH T (HΣH T + Σ ) -1 HΣ (26)
The marginal posterior distribution of any set of linear combinations of the state vector is given by

Eê|d ∼ N (Ee , EΣ E T ) (27)
where prediction operator E is of size n e by N . In this paper in particular, we are interested in restriction operators, i.e. Boolean rectangular matrices with a single non-zero element per line.

In practice, this amounts to predicting the components of the error vector corresponding to a particular subset of the FE DoFs, for instance all the DoFs corresponding to a subdomain Ω I ⊂ Ω, or to the boundary of this subdomain in the context of the Dirichlet hierarchical sub-modelling introduced previously.

Prior covariance structure

Great care is required when choosing the prior covariance. Computational tractability is of course vitally important, however it is key to note that the form of prior covariance will of course dictate how observations are propagated into the posterior distribution. If unrepresentative correlations are utilised across the domain it is entirely likely meaningless posterior distributions will be realised. In the context of structural modelling, a simple distance based prior may well suggest strong correlation between two neighbourhoods that, while physically close, are not mechanically connected and thus should not greatly influence one another. A simple example of this challenge is a "U" channel section, wherein the top edges of the section may be relatively close to one another but this proximity does not (necessarily) suggest strong correlation. What follows is a brief discussion of Matérn prior structures. Standard (kernel based) and random process based formulations are introduced, with the latter being utilised in the remainder of the work.

Kernel-based covariance. A standard choice is to use Kernels to define the covariance matrix, i.e.

Cov(c

i T ê, c j T ê) = k P i d , P j d δ i d - i d , j d - j d ( 28
)
2

The bayesian conditioning may also be expressed as a Kalman update

e = 0 -G(∆Q -H 0) Σ = (I -GH)Σ , (24)
where the Kalman gain is defined as

G = ΣH T (HΣH T + Σ ) -1 .
where c i is the i th canonical vector of R N and P k ∈ R d is the position of the FE vertex associated with the k th linear shape function. For instance, the Whittle-Matérn Kernel reads as

k (P i , P j ) = α 2 P i -P j 2 l ν K ν P i -P j 2 l ( 29
)
where K ν is the modified Bessel function of second kind of order ν, and α and l are hyperparameters of the Gaussian process3 .

Random-process-based covariance In this paper, we define the prior for uncertain state ê using a stochastic partial differential equation approach [START_REF] Lassi Roininen | Whittle-Matérn priors for Bayesian statistical inversion with applications in electrical impedance tomography[END_REF][START_REF] Sidén | Spatial 3D Mat\'ern priors for fast whole-brain fMRI analysis[END_REF]. The distribution of ê is implicitly defined through the solution of stochastic linear system

M + β 2 K ê = α √ MW (30)
which we rewrite as

Aê = b with A = M + β 2 K and b = α √ MW (31)
In the previous expressions, taken from the work of Roininen, Huttunen, and Lasanen [START_REF] Lassi Roininen | Whittle-Matérn priors for Bayesian statistical inversion with applications in electrical impedance tomography[END_REF], M and K are defined by M ij = Ω ψ i (x)ψ j (x) dΩ and K i,j = Ω ∇ψ j (x) • ∇ψ i (x)dΩ, respectively, and are (proportional to) the standard FE mass and stiffness matrices, respectively. W ∈ R N is a white noise vector with normally distributed independent components.

By applying the standard rules of linear operators applied to Gaussian random variables, we find the prior covariance for our error vector ê is

Σ = α 2 M + β 2 K -1 M M + β 2 K -1 (32)
Notice that mass matrix M is invertible, under weak assumptions on the quality of the FE discretisation. Moreover, both M and K are symmetric and positive. Hyperparameter β is a covariance length, which controls the smoothness of the random field, while hyperparameter α controls the overall amplitude of the process4 5 .

Parameter optimisation

Optimising the hyperparameters of the prior in a data-driven manner may be done by maximising the data-likelihood with respect to hyperparameter vector θ := (α β) T . In our case, the data log likelihood reads as [START_REF] Christopher | Gaussian Processes for Machine Learning[END_REF] log

(L θ (d) = - 1 2 d T (H T Σ(θ)H + Σ ) -1 d -log(Z(θ)) (36)
where

Z(θ) = (2π) M |H T Σ(θ)H + Σ |.
Providing that the number of observations remains small, this optimisation problem is tractable. The expression of the sensitivities of the datalikelihood can be found in [START_REF] Christopher | Gaussian Processes for Machine Learning[END_REF].

For now, this optimisation will be performed "by hand". We calibrate the overall size and length scale of the Gaussian process in order to match observations, by a procedure of trial-and-error. Fully automatising the algorithm is a perspective of this work but it not implemented here.

∀v ∈ U , (e, v) L 2 (Ω) + β 2 a (e, v) = α g, v (33)
where random linear functional g is such that

∀(v, w) ∈ U 2 , E ( g, v g, w ) = (v, w) L 2 (Ω) (34)
and e ∈ U × Ξ is the random continuous field of interest, and is defined over the set of all outcomes Ξ. 5 In paper [START_REF] Nguyen | Gaussian functional regression for state prediction using linear PDE models and observations[END_REF][START_REF] Cuong Nguyen | Functional Regression for State Prediction Using Linear PDE Models and Observations[END_REF], the authors propose, in another context, to make use of the following prior covariance:

Σ = α 2 K -1 M + β 2 K K -1 (35)
In our preliminary experiments with this prior, we failed to eliminate the spurious mesh dependencies that seem to appear when regularisation parameter β is non-vanishing. As a consequence, we did not manage to generate sufficient freedom to optimise this prior in a data-driven manner.

Computational bottleneck

It is now possible to compute the marginal distribution of the set of unseen linear combinations of components of ê that we are interested in. The posterior mean of Eê, is EΣH T (HΣH T + Σ ) -1 d. Product ΣH T appearing in innovation matrix HΣH T + Σ dominates the computational cost. It may be evaluated by considering the lines of observation operator H, one-by one, and leftapplying the covariance operator to each individual observation vector. This involves solving M linear problems of the type M + β 2 K . = .. This operation scales as o(M N 2 ). The inversion of (HΣH T + Σ ) -1 is of negligible cost as we wish to recover full error fields based on a limited number of observations, with typically M ≤ 20.

The computation of the posterior covariance is only tractable for low-rank prediction operators E. Indeed, the posterior covariance for the unseen quantities to be predicted is EΣE T -EΣH T (HΣH T +Σ ) -1 HΣE T . The second term can be computed by taking advantage of the fact that the number of observations is small. However, if prediction operator E has a large number of lines n e , the first term becomes intractable. In our case of hierarchical Dirichlet submodelling, E will extract all the DoFs of ∂Ω I .

We may generate realisations of Ee|d directly, without actually constructing the posterior covariance matrix. Such a sample s ∼ N (m , Σ ) is obtained as follows:

s = E x 0 + G d -Hx 0 ( 37
)
where x 0 ∼ N (0, Σ) is drawn from the prior distribution, d ∼ N (d, Σ ) is a randomly perturbed observation, and the Kalman gain is G = ΣH T (HΣH T + Σ ) -1 . Drawing a sample from the posterior distribution is dominated by the cost of solving two systems involving the matrix M + β 2 K on the left hand side (once to draw a sample from the prior distribution, and once to apply the Kalman gain). In this analysis, we consider that the innovation is already assembled and factorised to compute the posterior mean. This analysis of the numerical tractability of the Gaussian regression process is unfortunate news. Indeed, the goal of the previous developments was to replace the expensive computation of many observations of the error field via the adjoint method by a cheap interpolation between a few selected observations. However, we see here that the Gaussian process interpolation itself requires solving a series of FE systems, each as expensive to solve as an adjoint problem (M systems to build the innovation matrix, one more to compute the posterior mean and 2 per sample to evaluate the variability via Monte-Carlo). The method can only be of benefit if the cost of evaluating the posterior distribution is made orders of magnitude cheaper than that of actually computing the all QoIs directly, which we will do in the following section by proposing a dedicated low-rank approximation strategy.

Let us emphasise that the numerical difficulty is due to the dimension of the parameter space (as is often the case in state estimation, when using Kriging or Kalman filters to identify highdimensional spatial fields from sparse measurements), not the number of observations (as may appear when using Gaussian processes in machine learning, for instance to perform regression in relatively small feature spaces but with high volumes of data).

Results -Example Problem definition and application of full rank Bayesian State Estimation

Context and example problem definition

We beign by illustrating the Gaussian process enhanced GOEE methods developed above through a simple multiscale problem. We here consider a 2d (plane stress) continuum element "perforated L" plate, fixed along the lower edge and loaded in the vertical by a uniformly distributed load (see figure 1 (a)). Note that the "x" dof at the node at the lower left hand corner is constrained to prevent rigid body motion. To clarify the loading arrangement, the deformed (unity scale factor) and undeformed mesh is presented in figure 1 (b) (note that this is the primal problem solution, found by normal FE solution approaches). We consider the situation presented in figure 1 as our macro problem and introduce several potential micro models at one of two locations (model interfaces are shown in blue and red in figure 2) to complete the multiscale setting. The introduction of only a small number macro/micro boundaries clearly simplifies the problem, to the benefit of the error field recovery. In reality, it is highly likely that many potential interface will be possible over the global model, thereby greatly increasing the number of potential QoIs. In the present work, we limit selection of QoIs to the two boundaries highlighted in figure 2. In all cases, QoIs are centred on one of the boundary nodes using a Dirac delta function (furnishing us with point wise error estimates) and consider the resultant of the DoFs in a random direction. In this way, each QoI observation provides information on errors in both x and y. For the meshes presented in figure 2 and given the restrictions listed above, the "all duals" case (i.e. the situation in which all relevant adjoints are solved) amounts to solving an additional 51 problems. Coherency between meshes at the macro/micro model boundaries is maintained, thereby simplifying the propagation of errors through the scales (note however that this is done purely to simplify the implementation and is not a requirement). In the present work, error propagation is only performed using the "two cirle" micro model in the blue boundary position. Error fields are recovered over the entire macro model domain however, in the interest of space, repeated micro model implementations that do not add to reader comprehension are ommitted. In all cases, linearity in both material stiffness and geometry is assumed.

U Y =0 U X =0 X Y P (

Full Rank Bayesian Recovery of Goal Oriented Error Estimation field

We present here the recovery of error fields over the entire macro domain using the methods outlined above, wherein the stochastic partial differential equation representation of the Matérn prior is used (see equation 32), and consider how the number of observations (duals) influences the posterior distributions. In the following, we randomly select a number of duals from the 51 available (restrictions are outlined and justified in the previous section), optimise the hyperparameter α (see equation 32) by maximising the likelihood given by equation 36 (note β is taken as 1 in all cases, being a length scale on the order of the problem itself), evaluate the prior and posterior, and sample form the resultant to propagate error though the model scales. Given that β is assumed, the hyperparaemter optimisation becomes 1d. A starting value of 3.5 × 10 -4 is assumed for α in this optimisation (note that this value is preserved in realisations that make use of the prior but no observations) and the process amounts to a simple polynomial fitting approach (applied to evaluations of equation 36 for a wide range of α values around the initial estimate). The optimisation approach implemented here is clearly simplistic and far more sophisticated methods may be utilised for enhanced hyperparameter estimation, however given the scope of the present work the method is judged to be appropriate. It is important to note that some form of hyperparameter optimisation is required here, as the use of a single set of hyperparameters (regardless of number of observations) infers some form of prior knowledge that fixes the amplitude of the prior. We may sample from multivariate normal distributions (∼ N (µ, Σ)) using standard methods, such as the following. A covariance matrix Σ may be decomposed to the matrix A using Cholesky decomposition, such that

Σ = AA T (38)
A noise vector, z, may then be constructed whose components are independent and drawn from a standard normal. A sample, ζ is generated using the mean µ by

ζ = µ + Az ( 39
)
Note that, in the context of the present work, each sampling ζ is a realisation of the error field ê. Figure 3 presents 4 such samplings, all developed using the same z noise vector (shown in figure 3 (a)). By using the same z in all samplings, the influence of number of training duals (observations) can be visualised. Figure 3 (b) shows sampling performed on the prior, figure 3 (c) shows sampling performed on a 2 dual trained posterior, figure 3 (d) shows sampling performed on a 16 dual trained posterior, and figure 3 (e) shows sampling performed on a posterior trained by all (51) duals. Note that dual training sets are subsets of one another (that is to say, the 2 duals selected to train that particular posterior are also found in the 16 used to train that particular posterior), meaning features in the posterior resulting from particular observations should be preserved as the number of training cases is increased. Figures 3 (b)-(e) clearly show that, as the number of observations is increased, the posterior tends (as expected) to the case where all duals are considered (figure 3 (e)). We can also see the influence of the Dirichlet boundary conditions in all samples (clearly, we know these when generating the prior, hence variance in this neighbourhood is greatly reduced), and note that, as we move away from the RoI boundaries (and hence observations) we tend towards the prior distribution.

The influence of uncertainty at the RoI boundary (and the effect of number of QoI observations on this) can be further visualised by considering stress distributions over the entire the RoI. We consider the cases where 2, 16, and all available duals are used to train the posteriors. Monte-Carlo samplings are taken from each posterior and are used to calculate potential stress states in the micro model (as above). Distributions are then fitted to the results on an integration point by integration basis, allowing for the determination of the mean von Mises stress (M ean V M ) and its standard deviation (Std V M ) over the micro model. Results for the 2, 16, and all available duals training cases are shown in figures 4, 5, and 6, respectively. We may observe that, in these results, the number of training duals has little influence on the mean of the stress field. This is of course to be expected, as this field will be mainly governed by driving macro solution the macro solution (which is common to all realisations). Similar features are observed in all standard deviation fields (the largest deviations are located around stress concentrations), however the magnitude of the standard deviations decreases with an increasing number of training duals (the standard deviation of the 2 dual training cases is approximately 3 times greater than that of the all duals training case). This of course can be rationalised as a transition from the uncertainty associated with the prior (which is relatively large and is recovered in the low number of training duals case) to the uncertainty associated with the noise in our error observations (which is propagated in the all training dual case). Observations regarding the development of the mean and standard deviation of the von Mises stress fields can be further demonstrated by considering only the peak von Mises stresses in the micro model and noting how the resulting distributions vary with number of training duals. This is presented in figure 7, wherein we again see a preservation of the mean and a reduction in the standard deviation of the peak von Mises stress as we increase the number of training duals (from 2 to all available). Note that the standard deviation of the 32 training duals and all training dual distributions are near identical, suggesting that there is limited utility in calculating the additional QoIs.

Our discussions thus far have focused on a fixed set of training duals. That is to say, whenever we make a realisation for a 16 training dual case (for example), we have used the same 16 duals. While these we selected at random from the 51 available over the RoI boundary, one may question how representative the results are of all possible dual selections. Is it possible that the selected duals were particularly "good" (or "bad") by pure chance? This question is addressed in two ways here. First, we consider different random selections of training duals to develop the peak von Mises stress distributions, as shown in figure 8. While these distributions are distinct, they all possess the general properties noted above. We see, in all cases, a tendency from large amounts of uncertainty (associated with the prior) to relatively small amounts of uncertainty (associated with the noise in the error observations) as the number of training duals is increased. The Kullback-Leibler divergence (or relative entropy) can be used to formalise these observations. The Kullback-Leibler divergence between the approximate distribution p n1

Lẽ obtained using a prior covariance of rank n 1 and the approximate distribution p n2

Lẽ obtained using a prior covariance of rank n 2 , is

D KL (p n1 Lẽ ||p n2 Lẽ ) = 1 2 Tr Cov (p n2 Lẽ ) -1 Cov (p n1 Lẽ ) -(M + n e ) + ln |Cov (p n2 Lẽ ) | |Cov (p n1 Lẽ ) | (40)
In the context of the present work, we may use D KL to compare posteriors against the all dual training case. This is presented in figure 9. Note that, in figure 9, each evaluation of D KL is made using a posterior trained using a random selection of duals. As expected, as we increase the number of training duals, we see the spread in D KL reduce to the point where it represents only the noise in the observations. For small numbers of training duals we note an increasingly larger spread in D KL , suggesting a recovery of the prior distribution. Note that one may be "lucky" (or "unlunky") when randomly selecting duals and achieve an posterior that has a relatively low (or high) D KL . It is not possible, however, to ensure that this is the case without deliberate selection of the training duals based on some exterior selection pressure, which is outside the scope of the present work. 5 Low-rank Gaussian process by Galerkin Model Order Reduction

Stochastic reduced basis approximation

We now look to address the computational bottleneck identified in section 3.3 through the projection of the Gaussian process onto a low rank space. We will look for a reduced basis approximation of prior distribution of ê. The corresponding approximation reads as

ê ≈ ẽ = φγ (41)
where the n φ columns of φ are the deterministic reduced bases. Random vector γ may be optimally computed using Galerkin's principle (see e.g. [START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF][START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF][START_REF] Hoang | A fast, certified and "tuning free" two-field reduced basis method for the metamodelling of affinely-parametrised elasticity problems[END_REF])

∀ b ∈ R N , φ T (b -Aẽ) = 0 (42)
We find that ẽ = φ φ T Aφ -1

φ T b . (43)
Therefore, the approximate prior reads as

ẽ ∼ N 0, Σ := φΣ γ φ T (44)
where the covariance matrix for the n φ reduced random variables γ is

Σ γ = α 2 φ φ T (M + β 2 K)φ -1 φ T Mφ φ T (M + β 2 K)φ -1 φ T (45)
The approximate prior covariance Σ given in the expression above is rank-deficient. The systems of equations involved in computing the posterior mean and co-variance are now reduced to n φ equations, but the corresponding system matrices are fully populated.

Reduction of computational operations. The low-rank approximation can be taken advantage of when computing posterior distributions, typically by making use of the Sherman-Morrison formula. These computational tricks are classical and therefore not detailed in the paper (refer for instance to the usual implementation of the Ensemble Kalman Filter)

Projection subspaces

Karhunen-Loève transform.

If the basis vectors are made of all the eigenvectors of the full covariance matrix, i.e.

∀ i = 1 ... N, Kφ i = λ i Mφ i and ∀ j = i, φ T i φ j = δ ij , (46)
then the previous derivation is the Karhunen-Loève transform: the random coefficients are uncorrelated (and independent in the Gaussian case) and the covariance matrix Σ γ is diagonal. More precisely, we have

Σ = α 2 φΛ m (Λ m + β 2 Λ k ) -2 φ T (47)
where Λ m and Λ k are diagonal matrices with i th diagonal entry equal to φ T i Mφ i and φ T i Kφ i , respectively. The Karhunen-Loeve transform provides a natural method to truncate reduced basis expansions, but requires computing the spectrum of the discrete Laplacian operator. The associated computational cost is not acceptable in the present context.

Subspace generated by the set of precomputed dual solutions

What can we choose as good projection subspace? An obvious first choice is to reuse what is already at our disposal: the set of all adjoint states z i defined by

∀i ∈ 1, M , K T z i = Q i (48)
where Q T i is the i th line of observation operator H. These linear systems of equations have been solved to compute the GOEE corresponding to each component of data vector d. We will store concatenation Z = (z 1 z 2 ... z M ) ∈ R N ×M .

RBF enrichment and regularisation

Additional reduced basis functions are required to let the Gaussian process express uncertainty. Several options have been explored during the early phase of this work. In particular, we have successfully used coarse functions generated by Algebraic Multigrid Solvers [START_REF] Brezina | Adaptive smoothed aggregation (α SA) multigrid[END_REF] as projection vectors. In this paper, we report results related to the use of a set of analytically defined and numerically optimised readial basis functions (RBFs). This will be fully detailed in the following chapter. For now, we assume that these n b functions are represented in the FE space, and can be encoded by their nodal values, contained in matrix B ∈ R N ×n b .

Finally, the set of duals and RBFs are concatenanted as follows:

φ = (B Z) ∈ R N ×nΦ (49)
To avoid the potential (near-)linear dependencies of the columns of φ, an SVD of φ is computed, and φ is replaced by the matrix whose columns are the left-singular vectors of this decomposition that are associated non-vanishing singular values up to a numerical tolerance.

Accuracy measure

The discrepancy between the exact error vector and its low rank approximation as is a random variable defined by

∆e = ê -ẽ = A -1 I d -AB B T AB -1 B T b (50)
where we recall that b is a zero-mean Gaussian-distributed vector with covariance matrix α 2 M. Therefore, the prior distribution of discrepancies in errors in predicted QoIs Eê is given by

E∆e ∼ N 0, Σ ∆E := EΣ ∆ E T (51)
with

Σ ∆ = α 2 A -1 I d -AB B T AB -1 B T M I d -B B T AB -1 B T A T A -1 (52)
The computation of EΣ ∆ is not tractable, as we would need to compute A -1 E T , which is a matrix whose columns are made of all dual solutions, the computation of which we are precisely trying to avoid by using Gaussian process functional regression. However, choosing a selected few lines of E can provide valuable GOEE of the accuracy of the low-rank approximation to the full Gaussian process, which we will make use of later on.

6 Adaptive set of radial basis functions using stochastic adjoint-based goal oriented methodology and hierarchical clustering

Reduced basis vectors defined as finite element interpolants of RBFs

We propose to construct the reduced space as the span of a set of radial bases functions, projected in the FE space. To this end, let us define

∀ i ∈ 1, n b , ∀ x ∈ Ω, ζ i (x) = ω i/d exp - 1 2 x -c i/d 2 l i/d 2 e i-( i/d -1)d ( 53
)
ζ i (x) = η i/d (x) e i-( i/d -1)d (54)
and its FE interpolant

ζ i h (x) = I h ζ i (x) (55)
which may be represented by basis vector γ i ∈ R N , the component of which are such that

∀ i ∈ 1, n b , ∀ x ∈ Ω, ζ i h (x) = N j=1 ψ j (x)γ i j (56)
We assume in these developments that n b is a multiple of d. • For any i ∈ 1, n b /d , we define l i as the distance to the N th l closest point of set {c i } i∈ 1,n b /d , with N l ∈ R + * \{1}, in the sense of the standard euclidean distance. This technique ensures that the set of RBFs are overlapping to a satisfying degree.

Set {c i } i∈ 1,n b /d ∈ R d

• We define weight set {ω

i } i∈ 1,n b /d by ∀ i ∈ 1, n b /d , ω i (x) = 1 n b /d i=1 η i (x) (57)
• The set of centers C = {c i } i∈ 1,n b /d is constructed adaptively, using a goal oriented hierarchical clustering approach, as described in details in the next section.

To initialise the adaptation process, a standard coarse clustering of the vertices of FE tessellation T h is first performed. We use the agglomerative hierarchical clustering tree implemented in Matlab to partition the set of vertices of the mesh into n 0 b /d. Then, each element of set C 0 = {c 0 i } i∈ 1,n 0 b /d is computed as the centroid of the i th cluster of mesh vertices. Once the n b /d radial basis functions are constructed, we can build projection matrix B by concatenation:

B = (γ 1 γ 2 ... γ n b ) (58)

Stochastic dual weighted residual approach

We define the following error estimate

ν = M i=1 Var Q i T ∆e (59)
where {Q i } i∈ i M are an arbitrary set of extractors of scalar QoIs. As mentioned previously, we will use a random selection of the lines of prediction operator H.

We note that

Q i T ∆e = Q i T A -1 r = z i T r (60)
where i th dual z i is such that A T z i = Q i , and the random residual vector r := b -Aẽ has expression

r = I -Aφ φ T Aφ -1 φ T b =: ξb (61)
Therefore, r is a zero-mean Gaussian vector with covariance Cov(r) = ξMξ T . Hence, [START_REF] Zienkiewicz | The superconvergent patch recovery (spr) and adaptive finite element refinement[END_REF] reads as

ν = M i=1 z i T Cov(r)z i (62)
Equation ( 60) is usually coined Dual Weighted Residual in the error estimation literature [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF][START_REF] Oden | Estimation of Modeling Error in Computational Mechanics[END_REF], albeit usually employed in a deterministic setting. In order to make this expression small, an owing to Galerkin orthogonality, the suggestion is to use the adjoints {z i } i∈ i M as reduced basis vectors. They are, of course unavailable and we can only try and produce a similar effect by using selected RBFs.

Extending classical error estimation procedures, we use Galerkin orthogonality 6 to reveal that (62) reads as

ν = M i=1 z i -Π Φ z i T Cov(r) z i -Π Φ z i (64)
where

Π Φ = Φ Φ T Φ -1
Φ T is the norm two projector onto span(Φ), which we may further write as

ν = M i=1 w i T y i (65)
with w i := z i -Π Φ z i and y i := Cov(r)w i . Equation ( 65) can be seen as an extension of the classical dual weighted residual framework to discrete stochastic Galerkin approximations. As usual in dual weighed residual approaches, we may interpret the previous expression as a discrete error density, as

ν = M i=1 N j=1 w i j y i j = N j=1 M i=1 w i j y i j = N j=1 M i=1 |w i j y i j | . (66)
We may finally attempt to select the radial basis functions that approximately minimise ν, i.e. our chosen measure of the distance between exact and low-rank Gaussian priors for all QoIs. This can be done by cancelling large elements of set

M i=1 w i j y i j j∈ 1 N
, thereby making sure that for DoFs corresponding to large values of projected dual vectors {w i } i∈ 1 M , the associated components of "residual" vector {y i } i∈ 1 M are small (this is the central idea of the dual-weighted residual approach [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF][START_REF] Oden | Estimation of Modeling Error in Computational Mechanics[END_REF]). This implies refining the projection space in those areas, following Galerkin's principle.

Goal oriented hierarchical clustering

We propose to use the discrete GOEE density described previously to guide the hierarchical refinement of the agglomerative clustering algorithm described in section 6.1. This is done as follows. Given the current clustering of mesh vertices at iteration k of the construction process, associated with centroid set

C k = {c k i } i∈ 1,n k b /d , we compute, for every cluster index l ∈ 1 n l c q k l = n l c j=1 M i=1 w i,k I l (j) y i,k I l (j) (67)
where n l c is d times the number of vertices of cluster l, and I l is an index map between the clusterwise and global numbering of FE DoFs. In words, q k l is the contribution to the GOEE metric of cluster l at iteration k of the construction algorithm.

Then, we set cluster l for refinement, with

l = arg max l (q k l ) (68)
This cluster is subdivided into n r subclusters using K-means, the centroids of which will define new RBF centers, yielding the updated center set

C k+1 = {c k+1 i } i∈ 1,n k+1 b /d with n k+1 b = n k b + n r .

Rank selection by holdout/cross-validation

To make sure that the distribution of reduced-rank Gaussian process Eẽ is close enough to "truth" distribution of Eê, we will monitor the size of Σ ∆E through the greedy construction process. The RBF selection algorithm stops if assertion

1 M M i=1 Var Q i T ∆e ≤ ξ 1 M M i=1 Var Q i T ẽ (69)
with tolerance 0 < ξ 1, is satisfied. In the previous equation, the notation indicates that we use a new set of randomly selected lines of E to populate QoI extractors {Q i } i∈ i M . Indeed, the there exists a severe risk of over-fitting if using set {Q i } i∈ i M to build the reduced basis and monitor convergence. Unfortunately, this requires the computation of additional dual vectors. Alternatively, one could cross-validate the previous methodology, thereby avoiding this additional computational effort at the cost of a heavier implementation.

Results of Low Rank Approximation of Full Error Field Gaussian Process

We now implement the goal oriented hierarchical clustering approach described above to the example problem outlined in section 4.1. In the following, we consider a randomly selected 25 training duals in all cases and look to train a posterior based on a low rank approximation of the prior. RBFs define a projection space for the prior. We here use a hierarchical clustering approach to locate RBFs in the solution domain (RBFs are normalised to ensure partition of unity). A series to clusters are presented here to demonstrate the refining nature of the hierarchical approach, both at the macro and micro scales. Clusters for the example macro problem are shown in figure 10, with "Cluster 1" illustrating the initial cluster arrangement (due to k-medoids) and Clusters 3, 5, and 7 relating to the 3rd, 5th, and 7th refinement, respectively. Cluster refinement is achieved by evaluating the quantity q (as per equation 67, here referred to as the clustering density) for a given cluster arrangement (and associated RBF space), noting the cluster within which the weighted summation of q is maximised, and performing k-medoids on the selected cluster. In the results presented here, an initial 10 clusters (Cluster 1) are created to begin the hierarchical process, with subsequent refinements partitioning the selected cluster into 3. Example plots of the density q are presented in figure 11 for the clusters in figure 10. The example problem considered in the present work limits the location of QoIs to two RoI boundaries (see figure 2). As a result, the error in low rank QoI estimations (characterised by the density q) is concentrated in the upper right corner of the macro model. This can be clearly seen in figure 11, wherein there is not only the desired reduction in the magnitude of q with cluster refinement but also a concentration of q on the boundaries of clusters. For completeness example RBFs, resulting from the clusters in figure 10, are presented in figure 12.

We must now consider how, and if, the clustering based RBF refinement leads to convergence of the low rank approximation of the Gaussian process prior. We may first consider the entropy of the difference between our "exact" error field (furnished without model order reduction of the prior) and the RBF approximation of it. Recall that equation 52 gives the covariance of the discrepancy field. In lieu of a complete set of QoIs (E), we use the restricted set H which is made up of the same QoIs that inform our GOEE observations (see equation [START_REF] Liang | Multi-scale modeling of mechanical behavior of cured woven textile composites accounting for the influence of yarn angle variation[END_REF]. Through this substitution, we note that the calculation A -1 H T (used to evaluate the distribution in equation 51) becomes trivial, as it is equal to the solutions of the dual problems (z) we have already solved. The entropy of the difference between exact and low rank error fields is given by

1 2 ln det (2π exp (Σ ∆ H)) (70)
The evolution of this entropy term with cluster iteration can be seen figure 13. If the hierarchical clustering approach is effective we would expect a continued reduction in the discrepancy entropy, indicating a collapse of the discrepancy distribution as the low rank approximation approach the exact field. This feature is clearly demonstrated in figure 13. The near linear rate at which the entropy reduced with cluster iteration is somewhat surprising and suggests that, at least for the example problem, the improvements offered by cluster refinement are incremental despite the significant changes in q. Note that the entropy of the error field distribution itself is also presented in figure 13 for reference. We may also perform Monte-Carlo analyses on the micro models by sampling error correct DoF distributions that are generated by exact and low rank approximations of the Gaussian process terms (similar to the demonstrations performed in section 4.2). Peak von Mises stress may be extracted from the "2 circle" micro model (as used before) for the low rank approximations associated with each cluster refinement, as shown in figure 14. Note that as cluster refinement progresses we see a general broadening in peak stress distributions, resulting from the greater fidelity in error field approximation offered by higher dimensional projection spaces. Furthermore, we note a refinement of the distribution mean and standard deviation, again resulting from the great degree of variation captured in the refined cluster projections. In the example problem, it is clear that convergence of the sub spaces is oscillatory (as opposed to monotonic). It is unclear why this is the case, however we suggest that it is a symptom of the relatively coarse macro model mesh. Given the limited number of DoFs driving the micro model (and therefore the small number of possible QoIs) it is conceivable that, given some particular projection space, the propagation of certain QoIs over regions that in reality are not controlled by these quantities happens to result in preferable micro states. This feature of the low rank approximation is not explored here, however the convergent nature of the process should be noted. This satisfies the original motivation for the goal oriented projection.

Perspectives and Conclusions

A statistical learning approach is presented here for the estimation of full error fields in large FE models. The method is unobtrusive -it does not require any re-meshing operation and is computationally efficient, owing to the low-rank SPDE-based gaussian process specifically developed in this work. The stochastic PDE representation of the prior error field allows for posterior realisations to be generated using sparse matrix inversions, thereby ensuring the cost of generating such realisations, and that of tuning the GP hyperparameters using maximum marginal likelihood, is comparable to the cost of solving the FEA problem itself. It is important to note that the adoption of a stochastic PDE representation of the prior is of little help if done in isolation. It is the sparse approximation of this prior through a guided reduced basis projection that alleviates computational bottlenecks and enables the overall methodology to be computationally affordable.

Although attention is given to classical sub-modelling multiscale problems, the methods are applicable to many alternative cases. For example, consider a large domain where there exist many potential localisation features but the analyst does not know which is of the greatest concern. Consider also that, due to some limit on computational overhead, the mesh cannot be refined at all localisation features. The error field recovery approach outlined here would furnish the analyst with confidence intervals for QoIs at the localisation features, thereby indicating which should be investigated in "dive deeper" exercises. While it is highly unlikely for an analyst to be completely ignorant of potential RoIs, it is very likely in complex systems that many potential regions at the macro scale and features at the micro scale will be of interest. There is a need, therefore, to extract uncertainty information across the macro model such that this can be propagated down to various micro models. This observation drives the creation of approaches like the one developed in the present work. A strong constraint on computational cost is industry inspired and motivates the development of an efficient analysis method; a consideration often omitted in similar studies. Attention has also been limited to small strain linear (elastic) problems here, however the methods may be readily extended to more complicated non-linear conditions. GOEE concepts may, for example, be extended to non-linear problems by linearising constitutive matrices (as in the work of Cirak [START_REF] Cirak | A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem[END_REF] or that of Ghorashi and Rabczuk [START_REF] Sh | Goal-oriented error estimation and mesh adaptivity in 3d elastoplasticity problems[END_REF]). The integration of non-linearity at local levels is straightforward. Duals solutions (at the macro scale) can be updated with the iterative primal solution, with each leading to a modification of the recovered error field. Note that QoIs do not need to be the same between solution increments in the non-linear case. State estimation in space and time can provide inspiration, with dual solutions from previous increments being used to inform subsequent error field estimations [START_REF] Girolami | The statistical finite element method (statfem) for coherent synthesis of observation data and model predictions[END_REF]. Active learning can be implemented here to strategically define QoIs between increments. QoIs are defined arbitrarily in the present work, however future work should look to optimise the QoIs through adaptive learning/Greedy algorithms. Note that this can be accomplished in a simplistic way by characterising QoIs with the problem shape functions, however a more intriguing question lies in the full optimisation of QoIs in the problem domain.

Figure 1 :

 1 Figure 1: The example macroscopic system, showing (a) boundary conditions and loads, and (b) the deformed (red) and undeformed (black) mesh (with a unity deformation scaling factor).

Figure 2 :

 2 Figure 2: The example multiscale problem considered in the present work. Two potential RoI boundaries are defined (in red and blue), within which a multitude of micro models (that contain local features) can be placed, thereby motivating a description of error in driving DoFs on the RoI boundaries.

Figure 3 :

 3 Figure 3: A sampling from GOEE field distributions, showing convergence of the solutions with increasing numbers of training data sets (duals). Note the same noise vector, shown in (a), is used in all cases. Sub figures (b)-(e) show the stochastic differential equation determined Matérn prior, a posterior trained using 2 randomly selected duals, a posterior trained using 16 randomly selected duals, and a posterior trained using all duals on the RoI boundaries, respectively.

Figure 4 :

 4 Figure 4: Monte-Carlo error field propagation for the "2 hole" (position 2) microscopic model, showing (a) the mean field and (b) the standard deviation field. Results are based on a posterior trained using 2 randomly selected duals from the RoI boundaries.

Figure 5 :

 5 Figure 5: Monte-Carlo error field propagation for the "2 hole" (position 2) microscopic model, showing (a) the mean field and (b) the standard deviation field. Results are based on a posterior trained using 16 randomly selected duals from the RoI boundaries.

Figure 6 :

 6 Figure 6: Monte-Carlo error field propagation for the "2 hole" (position 2) microscopic model, showing (a) the mean field and (b) the standard deviation field. Results are based on a posterior trained using all duals from the RoI boundaries.

Figure 7 :

 7 Figure 7: Peak von Mises equivalent stresses (σ V M ), taken from a series of Monte-Carlo error field propagations for the "2 hole" (position 2) microscopic model. Microscopic models are driven by DoF fields modified by posterior distributions trained using varying numbers of dual solutions.

Figure 8 :

 8 Figure 8: Alternative peak von Mises equivalent stresses (σ V M ), taken from a series of Monte-Carlo error field propagations for the "2 hole" (position 2) microscopic model. Microscopic models are driven by DoF fields modified by posterior distributions trained using varying numbers of dual solutions. Note that, while selected duals are consistent in any one Monte-Carlo simulation, different randomly selected duals are used in subfigure (a)-(c).

Figure 9 :

 9 Figure 9: Relative entropy (Kullback-Leibler divergence) for Monte-Carlo developed posterior distributions (i.e. posterior distributions are developed using a random selection of duals from the RoI boundaries).

  n b /d is composed of the n b /d centers of the set of radial bases functions {η i } i∈ 1,n b /d ∈ R n b /d . Set {ω i } i∈ 1,n b /d ∈ R n b /d is a set of weights that allows the approximation to reproduce constant vectors, and {l i } i∈ 1,n b /d ∈ R n b /d is the set of length-scales associated with the RBFs. Sets {c i } i∈ 1,n b /d , {ω i } i∈ 1,n b /d and {l i } i∈ 1,n b /d are chosen as follows.

Figure 10 :

 10 Figure 10: Agglomerative hierarchical clustering development for the example macroscopic problem, using 10 initial clusters. The crosses indicate the centroids of each cluster of vertices. Clusters selected by the defined denisty quantity are refined into 3 sub-clusters using K-means at each cluster step.

Figure 11 :

 11 Figure 11: Clustering density (see equation 67) evolution during cluster refinement for the example macroscopic problem.

Figure 12 :

 12 Figure 12: Agglomerative hierarchical clustering development for the example macroscopic problem, using 10 initial clusters, showing RBF definition for each cluster.

Figure 13 :

 13 Figure13: Entropy of the discrepancy between exact error field on the union of all micro boundaries and its approximation using adaptive RBFs.

Figure 14 :

 14 Figure14: Peak von Mises equivalent stresses (σ V M ), taken from a series of Monte-Carlo error field propagations for the "2 hole" (position 2) microscopic model. Error fields that modify driving DoFs result from posterior distributions trained with 25 dual solutions. Projections are made using the cluster based RBF functions (with higher clustering numbers indicating subsequent clusterings that have made use of density driven refinement). A true "No RBFs" distribution is include for reference (this does not make use of any form of model order reduction in the calculation of the prior).

This can be seen as a maximum likelihood statement. The true error in QoI is a Gaussian with zero mean, centered on the ZZ estimate. The distance between the FE estimate of the QoI and the ZZ estimate then provides one data point from which the variance of this Gaussian distribution may be estimated, resulting in the methodology described previously.

The formulation described above corresponds to the introduction a continuous Gaussian process and its subsequent restriction to the vertices of the FE mesh

[START_REF] Marco | Multiscale modeling: A Bayesian perspective[END_REF] Let us indicate that the continuous counterpart of the prior introduced previously reads as

Galerkin orthogonality:∀i ∈ 1, n Φ , Φ i T r = 0(63)