A massively parallel implementation of multilevel Monte Carlo for finite element models

Available online at www.sciencedirect.com ScienceDirect Mathematics and Computers in Simulation 213 (2023) 18–39 www.elsevier.com/locate/matcom Original articles A massively parallel implementation of multilevel Monte Carlo for finite element models✩ Santiago Badiaa,b, Jerrad Hamptonb, Javier Principec,b,∗ aSchool of Mathematics, Monash University, Clayton, Victoria, 3800, Australia bCentre Internacional de Mètodes Numèrics en Enginyeria, Esteve Terrades 5, E-08860 Castelldefels, Spain cUniversitat Politècnica de Catalunya, Campus Diagonal Besòs, Av. Eduard Maristany 16, Edifici A (EEBE), 08019, Barcelona, Spain Received 7 July 2022; received in revised form 26 April 2023; accepted 18 May 2023 Available online 29 May 2023 Abstract The Multilevel Monte Carlo (MLMC) method has proven to be an effective variance-reduction statistical method for Uncertainty Quantification (UQ) in Partial Differential Equation (PDE) models, combining model computations at different levels to create an accurate estimate. Still, the computational complexity of the resulting method is extremely high, particularly for 3D models, which requires advanced algorithms for the efficient exploitation of High Performance Computing (HPC). In this article we present a new implementation of the MLMC in massively parallel computer architectures, exploiting parallelism within and between each level of the hierarchy. The numerical approximation of the PDE is performed using the finite element method but the algorithm is quite general and could be applied to other discretization methods. The two key ingredients of the implementation are a good processor partition scheme together with a good scheduling algorithm to assign work to different processors. We introduce a multiple partition of the set of processors that permits the simultaneous execution of different levels and we develop a dynamic scheduling algorithm to exploit it. The problem of finding the optimal scheduling of distributed tasks in a parallel computer is an NP-complete problem. We propose and analyze a new greedy scheduling algorithm to assign samples and we show that it is a 2-approximation, which is the best that may be expected under general assumptions. On top of this result we design a distributed memory implementation using the Message Passing Interface (MPI) standard. Finally we present a set of numerical experiments illustrating its scalability properties. ©2023TheAuthor(s).PublishedbyElsevierB.V.onbehalfofInternationalAssociationforMathematicsandComputersinSimulation (IMACS).ThisisanopenaccessarticleundertheCCBY-NC-NDlicense(http://creativecommons.org/licenses/by-nc-nd/4.0/). ✩This research was supported by the European Union’s Horizon 2020 research and innovation programme under the ExaQUte project, with grant agreement No 800898, the project RTI2018-096898-B-I00 from the “FEDER/Ministerio de Ciencia e Innovaci´ on - Agencia Estatal de Investigaci´ on” and the Australian Government through the Australian Research Council (project number DP210103092). The authors also acknowledge the Severo Ochoa Centre of Excellence (2019-2023), which financially supported this work under the grant CEX2018-000797-S funded by MCIN/AEI/10.13039/501100011033. The authors thankfully acknowledge the computer resources at MareNostrum and the technical support provided by Barcelona Supercomputing Center, Spain (IM-2019-3-0012). JH has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 712949 (TECNIOspring PLUS) and from the Agency for Business Competitiveness of the Government of Catalonia, Spain. ∗Corresponding author at: Universitat Polit` ecnica de Catalunya, Campus Diagonal Bes` os, Av. Eduard Maristany 16, Edifici A (EEBE), 08019, Barcelona, Spain. E-mail addresses: [email protected] (S. Badia), [email protected] (J. Hampton), [email protected] (J. Principe). https://doi.org/10.1016/j.matcom.2023.05.018 0378-4754/© 2023 The Author(s). Published by Elsevier B.V. on behalf of International Association for Mathematics and Computers in Simulation (IMACS). This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 Keywords: Multilevel Monte Carlo; Uncertainty quantification; Geometric uncertainty; Stochastic partial differential equations; Computational statistics; Parallel programming 1. Introduction UQ requires the solution of stochastic PDEs with random data. Some methods for solving stochastic PDEs, e.g. stochastic Galerkin [27,39], are based on a standard approximation in space like Finite Elements (FEs) or finite volumes, and different types of polynomial expansions in the stochastic direction [56]. Although generally powerful, these techniques typically suffer two significant drawbacks. The first one is being intrusive, i.e. a code that can be used to solve a deterministic problem must be modified to solve the stochastic analogue. The second one is the poor scaling in the dimension of the stochastic space, i.e. suffering from the so-called “curse of dimensionality” [2]. In contrast, sampling methods, i.e. Monte Carlo (MC) and its variants, have a convergence rate which is independent of the stochastic dimension and are non-intrusive. The only necessary assumption for such methods to converge to the exact statistics of the solution when the number of input samples tends to infinity is the existence of finite second moments, i.e. finite variance. This assumption is easily satisfied for many physical systems, making the method generally applicable to many problems of practical interest. They require the repeated evaluation of a deterministic model with different randomly generated inputs and can be implemented in a way that does not require intrusive modifications of the deterministic code. However, the number of samples required to achieve statistical convergence combined with the complexity of the computational model required to have enough spatial and/or temporal accuracy can make the MC method rather expensive. Indeed, the reduction of overall computational cost to achieve a desired error tolerance is the principal motivation for the development of its variants, like the MLMC method considered herein. Stochastic collocation methods [2] are also based on polynomial expansions but using orthogonal basis and appropriate quadratures. The final computation requires independent evaluations of the deterministic model, in the same way as sampling methods, thus sharing the non-intrusiveness property. Appropriate quadratures are required to mitigate the curse of dimensionality, which is a field of intense research [14,55]. Similarly, importance sampling methods [1,33,46], and the deterministic design of experiments [20] are also known to reduce the necessary number of model evaluations, particularly in the contexts of function approximation [40]. The algorithms discussed herein can be extended for parallelization of these “smart sampling” methods as well, but the required modifications are beyond the scope of this article. Focus on this important complexity reduction makes MLMC [28,29,38,41,42,45], and its adaptive variants [15, 16,24,30] arguably the most practical extensions of MC. These methods rely on different levels of computational effort for the deterministic model, e.g. a hierarchy of spatial or temporal meshes. These methods seek to combine samples on each level to benefit from the high accuracy of the expensive, higher levels, and the low computational cost of the less accurate, lower levels. A key aspect of the algorithm is the evaluation of how many samples should be used at each level. Sampling methods for UQ, including MLMC, belong to a general class of methods for outer-loop applications [44, Section 1.2], defined as computational applications in which there is an outer loop around a model that is called at each iteration to evaluate a function. This includes optimization and statistical inference apart from UQ. For example, in the case of optimization under uncertainty, a UQ problem is solved as an inner loop within each step of the outer loop of the optimization algorithm [19]. A big effort was performed during the 90 s for the development of the Dakota library [22,23] aimed to deal with this class of problems. Parallelization opportunities where classified into two areas, algorithmic and function evaluation1and it was already observed that the former requires very little or no inter-processor communication in contrast to the latter. Therefore, the first question posed was how to select the amount of parallelism used on each area, i.e. one has to decide between the assignment of processors to the parallelization of the model or the parallelization of the outer loop (with the concurrent execution of several models). The analysis in [22,23] shows that the most efficient choice is to use the minimum number of processors that permit to run the model, for the simple reason that some loss of efficiency occurs in strong scaling. The parallelization of MLMC has been considered in recent years [11,21,26,47–49,58] as multiple levels of spatial/temporal discretization introduce additional scheduling challenges to effectively utilize resources of a parallel 1The classification in [22,23] actually included four categories but two are enough for the current discussion about MLMC. 19 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 environment with minimal idling. One can conceptually exploit three levels of parallelism in MLMC. The first is the parallelization of the (deterministic) PDE solver itself, which is often dependent on the level, meaning that different levels will use different computational resources. The second is the parallelization of the different MC samples on each level of MLMC, which is the easier one given the independence between them and the large number required, especially at lower discretization levels. The third one is the parallelization between levels, which we found important for the optimality of the algorithm, as discussed below. The two main ingredients of the implementation are (i) a processor partitioning scheme and (ii) a scheduling strategy and the goal is to maximize efficiency. Scheduling strategies for MLMC sampling can be classified [21] in three different ways. First, they are classified according to the number of parallelization layers exploited simultaneously, inside samples, between samples and between levels. Second, they are classified as homogeneous, when all the processors are working to sample on the same discretization level with the same number of processors per sample, and heterogeneous otherwise. It is important to note that, according to this definition, heterogeneous scheduling includes many different possibilities, e.g. using different number of processors for running samples at the same level (which is considered in [21]) or running samples at different levels simultaneously (which is not considered in [21]). A final classification is between static scheduling, when the assignment of work is done before starting to sample, and dynamic scheduling, where assignment is made “on the fly”. It is also important to note that the term “dynamic” is used in [21] to refer to scheduling strategies where the work is assigned depending on a cost estimation made at the beginning of the calculation but before starting to sample. Several scheduling algorithms are studied in [21], mostly homogeneous, but also a couple of heterogeneous variants (with and without strong scaling), both static and dynamic, although three layer parallelism was not investigated. A scheduling with three layer parallelism was proposed in [48–50] with two variants, a static one in which the distribution is made based on work estimates given by the size of the discretization at each level and one in which the distribution of work depends also on the PDE coefficients (which are sample dependent) and is performed at execution time but before starting to sample. In these references, “adaptive load balancing” is used to refer to this second algorithm and it is explicitly stated that it is not “dynamic” see [49, page 111], although it is categorized in this way in [21]. The processor partitioning scheme in [48–50] divides available units into levels, each of them composed by samplers that contain several cores each, that is, the computational power is distributed into many levels and samples at a given time and assignment is made depending on the load balancing strategy just mentioned. For this scheme a 2-approximation is proved in [48,49] assuming strong scaling of the PDE solver, that is, the execution time is at most twice the optimal. The observed efficiency is improved by the adaptive load balancing algorithm although actual computational times are slower than in the static case [48, Figs. 3 and 4], which shows a very good strong and weak scaling up to 40k cores. A similar strategy was followed in [26], where an object oriented implementation of several sampling methods, including MC and Quasi Monte Carlo (QMC), either single or multilevel is presented. Another dynamic scheduling strategy for MLMC was recently proposed in [52]. It exploits general purpose scheduling algorithms in [51] which permit to deal with complex dependencies between tasks and this is exploited to extract parallelism between samples and levels. Although in principle it can be also used with parallel sampling, the mapping of tasks to processors and its effects on load balancing is not discussed. Moreover, parallelization at the sample level is not shown in the numerical examples. In this article we propose the first dynamic scheduling strategy with three layer parallelism that works on top of a multiple partition of the set of available processors. Given the number of levels of the MLMC algorithm and assuming a weak scaling of the PDE solver we define the number of processors per sample in terms of the number of processors used in the coarsest one, as in [21]. We then consider multiple partitions of the set of available processors by these numbers of processors, thus obtaining Ldifferent partitions, where Lis the maximum number of levels of the MLMC method. Together with this multiple partition we present a dynamic scheduling algorithm that assigns tasks to available resources prioritizing finer levels but allowing the concurrent execution of samples at different levels. This distinctive feature of our implementation admits a proof that the algorithm is a 2-approximation, with a very small idling actually occurring only at the end of the computation or at the end of an adaptive MLMC iteration. Besides, excellent scalability is obtained, both strong and weak. With some additional modifications discussed below, this scalability is observed even for the difficult case of models having low per-sample evaluation times. Therefore, the novelty of this approach includes 20 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 •a dynamic scheduling algorithm that exploits three layer parallelism: for each individual sample, across samples, and across levels; •the description of an MPI implementation based on a master–slave strategy with a multiple partition of slave processors; •a simple yet robust, dynamic batch sampling mechanism to reduce the communication between the master and the slaves; •a parallelization of the master for eliminating a bottleneck at extreme scales; •the numerical demonstration of the scalability of the implementation, also when stressed by short sampling times. The article is structured as follows. Section 2describes the MLMC method and an adaptive variant implemented herein and includes the description of the PDE we consider as an example. Section 3describes the scheduling in an abstract way and its 2-approximation property. Section 4describes the actual implementation, including the parallel partition strategy. Section 5presents several examples to demonstrate the weak and strong scalability of the implementation efficiency. 2. The multilevel Monte Carlo method The principal idea of MLMC is to exploit a hierarchy of model discretizations to reduce the overall computational cost by transferring the majority of the sampling cost to the cheaper models, and having an accuracy governed by the more expensive models. The method is appropriate for any hierarchy of models for which convergence of the Quantity of Interest (QoI) is known, although such knowledge is not generally necessary. 2.1. Model problem In this work we consider the following elliptic stochastic problem, although the algorithm permits consideration of a general class of PDE problems using parallel solvers. Given an oriented manifold M(ω)⊂Rdand its corresponding interior domain D(ω), find u(x, ω) such that −∇·(κ∇u)=fin D(ω),u=u0on M(ω),(1) where ω∈Ω, denotes the uncertainty, described by a complete probability space. Additionally, the stochastic coefficients, i.e. the diffusion κ=κ(x, ω), forcing term f=f(x, ω), and boundary condition u0=u0(x, ω) may be considered random fields too. We note that it is usual in the literature to assume these to be deterministic when stochastic domains are considered [13,17,18,34–36,43,57]. Let us assume that the realizations of M(ω) and D(ω) are bounded almost surely. We can define a bounded artificial domain Bthat contains all possible realizations of M(ω), i.e. D(ω)⊂B,∀ω∈Ω. We also assume that κand fare defined in B, independently of ω∈Ω. With the random solution of this problem at hand we aim to compute E(Q(u))where Qis a QoI, e.g. an integral of uon a sub-region or surface, or an evaluation of uat a specific point x, and Ethe expectation. The problem defined by (1) is well-posed under the assumption of uniform ellipticity [2,10,13], i.e. there exists κ0 such that κ(x)≥κ0,∀x∈B(and ∀ω∈Ωif a stochastic diffusion κ(x, ω) is considered). Under these assumptions, the bilinear form associated to the weak form of (1) is bounded and coercive and the Lax–Milgram lemma guarantees a solution for any ω∈Ωuniformly bounded by ∥f∥L2(B)and ∥u0∥H1/2(M(ω)) [13]. The discretization of (1) is constructed on top of a grid T, introducing a C0Lagrangian FE space. When the domain is stochastic we use the Aggregated Finite Element Method (AgFEM) method in which the grid Tis a shape regular partition of a background domain Bwhich is typically a bounding box. A discrete approximation is then constructed identifying cut cells of the background mesh and building a sub-triangulation on each of them. This construction permits the integration of the weak form of the problem and a judicious choice of the degrees of freedom can obtain a well posed-problem. The AgFEM method was introduced in [8], implemented in parallel in [54] and exploited to perform UQ in random domains in [3]. The reader is referred to these publications for further details. 21 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 2.2. The standard MLMC method The MLMC method distributes sampling on a hierarchy of discretizations to reduce the overall computational cost, with respect to that of sampling on the finest one, while keeping a similar accuracy. The hierarchy consists of L+1 meshes T0,T1,...,TLof sizes h0>h1>··· >hL. In particular, we consider hl=h0s−l(i.e. each mesh in the hierarchy is obtained by uniformly dividing each cell into sdsubcells) with s=2. Performing {Nℓ}L ℓ=0 simulations for different values of the random parameters ωon {Tℓ}L ℓ=0, the expectation is corrected using the whole hierarchy as E(QL)=E(Q0)+ L ∑ ℓ=1 E(Qℓ−Qℓ−1)≈ L ∑ ℓ=0 Yℓ:= ˜ QL, where Qℓ=Q(uℓ) is the approximation of the QoI computed using the FE solution uℓfor the discretization Tℓ and Yℓ=Qℓ−Qℓ−1for ℓ=1,...,Land Y0=Q0. The total error of this approximation is now given by e2(˜ QL)=E[(˜ QL−E(Q))2]=(E(QL−Q))2+V(˜ QL).(2) The first term, the (squared) discretization error, is reduced by enforcing that TLprovides a sufficiently accurate approximation. Assuming that it decays as chαfor some constant cand rate αthe maximum level required to make it smaller than, e.g. ε2/2, is L= ⌈1 αlogs(√2cε−1)⌉.(3) Theoretical estimates of αin terms of the regularity of the PDE and the properties of the discretization are available in some cases, e.g. for the linear finite element approximation of elliptic stochastic PDEs α=2, but cand αare generally unknown and can only be estimated a posteriori [24]. The second term in (2) is the statistical error, given by V(˜ QL)=∑L ℓ=0N−1 ℓV(Yℓ). Under general conditions V(Yℓ) tends to zero as ℓ→ ∞ and in this sense, the MLMC method can be understood as a variance reduction method. The total computational cost is given by CMLMC =∑L ℓ=0NℓCℓwhere Cℓis the complexity (computational cost) of evaluating Yℓ, i.e. the average cost to evaluate one sample at level ℓand another one at level ℓ−1. Its minimization allows the optimal number of samples to be taken on each level [28] to have a mean squared error smaller than ε2/2, as Nℓ=⌈2ε−2√V(Yℓ) Cℓ L ∑ i=0√V(Yi)Ci⌉.(4) Observe that the number of samples per level in (4) depends on V(˜ Qℓ), and Cℓand for some problems a theoretical estimate in terms of hℓis possible. These estimates, however, contain unknown constants which need to be determined at an initial screening phase. An alternative is to update them during the execution, as described in the next section. 2.3. An adaptive MLMC method A practical extension of MLMC is given by AMLMC [15,24,30] or a sophisticated variant named Continuation Multilevel Monte Carlo [16], which dynamically enforce the discretization and/or sampling error to be below a given tolerance. This is done by updating the maximum level Land the sample sizes {N0,...,NL}with (ideally) optimal values that reduce the MLMC error with minimal increase in computational cost. Since the actual expectation, variance, and cost are unknown, they are estimated from sample averages and variances, which we refer to as moment estimates. Extrapolation of these moments estimates are also needed when they are not available. The algorithm used here begins with a simple initial set of samples, which may be chosen by any number of considerations. It is presented in Section 4, consisting of a loop in which a set of samples are evaluated, convergence is checked and, if it is not reached, the maximum level Land the sample sizes {N0,...,NL}are updated. In order to update Lan estimation of the constants cand αin (3) is made using YLto estimate the discretization error [24]. 22 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 In order to update Nℓan estimation of variance and cost is required to evaluate (4). Variance is estimated by the sample variance, i.e. V(Yℓ)=1 Nℓ Nℓ ∑ i=1 (Yi ℓ−Yℓ)2. Cost estimates are computed for each level, by adding the total cpu-time used to compute each Yi ℓ, denoted here as Ci ℓand computing the average cost Cℓon that level. From these variance and cost estimates the number of samples needed to achieve the desired error are updated from (4) possibly with an additional fraction to ensure achieving the desired error without to avoid unnecessary outer loop iterations which reduce the efficiency of the algorithm. It is important to observe that each iteration of the adaptive algorithm, in which the number of samples and the maximum number of levels are updated, requires a reduction (to compute averages and sample variances). In practice, it results in a synchronization, as described in Section 4which makes the AMLMC more challenging to scale to large core counts than the standard MLMC. 3. Parallel scheduling In this section we present a new scheduling algorithm for the scheduling of sampling tasks required by the MLMC method in a parallel environment. As mentioned at the end of the previous section, the AMLMC algorithm iteratively updates the number of samples and levels, executing the MLMC algorithm at each iteration, see Remark 3.3. For a general introduction to the subject see [12,53]. In general (and imprecise) terms, the scheduling problem is the development of an algorithm for the assignment of tasks to processors to optimize some cost function, e.g. the total makespan (the maximum run-time). An important distinction to be made is between preemptive and nonpreemptive scheduling. In the first case, the scheduling algorithm is developed under the assumption that tasks can be preempted and restarted later on at no extra cost. We do not consider this possibility in our implementation and we therefore concentrate only on the second case. Some results are available under the assumption that each task is executed in one processor, as in [52]. Given a set of independent tasks T1,...,Tn, whose execution times are w1, . . . , wnrespectively, and a set of pprocessors, the makespan Wof a scheduling is bounded by the average work of each processor W=WT/pwhere WT= n ∑ i=1 wi, and the maximum work of a task is given by WLB =max {W,max i{wi}}. The bound W≥WLB is tight, the equality with the first argument of the maximum occurs when wiare all equal and with the second when p<n. This problem is NP-hard [12]. A simple greedy algorithm to solve this problem is to assign a task to the processor that has the least amount of work already assigned, regardless of the execution time of each task. In this case, it is easy to show that W<2Wopt, where Wopt is the makespan of the optimal scheduling [53], and it is said to be a 2-approximation. Moreover, the bound is tight, which is shown by considering the case of p2tasks of equal processing time followed by a single task whose processing time is ptimes larger. This example suggests that executing the longer tasks first reduces the makespan. This fact was exploited in [32] to develop a greedy algorithm in which tasks are ordered in decreasing execution time and then assigned to the first available processor, see algorithm 5.1.2 in [12]. The makespan bound for this algorithm is W≤(4 3−1 3p)Wopt, that is, a 4/3-approximation. The extension of this result to case of parallel execution of tasks is far from obvious. Besides, it is important to emphasize that the execution time of a given sample in the MLMC method is hardly known in advance. In some cases the execution time can be estimated before the calculation based on the PDE coefficients. For example, in [48], first order hyperbolic problems are approximated using explicit time integration methods with a time step determined by the (random) equation coefficients through the CFL stability condition. However, in general, tasks cannot be ordered by execution time beforehand and therefore a 2-approximation bound 23 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 is the best that can be generally expected. It is worth emphasizing that the actual makespan is much closer to the optimal in practice, see Section 5. Let us now move to the parallel case in which a task may require more than one processor. There are several factors to consider in the decision of how many processors are assigned to each task. In a general case [21], it depends on two factors, the idling that can be generated and the loss of efficiency in the strong scaling of the underlying solver. If no idling occurs, as in the algorithm we propose herein, and there are enough samples to assign work to all processors, the inevitable loss of efficiency in strong scaling of the solver implies that the best efficiency is always obtained assigning the minimum number of processors that permit to execute each sample, as already noticed in [22,23]. On the other hand, executing all tasks with the same number of processors would lead to larger execution times on finer levels of the MLMC hierarchy. Assuming a weakly scalable solver, and assigning a number of tasks proportional to the size of the mesh at each level gives constant mean execution time for each sample, with variations with the change of the random PDE coefficients. This assumption is made in [21], where the number of processors per task at level ℓof the MLMC algorithm is taken of the form 23ℓp0and p0is determined to maximize efficiency. We divide the tasks in charge of each sample in the MLMC algorithm into sets Tq= {Tq 1,...,Tq nq}that require qprocessors with q≤pfor fixed pand we denote by wq 1, . . . , wq nqthe execution (wall-clock) time of each of them. The problem of scheduling these tasks using pprocessors to minimize the makespan is also NP-complete. Therefore, we cannot expect to obtain a polynomial time algorithm, so we propose an algorithm to solve a restricted problem rather than the general one. We restrict the number of possible values of qto q0,q1,...,qMwhere q0|q1|...|qM (hereafter a|bmeans that ais an integer divisor of bwith both a,b∈N, as usual). Observe that the choice in [21] satisfies this assumption. We also assume that qM|p. This assumption of divisibility is not necessary in the scheduling algorithm, but permits to prove Theorem 3.1 presented below, see Remark 3.3. The total work is now WT=∑ q=q0,q1,...,qM q nq ∑ j=1 wq j. and WLB =max {W,max q,j{wq j}}. where, again W=WT/p. Under these restrictions (of the possible values of q0,q1,...,qM) we propose Algorithm 1, whose execution is illustrated in Fig. 1 and we prove that it is a 2-approximation in Theorem 3.1. Observe that because wallclock times of samples are unknown before their actual execution, implementing Algorithm 1requires dynamic programming [12]. A MPI implementation is described in the next section. Algorithm 1: Incremental greedy scheduling 1for q∈ {qM,qM−1,...,q1,q0}do 2for j=1,...,nqdo 3Assign Tq jto a set of qprocessors among those with less workload. 4end for 5end for Theorem 3.1. Given q0,q1,...,qMsatisfying q0|q1|...|qM|p, Algorithm 1is a 2-approximation. Proof. At the beginning the set of pprocessors is divided into groups p/qMof qMprocessors. Because ql−1|ql, once tasks Tql j,j=1,...,nqlhave been executed, each group of qlprocessors can be divided into ql/ql−1groups of ql−1tasks each, which are ready to execute tasks of the group Tql−1 j,j=1,...,nql−1. The previous statement is true for any l=M,...,0 and therefore, no processor idles until there are no more tasks to execute. Once this fact is established, the proof is similar to the case of single processor tasks given in [31] (see [53] for the simplified proof followed here). Consider the group of processors that finishes last and let wsbe the starting time of the last executed task and wlast its execution time. Because the algorithm assigns tasks to processors with less workload, all the rest of the processors are busy at time wsand therefore ws<W. Because wlast <maxi,j{wi j}we get ws+wlast <2WLB <2Wopt.□ 24 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 Fig. 1. An example of timeline execution with Algorithm 1with M=2 and q=2q2=4q1=8q0and n0=10 (green), n1=4 (red) and n2=2 (blue) samples. At t=t1, the n2samples requiring q2processors have been assigned and samples that require q1processors start to be assigned to the least loaded processors (those in the bottom half). At t=t3, the n1samples requiring q1processors have been assigned and samples that require q0processors start to be assigned. At t=t4there are no more samples to assign and some processors become idle. Idling time, occurring only at the end of the calculation, is shown in gray. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Remark 3.2. In the cases where the hypothesis of Theorem 3.1 are not satisfied it is still possible to apply the rationale behind Theorem 3.1. After all tasks Tql j,j=1,...,nqlhave been executed, the group of qlprocessors can be divided into ⌊ql/ql−1⌋subgroups to execute tasks Tql−1 j,j=1,...,nql−1and there will be rl=qlmod ql−1 processors left. If qm|rlfor some m<lthey can be assigned to execute tasks Tqm, otherwise they became idle. However, even if these rlprocessors can be assigned further work, it is possible that they idle before the rest of tasks have been executed. For instance, consider the case of q0=1,q1=4,q2=10. After execution of tasks T10 is finished two groups of 4 processors are created and tasks in T4start being executed while the remaining two processors execute tasks in T1. If all tasks in T1are executed before those in T4these two processors start idling. However, if n1≫n4, which is the case of the MLMC, a good balance can be expected in practice. Remark 3.3. It is important to keep in mind the hypothesis of independence of tasks made at the beginning of this section. Theorem 3.1 is only valid under this assumption. However, this is actually the case in the standard MLMC algorithm described in Section 2.2, apart from some simple post-processing at the end of the execution. In the case of the AMLMC algorithm described in Section 2.3 additional synchronization occurs at each iteration of the algorithm, introducing a dependency between tasks, i.e. those in the second iteration must be executed after those in the first. Therefore, the results of this section actually apply to the scheduling of tasks required to complete one iteration of the AMLMC which, in any case, represents a substantial amount of computational work. 4. A message passing implementation In this section we describe an implementation of the algorithm described in Section 3developed on top of the MPI standard. The two main ingredients of the implementation are a strategy for the partition of the set of processors into groups for parallel sampling, described in Section 4.1, and a dynamic scheduling algorithm, described in Section 4.2. We also include a modification that parallelizes the master coordinator that is designed for improving scalability at the most extreme scales, described in Section 4.3. 4.1. Processor partition strategy The whole strategy is based on a master–slave approach, with a master processor, the coordinator, in charge of the decisions required for scheduling MLMC tasks and the remaining slave processors in charge of sampling (also referred as slaves in the following). The communication between the master and slaves required to implement Algorithm 1is described in Section 4.2. We therefore assume we are given p+1 processors, one master and p slaves. 25 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 Fig. 2. A set of 33 processors and their family of partitions into groups of q2=16, q1=8 and q0=4 processors each. Here rank 0 is the master (not shown) and ranks 1 to 32 are slaves. They are first (concurrently) split into groups of 16 processors represented by green lines with roots 1 and 17. These groups are then split into groups represented by blue lines, whose roots are 1, 9, 17 and 25. Finally the groups are subsequently split into groups having 4 ranks each, which are represented in red (also used to signal their roots). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) As described in Section 3, slave processors will be required to work on different samples and, more importantly, coordinate with different slaves depending on the sample they are executing. Therefore, different partitions are required at different instants during the calculation. To satisfy this requirement we generate a family of partitions of the set of slaves. The total number of partitions, denoted by Mas in the previous section, is fixed during the calculation and levels of the MLMC algorithm are mapped to each partition. Therefore, the same index ℓis used in the following to denote the MLMC level and the partition in which samples of Yℓrun. When the AMLMC algorithm described in Section 2.3 is used, a value of Mmust be defined by the user such that M>Lduring the whole calculation and the number of samples for levels ℓ > Lare set to 0. The family of partitions is generated in a hierarchical manner, which requires q0<q1<··· <qM<p. The set of slave processors is first divided into ⌊p/qM⌋groups of qMand a group with rM=pmod qMprocessors. In a second step, each group of qMprocessors is divided into ⌊qM|qM−1⌋groups of qM−1and one group of rM−1=qM mod qM−1processors. The remaining group of rMprocessors is divided into ⌊rM|qM−1⌋(which can be 0) groups of qM−1processors and one group of size rM−1=rMmod qM−1. The process continues recursively until the final groups of q0processes are generated. The first processor of each group, which we will refer to as the root, will play an important role in the algorithm described in Section 4.2. To fix ideas, consider p+1=33 processors that will be used to sample with tasks requiring q2=16, q1=8 and q0=4 processors each, as illustrated in Fig. 2. The p=32 slave processors will be first split into 2 groups of 16 processors each. The second partition of the family will be constructed splitting each of them into 2 groups of size 8. Finally these 4 groups will be split into 2 sub-communicators having 4 processors each. Each slave processor will belong to 3 different groups. In this way, any slave processor may be used for the execution of a task requiring 4, 16 or 32 processors. On the other hand, consider p+1=31 processors that will be used to sample with tasks requiring q2=15, q1=6 and q0=3 processors each, as illustrated in Fig. 3. The p=30 slave processors will be first split into 2 groups of 15 processors each. The second partition of the family will be constructed splitting each group into 2 groups having 6 processors and one group having 3 processors. Finally, each group of the second partition will be split into one or two groups having 3 processors. 26 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 Fig. 4. Benchmarking speedup under weak scaling (left) and strong scaling (right) for different sampling times (short, medium and large), with and without parallel coordinator of the scheduling. Fig. 5. Benchmarking efficiency under weak scaling (left) and strong scaling (right) for different sampling times (short, medium and large), with and without parallel coordinator of the scheduling. Fig. 4 shows the speedup for the two scaling regimes given by (9), while Fig. 5 shows the efficiency. We note that the parallel benchmarking performs well for the Medium Time and Long Time models, but that the Short Time model provides a significant communication stress, and that in this case, the use of additional communication tasks improve scalability. We also consider two additional modifications, beyond the consideration of additional communicators and execution times. For these, we use only one communicator task, and µ=10−2s. First, is the variance of the execution time. We consider three cases: that of high σ, given by σH=0.5µ/√3≈0.289µ; that of medium σ, given by σM=0.2µ; and that of low σ, given by σL=0.01µ. These results are shown in Fig. 6. We note that changes to σprovide at most a marginal change in performance. Therefore, we can conclude that the variation of execution times among samples does not significantly affect the scalability properties of the algorithm, neither in a weak or strong scaling case. Additionally, we consider the use of batch sampling. Using a single master communicator, we consider the case when all batch sizes are of size 1 for µ=10−2s. These results are shown in Fig. 7. It can be seen that there is worse scalability at higher core counts without batch sampling, as the single master communicator cannot act fast enough to efficiently communicate each individual sample to the execution partitions. Therefore, we can conclude that the use of batching sampling can be useful to improve scalability for extreme core counts. We end this section with the analysis of the influence of the number of cores in the scalability. In Fig. 8 we compare the performance of two different sets of processors for sampling, namely, q=8,64,512 and q=9,81,729. As it can be observed in Fig. 8 the influence of this choice is minimal, as anticipated in Remark 3.2. Although this result could be surprising, we note that there is always a master processor in charge of the scheduling that introduces an oddness and therefore the set of processors is never partitioned in an optimal way in the architectures considered herein. For example, the leftmost point in all figures corresponds to 16 nodes with 48 33 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 Fig. 6. Benchmarking speedup under weak scaling (left) and strong scaling (right) for different variability in the sampling times. Fig. 7. Benchmarking speedup under weak scaling (left) and strong scaling (right) with and without batch scheduling. Fig. 8. Benchmarking speedup under weak scaling (left) and strong scaling (right) for different sets of processors. cores, that is, p+1=768 processors, one of which is used for the master. The partition of p=767 will result in some idling in both cases. In fact, in the first case it is possible to exploit 512 +3×64 +7×8=760 processors whereas in the second it is possible to exploit 729 +4×9=765 processors. In any case, this idling does not increase with the number of nodes and scalability is not affected. 5.2. 3D Poisson in a random domain Here we consider a discretized finite element model, where the number of tasks per model evaluation depends on the level. To this end, we describe the details of a 3D Poisson problem as in (1). 34 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 Here we assume the domain to be a stochastic “popcorn” shape as defined in [3], which provides a complex geometry that changes for each realized sample. This geometry is described by a relatively large and variable number of random variables, i.e. stochastic dimension is variable. It is composed of nsmall ellipses added on top of a large one, with n∈P(11), i.e. a Poisson distributed random variable with mean value 11, requiring a total of 5n+7 random variables to define the final shape. This gives a mean of 62 random variables and the extra complexity of changing between samples. Moreover, the diffusion term κ(x, ω), is a KLE where every point is log-normally distributed, and the associated normal distributed random variables are correlated by an exponential covariance function with correlation length 0.25. The normal random variables have 0 mean, and a standard deviation of 0.25, with the first 16 eigenfunctions in each dimension, used to construct the approximation with a total-order basis of order 16. Here, total-order 16 refers to a tensor product of polynomials in each dimension such that the sum of the orders of those polynomials is at most 16. This gives a total of 680 functions, which require an equal number of random variable coefficients. The boundary conditions are given by a sin|x|and the quantity of interest is the mean value of the solution in a square region inside the domain, see [3] for details. The linear system arising from the discretization of the PDE is solved using a parallel algebraic multigrid method [37] provided by HYPRE library [25] through the PETSC interface. We consider L=2 (or L0=1 in the case of AMLMC) levels of computation with q0=8,q1=64, and q2=512. Here we consider a uniform coarsest grid of 16 ×16 ×16 elements, and each subsequent grid has twice the number of elements in each dimension, up to 64 ×64 ×64 elements for ℓ=2. The AMLMC here has initial samples on all three levels as given by (8) and Table 1 for the appropriate Cfor {16,32,64,128,256,512,600} nodes. For strong scaling, we consider the same model parameters under strong scaling, with the sample multiplier Cfixed to 128, while for weak scaling Cis the number of computational nodes. Fig. 9 shows the speedup for the two scaling regimes given by (9). We note that there are no noticeable issues scaling up to the 600 nodes tested for this MLMC example. Finally, we utilize the AMLMC algorithm with this model, to test the effect of the synchronization points which that algorithm requires. For weak scaling we consider a set of desired tolerances for desired error, in this case of the volume of the associated geometric shape. Specifically, the error tolerances are chosen to be between 10−4and 1.617 ·10−5. These tolerances are chosen to allow the same initial sampling of the other 3D Poisson example, but such that those samples are insufficient to converge in the initial iteration. The strong scaling utilizes C=128, and an error tolerance corresponding of 3.535 ·10−5. We note that as the ultimate sample sizes for this experiment are variable, the performance is more variable. Of note, each execution leads to a different final achieved error. We consider a raw scaling, and also consider a scaling that mitigates this discrepancy, by scaling the results inversely proportional to the achieved error. That is, lower achieved errors will have used more samples, and thus more computational resources. We adjust the wall times to these values, so that the results are neutral with respect to the achieved error of the final AMLMC estimation. The results are presented in Fig. 10, and Fig. 11. We note that both the strong and weak scaling are significantly more efficient for smaller pfor this problem, but that this is not a result of significant idling, so that the scalability is reasonable as pincreases. 6. Conclusions In this work we present a new dynamic scheduling algorithm to compute samples in the MLMC method that fully exploits three levels of parallelism, within each sample, across samples, and across levels. The implementation of this algorithm based on the MPI standard described here is based on general primitives available in any version of the standard and therefore in any MPI distribution, i.e. it does not require any advanced feature like one-sided communication. We also demonstrate the performance of method when applied to MLMC, on a stress test of the scheduling algorithm on a benchmarking problem, as well as a 3D Poisson model solved on a complex, sample dependent domain via AMG. The implementation is a powerful, robust, and scalable means to perform UQ via MLMC or other similar statistical computations. The numerical examples presented herein show excellent scalability on the examples tested, when certain good scaling features are utilized. The algorithm implementation also performs well when utilized for AMLMC. Overall, our results show performance that competes with that demonstrated in 35 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 Fig. 9. 3D Poisson Popcorn speedup under weak scaling (left) and strong scaling (right). Fig. 10. 3D Poisson Popcorn AMLMC speedup under weak scaling (left) and strong scaling (right). Fig. 11. 3D Poisson Popcorn AMLMC efficiency under weak scaling (left) and strong scaling (right). the similar experiments in [21,52] (which only show results exploiting two layers of parallelism), showing high efficiency, and notable performance on the difficult case of rapidly executing models. The results in this work could also be extended by exploiting the techniques in [51] to account for the dependency between tasks in the context of parallel sampling. This is particularly relevant for the AMLMC method where some tasks have a dependency on other tasks that require more processors. Although this can be handled just by an iterative loop as described above, a dynamic scheduling with task dependency could potentially be more efficient. This development require to define appropriate priorities taking into account the amount of parallelism at sample level and is left for a future work. 36 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 References [1] Ben Adcock, Anders C. Hansen, Clarice Poon, Bogdan Roman, Breaking the coherence barrier: A new theory for compressed sensing, in: Forum of Mathematics, Sigma, Vol. 5, 2017. [2] Ivo Babuška, Fabio Nobile, Raúl Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM Rev. 52 (2) (2010) 317–355. [3] Santiago Badia, Jerrad Hampton, Javier Principe, Embedded multilevel Monte Carlo for uncertainty quantification in random domains, Int. J. Uncertain. Quantif. 11 (1) (2021) 119–142. [4] Santiago Badia, Alberto F. Martín, A tutorial-driven introduction to the parallel finite element library FEMPAR v1.0.0, Comput. Phys. Comm. 248 (2020) 107059. [5] Santiago Badia, Alberto F. Martín, Javier Principe, Implementation and scalability analysis of balancing domain decomposition methods, Arch. Comput. Methods Eng. 20 (3) (2013) 239–262. [6] Santiago Badia, Alberto F. Martín, Javier Principe, Multilevel balancing domain decomposition at extreme scales, SIAM J. Sci. Comput. 38 (1) (2016) C22–C52. [7] Santiago Badia, Alberto F. Martín, Javier Principe, FEMPAR: An object-oriented parallel finite element framework, Arch. Comput. Methods Eng. 25 (2) (2018a) 195–271. [8] Santiago Badia, Francesc Verdugo, Alberto F. Martín, The aggregated unfitted finite element method for elliptic problems, Comput. Methods Appl. Mech. Engrg. 336 (2018b) 533–553. [9] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil M. Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, Junchao Zhang, PETSc web page, 2023. [10] Andrea Barth, Christoph Schwab, Nathaniel Zollinger, Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic coefficients, Numer. Math. 119 (1) (2011) 123–161. [11] Niklas Baumgarten, Christian Wieners, The parallel finite element system M++ with integrated multilevel preconditioning and multilevel Monte Carlo methods, Comput. Math. Appl. 81 (2021) 391–406. [12] Jacek Blazewicz, Klaus H. Ecker, Erwin Pesch, Günter Schmidt, Malgorzata Sterna, Jan Weglarz, Handbook on Scheduling, Springer International Publishing, Cham, 2019. [13] Jehanzeb H. Chaudhry, Nathanial Burch, Donald Estep, Efficient distribution estimation and uncertainty quantification for elliptic problems on domains with stochastic boundaries, SIAM-ASA J. Uncertain. Quantif. 6 (3) (2018) 1127–1150. [14] Peng Chen, Sparse quadrature for high-dimensional integration with Gaussian measure, ESAIM Math. Model. Numer. Anal. 52 (2) (2018) 631–657. [15] K.A. Cliffe, M.B. Giles, R. Scheichl, A.L. Teckentrup, Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients, Comput. Vis. Sci. 14 (1) (2011) 3–15. [16] Nathan Collier, Abdul Lateef Haji-Ali, Fabio Nobile, Erik von Schwerin, Raúl Tempone, A continuation multilevel Monte Carlo algorithm, BIT Numer. Math. 55 (2) (2015) 399–432. [17] M. Dambrine, I. Greff, H. Harbrecht, B. Puig, Numerical solution of the Poisson equation on domains with a thin layer of random thickness, SIAM J. Numer. Anal. 54 (2) (2016) 921–941. [18] M. Dambrine, I. Greff, H. Harbrecht, B. Puig, Numerical solution of the homogeneous Neumann boundary value problem on domains with a thin layer of random thickness, J. Comput. Phys. 330 (2017) 943–959. [19] Subhayan De, Jerrad Hampton, Kurt Maute, Alireza Doostan, Topology optimization under uncertainty using a stochastic gradient-based approach, Struct. Multidiscip. Optim. 62 (5) (2020) 2255–2278. [20] Paul Diaz, Alireza Doostan, Jerrad Hampton, Sparse polynomial chaos expansions via compressed sensing and D-optimal design, Comput. Methods Appl. Mech. Engrg. 336 (2018) 640–666. [21] D. Drzisga, B. Gmeiner, U. Rüde, R. Scheichl, B. Wohlmuth, Scheduling massively parallel multigrid for multilevel Monte Carlo methods, SIAM J. Sci. Comput. 39 (5) (2017) S873–S897. [22] M. Eldred, W. Hart, Design and implementation of multilevel parallel optimization on the intel teraflops, in: 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, American Institute of Aeronautics and Astronautics, St. Louis, MO, 1998, pp. 44–54. [23] M. Eldred, W. Hart, B. Schimel, B. Waanders, Multilevel parallelism for optimization on MP computers - theory and experiment, in: 8th Symposium on Multidisciplinary Analysis and Optimization, American Institute of Aeronautics and Astronautics, Long Beach, CA, 2000. [24] Daniel Elfverson, Fredrik Hellman, Axel Malqvist, A multilevel Monte Carlo method for computing failure probabilities, SIAM-ASA J. Uncertain. Quantification 4 (1) (2016) 312–330. [25] Robert D. Falgout, Jim E. Jones, Ulrike Meier Yang, The design and implementation of hypre, a library of parallel high performance preconditioners, in: Are Magnus Bruaset, Aslak Tveito (Eds.), Numerical Solution of Partial Differential Equations on Parallel Computers, Springer Berlin Heidelberg, Berlin, Heidelberg, 2006, pp. 267–294. 37 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 [26] Robert N. Gantner, A generic C++ library for multilevel quasi-Monte Carlo, in: PASC 2016 - Proceedings of the Platform for Advanced Scientific Computing Conference, 2016, pp. 1–12. [27] Roger G. Ghanem, Pol D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag New York, 1991. [28] Michael B. Giles, Multilevel Monte Carlo path simulation, Oper. Res. 56 (3) (2008) 607–617. [29] Michael B. Giles, Christoph Reisinger, Stochastic finite differences and multilevel Monte Carlo for a class of SPDEs in finance, SIAM J. Financial Math. 3 (1) (2012) 572–592. [30] Michael B. Giles, Ben J. Waterhouse, Multilevel quasi-Monte Carlo path simulation, in: Advanced Financial Modelling, in: Radon Series Comp. Appl. Math, vol. 8, de Gruyter, Berlin, 2009, pp. 165–181. [31] R.L. Graham, Bounds for certain multiprocessing anomalies, Bell Syst. Tech. J. 45 (9) (1966) 1563–1581. [32] R.L. Graham, Bounds on multiprocessing timing anomalies, SIAM J. Appl. Math. 17 (2) (1969) 416–429. [33] Jerrad Hampton, Alireza Doostan, Compressive sampling of polynomial chaos expansions: Convergence analysis and sampling strategies, J. Comput. Phys. 280 (2015) 363–386. [34] Helmut Harbrecht, Jingzhi Li, First order second moment analysis for stochastic interface problems based on low-rank approximation, Math. Modelling Numer. Anal. 47 (5) (2013) 1533–1552. [35] H. Harbrecht, M. Peters, M. Siebenmorgen, Analysis of the domain mapping method for elliptic diffusion problems on random domains, Numer. Math. 134 (4) (2016) 823–856. [36] Helmut Harbrecht, Reinhold Schneider, Christoph Schwab, Sparse second moment analysis for elliptic problems in stochastic domains, Numer. Math. 109 (3) (2008) 385–414. [37] Van Emden Henson, Ulrike Meier Yang, Boomeramg: A parallel algebraic multigrid solver and preconditioner, Appl. Numer. Math. 41 (1) (2002) 155–177, Developments and Trends in Iterative Methods for Large Systems of Equations - in memorium Rudiger Weiss. [38] Ahmed Kebaier, Statistical romberg extrapolation: A new variance reduction method and applications to option pricing, Ann. Appl. Probab. 15 (4) (2005) 2681–2705. [39] Olivier Le Maitre, Omar M. Knio, Spectral Methods for Uncertainty Quantification: With Applications To Computational Fluid Dynamics, Springer, 2010, p. 536. [40] Nora Lüthen, Stefano Marelli, Bruno Sudret, Sparse polynomial chaos expansions: Literature survey and benchmark, SIAM-ASA J. Uncertain. Quantification 9 (2) (2021) 593–649. [41] S. Mishra, Ch Schwab, J. Šukys, Multi-level Monte Carlo finite volume methods for nonlinear systems of conservation laws in multi-dimensions, J. Comput. Phys. 231 (8) (2012a) 3365–3388. [42] S. Mishra, Ch Schwab, J. Šukys, Multilevel Monte Carlo finite volume methods for shallow water equations with uncertain topography in multi-dimensions, SIAM J. Sci. Comput. 34 (6) (2012b) B761–B784. [43] P. Surya Mohan, Prasanth B. Nair, Andy J. Keane, Stochastic projection schemes for deterministic linear elliptic partial differential equations on random domains, Internat. J. Numer. Methods Engrg. 85 (7) (2011) 874–895. [44] Benjamin Peherstorfer, Karen Willcox, Max Gunzburger, Survey of multifidelity methods in uncertainty propagation, inference, and optimization, SIAM Rev. 60 (3) (2018) 550–591. [45] M. Pisaroni, F. Nobile, P. Leyland, A continuation multi level Monte Carlo (C-MLMC) method for uncertainty quantification in compressible inviscid aerodynamics, Comput. Methods Appl. Mech. Engrg. 326 (2017) 20–50. [46] Holger Rauhut, Rachel Ward, Sparse Legendre expansions via ℓ1-minimization, J. Approx. Theory 164 (5) (2012) 517–533. [47] Nikolay Shegunov, Oleg Iliev, On dynamic parallelization of multilevel Monte Carlo algorithm, Cybern. Inf. Technol. 20 (6) (2020) 116–125. [48] Jonas Šukys, Adaptive Load Balancing for Massively Parallel Multi-Level Monte Carlo Solvers, in: Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 8384 LNCS, (PART 1) Springer, Berlin, Heidelberg, 2014a, pp. 47–56. [49] Jonas Šukys, Robust multi-level Monte Carlo finite volume methods for systems of hyperbolic conservation laws with random input data, Vol. 24 (Ph.D. thesis), 2014b, pp. 295–305. [50] Jonas Šukys, Siddhartha Mishra, Christoph Schwab, Static load balancing for Multi-Level Monte Carlo finite volume solvers, in: Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 7203 LNCS, (PART 1) Springer, Berlin, Heidelberg, 2012, pp. 245–254. [51] Enric Tejedor, Yolanda Becerra, Guillem Alomar, Anna Queralt, Rosa M. Badia, Jordi Torres, Toni Cortes, Jesús Labarta, Pycompss: Parallel computational workflows in Python, Int. J. High Perform. Comput. Appl. 31 (1) (2017) 66–82. [52] Riccardo Tosi, Ramon Amela, Rosa Badia, Riccardo Rossi, A parallel dynamic asynchronous framework for uncertainty quantification by hierarchical Monte Carlo algorithms, J. Sci. Comput. 89 (28) (2021). [53] Vijay V. Vazirani, Approximation Algorithms, Springer Berlin Heidelberg, Berlin, Heidelberg, 2003, pp. 423–469. [54] Francesc Verdugo, Alberto F. Martín, Santiago Badia, Distributed-memory parallelization of the aggregated unfitted finite element method, Comput. Methods Appl. Mech. Engrg. 357 (2019) 112583. [55] Dongbin Xiu, Jan S. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput. 27 (3) (2005) 1118–1139. [56] Dongbin Xiu, George Em Karniadakis, Modeling uncertainty in flow simulations via generalized polynomial chaos, J. Comput. Phys. 187 (1) (2003) 137–167. 38 S. Badia, J. Hampton and J. Principe Mathematics and Computers in Simulation 213 (2023) 18–39 [57] Dongbin Xiu, Daniel M. Tartakovsky, Numerical methods for differential equations in random domains, SIAM J. Sci. Comput. 28 (3) (2006) 1167–1185. [58] Petr Zakharov, Oleg Iliev, Jan Mohring, Nikolay Shegunov, Parallel Multilevel Monte Carlo Algorithms for Elliptic PDEs with Random Coefficients, in: Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 11958 LNCS, 2020, pp. 463–472. 39