J Tryoen † ‡

O Le Ma Ître

A Ern

ADAPTIVE ANISOTROPIC SPECTRAL STOCHASTIC METHODS FOR UNCERTAIN SCALAR CONSERVATION LAWS *

Keywords: uncertainty quantification, stochastic multiresolution, stochastic spectral method, adaptivity, Galerkin projection, conservation laws, hyperbolic systems AMS subject classifications. 60H35, 60H15, 65C20, 68U20

This paper deals with the design of adaptive anisotropic discretization schemes for conservation laws with stochastic parameters. A Finite Volume scheme is used for the deterministic discretization, while a piecewise polynomial representation is used at the stochastic level. The methodology is designed in the context of intrusive Galerkin projection methods with Roe-type solver. The adaptation aims at selecting the stochastic resolution level based on the local smoothness of the solution in the stochastic domain. In addition, the stochastic features of the solution greatly vary in the space and time so that the constructed stochastic approximation space depends on space and time. The dynamically evolving stochastic discretization uses a tree-structure representation that allows for the efficient implementation of the various operators needed to perform anisotropic multiresolution analysis. Efficiency of the overall adaptive scheme is assessed on the stochastic traffic equation with uncertain initial conditions and velocity leading to expansion waves and shocks that propagate with random velocities. Numerical tests highlight the computational savings achieved as well as the benefit of using anisotropic discretizations in view of dealing with problems involving a larger number of stochastic parameters.

1. Introduction. Stochastic spectral methods and so-called Chaos expansions provide effective tools to propagate parametric uncertainties in numerical models and have been applied successfully to different types of models. The determination of the stochastic spectral solution can be achieved by means of non-intrusive (sampling based) methods or a stochastic Galerkin projection. In this work, we consider the application of the Galerkin projection [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF][START_REF] Le Maître | Spectral Methods for Uncertainty Quantification[END_REF] to the resolution of scalar conservation laws involving uncertain data (such as model parameters, initial and boundary conditions) parametrized by N random variables with known distribution functions and taking values in a stochastic domain Ξ ⊂ R N . Previous applications of Galerkin projection to conservation laws include [START_REF] Lin | Predicting shock dynamics in the presence of uncertainties[END_REF][START_REF] Poette | Uncertainty quantification for systems of conservation laws[END_REF][START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF]. We recall that the Galerkin projection leads to a system of conservation laws governing the stochastic modes of the solution. One essential property of deterministic conservation laws is their capacity to develop non-smooth solutions in finite time, even for a smooth initial condition. This property is obviously present in the stochastic version of the model, with the additional characteristic that the singularity curves are also generally uncertain. We are interested in problems where a shock appears almost surely in finite time. In this case, since the shock speed and/or its location in space can be uncertain, the solution is discontinuous in space and in the stochastic domain. This feature calls for specific discretization techniques. In particular, we rely on Finite Volume schemes for space discretization and on piecewise polynomial discretizations in the stochastic domain [START_REF] Deb | Solution of stochastic partial differential equations using Galerkin finite element techniques[END_REF][START_REF] Le Maître | Multi-resolution analysis of Wiener-type uncertainty propagation schemes[END_REF][START_REF] Wan | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF]. In [START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF], we designed a Roe-type solver for the Galerkin system, and an entropy corrector was developed in [START_REF]Roe Solver with Entropy Corrector for Uncertain Nonlinear Hyperbolic Systems[END_REF].

The method proposed in [START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF][START_REF]Roe Solver with Entropy Corrector for Uncertain Nonlinear Hyperbolic Systems[END_REF], while able to deal with complex situations, remains computationally expensive since a very fine stochastic discretization is needed to properly represent the solution in the neighborhood of discontinuities. This observation calls for adaptive strategies. Since discontinuities are localized in space and evolve in time, we propose in this work stochastic representations depending on space and time, meaning that, at a given time, each Finite Volume cell supports its own stochastic discretization. Consequently, the overall discretization does not rely on a tensorization of stochastic and deterministic approximation spaces, a feature which, to our knowledge, constitutes an original contribution of the present work. The above methodology can be formulated within a multiresolution framework based on the concept of general binary trees to describe the discretization of the stochastic domain, similarly to previous work in the deterministic context [START_REF] Cohen | Fully adaptive multiresolution schemes for conservation laws[END_REF][START_REF] Cohen | Adaptive wavelet techniques in numerical simulation[END_REF]. Restriction and prediction operators are defined over general binary trees in the context of adaptive enrichment and coarsening procedures. For multidimensional stochastic domains, two key features are the use of binary trees and the derivation of (new) directional criteria for the anisotropy of the adaptive procedure, the computational effort being concentrated along the stochastic directions where singularities are observed.

The paper is organized as follows. In §2, we briefly recall the Galerkin projection of stochastic conservation laws and the Roe-type solver proposed in [START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF][START_REF]Roe Solver with Entropy Corrector for Uncertain Nonlinear Hyperbolic Systems[END_REF] in the non-adaptive context. Multiresolution analysis tools are then introduced in §3 to describe the stochastic discretization, while the adaptive Roe solver is presented in §4. Finally, simulations results are presented in §5 for the traffic equation in five stochastic dimensions, and conclusions are drawn in §6.

2. Galerkin Projection and stochastic Roe solver.

2.1. Stochastic scalar conservation laws. We are interested in stochastic nonlinear scalar conservation laws with uncertain input quantities parametrized for simplicity by N independent identically distributed random variables ξ := {ξ 1 . . . ξ N } uniformly distributed in Ξ := [0, 1] N . Let p ξ = 1 denote the density function of ξ ∈ Ξ and let L 2 (Ξ) be the space of second-order random variables defined on the probability space P ξ := (Ξ, B Ξ , p ξ ), where B Ξ is the Borel set of Ξ. The expectation operator in P ξ is denoted, for any random variable H defined on P ξ , by H := Ξ H(y)p ξ (y)dy. Let (x, t, ξ) ∈ Ω × [0, t f ] × Ξ, where Ω is the spatial domain and t f the simulation time. We seek for U (x, t, ξ) solving almost surely the following conservative problem

   ∂ ∂t U (x, t, ξ) + ∂ ∂x F (U (x, t, ξ); ξ) = 0, U (x, t = 0, ξ) = U IC (x, ξ), (2.1)
where F is the stochastic flux and U IC the initial condition. We assume that the problem (2.1) is well-posed and that for all (x, t) ∈ Ω × [0, t f ], U (x, t, •) ∈ L 2 (Ξ).

Stochastic discretization.

To approximate the solution in L 2 (Ξ), we need a stochastic discretization of the problem. This is obtained by considering a Hilbertian basis of random functionals in ξ spanning L 2 (Ξ) equipped with the inner product •, • ,

L 2 (Ξ) = span{Γ 1 (ξ), Γ 2 (ξ), . . . }, Γ α , Γ β = Ξ Γ α (ξ)Γ β (ξ)p ξ (ξ)dξ = δ αβ , (2.2)
where δ αβ is the Kronecker symbol. The discrete solution is sought in a finite dimensional subspace S M of dimension M constructed by truncating the Hilbertian basis:

S M = span{Γ 1 (ξ), Γ 2 (ξ), . . . , Γ M (ξ)} ⊂ L 2 (Ξ). (2.3)
We define the set M := {1 . . . M}. We seek an approximate solution U M ∈ S M in the form

U M (x, t, ξ) = α∈M u α (x, t)Γ α (ξ), (2.4)
where the deterministic fields u α (x, t) are called the stochastic modes of the solution.

