Peiyi Chen

Johann Guilleminot

Christian Soize

Concurrent Multiscale Simulations of Nonlinear Random Materials Using Probabilistic Learning

Keywords: Concurrent Methods, Nonlinear Elasticity, Probabilistic Learning, Random Media, Surrogates, Uncertainty Quantification

This work is concerned with the construction of statistical surrogates for concurrent multiscale modeling in structures comprising nonlinear random materials. The development of surrogates approximating a homogenization operator is a fairly classical topic that has been addressed through various methods, including polynomial-and deep-learning-based models. Such approaches, and their extensions to probabilistic settings, remain expensive and hard to deploy when the nonlinear upscaled quantities of interest exhibit large statistical variations (in the case of non-separated scales, for instance) and potential non-locality. The aim of this paper is to present a methodology that addresses this particular setting from the point of view of probabilistic learning. More specifically, we formulate the approximation problem using conditional statistics, and use probabilistic learning on manifolds to draw samples of the nonlinear constitutive model at mesoscale. Two applications, relevant to inverse problem solving and forward propagation, are presented in the context of nonlinear elasticity. We show that the framework enables accurate predictions (in probability law), despite the small amount of training data and the very high levels of nonlinearity and stochasticity in the considered system.

Introduction

Background

Concurrent nonlinear simulations involve the strong coupling between a macroscopic (or structural) formulation and a microscopic description capturing subscale details [START_REF] Ghosh | Two scale analysis of heterogeneous elastic-plastic materials with asymptotic homogenization and voronoi cell finite element model[END_REF][START_REF] Smit | Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi-level finite element modeling[END_REF][START_REF] Feyel | Multiscale fe2 elastoviscoplastic analysis of composite structures[END_REF][START_REF] Feyel | Fe2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials[END_REF][START_REF] Terada | A class of general algorithms for multi-scale analyses of heterogeneous media[END_REF][START_REF] Kouznetsova | Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme[END_REF][START_REF] Geers | Multiscale modeling of microstructure-property relations[END_REF][START_REF] Raju | A review of the fe 2 method for composites[END_REF]. One popular approach is the so-called FE 2 method [START_REF] Feyel | Multiscale fe2 elastoviscoplastic analysis of composite structures[END_REF][START_REF] Feyel | Fe2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials[END_REF], where information (in the form of a deformation gradient and any adapted stress variable) is transferred back-and-forth between quadrature points at the macroscale and statistical or representative volume elements (depending on whether the separation of scales exists or not). While versatile and powerful, such frameworks require significant computational resources that often surpass the capabilities of intermediate-power computers, especially when the underlying behavior is highly nonlinear. In this context, the development of surrogate models for largescale systems has become a very active research domain and has generated a substantial body of literature. Various methodologies have been proposed to address this challenge, including (in a non-exhaustive manner) the development of deterministic representations [START_REF] Yvonnet | The reduced model multiscale method (r3m) for the non-linear homogenization of hyperelastic media at finite strains[END_REF][START_REF] Monteiro | Computational homogenization for nonlinear conduction in heterogeneous materials using model reduction[END_REF][START_REF] Yvonnet | Numerically explicit potentials for the homogenization of nonlinear elastic heterogeneous materials[END_REF][START_REF] Temizer | An adaptive method for homogenization in orthotropic nonlinear elasticity[END_REF][START_REF] Temizer | A numerical method for homogenization in non-linear elasticity[END_REF][START_REF] Fritzen | Reduced basis homogenization of viscoelastic composites[END_REF][START_REF] Bhattacharjee | A nonlinear manifold-based reduced order model for multiscale analysis of heterogeneous hyperelastic materials[END_REF][START_REF] Yvonnet | Computational homogenization method and reduced database model for hyperelastic heterogeneous structures[END_REF][START_REF] Fritzen | Two-stage data-driven homogenization for nonlinear solids using a reduced order model[END_REF][START_REF] Liu | Self-consistent clustering analysis: an efficient multi-scale scheme for inelastic heterogeneous materials[END_REF][START_REF] Han | Efficient multiscale modeling for woven composites based on self-consistent clustering analysis[END_REF][START_REF] Feng | A concurrent multiscale framework based on self-consistent clustering analysis for cylinder structure under uniaxial loading condition[END_REF] and probabilistic/statistical-based approaches [START_REF] Efendiev | Multiscale finite element methods for nonlinear problems and their applications[END_REF][START_REF] Efendiev | Multiscale finite element methods: theory and applications[END_REF][START_REF] Weinan | The heterognous multiscale methods[END_REF][START_REF] Weinan | Heterogeneous multiscale methods: a review[END_REF][START_REF] Abdulle | Adaptive reduced basis finite element heterogeneous multiscale method[END_REF][START_REF] Abdulle | The heterogeneous multiscale method[END_REF], and more recently, the integration of machine learning (ML) tools, both with and without probabilistic/statistical formulations; see, e.g., [START_REF] Le | Computational homogenization of nonlinear elastic materials using neural networks[END_REF][START_REF] Rao | Three-dimensional convolutional neural network (3d-cnn) for heterogeneous material homogenization[END_REF][START_REF] Lu | A data-driven computational homogenization method based on neural networks for the nonlinear anisotropic electrical response of graphene/polymer nanocomposites[END_REF][START_REF] Leung | Nh-pinn: Neural homogenization-based physics-informed neural network for multiscale problems[END_REF][START_REF] Nguyen-Thanh | A surrogate model for computational homogenization of elastostatics at finite strain using the hdmr-based neural network approximator[END_REF][START_REF] Nguyen-Thanh | A surrogate model for computational homogenization of elastostatics at finite strain using high-dimensional model representation-based neural network[END_REF][START_REF] Feng | Finite strain fe2 analysis with data-driven homogenization using deep neural networks[END_REF][START_REF] Han | A neural network approach for homogenization of multiscale problems[END_REF][START_REF] Wessels | Computational homogenization using convolutional neural networks[END_REF][START_REF] Black | Deep neural networks for parameterized homogenization in concurrent multiscale structural optimization[END_REF].

In most of the above contributions, the multiscale surrogate is built either through polynomial approximations or deep learning models, and applied locally at each point of the coarse scale discretization. In these settings, the intrinsic randomness -and potential non-locality -induced by random media with non-separated scales is very challenging to capture due to representation limitations. In this work, we explore an alternative path to address this problem and seek to construct a statistical surrogate model where the forward map of interest (specifically, the non-local constitutive model) is approximated using statistical conditioning. Instead of calibrating a regression model between the input (e.g., the deformation gradient) and the output (say, a stress measure), we aim to directly generate samples from the input-output joint probability measure, and to estimate quantities of interest through conditional means. This viewpoint requires the use of a generative model capable of accurately capturing measure concentration and the (unknown) geometry of the support of the measure in the small data limit -a task that remains particularly challenging for strongly non-Gaussian distributions in high dimensions. We note that the construction of generative models is a vibrant topic across many scientific communities, and providing an extensive review on existing techniques is beyond the scope of this paper. In the present study, we employ probabilistic learning on manifolds (PLoM) [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF] to perform this task. The choice of this technique is motivated by (i) its capability to sample the probability measure defined by the training dataset and in particular, to respect measure concentration and support information (as demonstrated in [START_REF] Farhat | Feasible probabilistic learning method for model-form uncertainty quantification in vibration analysis[END_REF][START_REF] Ghanem | Design optimization of a scramjet under uncertainty using probabilistic learning on manifolds[END_REF][START_REF] Arnst | Computation of sobol indices in global sensitivity analysis from small data sets by probabilistic learning on manifolds[END_REF][START_REF] Capiez-Lernout | Nonlinear stochastic dynamics of detuned bladed disks with uncertain mistuning and detuning optimization using a probabilistic machine learning tool[END_REF][START_REF] Ghanem | Probabilistic learning and updating of a digital twin for composite material systems[END_REF][START_REF] Safta | Trajectory design via unsupervised probabilistic learning on optimal manifolds[END_REF][START_REF] Sinha | Aeroacoustic liner impedance metamodel from simulation and experimental data using probabilistic learning[END_REF][START_REF] Zhong | Surrogate modeling of structural seismic response using Probabilistic Learning on Manifolds[END_REF][START_REF] Ezvan | Updating an uncertain and expensive computational model in structural dynamics based on one single target frf using a probabilistic learning tool[END_REF][START_REF] Almeida | A probabilistic learning approach applied to the optimization of wake steering in wind farms[END_REF][START_REF] Zhong | Surrogate modeling of structural seismic response using probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic-learning-based stochastic surrogate model from small incomplete datasets for nonlinear dynamical systems[END_REF]), (ii) relative ease of implementation, and (iii) its reliance on low-dimensional, interpretable parameterization. Our main contributions are as follows. First, we formulate the approximation of the non-local homogenized response in nonlinear elasticity as a learning problem. Second, we perform extensive numerical studies and address the validation of the framework under two scenarios relevant to inverse problem solving and forward propagation. In the former case, the approach can be used, for instance, to calibrate hyperparameters in the material model at fine scale, integrating data at the coarse scale. The latter case represents the classical surrogate setting with aleatoric uncertainties induced by subscale randomness (without separation of scales). Notice that while the proposed developments are derived in the context of nonlinear elasticity, they remain applicable to other classes of constitutive models -at the expense of adapting the mechanistic parameterization.