2.3. Galerkin system. Plugging U M into (2.1) and requiring the residual to be orthogonal to the subspace S M , we obtain the Galerkin system which couples all the stochastic modes in the form

   ∂ ∂t u(x, t) + ∂ ∂x f (u(x, t)) = 0, u(x, t = 0) = u IC (x), (2.5)
where u(x, t) = (u α (x, t)) α∈M ∈ R M and f (u(x, t)) = (f α (u)) α∈M ∈ R M are the vector of the stochastic modes and the Galerkin flux vector with

f α (u) := F U M ; • , Γ α , (2.6)
and

u IC = U IC , Γ α α∈M . The Galerkin Jacobian matrix ∇ u f of order M is (∇ u f (u)) α,β∈M = ∇ U F (U M ; •), Γ α Γ β α,β∈M . (2.7)
The hyperbolicity of the Galerkin system (2.5), that is, the R-diagonalizability of the Galerkin Jacobian matrix ∇ u f , has been extensively studied in [START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF] (in fact, in the more general context of systems of conservation laws). In particular, for scalar problems as considered in this work, the Galerkin system (2.5) is proven to be hyperbolic.

2.4

. Stochastic Roe solver. The Galerkin system (2.5) is discretized in space and time using a Finite Volume method in the form

u n+1 i = u n i - ∆t n ∆x ϕ(u n i , u n i+1 ) -ϕ(u n i-1 , u n i ) , (2.8)
where u n i is an approximation to the cell-average in the spatial domain of the solution u in the cell of center x i := i∆x with width ∆x at the discrete time t n , ∆t n is the n-th time step, and ϕ(•, •) is the first-order Galerkin numerical flux chosen in the form

ϕ(u L , u R ) = f (u L ) + f (u R ) 2 -|a Roe LR | u R -u L 2 .
(2.9)

Following [START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF], a Roe LR is the Roe-linearized Galerkin Jacobian matrix defined by

a Roe LR := ∇ U F (U Roe LR ; •), Γ α Γ β α,β∈M , (2.10)
where U Roe LR ∈ L 2 (Ξ) is the stochastic Roe state reconstructed from u L and u R . To avoid the expensive spectral decomposition of the Roe Galerkin Jacobian matrix a Roe LR when computing its absolute value, we proposed to approximate |a Roe LR | by a low degree polynomial transformation applied to a Roe LR , constructed using an approximation of its spectrum (see [START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF] for details). Finally, the time-step ∆t n is computed using a CFLtype condition in the form ∆t n = CFL × ∆x/Λ where Λ is the maximum approximate eigenvalue of a Roe LR over all LR interfaces and CFL denotes a user-dependent positive parameter ≤ 1.

As motivated in the introduction, since shock velocities and shock location can be uncertain, the solution is not smooth in the stochastic domain. Consequently, as in [START_REF] Deb | Solution of stochastic partial differential equations using Galerkin finite element techniques[END_REF][START_REF] Le Maître | Multi-resolution analysis of Wiener-type uncertainty propagation schemes[END_REF][START_REF] Wan | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF], we rely on piecewise polynomial approximations for the stochastic discretization. Numerical tests presented in [START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF][START_REF]Roe Solver with Entropy Corrector for Uncertain Nonlinear Hyperbolic Systems[END_REF] on stochastic Burgers and Euler equations using uniform isotropic partitions of the stochastic domain have demonstrated the robustness and accuracy of the above stochastic Roe solver. We also mention that this solver (as any Roe solver) requires a non-linear entropy corrector in presence of sonic points ; such a corrector was designed in [START_REF]Roe Solver with Entropy Corrector for Uncertain Nonlinear Hyperbolic Systems[END_REF].

Stochastic discretization.

3.1. Binary trees. Binary trees provide a convenient representation of nonuniform, anisotropic partitions of the stochastic domain Ξ = [0, 1] N . In a binary tree T, every node has either zero or two children and every node, except the root node denoted by n 0 , has a unique parent. Nodes are collected in the set N (T). Nodes with no children are called leaves and are collected in the set L(T), while nodes with two children are collected in the set N (T) := N (T) \ L(T). The two children of a node n ∈ N (T) are called "left" and "right children" (and also sisters) and are denoted by c -(n) and c + (n). The parent of a node n ∈ N (T) \ {n 0 } is denoted by p(n).

To each node n ∈ N (T), we assign a support S(n) ⊂ Ξ constructed as follows. The supports have the tensor-product form

S(n) = [x - n,1 , x + n,1 ] × • • • × [x - n,N , x + n,N ]
. We set S(n 0 ) = Ξ. The supports of the other nodes are defined recursively by a dyadic partition of the support of the parent node. To this purpose, in the multidimensional case (N > 1), to each node n ∈ N (T), we first assign an indicator, denoted by d(n) ∈ {1 . . . N}, of the direction along which the dyadic partition of its support S(n) is performed. Then, the support of the left and right children are respectively In practice, we consider binary trees T with a fixed maximum number of successive partitions allowed in each direction d ∈ {1 . . . N}. This quantity is called the resolution level and is denoted by Nr. As a result, there holds, for all n ∈ N (T), |S(n)| d ≥ 2 -Nr . A particular case of binary tree is the complete binary tree where |S(l)| d = 2 -Nr for all the leaves and all directions. Thus, there are 2 NNr leaves in a complete binary tree. Such trees, which are associated with uniform isotropic partitions of Ξ, were (implicitly) considered in [START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF][START_REF]Roe Solver with Entropy Corrector for Uncertain Nonlinear Hyperbolic Systems[END_REF] for N = 1 and N = 2. Clearly, for large values of N, the resulting uniform partitions of Ξ are not tractable numerically, so that we rely here on the coarser and anisotropic dyadic partition of Ξ given by (3.1) for a binary tree T that is not complete. An example of one-dimensional complete and incomplete binary trees and their corresponding partitions are shown in the top panel of Figure 3.1, while a two-dimensional incomplete binary tree and its corresponding partition are shown in the bottom panel of Figure 3.1. There is an essential difference between one-dimensional and multidimensional binary trees: for N > 1, there are in general more than one tree with the same set of leaves, i.e., yielding the same partition of Ξ. This is illustrated in Figure 3.2 for N = 2. Consequently, we say that two trees T and T ′ are equivalent if they share the same set of leaves,

S(c -(n)) = [x - n,1 , x + n,1 ] × • • • × [x - n,d , (x - n,d + x + n,d )/2] × • • • × [x - n,N , x + N ] and S(c + (n)) = [x - n,1 , x + n,1 ] × • • • × [(x - n,d + x + n,d )/2, x + n,d ] × • • • × [x - n,N , x + n,
T ≡ T ′ ⇔ L(T) = L(T ′ ). (3.2)
The notion of equivalent trees is needed in the coarsening and enrichment procedures of §4.2.

Stochastic approximation spaces.

Let S(T) be the stochastic approximation space associated with a tree T, spanned by piecewise polynomials on the partition of Ξ given by (3.1):

S(T) = {U T : Ξ → R; ∀l ∈ L(T), U T | S(l) ∈ P N No [ξ]}, (3.3)
where P N No [ξ] is the N-dimensional polynomial space with (partial or total) order less than or equal to No. Let P denote the dimension of the polynomial space P N No [ξ] and set P := {1 . . . P}. The space S(T) has dimension dim(S(T)) = card(L(T)) × P.

(3.4)

Any function U T ∈ S(T) can be written as

U T (ξ) = l∈L(T) α∈P u l α Φ l α (ξ) , (3.5)
where the deterministic coefficients u l α ∈ R are called the Stochastic Element (SE) coefficients of U T in S(T) and the family {Φ l α } l∈L(T),α∈P forms the SE basis. Letting {Φ * α } α∈{1...No+1} be the set of normalized Legendre polynomials on [0, 1] with degree less than or equal to No, let {Φ n0 α } α∈P be the basis of P N No [ξ] obtained by (partial or full) tensorization of the polynomials Φ * α . Then, for all l ∈ L(T) and all α ∈ P, the polynomials Φ l α have for expression

Φ l α (ξ) = |S(l)| -1/2 Φ n0 α (M l (ξ)), ξ ∈ S(l), 0, otherwise. (3.6)
The SE basis is orthonormal, that is,

Φ l α , Φ l ′ β = δ l,l ′ δ α,β , ∀l, l ′ ∈ L(T), ∀α, β ∈ P. (3.7)
The SE expansion (3.5) can be rewritten

U T (ξ) = j∈E(T) u j Φ j (ξ), (3.8)
where E(T) is the multi-index set of the SE expansion whose cardinality is equal to the dimension of S(T) given by (3.4). Alternatively, any function U T ∈ S(T) can be expressed in terms of hierarchical details over the nodes n ∈ N (T) as

U T (ξ) = α∈P u n0 α Φ n0 α (ξ) + n∈ b N (T) α∈P ũn α Ψ n,d(n) α (ξ) . (3.9)
The functions Ψ n,d(n) α

, called multi-wavelets (MW), are supported in S(n). They can be conveniently defined from N sets of P mother functions {Ψ n0,d α } α∈P , d ∈ {1 . . . N}, that are piecewise polynomials on Ξ such that

span α {Φ c - d (n0) α , Φ c + d (n0) α } = span α {Φ n0 α } ⊥ ⊕ span α {Ψ n0,d α }, (3.10)
where c ± d (n 0 ) result from the dyadic partition of n 0 in the direction d, see [START_REF] Le Maître | Multi-resolution analysis of Wiener-type uncertainty propagation schemes[END_REF][START_REF] Tryoen | Adaptive stochastic Galerkin methods for parametric uncertainty propagation in hyperbolic systems[END_REF] for more details (see also [START_REF] Alpert | A class of bases in L 2 for the sparse representation of integral operators[END_REF] in 1D). Then, for all n ∈ N (T) and all α ∈ P, the piecewise polynomials Ψ n,d(n) α have for expression

Ψ n,d(n) α (ξ) = |S(n)| -1/2 Ψ n0,d(n) α (M n (ξ)), ξ ∈ S(n), 0, otherwise, (3.11)
so that, owing to (3.10), for all n ∈ N (T),

span α {Φ c - d (n) α , Φ c + d (n) α } = span α {Φ n α } ⊥ ⊕ span α {Ψ n,d(n) α }. (3.12)
Finally, the MW expansion (3.9) can be rewritten

U T (ξ) = j∈W(T) ũj Ψ j (ξ), (3.13)
where W(T) is the multi-index set of the MW expansion.

Restriction and prediction operators.

Restrictions and predictions operators are useful tools in the adaptive context. Let T 1 and T 2 be two binary trees. We say that

T 1 ⊂ T 2 if ∀l 2 ∈ L(T 2 ), ∃!l 1 ∈ L(T 1 ) s.t. S(l 1 ) ⊂ S(l 2 ). (3.14)
Clearly, if T 1 ⊂ T 2 , then S(T 1 ) ⊂ S(T 2 ).

Restriction operator.

Let T 1 and T 2 be two binary trees such that T 1 ⊂ T 2 . Given U T2 ∈ S(T 2 ), we define the restriction of U T2 to S(T 1 ), denoted R ↓T1 U T2 , as the orthogonal L 2 (Ξ)-projection of U T2 onto S(T 1 ), i.e., (U T2 -R ↓T1 U T2 ) ⊥ S(T 1 ). In terms of MW coefficients, the restriction operation is straightforward. Letting ũn α be the MW coefficients of U T2 and using the orthonormality of the MW basis yields, for all n ∈ N (T 1 ) and all α ∈ P,

R ↓T1 U T2 n α = ũn α . (3.15)
The computation of the SE coefficients of the restriction is not as immediate. Assuming that the SE expansion of U T2 is known, we construct a sequence of trees

T (i) such that T 2 = T (0) ⊃ • • • ⊃ T (i) ⊃ • • • ⊃ T (l) = T 1
, where two consecutive trees differs from one generation only, i.e., a leaf of T (i+1) is either a leaf or a node with leaf children in T (i) . Therefore, the transition from T (i) to T (i+1) consists in removing pairs of sister leaves. The process is illustrated in the left part of Figure 3.3 for the removal of a single pair of sister leaves. Focusing on the removal of a (left-right ordered) pair of sister leaves {l -, l + }, the SE coefficients of the restriction of U T (i) associated with the

new leaf l = p(l -) = p(l + ) ∈ L(T (i+1) ) in direction d(l) are u l α = β∈P R -,d(l) α,β u l - β + R +,d(l) α,β u l + β , (3.16)
where, for all d ∈ {1 . . . N}, the transition matrices R ±,d of order P have entries given by R

±,d α,β = Φ n0 α , Φ c ± d (n0) β .
Restriction Prediction

T (i) T (i+1) l l + - + - l = p(l ) = p(l ) T (i+1) T (i) c (n) - n n c (n) + Fig. 3.3.
Schematic representation of the elementary restriction (left) and prediction (right) operators through the removal and creation respectively of the (leaves) children of a node.

Prediction operator.

Let T 1 and T 2 be two binary trees such that T 1 ⊂ T 2 . The prediction operation consists in extending U T1 ∈ S(T 1 ) to the larger stochastic space S(T 2 ). We denote by P ↑T2 U T1 this prediction. Different predictions can be used (see [START_REF] Cohen | Fully adaptive multiresolution schemes for conservation laws[END_REF][START_REF] Cohen | Adaptive wavelet techniques in numerical simulation[END_REF]); here we have considered the simplest one, where no information is generated by the prediction. As for the restriction operation, the MW expansion of the prediction is immediately obtained from the MW coefficients of U T1 . We obtain, for all n ∈ N (T 2 ) and all α ∈ P,

P ↑T2 U T1 n α = ũn α , n ∈ N (T 2 ), 0, otherwise. (3.17)
For the SE coefficients of the prediction, we can again proceed iteratively using a series of increasing intermediate trees, differing by only one generation from one to the other. This time, the elementary operation consists in adding children to some leaf of the current tree in a chosen direction d, as illustrated in right part of Figure 3.3.

The SE coefficients associated to the new leaves of a node n are given by

u c -(n) α = β∈P R -,d α,β u n β , u c + (n) α = β∈P R +,d α,β u n β , (3.18)
with the same transition coefficients as those used in (3.16). For two trees T 1 ⊂ T 2 , we observe that R ↓T1 • P ↑T2 = I T1 , while in general P ↑T2 • R ↓T1 = I T2 (I denoting the identity).

4. Adaptive stochastic Roe solver. Singularities in the solution of stochastic conservation laws remain localized both in the deterministic domain Ω × [0, t f ] for each ξ ∈ Ξ, and in the stochastic domain Ξ for each (x, t) ∈ Ω × [0, t f ]. In other words, the solution is almost everywhere smooth on Ω × [0, t f ] × Ξ. This observation strongly advocates for the use of adaptive strategies where the computational effort is concentrated along the singularity curves, while coarser discretizations are used where the solution is smooth. In what follows, we only consider adaptation of the stochastic discretization by relying on a fixed spatial mesh (the time-step being adapted to satisfy a global CFL condition). However, the adapted stochastic discretization depends on the spatial variable x and the time t. In the context of the stochastic discretization framework introduced in the previous section, it amounts to an indexation with both x and t of the trees T defining the stochastic approximation space S(T). Specifically, we now denote T n i the tree associated with the i-th cell of the spatial mesh at the discrete time t n , such that the approximate stochastic solution on the i-th cell at time t n has for expansions

U n i (ξ) = j∈E(T n i ) (u n i ) j Φ j (ξ) = j∈W(T n i ) ( u n i ) j Ψ j (ξ) ∈ S(T n i ). (4.1)
4.1. Algorithm. In this section, we present an overview of the adaptive stochastic Roe-type solver used for the time-integration of the stochastic conservation law. The structure of the solver is outlined in Algorithm 4.1. The algorithm starts with the definition of the initial data consisting, for each spatial cell, of a tree T 0 i and the approximation of the cell-averaged initial condition yielding U 0 i ∈ S(T 0 i ). The accuracy parameter η (to be used in §4.2) and the resolution level Nr are also required before proceeding with the time-iterations that constitute the core of the algorithm. A timeiteration consists in four main steps: an enrichment of the current stochastic approximation spaces, the computation of the fluxes at the interfaces, the time-advancement of the solution, and finally the coarsening of the underlying trees. We remark that in Algorithm 4.1, the time integration and coarsening steps have been distinguished for clarity, but in fact the coarsening can be applied to each spatial cell immediately after the its time-advancement, resulting in a more efficient implementation. We briefly outline the role of the different procedures appearing in Algorithm 4.1. Construct T 0 i and U 0 i 4. end for 5. n ← 0, t ← 0 6. while t < t f do 7.

Step I: enrichment

8.
for all cells i of the spatial mesh do 9.

(

T n i , U n i ) ← Enrich (T n i , U n i , η, Nr)
[Enrich the approximation spaces] 10.

end for 11.

Step II: computation of the fluxes

12.
for all interfaces i -1/2 of the spatial mesh do 13.

T i-1/2 ← Union T n i-1 , T n i [Construct flux space] 14. U L ← Predict(T i-1/2 , U n i-1 ), U R ← Predict(T i-1/2 , U n i ) [Predict states] 15. Φ i-1/2 ← 0 [Initialize flux] 16.
for all l ∈ L(T i-1/2 ) do 17.

Φ i-1/2 ← Φ i-1/2 + Roeflux(l, U L , U R ) [Add leaf's contribution to the flux] 18.
end for

19.
end for 20.

Select the time-step ∆t n 21.

Step III: time integration 22.

for all cells i of the spatial mesh do 23.

T * i ← Union T i-1/2 , T i+1/2 [Construct integration space] 24. δΦ ← Predict(T * i , Φ i-1/2 ) -Predict(T * i , Φ i+1/2 ) [Compute flux difference] 25. U n+1 i ← Predict(T * i , U n i ) + ∆t n ∆x δΦ [Integrate in time] 26.
end for 27.

Step IV:

Coarsening 28.
for all cells i of the spatial mesh do 29.

(T n+1 i , U n+1 i ) ← Coarsen(T * i , U n+1 i , η, Nr) [Threshold the solution] 30.
end for

31. t ← t + ∆t n , n ← n + 1 [Increment time]
32. end while

Step I. The purpose of this step is to enrich the stochastic approximation spaces in order to anticipate the emergence during the time-step of additional local stochastic features in the solution requiring more resolution. The procedure Enrich extends the current tree by refining some of its leaves using one of the two enrichment strategies described in §4.2.

Step II. In this step, the numerical fluxes at all cell interfaces are evaluated. The first procedure, Union(T 1 , T 2 ), constructs the minimal tree encompassing both T 1 and T 2 . This union is needed because in general T n i-1 = T n i , for two neighboring cells i-1 and i sharing an interface. Since the solutions over the two cells are not defined with respect to the same stochastic basis, we first construct a common stochastic approximation space defined by the union of the two cells trees. The formal definition of the union of trees is as follows. Given two generic trees T 1 and T 2 , we define their union-tree T 1∪2 := T 1 ∪ T 2 as (one of) the minimal tree(s) (in terms of number of leaves) such that, for all l ∈ L(T 1∪2 ),

∃!l 1 ∈ L(T 1 ), ∃!l 2 ∈ L(T 2 ), S(l) = S(l 1 ) ∩ S(l 2 ). (4.2)
The union of two trees is illustrated in Figure 4.1 for N = 2. The union-tree is not unique whenever N > 1, since different minimal trees can be constructed to satisfy (4.2), but these union-trees are equivalent and yield the same stochastic space S(T 1∪2 ). Moreover, S(T 1 ) and S(T 2 ) are subspaces of S(T 1∪2 ) since T 1 ⊂ T 1∪2 and T 2 ⊂ T 1∪2 . Thus, we can compute the SE coefficients in S(T 1∪2 ) of the stochastic quantities U T1 and U T2 by means of the prediction operators introduced previously. The procedure Predict in Algorithm 4.1 precisely implements the recursive prediction introduced in §3.3. The solutions of the left and right cells being now defined in S(T i-1/2 ), the numerical flux Φ i-1/2 at the interface can be computed. Following [START_REF] Tryoen | Intrusive Projection Methods with Upwinding for Uncertain Nonlinear Hyperbolic Systems[END_REF], the flux is computed in terms of its SE expansion coefficients, and we exploit the decoupling of the flux evaluation for distinct leaves to perform the calculation in parallel. The procedure Roeflux implements the computation of the local numerical flux at the level of a leaf l. We have briefly described in §2.4 the computation of the expansion coefficients of the numerical flux in a generic stochastic basis, given the expansion coefficients of the left and right states (u L and u R ) in this basis. We can readily reuse this method with the basis functions {Φ n0 α } α∈P , relying on the affine maps M l : S(l) → Ξ and a scaling of the SE coefficients. In addition, although this is not made explicit in Algorithm 4.1, Roeflux also returns the maximal velocity over the leaf, to be used for the CFL-based selection of the time-step ∆t n .

Step III. In this step, we advance in time the solution of each cell by integration of the difference in the Galerkin fluxes of its interfaces according to (2.8). The procedures Union and Predict are used to define all the relevant quantities in a common stochastic space.

Step IV. The purpose of this step is to control, in each cell of the stochastic mesh, the growth of its tree resulting from the enrichment and union operations, by performing a coarsening of the tree. The procedure Coarsen, applied to each cell (possibly in parallel), performs this coarsening by means of thresholding, where nonsignificant details in the solution are removed to achieve minimal resolution for a target accuracy fixed by η. Details on the procedure Coarsen are provided in §4.2.

Adaptivity.

In this section, we detail the essential adaptivity tools needed for the control of the local stochastic resolution, with the objective to efficiently reduce the complexity of the computations. There are two essential tools, the coarsening and enrichment procedures. Thresholding error. Let η > 0 be the tolerance fixed in Algorithm 4.1 and recall that Nr denotes the resolution level. Let ũn α denote the MW expansion coefficients of U T , see (3.9). We define D(η, Nr) as the subset of N (T) such that

D(η, Nr) := n ∈ N (T); ũn ℓ 2 ≤ 2 -|n|/2 (NNr) -1/2 η , (4.3)
where ũn := (ũ n α ) α∈P and ũn 2 ℓ 2 = α∈P (ũ n α )

2 . The motivation for (4.3) is that, letting U T be the thresholded approximation of U T obtained by omitting in the second sum of (3.9) the nodes n ∈ D(η, Nr), there holds

Û T -U T 2 L 2 (Ξ) = n∈D(η,Nr) ũn 2 ℓ 2 ≤ n∈D(η,Nr) 2 -|n| (NNr) -1 η 2 ≤ η 2 , (4.4
)

since n∈D(η,Nr) 2 -|n| = NNr-1 j=0 #{n ∈ D(η, Nr); |n| = j}2 -j ≤ NNr-1 j=0
1 = NNr. Coarsening procedure. Two points deserve particular attention. The first one is that N (T) \ D(η, Nr) has not a binary tree structure in general, so that a procedure is needed to maintain this structure when removing nodes of T. Here, we choose a conservative approach where the resulting subtree T -may still contain some nodes in the set D(η, Nr). Specifically, we construct a sequence of imbricated trees, obtained through the removal of pairs of sister leaves from one tree to the next: a couple of sister leaves having node n for parent is removed if n ∈ D(η, Nr). The coarsening sequence is stopped whenever no couple of sister leaves can be removed, and this yields the desired subtree T -. The second point is that the above algorithm only generates trees such that, along the sequence, the successive (coarser and coarser) partitions of Ξ follow, in backward order, the partition directions d(n) prescribed by T. This is unsatisfying because for N > 1, there are many trees equivalent to T, and we would like the coarsened tree to be independent of any particular choice in this equivalence class. To avoid arbitrariness, the trees of the sequence are periodically substituted by equivalent ones, generated by searching in the current tree the pattern of a node n whose children c -(n) and c + (n) are not leaves and are subsequently partitioned along the same direction d(c

+ (n)) = d(c -(n)) which differs from d(n); when such a pattern is found, partition directions are exchanged, d(n) ↔ d(c -(n)) = d(c + (n))
, together with the corresponding permutation of the descendants of the children nodes. This operation, illustrated in Figure 4.2, is applied periodically and randomly along the coarsening procedure.

c (n) + c (n) + c (n) - c (n) - Tc Td Ta Tc Tb Td n n
Tb Ta Fig. 4.2. Illustration of the elementary operation to generate equivalent trees: the pattern of a node with its children divided along the same direction (left) is replaced by the same pattern but with an exchange of the partition directions (right) plus the corresponding permutation of the descendants of the children.

Anisotropic enrichment.

Let T be a binary tree and let U T ∈ S(T). The purpose of the enrichment is to increase the dimension of S(T), by adding descendents to some of its leaves. Enrichment of the stochastic space is made necessary by the possible emergence in time of new features in the stochastic solution, such as shocks, that require more resolution. The underlying assumption is that the time-step is small enough so that a tree containing significant details at both times t n and t n+1 can be constructed from the details at t n , see [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF]. In what follows, the enrichment is limited to one partition along each dimension at most.

The simplest enrichment procedure consists in systematically partitioning all the leaves l ∈ L(T) once for all d ∈ {1 . . . N} provided |S(l)| d > 2 -Nr . This procedure generates a tree T + that typically has 2 N card(L(T)) leaves, which is only practical when N is small. More economical strategies are based on the analysis of the MW coefficients in U T to decide which leaves of T need be partitioned and along which direction (see for instance [START_REF] Le Maître | Multi-resolution analysis of Wiener-type uncertainty propagation schemes[END_REF][START_REF] Le Maître | Multi-resolution-analysis scheme for uncertainty quantification in chemical systems[END_REF]). We derive below two new directional enrichment criteria in the context of N-dimensional binary trees. MultiD enrichment criterion. Classically, the theoretical decay rate of the MW coefficients with resolution level is used to decide the partition of a leaf from the norm of MW coefficients of its parent (see for instance [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF][START_REF] Cohen | Fully adaptive multiresolution schemes for conservation laws[END_REF] in the deterministic case).

We first recall some background in the 1D case (N = 1). Let U ∈ L 2 (Ξ) with Ξ = [0, 1]. Let T 1D be a 1D binary tree and let U T1D be the L 2 (Ξ)-orthogonal projection of U onto S(T 1D ). Let ũn α denote the MW coefficients of U T1D . Then, if U is locally smooth enough, the magnitude of the MW coefficients ũn α of a generic node n ∈ N (T 1D ) can be bounded as

|ũ n α | = inf P ∈PNo[ξ] | (U -P ), Ψ n α | ≤ C|S(n)| No+1 U H No+1 (S(n)) , (4.5)
where H No+1 (S(n)) is the usual Sobolev space of order (No + 1) on S(n). Recalling that |S(n)| = 2 -|n| , the bound (4.5) shows that the norm of the MW coefficients decays roughly as O(2 -|n|(No+1) ) for smooth U . Therefore, the norm of the (unknown) MW coefficients of a leaf l ∈ L(T 1D ) can be estimated from the norm of the (known) MW coefficients of its parent as ũl ℓ 2 ∼ 2 -(No+1) ũp(l) ℓ 2 . This estimate can, in turn, be used to derive an enrichment criterion; specifically, a leaf l is partitioned if the estimation for ũl ℓ 2 does not match the thresholding criterion (4.3), that is, if

ũp(l) ℓ 2 ≥ 2 No+1 2 -|l|/2 Nr -1/2 η and |S(l)| > 2 -Nr . (4.6)
The extension to N > 1 of the enrichment criterion (4.6) is not straightforward in the context of binary trees. Indeed, the MW coefficients associated with a node n carry an information essentially related to the splitting direction d(n). Thus, for a leaf l ∈ L(T), they cannot be used for an enrichment criterion in a direction d = d(p(l)). To address this issue, we define, for any leaf l ∈ L(T) and any direction d ∈ {1 . . . N}, its virtual parent p d (l) as the (virtual) node that would have l as a child after a dyadic partition along the d-th direction. Consistently, s d (l) denotes the virtual sister of l along direction d. Note that p d (l) ∈ N (T) only for d = d(p(l)); moreover, in general s d (l) / ∈ N (T). These definitions are illustrated in Figure 4.3 which shows for N = 2 the partition associated with a tree T (left plot), and the virtual sisters of two leaves. The SE coefficients of the virtual sisters,

u s d (l) α := U T , Φ s d (l) α , α ∈ P, (4.7)
are efficiently computed by exploiting the binary structure of T and relying on the elementary restriction and prediction operators defined in §3.3. Without going into too many details, let us mention that the computation of the SE coefficients of s d (l) Illustration of the virtual sisters of a leaf l of a tree T whose partition is shown in the left plot. In the center plot, the leaf l is hatched diagonally in blue and its two virtual sisters for d = 1 and 2 (hatched horizontally and vertically respectively) are leaves of T, both being c + (p d (l)). In the right plot, a different leaf l is considered (still hatched diagonally in blue) with virtual sisters s d (l) which for d = 1 (hatched horizontally) is a node of T but not a leaf, and which for d = 2 (hatched vertically) is not a node of T.

✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✁ ✂ ✁✂
amounts to i) finding the subset of leaves in L(T) whose supports overlap with S(s d (l)), ii) constructing the subtree having for leaves this subset, and iii) restricting the solution over this subtree up to s d (see Appendix B). In practice, one can reuse the restriction operator defined in §3.3 to compute the usual details in the {Ψ n,d α } α∈P basis for a chosen direction d.