This paper is organized as follows. The multiscale mechanistic framework is first introduced in Section 2. The deterministic scale-coupling problems (and their stochastic counterparts) are presented, together with the stochastic model enabling the representation of material randomness at mesoscale. Section 3 provides a comprehensive overview on the probabilistic learning framework, including both theoretical and algorithmic aspects. In Section 4, the proposed framework is applied in the context of finite elasticity. The two aforementioned scenarios are specifically introduced to assess the robustness of the method (in probability law). Concluding comments are finally provided in Section 5.

Main Notation (i) Conventions for variables.

A lower-case Latin or Greek letter, such as x or η, is a deterministic real variable. A boldface lower-case Latin or Greek letter, such as x or η, is a deterministic vector. An upper-case Latin or Greek letter, such as X or Ξ, is a real-valued random variable. A boldface upper-case Latin letter, such as X, is a vector-valued random variable. A lower-or upper-case Latin letter between brackets, such as [x] or [X], is a deterministic matrix. A boldface upper-case letter between brackets, such as [X], is a matrix-valued random variable.

(ii) Probability space, random variable, probability measure, and probability density function. For any finite integer m ≥ 1, the Euclidean space R m is equipped with the σ-algebra B R m . If Y is a R m -valued random variable defined on the probability space (Θ,

T , P), Y is a mapping θ → Y(θ) from Θ into R m , measurable from (Θ, T ) into (R m , B R m ),
n . x = (x 1 , . . . , x n ): point in R n . ⟨x, y⟩ = x 1 y 1 + . . . + x n y n : inner product in R n . ∥ x ∥: norm in R n such that ∥ x ∥ = ⟨x, x⟩. [x] T : transpose of matrix [x]. ∥ [x] ∥: Frobenius norm of matrix [x]. δ kk ′ : Kronecker's symbol.

Description of the Mechanistic Framework

Definition of the Structural Problem

Let Ω str be an open bounded domain in R d (here, d = 2 without loss of methodological generality) representing the reference configuration for the structure of interest, and denote by ∂Ω str the boundary of Ω str . For any material point x ∈ Ω str , the spatial point x φ in the deformed configuration Ω φ str is given by x φ = φ(x), where φ is the deformation map. To make the presentation concrete, we assume that the material (at fine scale) is hyperelastic, compressible and isotropic. For any x ∈ Ω str , the deformation gradient [F] is a second-order tensor defined as

[F] = [∇ x x φ ]. The right Cauchy-Green deformation tensor is defined as [C] = [F] T [F],
and the Green-Lagrange strain tensor defined as

[E GL ] = 1 2 ([C] -[I]
). For the sake of simplicity, we consider a Saint Venant-Kirchhoff model, with a strain energy density function given by

ψ([E GL ]) = λ 2 [tr([E GL ])] 2 + µ tr([E GL ] 2 ) , (1)
where λ and µ are the Lamé parameters (see Section 2.3.1 for details). It is well-known that the above Saint Venant-Kirchhoff material is not polyconvex (see, e.g., Section 4.3 in [START_REF] Ciarlet | Mathematical Elasticity: Three-Dimensional Elasticity[END_REF]). The use of this model may thus lead to poor numerical stability and pathological behaviors in general. Such issues were not observed in the applications presented in this paper, given the multiscale surrogate modeling context. In particular, the boundary conditions inherited from the structural boundary value problem did not generate asymptotic behavior. The results supporting the relevance of the proposed methodology (and more specifically, the ability to approximate the homogenized constitutive model) are therefore not expected to be fundamentally affected by this choice. Note also that the proposed approach can accommodate other types of constitutive behaviors, and that the above choice pertaining to the strain energy density function is not expected to impact the methodological results presented in this research. In a general setting, the strong form (resulting from the balance of linear momentum) of the boundary value problem (BVP) in the reference configuration is stated as [START_REF] Ciarlet | Mathematical Elasticity: Three-Dimensional Elasticity[END_REF] Div

[P(x)] + b(x) = 0 , ∀ x ∈ Ω str , (2)
u(x) = u(x) , ∀ x ∈ ∂Ω D str , (3)
[P(x)] • n(x) = t(x) , ∀ x ∈ ∂Ω N str , (4)
where Div denotes the divergence operator in the reference configuration, [P] is the first Piola-Kirchhoff stress tensor defined as

[P] = ∂ψ([F]) ∂[F] , (5)
the vector b is the body force, n is unit vector normal to the boundary in the reference configuration, u and t are given smooth vector fields on the Dirichlet and Neumann boundaries, denoted by ∂Ω D str and ∂Ω N str respectively. The solution to the above problem is classically sought (in an appropriate function space) as a stationary point of the following energy functional [START_REF] Ciarlet | Mathematical Elasticity: Three-Dimensional Elasticity[END_REF][START_REF] Wriggers | Nonlinear finite element methods[END_REF][START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF]:

Π(φ) = B ψ([F]) dV - B b • φ dV - ∂B N t • φ dA . (6)
In this work, we apply a Dirichlet boundary condition u on ∂Ω str (i.e., no traction is applied, and the body force is neglected).

Definition of the Macroscopic Problem in the Context of Concurrent Multiscale Approaches

Let Ω mac ⊂ Ω str denote the reference configuration for the subdomain where the surrogate must be constructed, and denote by ∂Ω mac its boundary (see Fig. 1). Let u mac be the restriction of the solution to the structural problem (defined in the previous section) to the boundary ∂Ω mac . The strong form of the boundary value problem in the reference configuration of Ω mac is stated as

Div [P mac (x)] = 0 , ∀ x ∈ Ω mac , (7) u(x) = u mac (x) , ∀ x ∈ ∂Ω mac , (8)
where [P mac ] is the first Piola-Kirchhoff stress tensor at macroscale. In order to define the multiscale setting, we consider a statistical volume element Ω mes (x) located at point x ∈ Ω mac , and denote by ∂Ω mes the boundary of Ω mes (as shown in Fig. 1). Given a finite element discretization of Ω str and at a given iteration in the nonlinear (Newton-Raphson) solver, the concurrent method proceeds by estimating the deformation gradient [F mac (x q )] at any quadrature point x q in Ω mac , and by evaluating the apparent first Piola-Kirchhoff stress tensor defined as

[P mac (x q )] = ∂ψ mac ([F mac (x q )]; Ω mes (x q )) ∂[F mac (x q )] , (9)
where ψ mac (•; Ω mes (x q )) is the apparent strain energy density function associated with the mesoscopic domain Ω mes (x q ), using localization (through [F mac (x q )]) and homogenization (via ψ mac (•; Ω mes (x q ))). Note that as previously pointed out, scale separation is not enforced and thus, all quantities obtained by upscaling are termed apparent, following the convention introduced by Huet [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF] (see also [START_REF] Ostoja-Starzewski | Microstructural Randomness and Scaling in Mechanics of Materials[END_REF]). In order to compute [P mac ] at each quadrature point x q ∈ Ω mac , we use the FE 2 method [START_REF] Feyel | Fe2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials[END_REF][START_REF] Feyel | A multilevel finite element method (fe2) to describe the response of highly non-linear structures using generalized continua[END_REF] and solve the boundary value problem defined as Div [P mes (x)] = 0 , ∀ x ∈ Ω mes (x q ) , (10)

u(x) = ([F mac (x q )] -[I])x , ∀ x ∈ ∂Ω mes (x q ) , (11)
in the reference configuration of Ω mes (x q ), where [F mac (x q )] is the deformation gradient inherited from the macroscale boundary value problem at x q . The apparent first Piola-Kirchhoff stress tensor at x q is then evaluated as [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF] [P mac (x q )] = 1 |Ω mes (x q )| Ω mes (x q )

[P mes (x)] dx .

As we will explain in the next section, the pairs of associated deformation gradients and first Piola-Kirchhoff stress tensors at all quadrature points in the macroscopic domain Ω mac are then used to compute the pairs of associated right Cauchy-Green deformation tensors and second Piola-Kirchhoff stress tensors, which are the (objective) mechanistic variables considered in the probabilistic learning process introduced in Section 3. Note that while the deformation gradient and the first Piola-Kirchhoff stress tensor could also be used in the learning approach, the choice of the right Cauchy-Green deformation tensor and second Piola-Kirchhoff stress tensor as quantities of interest leads to smaller dimensions since both tensors are symmetric. It should also be pointed out that preserving mechanical variables over the entire domain (namely, Ω mac ) enables the consideration of a nonlocal apparent constitutive model, as opposed to the calibration of a surrogate at one particular point in the domain (which is more relevant to local constitutive models).