We now return to the design of a multiD enrichment criterion. A natural extension of (4.6) is that a leaf l is partitioned in the direction d if

ũp d (l) ℓ 2 ≥ diam(S(p d (l))) diam(S(l)) No+1 2 -|l|/2 (NNr) -1/2 η and |S(l)| d > 2 -Nr . (4.8)
This criterion is motivated by the following multiD extension of the bound (4.5) for the magnitude of the MW coefficients ũn α in the direction d for a generic node n,

|ũ n α | = inf P ∈P N No [ξ] (U -P ), Ψ n,d α ≤ Cdiam(S(n)) No+1 U H No+1 (S(n)) . (4.9)
Directional enrichment criterion. We want to improve the criterion (4.8) since the isotropic factor diam(S(p d (l)))/diam(S(l)) can take the value 1 in the context of anisotropic refinement. Instead, we would like to devise a criterion with the factor 2 No+1 , since this will lead to smaller enriched trees. To this purpose, we derive an alternative criterion that is fully directional. For any direction d ∈ {1 . . . N} and any node n ∈ T, we define the directional detail coefficients ūn,d β∈{1...No+1} through ūn,d

β := U, Ψn,d β , Ψn,d β (ξ) =    |S(n)| -1/2 Ψ * β ξ d -x - n,d x + n,d -x - n,d , ξ ∈ S(n),
|ū n,d β | = inf P (ξ)∈PNo[ξ d ] (U -P ), Ψn,d β ≤ C|S(n)| No+1 d U L 2 (S ∼d (n),H No+1 (S d (n))) .
(4.11)