Description of Material Uncertainties

Definition of the Stochastic Model

In this section, we detail the construction of the stochastic model for the strain energy density function defined by Eq. [START_REF] Ghosh | Two scale analysis of heterogeneous elastic-plastic materials with asymptotic homogenization and voronoi cell finite element model[END_REF]. Given the scope of this work, which is focused on the learning perspective rather than stochastic modeling, the hyperelastic model is randomized by defining the first Lamé parameter, λ, as a random field. Models enabling the randomization of all parameters in various classes of strain energy density functions can be found in the references provided after Eq. [START_REF] Bhattacharjee | A nonlinear manifold-based reduced order model for multiscale analysis of heterogeneous hyperelastic materials[END_REF], and in [START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF][START_REF] Staber | Stochastic modeling and generation of random fields of elasticity tensors: A unified information-theoretic approach[END_REF] for linear elasticity (for all symmetry classes). We also note that results published elsewhere reporting on the (first-order) marginal cross-correlation of elastic moduli suggest that one latent random field may be sufficient to induce multiscale-informed stochasticity in the isotropic case (depending on the material under consideration; see, e.g., [START_REF] Hun | Stochastic multiscale modeling of crack propagation in random heterogeneous media[END_REF] for a reinforced composite material).

The first Lamé parameter random field is denoted by {λ(x), x ∈ Ω str }, is defined on probability space (Θ, T , P), and takes values in R + \{0} . In this paper, we define {λ(x), x ∈ Ω str } as

λ(x) = H (Ξ(x)) , ∀x ∈ Ω str , (13)
where H denotes a so-called transport map, constructed to enforce admissibility (in the almost sure sense), and {Ξ(x), x ∈ R d } is a centered homogeneous Gaussian random field. This Gaussian field is completely defined by its correlation function (x, x ′ ) → ρ(x, x ′ ) = E{Ξ(x)Ξ(x ′ )}, which is taken as

ρ(x, x ′ ) = d i=1 exp - x i -x ′ i ℓ c 2 , ∀(x, x ′ ) ∈ R d × R d , (14)
for the sake of illustration, with ℓ c a model parameter such that

+∞ 0 exp -(τ/ℓ c ) 2 dτ = L c
, where L c is the spatial correlation length of the Gaussian random field (which is assumed to be independent of the direction, for simplicity). Following the methodology introduced in [START_REF] Soize | Non-gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF] in the context of anisotropic linear elasticity, the transport map is constructed using information theory and the principle of maximum entropy [START_REF] Shannon | A mathematical theory of communication[END_REF][START_REF] Jaynes | Information theory and statistical mechanics i[END_REF][START_REF] Jaynes | Information theory and statistical mechanics ii[END_REF]; see [START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF] for an introduction to concepts and methodologies, as well as [START_REF] Guilleminot | 12 -modeling non-gaussian random fields of material properties in multiscale mechanics of materials[END_REF] for specific results in (linear and nonlinear) mechanics of materials. Specifically, H is defined by imposing that

λ(x) = H (Ξ(x)) ∼ P ME , ∀x ∈ Ω str , (15)
where P ME is the probability measure induced by entropy maximization under constraints. General methodologies and information-theoretic results for a large class of models in nonlinear elasticity can be found in [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF] and [START_REF] Staber | Stochastic modeling of the Ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF], for the cases of isotropic incompressible and compressible materials, respectively. Extensions to spatially-dependent anisotropic hyperelastic models can be found in [START_REF] Staber | A random field model for anisotropic strain energy functions and its application for uncertainty quantification in vascular mechanics[END_REF][START_REF] Chen | Spatially-dependent material uncertainties in anisotropic nonlinear elasticity: Stochastic modeling, identification, and propagation[END_REF], and applications including calibration and validation using experimental data are available in [START_REF] Staber | Stochastic hyperelastic constitutive laws and identification procedure for soft biological tissues with intrinsic variability[END_REF][START_REF] Staber | Stochastic modeling and identification of a hyperelastic constitutive model for laminated composites[END_REF][START_REF] Chen | Spatially-dependent material uncertainties in anisotropic nonlinear elasticity: Stochastic modeling, identification, and propagation[END_REF] (see also [START_REF] Mihai | Stochastic Elasticity: A Nondeterministic Approach to the Nonlinear Field Theory[END_REF] and the references therein for a review of applications to canonical mechanics problems). Since λ corresponds to an elasticity parameter, results obtained in the context of stochastic linearized elasticity can also be invoked. Accounting for the positiveness constraint, as well as for the existence of second-order moments for the linearized elasticity tensor and its inverse [START_REF] Soize | Non-gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF][START_REF] Guilleminot | Encyclopedia of Continuum Mechanics[END_REF], it can be shown that P ME corresponds to a Gamma distribution. Denoting by λ and δ λ the mean and coefficient of variation of λ, it follows that

H = F -1 G(δ -2 λ ,λδ -2 λ ) • F N(0,1) , (16)
where F -1 G is the inverse cumulative distribution function of the Gamma distribution with shape and scale parameters given by δ -2 λ and λδ -2 λ , respectively, and F N(0,1) is the cumulative distribution function of the standard Gaussian distribution. Notice that these hyperparameters can be made spatially-dependent to improve expressiveness in the model: this sophistication is, however, irrelevant for the objectives pursued in this paper.

In the applications presented below, the underlying Gaussian random field {Ξ(x), x ∈ R d } is sampled using a truncated Karhunen-Loève expansion, with an order of truncation determined such that the L 2 error falls below a given threshold (chosen as 1e-2). Samples for both the Gaussian and non-Gaussian random fields are shown in Fig. 2, for λ = 40 000, δ λ = 0.2, and L c = 0.3.

Definition of the Stochastic Boundary Value Problems

Considering the Lamé parameter random field defined in Section 2.3.1 in the BVPs introduced in Sections 2.1 and 2.2 leads to the definition of stochastic boundary value problems (SBVPs), which are briefly described below for the sake of readability. All equalities below hold in the almost sure sense.

Following the retained modeling setup, the structural stochastic boundary value problem is given by Div

[P(x)] = 0 , ∀ x ∈ Ω str , (17)
U(x) = u(x) , ∀ x ∈ ∂Ω D str , (18)
where [P] is the stochastic Piola-Kirchhoff stress tensor (arising from the randomization of the strain energy density function via λ), {U(x), x ∈ Ω str } is the displacement solution random field and x → u(x) is the known deterministic field introduced in Section 2.1. Similarly, the macroscopic SBVP on the domain of interest Ω mac (where the statistical surrogate is built) writes

Div [P mac (x)] = 0 , ∀ x ∈ Ω mac , (19)
U(x) = U mac (x) , ∀ x ∈ ∂Ω mac , (20)
where {U mac (x), x ∈ ∂Ω mac } is now a random field with values in R d , due to the fact that the background medium is stochastic. Finally, the SBVP considered in the concurrent approach (for any subdomain Ω mes (x q ) centered at quadrature point x q in Ω mac ) is

Div [P mes (x)] = 0 , ∀ x ∈ Ω mes (x q ) , (21)
U(x) = ([F mac (x q )] -[I])x , ∀ x ∈ ∂Ω mes (x q ) , (22)
where [F mac (x q )] is the stochastic deformation gradient at x q , defined through localization.

In this work, we consider the strong stochastic solutions to the weak formulations (using the Galerkin method and a finite element discretization) of the above SBVPs. The Monte Carlo approach is chosen as the stochastic solver. The pairs of right Cauchy-Green deformation tensors and second Piola-Kirchhoff stress tensors (denoted by

{[C mac (x q )], [S mac (x q )]} Q mac q=1
) are collected at all quadrature points, for all samples of the Lamé parameter random field, to constitute the dataset for the probabilistic learning procedure introduced in Section 3.

Implementation Verification for the FE 2 Method (With Deterministic Background Media)

For the sake of illustration, we consider the structural domain Ω str = [0, 3] 2 , and define the macroscopic domain of interest as Ω mac = [1, 2] 2 . The mesoscopic domain Ω mes (x q ) at quadrature point x q is defined by a characteristic length L Ω mes = 0.025. Spatial discretization is realized using Q4 elements at all mesh resolutions. The number of elements per direction is 15 in Ω str , 5 in Ω mac , and 5 in Ω mes . The Dirichlet boundary conditions applied on the boundary of Ω str are given by

u(x) = 0 , x 1 = 0 , ∀x 2 ∈ [0, 3] , (23)
u(x) = 0.1 0 , x 1 = 3 , ∀x 2 ∈ [0, 3] . ( 24
)
As described in the previous section, the solution vector u in Ω str along the boundary ∂Ω mac of Ω mac is applied as the Dirichlet boundary condition for the multi-scale problem. Implementation was performed within the MOOSE finite element framework [START_REF] Permann | MOOSE: Enabling massively parallel multiphysics simulation[END_REF]. A convergence study on the solution on Ω str was first conducted. A manufactured displacement field taken as u MMS (x) = (0.01 sin(y), 0.01 sin(x)) T is considered, with material parameters given by λ = 40 000 [kg/cm 2 ] and µ = 10 000 [kg/cm 2 ] (these values, taken from [START_REF] Picinbono | Non-linear anisotropic elasticity for real-time surgery simulation[END_REF], correspond to a soft biological tissue, modeled as a Saint Venant-Kirchhoff material). Dirichlet boundary conditions in accordance with the above solution are prescribed on all boundaries. A body force is defined such that the manufactured solution corresponds to the nonlinear boundary value problem defined in Section 2.1. Regarding numerical solving, a standard Newton-Raphson solver was used with a maximum number of nonlinear iterations set to 25, with a relative tolerance taken as 1e-10, and an absolute tolerance given by 1e-12. The convergence order is measured by the L 2 -norm of the difference between the approximation u h and the reference solution u ref within the domain [0, 3] 2 . Standard h-convergence is observed, as illustrated in Fig. 3. Next, the implementation of the FE 2 method was verified by comparing the normalized L 2 -norm error between the solution vector (at all nodes) on Ω str without multiscale coupling, and the solution vector (at all nodes) obtained by using the FE 2 method in the subdomain Ω mac = [1, 2] 2 . Fig. 4 shows the first sample of the material random field λ (for λ = 40 000, δ λ = 0.2, and L c = 0.3), as well as solutions to the structural and macroscale problems. In order to perform a statistical analysis on the error, 500 independent samples of λ were generated. The mean of the normalized L 2 error is 4.69e-4, and the coefficient of variation is 0.22. The probability density function of the normalized L 2 norm error is shown in Fig. 5. These results indicate proper implementation of the concurrent multiscale method, which is used to build the dataset for the probabilistic learning technique (introduced in the next section).

Overview of the Probabilistic Learning on Manifolds (PLoM) Algorithm and its Parameterization

In this section, we provide a concise overview of the PLoM algorithm. The reason for providing this review is twofold. First, we aim to assist readers in analyzing and comprehending the underlying parameterization, the values chosen for the parameters, and the results pertaining to algorithmic control and convergence. Second, there is no published paper that summarizes all the ingredients of the PLoM approach that are used in this work. Early developments addressing, for instance, the quantification of probability measure concentration and the estimation of the smoothing parameter in the calculation of the diffusion maps basis and its truncation order, are disseminated in a series of papers (see below). On the other hand, some presented results are new, including the expression of the relaxation parameters as a function of the iteration number (in the implementation of the learning algorithm), and the expression of the drift matrix associated with the normalization condition.

The PLoM approach [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF] starts with the consideration of a training dataset D d , comprising a relatively small number n d of points generated from an underlying stochastic manifold associated with a R n -valued random variable X = (Q, W), defined on a probability space (Θ, T , P). Here, Q represents the quantity of interest and is a R n q -valued random variable, while W denotes the control parameter and is a R n w -valued random variable. The total dimension is n = n q + n w . Another R n u -valued random variable U, defined on (Θ, T , P), is also considered as an uncontrolled parameter. The random variable Q is assumed to be expressed as Q = f(U, W), where the measurable mapping f is not explicitly known. The joint probability distribution P W,U (dw, du) of W and U is assumed to be given. The non-Gaussian probability measure P X (x) = P Q,W (dq, dw) of X = (Q, W) is concentrated in a region of R n , for which the only available information is the cloud of points in the training dataset D d . The PLoM method enables the generation of the learned dataset D ar for X, consisting of n MC ≫ n d points (learned realizations) generated by the non-Gaussian probability measure estimated using the training dataset. The preservation of the probability measure concentration is guaranteed by the utilization of a diffusion-maps basis, which enriches the available information from the training dataset. Utilizing the learned dataset D ar , PLoM enables the computation of conditional statistics, such as w → P O|W (do|W = w), on C w . Here, O = ξ(Q), where ξ is a measurable mapping from R n q into R n o , allowing for the construction of statistical surrogate models (metamodels) within a probabilistic framework. The formulas for the computation of conditional mathematical expectations, conditional probability density functions, and conditional cumulative distribution functions, given any w 0 in C w , are given in Appendix A.

The training dataset D d comprises n d independent realizations

x j d = (q j d , w j d ) in R n = R n q × R n w for j ∈ {1, . . . , n d } of random variable X = (Q, W), in which q j d = f(u j d , w j d ).
The PLoM method allows for generating the learned dataset D ar , consisting of n ar ≫ n d learned realizations {x ℓ ar , ℓ = 1, . . . , n ar } of random vector X. Once the learned dataset is constructed, the learned realizations for Q and W can be extracted as (q ℓ ar , w ℓ ar ) = x ℓ ar for ℓ = 1, . . . , n ar .

Construction of a Reduced Representation

The n d independent realizations {x j d , j = 1, . . . , n d } are represented by the matrix [

x d ] = [x 1 d . . . x n d d ] in M n,n d . Let [X] = [X 1 , . . . , X n d
] be the random matrix with values in M n,n d , where its columns are n d independent copies of random vector X. Utilizing Principal Component Analysis (PCA) of X, random matrix [X] is written as,

[X] = [x] + [φ] [µ] 1/2 [H] , (25)
where

[H] = [H 1 , . . . , H n d ] is a M ν,n d -valued random matrix (ν ≤ n),
] = [η 1 d . . . η n d d ] ∈ M ν,n d of [H] is computed by [η d ] = [µ] -1/2 [φ] T ([x d ] -[x]
). The value ν is classically calculated in order that the L 2 -error function ν → err X (ν), defined by

err X (ν) = 1 - ν α=1 µ α E{∥X∥ 2 } , (26)
be smaller than ε PCA . If ν < n -1, statistical reduction occurs.

Probability Measure of H

The probability measure P H of H has to be estimated in order to construct the probability measure of random matrix [H] used in the PLoM methodology. Let P H (dη) = p H (η) dη be the probability measure on R ν of H, whose probability density function η → p H (η) on R ν is estimated by using the Gaussian kernel-density estimation (KDE) with the training dataset D train (η) = {η j , j = 1, . . . , n d },

p H (η) = 1 n d n d j=1 1 ( √ 2π ŝ ) ν exp - 1 2 ŝ2 ∥ ŝ s ν η j -η ∥ 2 , ∀η ∈ R ν . (27)
In these equations, ŝν is a modification of the standard Silverman bandwidth s ν , defined by

ŝν = s ν s 2 ν + n d -1 n d , s ν = 4 n d (2 + ν) 1/(ν+4)
.

With such a modification, the normalization of H is preserved for any value of n d , that is to say,

E{H} = R ν η p H (η) dη = 1 2 ŝ2 η = 0 ν , E{H ⊗ H} = R ν (η ⊗ η) p H (η) dη = ŝ2 [I ν ] + ŝ2 s 2 ν (n d -1) n d [ C H ] = [I ν ] ,
where η ∈ R ν and [ C H ] ∈ M + ν are the estimates of the mean value and the covariance matrix of H, performed with D train (η). Theorem 3.1 in [START_REF] Soize | Probabilistic learning on manifolds[END_REF] proves that, for all η fixed in R ν , Eq. ( 27) is a consistent estimation of the sequence {p H } n d for n d → +∞.

Development of a Reduced-Order Basis Using Diffusion Maps

To preserve the concentration of the learned realizations in the region where the points of the training dataset are concentrated, PLoM relies on an algebraic basis in vector space R n d , constructed using the diffusion-maps basis [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF]. Let [K] and [b] be matrices such that, for all i and j in {1, . . . , The eigenvector ψ 1 associated with λ 1 = 1 is a constant vector. For a given integer κ ≥ 0, the diffusion-maps basis {g 1 , . . . , g α , . . . , g n d } forms a vector basis of R n d defined by g α = λ κ α ψ α . The reduced-order diffusion-maps basis of order m is defined, for a given integer m, as the set {g 1 , . . . , g m }, represented by the matrix

n d }, [K] i j = exp{-(4 ε DM ) -1 ∥η i d -η j d ∥ 2 } and [b] i j = δ i j b i with b i = n d j=1 [K] i j ,
[g m ] = [g 1 . . . g m ] ∈ M n d ,m with g α = (g α 1 , . . . , g α n d
) and [g m ] jα = g α j . This basis depends on two parameters, ε DM and m, which need to be identified. It is proven in [START_REF] Soize | Probabilistic learning on manifolds[END_REF], that the PLoM method does not depend on κ, which can therefore be chosen to 0. We aim to determine the optimal value m opt ≤ n d for m and the smallest value ε opt > 0 for ε DM such that (see [START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF])