Proceeding as previously, the enrichment criterion states that a leaf l is partitioned along direction d if

ūp d (l) ℓ 2 ≥ 2 No+1 2 -|l|/2 (NNr) -1/2 η and |S(l)| d > 2 -Nr . (4.
12)

The details norm associated with the basis { Ψn,d β } β∈P can be obtained directly from the vector of MW coefficients ũn,d by averaging it in all but the d-th direction.

5.

Results. The effectiveness of the proposed method is assessed on the so-called traffic equation, where the normalized density U of vehicles is governed by (2.1) with

F (U ; ξ) = A(ξ)U (ξ)(1 -U (ξ))
, where A(ξ) is almost surely positive and represents the reference velocity.

Problem definition.

We consider a spatial domain Ω = [0, 1] with periodic boundary conditions. The considered uncertainties are on the initial (1-periodic) data U IC (x, ξ) and on the velocity A(ξ). Specifically, the uncertain initial condition consists of four piecewise constant uncertain states in x, parametrized using four independent random variables ξ 1 , ξ 2 , ξ 3 , and ξ 4 , with uniform distributions in [0, 1]:

U IC (x, ξ) = U (ξ 1 ) -U -(ξ 2 )I [0.1,0.3] + U + (ξ 3 )I [0.3,0.5] -U -(ξ 4 )I [0.5,0.7] , (5.1)
where

U (ξ 1 ) = 0.25 + 0.01ξ 1 ∼ U[0.25, 0.26], U -(ξ d ) = 0.2 + 0.015ξ d ∼ U[0.2, 0.215],
d ∈ {2, 4}, and U + (ξ 3 ) = 0.1 + 0.015ξ 3 ∼ U[0.1, 0.115] (here, I Z denotes the characteristic function of the set Z). In addition, the uncertain velocity is parametrized by one random variable ξ 5 with uniform distribution in [0, 1]:

A(ξ 5 ) = 1 + 0.05(2ξ 5 -1) ∼ U[0.95, 1.05]. (5.2)
This problem has therefore five stochastic dimensions (N = 5). The space P N No [ξ] is spanned by the partially tensorized Legendre polynomials with degree ≤ No, so that P = (N+No)! N!No! . The resolution level is set to Nr = 6. 5.2. Numerical results. We fix η = 10 -4 , No = 3, and Nc = 200 uniform cells for the spatial discretization. The initial condition is illustrated in the left panel of Figure 5.1, where 20 realizations of U IC (ξ) are shown. In the middle and right panels of Figure 5.1, we show 20 realizations of the solution at times t = 0.4 and 0.9, reconstructed from the adaptive expansion and corresponding to the realizations of the initial data from the left panel. We observe the generation of two expansion waves from x = 0.1 and x = 0.5, and of two shock waves from x = 0.3 and x = 0.7. As time evolves, the first expansion wave reaches the first shock, while the second expansion waves reaches the second shock. Because of uncertainties in the wave velocities, the instants where the waves catch up are uncertain. The dynamics and the impact of uncertainties can be better appreciated on the space-time diagram of the solution expectation and standard deviation plotted for t ∈ [0, 1] in the left and right panels of Figure 5.2. We observe the smooth nature of the solution expectation and the steep variations in the solution standard deviation, with maxima reached along the paths of the two shocks.

The respective contribution of each input parameter on the variance of the solution is represented by the first-order sensitivity indices S d in Figure 5.3 (see Appendix B for the definition and computation of the sensitivity indices). Before the merging of the expansion and shock waves (t < 0.4), significant values are observed for S 1 , S 2 , S 3 , and S 4 over portions of the computational domain corresponding to the three dependence cones between the waves, where the solution takes one of the three initial uncertain states. The portions of the spatial domain where S 1-4 take significant values reduces as time increases, indicating the emergence of more and more interactions between the random parameters. On the contrary, because ξ 5 parametrizes the uncertain velocity A, the significant values of S 5 appears along paths of the different waves and affect a portion of the spatial domain that increases with time. The emergence of interactions between parameters can be appreciated from the most right panel of Figure 5.3, where the quantity 1 -N d=1 S d , i.e., the fraction of the variance due to higher-order sensitivity indices, is plotted. This figure shows that interactions primarily take place along the shock paths. We also present the total sensitivity indices T d which measure the total sensitivity of the solution with respect to the parameter ξ d . These total sensitivity indices are displayed in Figure 5.4 as a function of x at the same times as in Figure 5.1. We recall that T d ≤ 1, while d T d > 1 in general. We observe that T 2 and T 3 (resp. T 4 ) take significant values over supports that are compact in the neighborhood of the first (resp. second) shock wave, and that their magnitude tends to decay in time. On the contrary, the portion of the spatial domain where T 5 reaches a value close to 1 becomes larger as time increases, indicating the extension of the domain of influence of the uncertainty in A. For instance, for t = 0.9, the set {T 5 ≈ 0} is included in x ∈ [0.4, 0.5], that is, to the only remaining part of the domain where the stochastic solution is spatially constant (see the right plot of Figure 5.1). Finally, the dynamics of T 1 , which is related to an uncertainty in the initial data that is non-local, is much more complex. Specifically, T 1 continues to be significant in areas where the stochastic solution is piecewise constant in space and along the shocks, while in rarefaction waves T 1 becomes quickly insignificant.

Refinement analysis.

We have tested the two enrichment criteria (multiD (4.8) and directional (4.12)) using different values for η and No. The results of these experiments indicate that, at fixed η and No, the multiD criterion leads to more refined stochastic discretizations. However, the finer stochastic discretizations resulting from the multiD criterion only achieve a marginal reduction of the approximation error (as measured by the stochastic approximation error ǫ sd defined by (5.3) below) compared to the directional criterion. This is illustrated in Figure 5.5 where we report the evolution in time of the total number of SE for the two enrichment criteria, different values of η, and No = 3. The right most plot shows the corresponding of ǫ sd as a function to the total number of SE at t = 0.5. Because the two enrichment criteria have similar computational complexity for their evaluation, the results presented in what follows all use the directional criterion (4.12). The left plot in Figure 5.6 displays the space-time diagram of the averaged depths of the trees measured as log 2 card(L(T n i )) for each (x i , t n ) and with η = 10 -4 and No = 3. This plot shows the adaptation of the stochastic resolution in space and time to the local stochastic smoothness; as expected, a finer stochastic discretization along the path of the shock waves is necessary, while a coarser discretization suffices in the expansion waves and in the regions where the solution is spatially constant. The right plot in Figure 5.6 shows the time evolution of the total number of leafs in the stochastic discretization. We observe a monotonic increase in the number of leafs, with higher rates when additional wave interactions occur and, subsequently, with a roughly constant rate since the stochastic shocks, which dom- Because ξ 1 parametrizes the uncertain initial condition on the whole domain, this variable affects the velocity of the two shock waves, so that the discretization is finer in the neighborhood of the two shocks. Then, ξ 2 and ξ 3 (resp. ξ 4 ) affect the velocity of the first shock wave (resp. the second), so that the discretization is finer in the neighborhood of the first (resp. the second) shock. Finally, ξ 5 , which parametrizes the velocity A and therefore affects the velocity of the two shocks, is observed to be the most influent parameter, so that the trees are deeper in the fifth direction ; this fact explains the high values of the aspect ratio near the shocks. 5.4. Convergence and computational time analysis. The convergence of the adaptive stochastic method is numerically investigated in a first series of experiments. We fix the number of spatial cells to Nc = 200, and compute the solution at t n = 0.5 for different values of η and No. We characterize the approximation error in the semi-discrete solution by the following measure