1 = λ 1 > λ 2 (ε opt ) ≃ . . . ≃ λ m opt (ε opt ) ≫ λ m opt +1 (ε opt ) ≥ . . . ≥ λ n d (ε opt ) > 0 , (28)
with an amplitude jump equal to an order of magnitude (a factor 10, as demonstrated in [START_REF] Soize | Probabilistic learning on manifolds[END_REF]) between λ m opt (ε opt ) and λ m opt +1 (ε opt ). A more detailed analysis leads to the following algorithm for estimating ε opt and m opt . Let ε DM → Jump(ε DM ) be the function on ]0, +∞[ defined by

Jump(ε DM ) = λ m opt +1 (ε DM )/λ 2 (ε DM ) . ( 29
)
The algorithm is as follows: set the value of m to m opt = ν + 1 and identify the smallest possible value ε opt of ε DM such that Jump(ε opt ) ≤ 0.1 and Eq. ( 28) is satisfied. ar , ℓ ′ = 1, . . . , n ar } of H. These learned realizations enable the estimation of converged conditional statistics, which are then utilized to construct statistical surrogate models related to X, and subsequently, to (Q, W).

Quantifying the Probability Measure Concentration of Random Matrix [H m opt ]

For m ≤ n d , the probability measure concentration of random matrix [H m ] is defined (see [START_REF] Soize | Probabilistic learning on manifolds (PLoM) with partition[END_REF]) by:

d 2 n d (m) = E{∥[H m ] -[η d ]∥ 2 }/∥[η d ]∥ 2 . ( 30
)
Let M opt = {m opt , m opt + 1, . . . , n d }, where m opt represents the optimal value of m as defined earlier. Theorem 7.8 of [START_REF] Soize | Probabilistic learning on manifolds[END_REF] shows that min

m∈M opt d 2 n d (m) ≤ 1 + m opt /(n d -1) < d 2 n d (n d ),
indicating that the PLoM method, for m = m opt and [g m opt ], is a better method than the standard one corresponding to

d 2 n d (n d ) = 1 + n d /(n d -1) ≃ 2. Using the n MC realizations {[η ℓ ar ], ℓ = 1, . . . , n MC } of [H m opt ], we have the estimate d 2 n d (m opt ) ≃ (1/n MC ) n MC ℓ=1 {∥[η ℓ ar ] -[η d ]∥ 2 }/∥[η d ]∥ 2 .
3.6. Generation of Learned Realizations {η ℓ ′ ar , ℓ ′ = 1, . . . , n ar } for the Random Vector H The MCMC generator is detailed in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF]. Let {([Z(t)], [Y(t)]), t ∈ R + } be the unique asymptotic (as t → +∞) stationary diffusion stochastic process with values in M ν,m opt × M ν,m opt , representing the following reduced-order ISDE (stochastic nonlinear second-order dissipative Hamiltonian dynamic system), for t > 0,

d[Z(t)] = [Y(t)] dt , d[Y(t)] = [L([Z(t)])] dt - 1 2 f 0 [Y(t)] dt + f 0 [dW wien (t)] , with [Z(0)] = [η d ] [a] and [Y(0)] = [N] [a], in which [a] = [g m opt ] ([g m opt ] T [g m opt ]) -1 ∈ M n d ,m opt . (i) [L([Z(t)])] = [L([Z(t)] [g m opt ] T )] [a] is a random matrix with values in M ν,m opt . For all [u] = [u 1 . . . u n d ] in M ν,n d with u j = (u j 1 , . . . , u j ν ) in R ν , the matrix [L([u]
)] in M ν,n d is defined, for all k = 1, . . . , ν and for all j = 1, . . . , n d , by (v) The free parameter f 0 , satisfying 0 < f 0 < 4/ ŝν , allows for the control of the dissipation term in the nonlinear second-order dynamic system (a dissipative Hamiltonian system) to quickly damp the transient effects induced by the initial conditions. A commonly used value is f 0 = 4 (noting that ŝν < 1). Consequently, the ISDE is solved over the interval ]0, T ], where T depends on f 0 and represents the smallest integration final time allowing [Z m opt ] to be chosen as [Z(T )] while being in the stationary regime. (vi) The learned realizations {x ℓ ′ ar , ℓ ′ = 1, . . . , n ar } of random vector X are then obtained by reshaping the realizations

[L([u])] k j = 1 p(u j ) {∇ u j p(u j )} k , (31)
p(u j ) = 1 n d n d j ′ =1 exp{- 1 2 ŝ 2 ν ∥ ŝν s ν η j ′ -u j ∥ 2 } , ∇ u j p(u j ) = 1 ŝ 2 ν n d n d j ′ =1 ( ŝν s ν η j ′ -u j ) exp{- 1 2 ŝ 2 ν ∥ ŝν s ν η j ′ -u j ∥ 2 } . (ii) [W wien (t)] = [W wien (t)] [a] where {[W wien (t)], t ∈ R + } is the M ν,n d -
{[x ℓ ar ] = [x] + [φ] [µ] 1/2 [η ℓ
ar ] , ℓ = 1, . . . , n MC } (see Eq. ( 25)) with n ar = n MC × n d .

Preserving Normalization Conditions Through Constraints on Second-Order Moments of H

In general, the mean value of H estimated using the n ar learned realizations {η ℓ ′ ar , ℓ ′ = 1, . . . , n ar }, is sufficiently close to zero. Similarly, the estimate of the covariance matrix of H is also sufficiently close to the identity matrix. However, there are instances where the mean value may not be sufficiently small, and the diagonal entries of the estimated covariance matrix can fall below 1. The normalization conditions can be reestablished during the learning algorithm described in Section 3.6 by imposing, for all k = 1, . . . , ν, the constraints E{H k } = 0 and E{(H k ) 2 } = 1. These constraints can be globally rewritten as (i) Methodology for imposing the constraint in the learning algorithm. We apply the Kullback-Leibler minimum cross-entropy principle as proposed in [START_REF] Soize | Physics-constrained non-Gaussian probabilistic learning on manifolds[END_REF][START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF]. Let H c be the R ν -valued random variable that satisfies the constraint defined by Eq. ( 31), expressed as E{h(H c )} = b. The learned probability measure P H c (dη) = p H c (η) dη, represented by a density p H c on R ν , which satisfies the constraint and which is closest to p H defined by Eq. ( 27), is the solution of the following optimization problem,

E{h(H)} = b on R n c , (32)
p H c = arg min p∈C ad,p R ν p(η) log p(η) p H (η) dη , (33)
in which the admissible set C ad,p is defined by

C ad,p = η → p(η) : R ν → R + , R ν p(η) dη = 1 , R ν h(η) p(η) dη = b . ( 34
)
It is proven that there exists a unique solution to the optimization problem defined by Eqs. ( 33) and [START_REF] Han | A neural network approach for homogenization of multiscale problems[END_REF], which is reformulated using Lagrange multipliers to account for the constraints in the admissible set (refer to Proposition 1 in [START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF] for the construction of the probability measure of H c and the proof of its existence and uniqueness).

(ii) Learning algorithm implementation. To take into account the constraints in the learning algorithm defined in Section 3.6, Eq. ( 31) is replaced by the following one,

[L λ ([u])] k j = 1 p(u j ) {∇ u j p(u j )} k -λ k -2 λ k+ν u j k .
in which the Lagrange multiplier λ ∈ R n c , associated with the constraint defined by Eq. ( 32), is calculated using an iteration algorithm (see [START_REF] Soize | Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints[END_REF]). At each iteration i, the value of λ is denoted by λ i and the corresponding random vector H c is denoted by H λ i . The value λ i+1 is computed as a function of λ i by

λ i+1 = λ i -α i [Γ ′′ (λ i )] -1 Γ ′ (λ i ) , i ≥ 0 , λ 0 = 0 n c , in which Γ ′ (λ i ) = b -E{h(H λ i )} and [Γ ′′ (λ i )] = [cov{h(H λ i )}
] (the covariance matrix). The positive coefficient α i is a relaxation parameter (less than 1) that is introduced for controlling the convergence of the iteration algorithm. For given i 2 ≥ 2, for given β 1 and