ǫ 2 sd = ∆x Nc i=1 Ξ U n i (ξ) -U n ex,i (ξ) 2 dξ, (5.3)
where U n ex,i denotes the exact stochastic semi-discrete solution. This error measure is approximated by means of a Monte Carlo simulation, consisting in a) sampling uniformly ξ ∈ Ξ, b) solving the corresponding discrete deterministic problems with a deterministic Roe solver, c) computing the difference with the computed adapted solution, and d) averaging over samples to get an empirical estimate of ǫ sd . In practice, 10,000 MC samples suffice to obtain a well converged error measure. In all these experiments, a fixed time-step ∆t = 1/200 is used. Figure 5.8 shows the decay of ǫ 2 sd when the tolerance η in the adaptive algorithm is decreased. The different curves correspond to polynomial orders No ∈ {2 . . . 5}. The left plot depicts the error measure as a function of the total number of elements (leaves) in the adaptive stochastic discretization at t n = 0.5, namely the sum over all cells i of card(L(T n i )). The convergence of the semi-discrete solution as η is lowered is first observed for all polynomial orders tested. In addition, the higher No, the lower the error and the faster the convergence rate, owing to richer approximation spaces for equivalent number of stochastic elements. However, if the error measure ǫ sd is plotted as a function of the total number of degrees of freedom in the stochastic approximation spaces, i.e., the total number of leaves times P, as shown in the right plot of Figure 5.8, we observe that for low resolution (largest η), low polynomial orders are more efficient than larger ones. On the contrary, for highly resolved computations (lowest values of η), high polynomial orders achieve a more accurate approximation for a lower number of degrees of freedom. This behavior is explained by the diffusivity of the Roe solver that (slightly) smoothes the solution at the stochastic level, even if steep dependences with respect to ξ are observed in particular along the shocks paths.

To complete the analysis of the adaptive method, we briefly discuss its computational efficiency. Our main purpose is to demonstrate that the overhead due to the adapted stochastic discretization in space and time is limited. We first observe that our Roe scheme being explicit in time, parallelization by spatial domain decomposition is immediate. However, a dynamical partition of the spatial domain would be necessary to balance the computational loads as the stochastic discretization evolves in time. A second level of parallelism concerns the flux computations using the procedure Roeflux, which can be applied to different sets of leaves in parallel. Therefore, we are mostly concerned with the efficiency of the procedures Union, Predict, Enrich, and Coarsen, in particular the scaling of their computational times with the size of the stochastic discretization. The two plots in Figure 5.9 report the CPU times (in arbitrary units) for the advancement of the solution over a time step using the discretization parameters No = 2, η = 10 -3 , and No = 3, η = 10 -4 respectively. The CPU times are given as a function of the total number of leaves involved in the flux evaluation (number of calls to Roeflux). These numerical experiments show that, owing to the representation of the stochastic approximation spaces using binary tree structures, an asymptotically linear computational time in the number of leaves is achieved. The contributions of the different parts of the algorithms are also detailed. For the two discretization parameters, the most time consuming part of the algorithm is the flux evaluation which significantly dominates the computational times for enrichment and coarsening.

Conclusion.

We have proposed an adaptive anisotropic strategy in the context of multiresolution analysis for uncertain conservation laws with a locally refined stochastic approximation space depending on space and time. The binary tree structure used to represent the stochastic discretization permits efficient implementation of the operators needed to deal with anisotropic adaptation, in particular the enrichment and coarsening procedures, but also for the post-processing of the numerical solution to determine complex information such as sensitivity indices. Two anisotropic criteria have been derived to decide for the enrichment along the different stochastic dimensions of the problem. The present results illustrate the ability of the method to deal with multidimensional stochastic nonlinear scalar conservation laws including shocks with significant computational savings owing to the adapted anisotropic discretization. Future work can aim at further improvements of the present adaptive strategy, in particular, by considering adaptive spatial discretization or higher-order numerical fluxes. If both S d and T d are small, ξ d is deemed non-influent. If S d is small and T d is not small, then ξ d is influent only through its interaction with other random parameters. Consider now a binary tree T and a functional U T ∈ S(T). The tree data structure and local projection restriction operators can be exploited to efficiently compute the marginals of U T . For instance, to compute U T \ξ v for some set v ⊆ {1 . . . N}, we construct a N-dimensional tree T v which is sufficiently rich to exactly represent U T \ξ v . Because T v T in general, we first assemble the set of leaves centers {x l , l ∈ L(T)} such that, for all l ∈ L(T), Finally, the projection of U T in S(T v ) is marginalized, locally over each leaf l ∈ L(T v ), along every direction d ∈ ∼v. In fact, for any v ′ ⊂ v, the procedure can be applied recursively to construct T v ′ from T v and to compute U T \ξ v ′ from U T \ξ v .

  N ] where d = d(n). This construction leads to a partition of the stochastic domain Ξ in the form Ξ = l∈L(T) S(l). (3.1) For a node n ∈ N (T), its depth |n| is defined as the number of generations it takes to reach n from the root node n 0 . It is readily seen that the support of node n has measure |S(n)| := 2 -|n| . We also define the measure of S(n) in direction d as |S(n)| d := x + n,dx - n,d , its diameter as diam(S(n)) := max d |S(n)| d , and its volume in all directions except d as |S(n)| ∼d := |S(n)|/|S(n)| d . Finally, for any node n ∈ N (T), M n denotes the affine map from S(n) onto the reference stochastic domain Ξ.

Fig. 3 . 1 .

 31 Fig. 3.1. Top : Complete binary tree (left) and incomplete binary tree (right) for N = 1; the corresponding partitions of Ξ = [0, 1] are shown below the trees. Bottom : Multidimensional binary tree for N = 2 (left). Dash (resp. full) segments represent a partition along the first (resp. second) direction. Corresponding partition of Ξ = [0, 1] 2 (right).

Fig. 3 . 2 .

 32 Fig. 3.2. Example of two equivalent trees for N = 2. The solid (resp. dash) segments represent a partition along the first (resp. second) direction. The partition of Ξ is shown at the center.

Algorithm 1

 1 Structure of the adaptive stochastic Roe-type solver 1. Set η > 0 and Nr > 0 [Set precision and resolution level] 2. for all cells i of the spatial mesh do 3.

Fig. 4 . 1 .

 41 Fig. 4.1. Example of union of two trees for N = 2: T 1 (left), T 2 (center) and T 1∪2 = T 1 ∪ T 2 (right). For each tree the corresponding partition is also shown.

4. 2 . 1 .

 21 Coarsening. Let T be a binary tree and let U T ∈ S(T). The coarsening procedure aims at constructing a subtree T -⊂ T (or, equivalently, a stochastic approximation subspace S(T -) ⊂ S(T)) through a thresholding of the MW expansion coefficients of U T .