β 2 such that 0 < β 1 < β 2 ≤ 1, α i is defined, for i ≤ i 2 , by α i = β 1 +(β 2 -β 1 )(i-1)/(i 2 -

Application

Notation and Definition of Case Studies

In order to deploy and assess the performance of the PLoM approach, two scenarios are introduced as follows.

• In the first scenario, relevant to inverse problem solving, the stochastic boundary conditions on Ω mac and the hyperparameters in the random field model (see Section 2.3) are used as control parameters. Let W disp be the random vector corresponding to the discretization of the random field {U mac (x), x ∈ ∂Ω mac } (see Section 2.3.2), and let W hyp be the random vector associated with the randomization of the mean value, the coefficient of variation, and the correlation length of λ. The control variables are W hyp and W disp .

• In the second scenario, related to propagation, the hyperparameters for the random field model are set to their mean values and are not considered as control variables. The latter only comprise the boundary displacements on Ω mac , gathered in W disp .

In the context of statistical surrogate modeling, our goal is then to estimate conditional distributions for the mechanistic variables in the homogenized constitutive model (namely, the right Cauchy-Green deformation tensor and the second Piola-Kirchhoff stress tensor), at any quadrature point of interest in the region where concurrent multiscale coupling is considered, given specific values of the control parameters. For the first scenario, this would enable, through the formulation of an ad hoc statistical inverse problem, the identification of the hyperparameters in the material random field model, assuming that (e.g., experimental) data on the above tensors are available. In the second scenario, this would allow for the estimation of the stress and strain variables (defined as conditional means) for given values of the (Dirichlet) boundary conditions on the macroscopic domain -a statistical extension to standard surrogate modeling where the mapping [C mac (x q , θ)] → [S mac (x q , θ)], 1 ≤ q ≤ Q mac would be approximated using, e.g., regression.

For the sake of illustration, non-informative prior models are chosen in the first scenario. More specifically, the mean, coefficient of variation, and the correlation lengths are assumed to be statistically independent and uniformly distributed on [32000, 48000], [0.1, 0.3], and [0.2, 0.4], respectively. The truncation order in the Karhunen-Loève expansion is computed for each sample of the correlation length, using the same threshold (1e-2).

The relationships between the notations used for the PLoM algorithm summarized in Section 3 and the notations of mechanical quantities are as follows. We recall that, for the n p = 100 Gauss points referred to by the set of indices I = {i = 1, ..., n p }, {C i j , j = 1, 2, 3} represents the 3 components of the random Cauchy-Green tensor at index point i ∈ I, and {S i j , j = 1, 2, 3} represents the 3 components of the random second Piola-Kirchhoff stress tensor. We will also use the notations C and S for the reshaping of these two random tensors, which are then random vectors with values in R 3n p . Finally, plots of probability density functions are obtained using (non-parametric) kernel density estimators, while the conditional probability density functions are estimated using Eq. ( 31).

(i) Quantity of interest. The R n -valued random variable Q is defined by Q = (S, C) with n q = 2 × (3 n p ) = 600.
(ii) Control parameter. Following the previous discussion, let W disp be the R n w,disp -valued random variable for which the n w,disp = 48 components are the discretized displacements on the boundary. Let W hyp = (W hyp,1 , W hyp,2 , W hyp,3 ) be the hyperparameters that control the prior probability model of the random medium (defined in Section 2.3). For the construction of the training and learned datasets, the definition of the R n w -control parameter W depends on the scenario.

• Scenario 1: the training and learned datasets are constructed with W = (W hyp , W disp ) and n w = 3 + n w,disp = 3 + 48 = 51. The conditional statistics are constructed for given

W hyp = w hyp,o ∈ R 3 .
• Scenario 2: the training and learned datasets are constructed with W = W disp and n w = n w,disp = 48. The value of W hyp is fixed to the statistical mean value w hyp = (w hyp,1 , w hyp,2 , w hyp,3 ) of the prior probability model of W hyp . The conditional statistics are constructed for given W disp = w disp,o ∈ R 48 .

The random variable U (uncontrolled parameter) corresponds to the stochastic germs in the Karhunen-Loève expansion of the underlying Gaussian random field, which are statistically independent normalized Gaussian random variables.

Conditional Statistics

For the validation of the proposed methodology, the conditional statistics estimated with the learned dataset and those estimated with the validation dataset (considered as a reference) will be compared. Below, we define the considered conditional statistics, including the conditional probability density functions and the conditional mean functions (given the control parameter).

• Scenario 1: for all i ∈ I, for j ∈ {1, 2, 3}, for w o given in R 3 , and for all q in R, we consider

q → p C i j | W hyp (q | w o ) , q → p S i j | W hyp (q | w o ) , i → E{C i j | W hyp = w o } = R q p C i j | W hyp (q | w o ) dq and i → E{S i j | W hyp = w o } = R q p S i j | W hyp (q | w o ) dq .
• Scenario 2: for all i ∈ I, for j ∈ {1, 2, 3}, for w o given in R 48 , and for all q in R, we consider In this section, we define the parameter values used by the PLoM algorithms (as summarized in Section 3), and we present the convergence analysis for both scenarios. Notations are those introduced in Section 3. To simplify referencing with respect to each scenario, the first provided value corresponds to Scenario 1, while the second value corresponds to Scenario 2. For instance, "n c = 280 and 272" means that n c = 280 for Scenario 1 and n c = 272 for Scenario 2. When a single value is given, it applies to both scenarios. For example "n d = 500" means that n d = 500 for Scenario 1 and Scenario 2.

q → p C i j | W disp (q | w o ) , q → p S i j | W disp (q | w o ) , i → E{C i j | W disp = w o } = R q p C i j | W disp (q | w o ) dq , and i → E{S i j | W disp = w o } = R q p S i j | W disp (q | w o ) dq .
(i) Values of the general parameters. The total dimension of X = (Q, W) is n = n q + n w = 651 and 648. The number of points in the training dataset D d is n d = 500.

(ii) Reduced representation and diffusion-maps basis. Fig. 6a displays the graph of the function ν → err X (ν) defined by Eq. ( 26). For Scenario 1, convergence of the representation is achieved at ν = 140, resulting in an error of err X (140) = 2.85 × 10 -4 . In Scenario 2, convergence occurs at ν = 136, with an error of err X (136) = 2.99 × 10 -4 . Regarding the computation of the diffusion-maps basis introduced in Section 3.3, the optimal value of the smoothing parameter ε DM is determined as ε opt = 454 and 278, corresponding to the optimal value m opt = 141 and 137 of parameter m. (iii) Generating the learned realization. The learned realizations are generated as explained in Section 3.6, incorporating the constraints outlined in Section 3.7. The free parameter f 0 is chosen as 4, and the integration step ∆t of the Störmer-Verlet scheme is 0.21. Each realization [z ℓ α ] represents the ℓth realization of [Z(T )] for T = 30 × ∆t (due to the damping controlled by f 0 = 4, T is a time at which the stationary response is reached). We have chosen n MC = 100, resulting in n ar = n MC × n d = 50 000.

(iv) Constraints on second-order moments of H. For Scenario 1, the number of constraints is n c = 280, and the relaxation parameter α i is defined by β 1 = 0.01, i 2 = 10, and β 2 = 0.1. For Scenario 2, n c = 272, and β 1 = 0.01, i 2 = 20, and β 2 = 0.5. The convergence of the iterative algorithm presented in Section 3.7-(ii), to take into account the constraints on second-order moments of H in PLoM, as a function of iteration number i, is studied in analyzing the graph of the error function i → err(i) defined by Eq. ( 35) (see Fig. 7a) and the graph of the function i → ∥λ i ∥ (see Fig. 7b). A very good convergence is observed, with err(2000) = 10 -3 (Scenario 1) and err(109) = 3.2×10 -4 (Scenario 2). For both scenarios and for all k in {1, . . . , ν}, at convergence, it holds that |E{H c,k }| ≤ 10 -7 and 0.999 ≤ E{H 2 c,k } ≤ 1. (v) Concentration of the learned probability measure. As explained in Section 3.5, the quality of the PLoM algorithm is assessed by examining the value of d 2 n d (m opt ), defined by Eq. [START_REF] Leung | Nh-pinn: Neural homogenization-based physics-informed neural network for multiscale problems[END_REF]. At convergence, we obtain d 2 n d (m opt ) = 0.17 and 0.16, indicating excellent quality of PLoM to preserve the concentration of the learned probability measure.

(vi) Illustration of the learned pdf of components of H and of the clouds of the learned points. In H c , the subscript c is removed to simplify the writing.

Validation Analysis

In this section, we present a validation of the proposed methodology. This methodology is based on the construction of conditional statistics (statistical surrogate model) defined in Section 4.2, which are estimated with the n ar points of the learned dataset, generated with the PLoM algorithm. The estimates of these conditional statistics are converged because n ar is large. For validation purposes, these conditional statistics must be compared with a reference. This reference can only be obtained by constructing a validation dataset with n v points generated with the nonlinear computational model used to construct the n d points of the training dataset. Ideally, n v ≃ n ar and the construction of a validation dataset can be achieved for both scenarios. Here, due to limitations in computational resources, we only consider the construction of a new validation dataset for Scenario 2. This validation dataset D v comprises n v = 20 000 independent realizations (q ℓ v , w ℓ v ) in R n R n q R n with n q = 600 and n w = 48, for ℓ ∈ {1, . . , n v } of random variable (Q, W) (note that n v < n ar in this case). The constitution of this dataset requires solving 2 040 000 nonlinear finite-element computations overall. The following settings are then considered.

• For Scenario 1: as explained in Section 4. (i) Validation for Scenario 1. Concerning the conditional probability density functions, we select the first components at point #60, corresponding to the random variables C 60,1 and S 60,1 relative to the random tensors C and S, as quantities of interest. This point is relatively central in the macroscopic domain, and the first components at this location present significant fluctuations -hence making this choice a reasonable one from a validation standpoint. Fig. 11 displays the graphs of the conditional probability density functions q → p C 60,1 |W hyp (q|w o ) of random variable C 60,1 and q → p S 60,1 |W hyp (q|w o ) of random variable S 60,1 given W hyp = w 0 , estimated with the learned dataset and the validation dataset. It can be seen that the width of the supports, which control the variances, is well predicted. For the three components indexed by j ∈ {1, 2, 3}, Fig. 12 displays the graphs of the conditional mathematical expectations i → E{C i, j |W hyp = w o } for the family {C i, j , i ∈ I} of random variables, while Fig. 13 displays the graphs of i → E{S i, j |W hyp = w o } for the family {S i, j , i ∈ I}. While very large variations are observed over the quadrature points, the accuracy in the predictions of the conditional expectations remains remarkably good given the small number of training data. It should be noted that the graphs in Figs. 12b and 13b differ only from a constant. for random variable 60,1 given W disp = w 0,k , estimated the learned dataset and the dataset. Similarly, Fig. 15 displays the graphs of q → p S 60,1 |W disp (q|w o,k ) for random variable S 60,1 given W disp = w 0,k . Again, it is seen that the supports of the (marginal) distributions are well predicted, and the validation results are reasonably good.

For j ∈ {1, 2, 3} and the two values of the control parameter (k = 1 and k = 2), Figs. 16 and 17 display the graphs of the conditional mathematical expectations i → E{C i, j |W disp = w o,k } for the family {C i, j , i ∈ I} of random variables, while Figs. 18 and 19 display the graphs of i → E{S i, j |W hyp = w o,k } for the family {S i, j , i ∈ I}. Good accuracy is observed over all quadrature points, which demonstrates the capability of the framework to capture the non-smooth large variations generated by the localization in the nonlinear stochastic boundary value problems.

Remarks on CPU Time

For the probabilistic learning algorithm, the total CPU time is due to the construction of the learned dataset, the numerical cost of conditional statistics being completely negligible. For the construction of learned datasets, the CPU time is directly proportional to the number of iterations required to satisfy the normalization constraints. For the first scenario, 2, 000 are used and the total CPU time is 4 hr and 26 min. For the second scenario, only 109 iterations are necessary to reach excellent convergence and the total CPU time is 19 min. It should be noted that for scenario 1, the convergence tolerance could have been reduced without significantly penalizing the quality of the results, and this would have significantly reduced the number of iterations and therefore CPU time.

Regarding the generation of the training dataset using nonlinear finite element simulations and the FE 2 approach, the average CPU time to complete one simulation is 449.80 s, and one multiscale simulation at any Gauss point (and final increment) takes from 0.5 to 5.401 s (depending on both the sample of the local material properties and the applied boundary conditions). Assuming that conditional statistics would be computed using 50,000 samples, the CPU time associated with a brute force approach (with sequential execution) would then be ≈ 50, 000 × 450 = 6, 250 hr. While this time can be reduced using distributed computing, it remains much larger than the cost induced by the learning method.

Conclusion

In this work, we have introduced a statistical surrogate model for concurrent multiscale simulations involving nonlinear materials with non-separated scales. The methodology combines probabilistic learning on manifolds, a generative model that allows for measure concentration and support information to be accurately captured, with the use of conditional statistics to approximate the mapping between apparent strain and stress variables -namely, the right Cauchy-Green tensor and the second Piola-Kirchhoff stress tensor -at a finite set of points (defining a subregion of interest where concurrent coupling must be deployed). As opposed to standard techniques relying on, e.g., polynomial or neural network surrogates, the proposed approach (1) can readily accommodate the (aleatoric) randomness raised by the multiscale setting, (2) enables the seamless integration of nonlocal interactions through the consideration of joint distributions, and (3) can perform efficiently in the small data regime. Two applications, relevant to inverse problem solving and forward propagation, were presented in the context of nonlinear elasticity. In the first case, the hyperparameters for the prior model defining the random media and boundary conditions at mesoscale are considered as control parameters. This setting can be used, for instance, to identify the hyperparameters when experimental observations are available. In the second application, only mesoscopic boundary displacements are used as control variables. In this case, the prior model at fine scale is fixed, and the effect of material randomness can be quantified. In both cases, large spatial variations, large statistical fluctuations, and strong non-Gaussianity are observed. It was shown that despite these challenges, the framework remains capable of delivering reasonably accurate estimations, even with a fairly limited amount of training data. The information contained in the latter is a critical aspect in the methodology: this information must be rich enough to learn the probability measure and discover the geometry of the manifold defining its support.

  and Y(θ) is a realization (sample) of Y for θ ∈ Θ. The probability distribution of Y is the probability measure P Y (dy) on the measurable set (R m , B R m ) (we will simply say on R m ). The Lebesgue measure on R m is denoted by dy and P Y (dy) = p Y (y) dy, with p Y the probability density function (pdf) on R m of P Y (dy) with respect to dy. Finally, E denotes the mathematical expectation operator. (iii) Algebraic notations. R: set of all the real numbers. R n : Euclidean vector space on R of dimension n. M n,m : set of all the (n × m) real matrices. M n : set of all the square (n × n) real matrices. M + n : set of all the positive-definite symmetric (n × n) real matrices. [I n ]: identity matrix in M

Figure 1 :

 1 Figure 1: Definition of scales in the concurrent multiscale simulations. Fluctuations in the first Lamé parameter λ are introduced at fine scale.

  (a) Sample of {Ξ(x), x ∈ Ω str }. (b) Sample of {λ(x), x ∈ Ω str } computed from (a).

Figure 2 :

 2 Figure 2: Realizations of the underlying Gaussian field (left) and material parameter random field (right).

Figure 3 :

 3 Figure 3: Convergence of the L 2 error (h-refinement) for the reference solution.

Figure 4 :

 4 Figure 4: Sample of the material parameter λ (left) and associated displacement magnitudes for the solutions on Ω str (middle) and Ω mac (right).

Figure 5 :

 5 Figure 5: Probability density function of the normalized L 2 -norm error (estimated with 500 samples).

  and [µ] is the (ν × ν) diagonal matrix of the ν positive eigenvalues of the empirical estimate of the covariance matrix of X. The (n × ν) matrix [φ] consists of the associated eigenvectors such [φ] T [φ] = [I ν ]. The matrix [x] in M n,n d has identical columns, each being equal to the empirical estimate x ∈ R n of the mean value of random vector X. The columns of [H] are n d independent copies of a random vector H with values in R ν , satisfying the normalization conditions, E{H} = 0 ν and E{H ⊗ H} = [I ν ]. The realization [η d

  where ε DM > 0 is a smoothing parameter. Let [P] = [b] -1 [K] be the matrix in M n d , with positive entries, satisfying n d j=1 [P] i j = 1 for all i = 1, . . . , n d . Matrix [P] can be regarded as the transition matrix of a Markov chain that represents the probability of transition in one step. The eigenvalues λ 1 , . . . , λ n d and the associated eigenvectors ψ 1 , . . . , ψ n d of the right-eigenvalue problem [P] ψ α = λ α ψ α satisfy 1 = λ 1 > λ 2 ≥ . . . ≥ λ n d and are computed by solving the generalized eigenvalue problem [K] ψ α = λ α [b] ψ α with the normalization condition ⟨[b] ψ α , ψ β ⟩ = δ αβ .

3. 4 .

 4 Reduced-Order Representation of Random Matrices [H ] and [X ] to Preserve Probability Measure ConcentrationThe diffusion-maps vectors g 1 , . . . , g m ∈ R n d span a subspace of R n d that characterizes, for the optimal values m opt and ε opt of m and ε DM , the local geometry structure of dataset {η j d , j = 1, . . . ,n d }. The PLoM method introduces the M ν,n d -valued random matrix [H m ] = [Z m ] [g m ] T with m ≤ n d , corresponding to a data-reduction representation of random matrix [H], in which [Z m ] is a M ν,m -valued random matrix.The MCMC generator of random matrix [Z m ] is chosen from the class of Hamiltonian Monte Carlo methods, explicitly described in[START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF], and mathematically detailed in Theorem 6.3 of[START_REF] Soize | Probabilistic learning on manifolds[END_REF]. For generating the learned dataset, the best probability measure of [ H m ] is obtained for m = m opt and by using the previously defined basis [g m opt ]. For these optimal quantities m opt and [g m opt ], the generator allows for computing n MC realizations {[z ℓ ar ], ℓ = 1, . . . , n MC } of [Z m opt ] and therefore, for deducing the n MC realizations {[η ℓ ar ], ℓ = 1, . . . , n MC } of [H m opt ]. The reshaping of matrix [η ℓ ar ] ∈ M ν,n d allows for obtaining n ar = n MC × n d learned realizations {η ℓ ′

  in which n c = 2ν. Here, the function h = (h 1 , . . . , h n c ) and the vector b = (b 1 , . . . , b n c ) are defined such that, for all k in {1, . . . , ν}, we have h k (H) = H k , h k+ν (H) = (H k ) 2 , b k = 0, and b k+ν = 1.

4. 3 .

 3 Parameter Values and Convergence Analysis of PLoM Algorithms for Scenarios 1 and 2

  Fig. 6b shows the graph of the function α → λ α representing the eigenvalues of the transition matrix [P].

Figure 6 :

 6 Figure 6: Convergence of the PCA reduced representation (a). Eigenvalues of the transition matrix [P] for the diffusion-maps basis (b). Scenario 1 (thin line), Scenario 2 (thick line).

Figure 7 :

 7 Figure 7: Convergence analysis of the iterative algorithm to take into account the constraints on second-order moments of H in PLoM as a function of iteration number i, presented in log-log scales. Scenario 1 (thin line), Scenario 2 (thick line).

Fig. 8 (

 8 Scenario 1) and Fig.9(Scenario 2) depict the probability density functions of components H 1 , H 30 , and H 70 of H, estimated using the training dataset and the learned dataset. These figures exhibit a good coherence. It should be noted that the convergence of these estimates is not the same, as there are n d = 500 points in the training dataset and n ar = 50 000 points in the learned dataset. Fig.10a(Scenario 1) and Fig.10b(Scenario 2) display the clouds of the learned points and the training points of H, associated with components H 1 , H 30 , and H 70 . These results illustrate the preservation of the concentration of the probability measure.

Figure 8 :

 8 Figure 8: Scenario 1: probability density functions of components H 1 , H 30 , and H 70 of H, estimated with the training dataset (dashed line) and with the learned dataset (solid line).

Figure 9 :

 9 Figure 9: Scenario 2: probability density functions of components H 1 , H 30 , and H 70 of H, estimated with the training dataset (dashed line) and with the learned dataset (solid line).

(a) Scenario 1 :Figure 10 :

 110 Figure 10: Clouds of the learned points and the training points of H, associated with components H 1 , H 30 , and H 70 .

  1, we have W = (W hyp , W disp ) and n w = 3 + 48 = 51. To construct the reference conditional statistics with the validation dataset D v , we can therefore consider a single value n o = 1, w o = w hyp ∈ R 3 of the random control parameter W hyp . The realizations {w ℓ v } ℓ≥1 are not useful, and only the realizations {q ℓ v } ℓ≥1 are used. • For Scenario 2: for the validation of the conditional statistics, the control variable is W = W disp with n w = n w,disp = 48. To construct the reference conditional statistics with the validation dataset D v , we introduce n o values, {w o,k ∈ R n w , k = 1, . . . , n o }, of the random control parameter W with values in R n w . We have randomly drawn, with a uniform distribution, the n o vectors {w o,k , k = 1, . . . , n o } from the set {w ℓ v , ℓ = 1, . . . , n v }. Due to space limitations, we consider n o = 2: these two vectors correspond to the realizations #1 882 and #19 502 in the set {w ℓ v , ℓ = 1, . . . , 20 000}. The validation results are presented below for each scenario.

Figure 11 :

 11 Figure 11: Validation of Scenario 1: Graph of the conditional pdf given W hyp = w 0 of component 1 of the random tensors C and S at point 60. Learned dataset (solid line), validation dataset (dashed line).

Figure 12 :

 12 Figure 12: Validation of Scenario 1: Graph of the conditional mathematical expectation as function of points i given W hyp = w 0 for the three components of the random tensor C. Learned dataset (solid line), validation dataset (dashed line).

(ii) Validation for Scenario 2 .Figure 13 :

 213 Figure 13: Validation of Scenario 1: Graph of the conditional mathematical expectation as function of points i given W hyp = w 0 for the three components of the random tensor S. Learned dataset (solid line), validation dataset (dashed line).

Figure 14 :

 14 Figure 14: Validation of Scenario 2: Graph of the conditional pdf given W disp = w 0,1 or w 0,2 of component 1 of the random tensor C at point 60. Learned dataset (solid line), validation dataset (dashed line).

Figure 15 :

 15 Figure 15: Validation of Scenario 2: Graph of the conditional pdf given W disp = w 0,1 or w 0,2 of component 1 of the random tensor S at point 60. Learned dataset (solid line), validation dataset (dashed line).

Figure 16 :

 16 Figure 16: Validation of Scenario 2: Graph of the conditional mathematical expectation as function of points i given W disp = w 0,1 for the three components of the random tensor C. Learned dataset (solid line), validation dataset (dashed line).

Figure 17 :

 17 Figure 17: Validation of Scenario 2: Graph of the conditional mathematical expectation as function of points i given W disp = w 0,2 for the three components of the random tensor C. Learned dataset (solid line), validation dataset (dashed line).

Figure 18 :

 18 Figure 18: Validation of Scenario 2: Graph of the conditional mathematical expectation as function of points i given W disp = w 0,1 for the three components of the random tensor S. Learned dataset (solid line), validation dataset (dashed line).

Figure 19 :

 19 Figure 19: Validation of Scenario 2: Graph of the conditional mathematical expectation as function of points i given W disp = w 0,2 for the three components of the random tensor S. Learned dataset (solid line), validation dataset (dashed line).

  T . It should be noted that the calculation of the realizations [Z m opt ] is done in parallel computation, each realization [z ℓ ar ] of [Z m opt ] being calculated on a "worker" associated with a realization [W wien,ℓ ] of the Wiener process [W wien ].

valued normalized Wiener process. (iii)

[N]

is the M ν,n d -valued normalized Gaussian random matrix that is independent of process [W wien ]. (iv) We then have [Z m opt ] = lim t→+∞ [Z(t)] in probability distribution. The Störmer-Verlet scheme is used for solving the reduced-order ISDE (see

[START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF]

), which allows for generating the learned realizations, [z 1 ar ], . . . , [z n MC ar ], and then, generating the learned realizations [η 1 ar ], . . . , [η n MC ar ] such that [η ℓ ar ] = [z ℓ ar ] [g m opt ]

Acknowledgments

The work of the P. C. and J. G. was supported by the National Science Foundation, United States, under awards DGE-2022040 and CMMI-1942928.

Appendix A. Formulas for Conditional Statistics

In this Appendix, we use the following notation. The real-valued random variable Q denotes any component of the R n q -valued random variable Q, and the real variable q represents the corresponding component of the vector q in R n q . Let Q and W = ( W 1 , . . . , W n w ) be the normalized random variables defined by

where q, w k , and σ Q , σ W k are the mean values and standard deviations of the random variables Q and W k . These values are estimated using empirical statistical estimators based on the learned realizations {(q ℓ ar , w ℓ ar ), ℓ = 1, . . . , n ar }. For any value w 0 = (w 0,1 , . . . , w 0,n w ) of the control parameter given in C w ⊂ R n w , we define the vector

The Gaussian KDE estimation of the joint probability distribution of Q and W, with respect to d q d w, is written as

In Eq. (A.3), s is the Silverman bandwidth defined by

and, for ℓ = 1, . . . , n ar , q ℓ ar and w ℓ ar,k are defined by

From Eq. (A.3), the following formulas for conditional statistics are derived.

The conditional probability density function p Q|W (q | w 0 ) with respect to dq of Q given W = w 0 in C w is defined as

, q = (qq)/σ Q . (A.7)

n ar ℓ=1 exp(-1 2s 2 ∥ w 0w ℓ ∥ 2 )

, q * = (q *q)/σ Q , (A.8) with

erf( 1 √ 2 s ( q *q ℓ ar )) , erf(y) = 2 √ π y 0 e -t 2 dt . (A.9)