Fig. 4 .

 4 Fig.4.3. Illustration of the virtual sisters of a leaf l of a tree T whose partition is shown in the left plot. In the center plot, the leaf l is hatched diagonally in blue and its two virtual sisters for d = 1 and 2 (hatched horizontally and vertically respectively) are leaves of T, both being c + (p d (l)). In the right plot, a different leaf l is considered (still hatched diagonally in blue) with virtual sisters s d (l) which for d = 1 (hatched horizontally) is a node of T but not a leaf, and which for d = 2 (hatched vertically) is not a node of T.

  β } β∈{1...No+1} is the set of 1D wavelet functions defined on [0, 1]. The vector of coefficients ūn,d measures details in U at the scale |S(n)| d , in direction d only, by averaging-out any variability in U along the other directions. The estimate for the directional details magnitude is now (see Appendix A)

Fig. 5 . 1 .

 51 Fig. 5.1. Stochastic traffic equation: sample set of 20 realizations of the initial condition (left) and computed solution at t = 0.4 (middle) and t = 0.9 (right).

Fig. 5 . 2 .

 52 Fig. 5.2. Space-time diagrams of the solution expectation (left) and standard deviation (right).

Fig. 5 . 3 .

 53 Fig. 5.3. Space-time diagrams (x, t) ∈ [0, 1] × [0, 1] of the first-order sensitivity indices and the contribution of sensitivity indices of higher order.

5 Fig. 5 . 4 .

 554 Fig. 5.4. Total sensitivity indices as a function of x ∈ [0, 1] at t = 0.4 (left) and t = 0.9 (right).

Fig. 5 . 5 .

 55 Fig. 5.5. Comparison of the two enrichment criteria for No = 3 and different values of η as indicated. Evolution in time of the total number of stochastic elements in the discretization for the multiD criterion (4.8) (left plot) and the directional criterion (4.12) (center plot). Right plot: corresponding error measures ǫ sd at t = 0.5 as a function of the total number of SE for the two enrichment criteria.

Fig. 5 . 6 .

 56 Fig. 5.6. Space-time diagrams of the averaged depth of local trees in log 2 scale (left) and evolution in time of the total number of stochastic elements (right).

Fig. 5 . 7 .

 57 Fig. 5.7. Space-time diagrams (x, t) ∈ [0, 1] × [0, 1] of the averaged directional depths and of the aspect ratio.

Fig. 5 . 8 .

 58 Fig. 5.8. Convergence of the semi-discrete error ǫ sd at time t n = 0.5 for different values of η ∈ [10 -2 , 10 -5 ] and different polynomial orders No ∈ {2 . . . 5}. The left plot reports the error as a function of the total number of stochastic elements, while the right plot shows the error as a function of the total number of degrees of freedom in the stochastic approximation space.

Fig. 5 . 9 .

 59 Fig. 5.9. Dependence of the CPU time (per time iteration) on the stochastic discretization measured by the total number of leaves; left: No = 2 and η = 10 -3 ; right: No = 3 and η = 10 -4 . The contributions of the various steps of the adaptive algorithm are also shown.

2 HU 2 dξU L 2 (

 222 Appendix A. Derivation of the directional indicator. Let d ∈ {1 . . . N} and let ∼d denote all the directions except d. Let n be a node of a binary tree T. Let U ∈ L 2 (Ξ). We recall that Ψn,dβ , β ∈ {1 . . . No + 1}, is a function of ξ d only such that Ψn,d β L 2 (Ξ) = 1. Therefore, |ū n,d β | = inf P ∈PNo[ξ d ] U -P, Ψn,d β = inf P ∈PNo[ξ d ] S(n) (U (ξ ∼d , ξ d ) -P (ξ d )) Ψn,d β (ξ d )dξ = inf P ∈PNo[ξ d ] |S(n)| ∼d S d (n) Ū n ∼d (ξ d ) -P (ξ d ) Ψn,d β (ξ d )dξ ≤ C|S(n)| ∼d |S(n)| No+1 d Ū n ∼d H No+1 (S d (n)) H No+1 (S d (n)) ,(A.1)whereŪ n ∼d (ξ d ) = |S(n)| -1 ∼d S ∼d (n) U (ξ ∼d , ξ d )dξ∼d is the marginalization of U (ξ) over the support S(n) in all the directions ∼d. Furthermore (omitting the reference to the node n),Ū∼d No+1 (S d ) = (ξ ∼d , ξ d )dξ ∼d 2 d U (ξ ∼d , ξ d )dξ ∼d S ∼d ,H No+1 (S d )) , (A.2) with anisotropic Sobolev norm U 2 L 2 (S ∼d ,H No+1 (S d )) = S ∼d U (ξ ∼d , •) 2 H No+1 (S d) dξ ∼d . Combining (A.1) with (A.2) yields (4.11). Appendix B. Sensitivity indices. Let U ∈ L 2 (Ξ) with Ξ = [0, 1] N . The random functional U admits a unique hierarchical orthogonal decomposition, called Hoeffding or Sobol decomposition, of the form [4, 8, 15]U (ξ) = v⊆{1...N} U v (ξ v ), (B.1)where v is a (possibly empty) set of ordered integers v = (v 1 . . . v m ) with m = card(v) =: |v|, ξ v = (ξ v1 . . . ξ vm ), and such that[0,1] U v (ξ v )dξ d = 0 ∀d ∈ v, ∀v ⊆ {1 . . . N} (B.2) Ξ U v (ξ v )U v ′ (ξ v ′ )dξ = 0 ∀v, v ′ ⊆ {1 . . . N}, v = v ′ . (B.3)Each function U v in the decomposition (B.1) can be explicitly expressed in terms of marginals, namelyU ∅ = U , U {d} (ξ d ) = U \ξ d -U , U v (ξ v ) = U \ξ vv ′ ⊂v U v ′ , |v| ≥ 2, (B.4)whereU \ξ v = [0,1] |v| U (ξ ∼v , ξ v )dξ vdenotes the marginalization of U over Ξ with respect to all variables ξ d , d ∈ v, and ∼v := {1 . . . N} \ v. The Sobol decomposition is very useful for sensitivity analysis. Denoting Var(U ) the variance of U , we infer, owing the orthogonality property (B.3), Var(U ) = v⊆{1...N} Var(U v ). (B.5) Of particular importance [9] in characterizing the relative influence of the uncertain parameters ξ d on the variability of U are the first-order sensitivity indices S d and total sensitivity indices T d defined respectively by

7 )

 7 We then build T v as the minimal tree such that the interior of the support of each leaf contains one point xl . The procedure is schematically illustrated in Figure B.1. This construction yields |S(l)| d = 1, for all l ∈ L(T v ) and d ∈ ∼v. The next step consists in projecting U T ∈ S(T) onto S(T v ), an operation which amounts to performing restrictions along the directions d / ∈ v and predictions (injection) in directions d ∈ v.

Fig. B. 1 .

 1 Fig. B.1. Construction principle for T {1} from a tree T and N = 2. The initial tree T is shown in the left plot, with the leaves centers plotted as circles. The right plot shows the set of distinct leaves centers {x l }, as triangles, and the resulting tree T {1} .

* This work is partially supported by GNR MoMaS (ANDRA, BRGM, CEA, EdF, IRSN, PACEN-CNRS). O.P. Le Maître is partially supported by the French National Research Agency (Grants ANR-08-JCJC-0022 & ANR-2010-Blan-0904).


Jean-François Casta

Olivier J Ramond

Stephen W Lin

Value

Value Relevance of Summary Accounting Income Measures: Evidence from Major European Capital Markets

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.