Optimal fidelity multi-level Monte Carlo for quantification of uncertainty in simulations of cloud cavitation collapse

OPTIMAL FIDELITY MUL TI-LEVEL MONTE CARLO F OR QUANTIFICA TION OF UNCER T AINTY IN SIMULA TIONS OF CLOUD CA VIT A TION COLLAPSE ∗ JONAS ˇ SUKYS † , URSULA RASTHOFER ‡ , F ABIAN WERMELINGER ‡ , P ANAGIOTIS HADJIDOUKAS ‡ , AND PETROS KOUMOUTSAK OS ‡ § Abstract. W e quantify uncertainties in the location and magnitude of extreme pressure sp ots revealed from large scale multi-phase ﬂow sim ulations of cloud ca vitation collapse. W e examine clouds containing 500 ca vities and quantify uncertain ties related to their initial spatial arrangemen t. The resulting 2000-dimensional space is sampled using a non-intrusiv e and computationally eﬃcient Multi-Level Mon te Carlo (MLMC) methodology . W e introduce no vel optimal control v ariate coeﬃ- cients to enhance the v ariance reduction in MLMC. The prop osed optimal ﬁdelity MLMC leads to more than tw o orders of magnitude speedup when compared to standard Monte Carlo metho ds. W e identify large uncertainties in the location and magnitude of the p eak pressure pulse and presen t its statistical correlations and joint probability density functions with the geometrical characteristics of the cloud. Characteristic prop erties of spatial cloud structure are iden tiﬁed as potential causes of signiﬁcant uncertainties in exerted collapse pressures. Key w ords. compressible multi-phase ﬂow, high performance computing, diﬀuse interface method, bubble collapse, cloud ca vitation, uncertain ty quan tiﬁcation, multi-lev el Mon te Carlo, opti- mal control v ariates, fault tolerance AMS sub ject classiﬁcations. 76B10, 68W10, 65C05, 65C60, 68M15. 1. Introduction. Cloud ca vitation collapse pertains to the inception of m ultiple gas ca vities in a liquid, and their subsequent rapid collapse driven by an increase of the am bien t pressure. Shock w av es emanate from the cavities with pressure p eaks up to tw o orders of magnitude larger than the ambien t pressure [ 46 , 12 , 7 ]. When suc h shock w av es impact on solid w alls, they ma y cause material erosion, considerably reducing the p erformance and longevit y of turb o-mac hinery and fuel injection engines [ 65 , 38 ]. On the other hand, the destructive potential of cavitation can b e harnessed for non-in v asive biomedical applications [ 17 , 8 ] and eﬃcient tra v el in military sub- surface applications [ 5 ]. Prev alen t conﬁgurations of cavitation often inv olv e clouds of h undreds or thousands of bubbles [ 7 ]. The cloud collapse entails a location depen- den t, often non-spherical, collapse of eac h ca vity that progresses from the surface to the cen ter of the cloud. Pressure wa v es emanating from collapsing cavities lo cated near the surface of the cloud act as ampliﬁers to the subsequen t collapses at the cen- ter of the cloud. The in teraction of these pressure wa v es increases the destructive p oten tial as compared to the single bubble case. Ca vitation, in particular as it o c- curs in realistic conditions, presents a formidable challenge to exp erimental [ 6 , 10 , 79 ] and computational [ 1 , 72 ] studies due to its geometric complexit y and m ultitude of spatio-temp oral scales. Blake et al. [ 4 ] studied the single bubble asymmetric collapse ∗ Submitted to the editors 9 May 2017. F unding: INCITE pro ject CloudPredict; Argonne Leadership Computing F acility , DOE project DE-AC02-06CH11357; PRA CE: JFZ and CINECA, pro jects PRA091 and Pra09 2376; CSCS pro ject ID s500. † Computational Science and Engineering Laboratory , ETH Z¨ urich, Switzerland. Current address: Eaw ag, Swiss F ederal Institute of Aquatic Science and T echnology , Switzerland. ( jonas.sukys@eaw ag.ch , science.sukys.ch ). ‡ Computational Science and Engineering Laboratory , ETH Z¨ urich, Switzerland ( urasthofer,fabianw,phadjido,petros@mavt.ethz.c h , cse- lab.ethz.ch ). § Corresponding author. 1 2 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS using a b oundary integral method. Johnsen and Colonius [ 34 ] in v estigated the p oten- tial damage of single collapsing bubbles in b oth spherical and asymmetric regime for a range of pulse peak pressures in sho c k-induced collapse. Lauer et al. [ 37 ] s tudied collapses of arrays of cavities under sho c k wa ves using the sharp interface technique of Hu et al. [ 31 ]. Recen t numerical inv estigation of cloud cavitation inv olved a cluster of 125 v ap or bubbles inside a pressurized liquid at 40 bar [ 18 , 1 ], and a cluster of 2500 gas bubbles with ambien t liquid pressure of 100 bar [ 61 ]. Large scale numerical sim ulations of cloud cavitation collapse considered clouds containing 50’000 bubbles [ 28 ]. Visualizations of suc h a collapsing cloud and the resulting fo cused pressure p eak at the center are reproduced in Figure 1 . Ho wev er the computational demands of these simulations do not allo w for further parametric studies. Fig. 1 . Surfac e view of a c ol lapsing cloud c ontaining 50’000 gas cavities (left) and pr essur e p e ak at the c enter of the cloud (right). Outer bubbles evolve into c ap-like shap es, forming an array of inwar ds dir e cte d velo city micro-jets. The r esulting focuse d pr essur e p e ak ampliﬁcations of two or ders of magnitude ar e exerte d at the cloud center [ 28 ]. A challenge in mo deling and quantifying cloud cavitation collapse is the dep en- dence of critical Quantities of In terest (QoIs), such as p eak pressure or collapse time, on a particular (random) cloud conﬁguration [ 72 ] (see also Figure 2 ). The systematic 0.0 0.2 0.4 0.6 0.8 1.0 T 0 5 10 15 20 25 Maximum Pressure cloud #1 cloud #2 Fig. 2 . Evolution of the normalize d p e ak pr essur e versus normalize d time in fr e e-ﬁeld c ol lapse of two spheric al clouds c ontaining 1’000 bubbles. Two samples are dr awn fr om the same uniform dis- tribution (for the p ositions of the c avities) and lo g-normal distribution (for the r adii of the c avities). OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 3 study of such dep endencies can b e addressed through an Uncertain ty Quan tiﬁcation (UQ) framew ork, recen tly applied in [ 16 , 5 ]. In [ 42 , 43 , 44 ], a mathematical framew ork w as provided for uncertain solutions of hyperb olic equations. P opular probability space discretization methods include generalized P olynomial Chaos (gPC) tec hniques (please see [ 13 , 39 , 74 , 58 , 75 , 24 ] and references therein). An alternative class of metho ds for quantifying uncertaint y in PDEs are the sto chastic collo cation meth- o ds [ 78 , 40 , 77 ]. How ever, the lack of regularity of the solution with resp ect to the sto c hastic v ariables imp edes eﬃcient p erformance of b oth the sto chastic Galerkin as w ell as the sto chastic collo cation methods, in particular for high-dimensional param- eter spaces. Here, we propose the dev elopment and implementation of non-intrusiv e Mon te Carlo (MC) methods for UQ of cloud cavitation collapse. In MC metho ds, the gov erning equations are solved for a sequence of randomly generated samples, whic h are combined to ascertain statistical information. Ho wev er, the robustness of MC methods with resp ect to solution regularity comes at the price of a low er- ror con v ergence rate regarding the n um b er of samples. Drawbac ks of the mentioned n umerical uncertaint y quantiﬁcation metho ds inspired the developmen t of v arious m ulti-ﬁdelity metho ds, suc h as Multi-Level Monte Carlo (MLMC) [ 21 ], m ulti-ﬁdelity recursiv e co-kriging and Gaussian–Marko v random ﬁeld [ 56 ]. F urther developmen ts include m ulti-ﬁdelit y Gaussian Process Regression (GPR) based on co-kriging [ 51 ], and purely data-driven algorithms for linear equations using Gaussian pro cess priors to completely circum ven t the use of numerical discretization sc hemes [ 60 ]. MLMC metho ds were introduced b y Heinrich for numerical quadrature [ 29 ], then pioneered b y Giles for Itˆ o SPDEs [ 21 ], and hav e b een lately applied to v arious stochastic PDEs [ 3 , 14 , 23 , 49 ]. The MLMC algorithm was also extended to hyperb olic conserv ation la ws and to massively parallel simulations of the random m ulti-dimensional Euler, magneto-h ydro dynamics (MHD) and shallo w water equations [ 42 , 43 , 44 , 45 , 70 , 67 ]. Subsequen t MLMC impro vemen ts include Ba yesian inference for fusing knowledge on empirical statistical estimates and deterministic con vergence rates [ 15 ], an alternative m ulti-ﬁdelity setting where sample recycling with optimal control v ariate co eﬃcients are emplo yed [ 54 , 55 ], and m ulti-level Marko v Chain Monte Carlo p osterior sampling in Bay esian inference problems [ 19 , 20 ]. The ensem ble uncertain ty of the clouds is parametrized b y means of probabilit y distributions of ca vity radii, p ositions and initial in ternal pressures. Our goal is to p erform sim ulations of cloud ca vitation collapse with unpreceden ted n um b er of inter- acting cavities with full UQ on QoIs (p eak pressure magnitudes, lo cations and cloud collapse times) in terms of means, v ariances, conﬁdence in terv als and even proba- bilit y densit y functions. T o the best of our knowledge, such extensive UQ analysis of (even small) clouds has not b een accomplished b efore. W e propose an improv ed non-in trusive MLMC metho d with nov el optimal control v ariate co eﬃcients coupled with the state-of-the-art n umerical solv er for cloud ca vitation collapse sim ulations and pro vide robust statistical estimates on relev an t quantities of in terest, emphasizing the imp ortance of UQ in such problems. The paper is structured as follo ws: section 2 in tro duces the go verning m ulti-phase equations and the ﬁnite volume metho d used for their solution. It also presents the optimal con trol v ariate MLMC metho d for the statistical sampling of QoIs. Numerical exp erimen ts quan tifying suc h uncertain ties and iden tifying their relations to the geo- metric properties of the cloud b y means of join t probability densit y function estimates are rep orted in section 3 . W e summarize our ﬁndings in section 4 . 4 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS 2. Computational metho ds. The go verning system of equations for multi- phase ﬂows are discretized using a Finite V olume Metho d (FVM) that is eﬃcien tly implemen ted so as to take adv an tage of sup ercomputer architectures. The sampling needed to estimate statistical properties from an ensemble of ev olved QoIs is performed b y a nov el Optimal Fidelity MLMC (OF-MLMC) metho d. W e introduce a Python implemen tation PyMLMC [ 69 ] of OF-MLMC and its embedding in an op en source uncertain ty quan tiﬁcation framework [ 26 ]. 2.1. Gov erning equations. The dynamics of ca vitating and collapsing clouds of bubbles are gov erned by the compressibilit y of the ﬂow, with viscous dissipation and capillary eﬀects taking place at orders of magnitude slo wer time scales. Hence, w e describ e cloud cavitation collapse through the Euler equations for inviscid, compress- ible, multi-phase ﬂows. The system of equations, deriv ed from the Baer-Nunziato mo del [ 2 ], describes the evolution of density , momentum and total energy of the m ulti-phase ﬂow in the domain D ⊂ R 3 as [ 50 , 57 ]: ∂ α 1 ρ 1 ∂ t + ∇ · ( α 1 ρ 1 u ) = 0 , (1) ∂ α 2 ρ 2 ∂ t + ∇ · ( α 2 ρ 2 u ) = 0 , (2) ∂ ( ρ u ) ∂ t + ∇ · ( ρ u ⊗ u + p I ) = 0 , (3) ∂ E ∂ t + ∇ · (( E + p ) u ) = 0 , (4) ∂ α 2 ∂ t + u · ∇ α 2 = K ∇ · u , (5) where (6) K = α 1 α 2 ( ρ 1 c 2 1 − ρ 2 c 2 2 ) α 1 ρ 2 c 2 2 + α 2 ρ 1 c 2 1 . All quantities, unless otherwise stated, depend on spatial v ariable x ∈ D and time v ariable t ≥ 0. This system comprises tw o mass conserv ation equations, one for eac h phase, conserv ation equations for momen tum and total energy in single- or mixture-ﬂuid formulation as w ell as a transp ort equation for the volume fraction of one of the tw o phases with source/sink term on the righ t-hand side. In ( 1 )–( 5 ), u denotes the v elo city , p the pressure, I the iden tit y tensor, ρ the (mixture) density , E the (mixture) total energy E = ρe + 1 / 2 ρ ( u · u ), where e is the (mixture) sp eciﬁc in ternal energy . Moreov er, ρ k , α k and c k with k ∈ { 1 , 2 } are densit y , volume fraction and speed of sound of the tw o phases. F or the mixture quan tities, the following addi- tional relations hold: α 1 + α 2 = 1, ρ = α 1 ρ 1 + α 2 ρ 2 , and ρe = α 1 ρ 1 e 1 + α 2 ρ 2 e 2 . W e do not account for mass transfer betw een diﬀeren t phases (ev ap oration or condensation), so that the ab ov e equations describ e a multi-component ﬂow. The source term in ( 5 ) for homogeneous mixtures [ 35 ] describ es the reduction of the gas volume fraction in a mixture of gas and liquid when a compression wa ve trav els across the mixing region and vice versa for an expansion w av e. F or a more detailed analysis on the p ositive inﬂuence of this term on the accuracy of the mo del equations, we refer to [ 61 ]. The equation system is closed by an appropriate equation of state for eac h of the phases. W e employ the stiﬀened equation of state (see [ 41 ] for a review) to capture liquids and gases. This enables a simple, analytic approximation to arbitrary ﬂuids OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 5 and is expressed by (7) p = ( γ k − 1) ρ k e k − γ k p c ,k , where isobaric closure is assumed [ 57 ]. Parameters γ k and p c ,k dep end on the material. F or p c ,k = 0, the equation of state for ideal gases is reco vered. F or sim ulations in this man uscript, γ 1 = 6 . 59 and p c , 1 = 4 . 049 · 10 8 P a are used for water and γ 2 = 1 . 4 and p c , 2 = 0 P a for air. 2.2. Numerical metho d. The go verning system ( 1 )–( 5 ) can be recast in to the quasi–conserv ative form ∂ Q ∂ t + ∂ F ( Q ) ∂ x + ∂ G ( Q ) ∂ y + ∂ H ( Q ) ∂ z = R ( Q ) , (8) where Q = ( α 1 ρ 1 , α 2 ρ 2 , ρ u , E , α 2 ) T is the v ector of conserv ed (except for α 2 whic h has a non-zero source term) v ariables and F ( Q ), G ( Q ), H ( Q ) are v ectors of ﬂux functions F ( Q ) =           α 1 ρ 1 u x α 2 ρ 2 u x ρu 2 x + p ρu y u x ρu z u x ( E + p ) u x α 2 u x           , G ( Q ) =           α 1 ρ 1 u y α 2 ρ 2 u y ρu x u y ρu 2 y + p ρu z u y ( E + p ) u y α 2 u y           , H ( Q ) =           α 1 ρ 1 u z α 2 ρ 2 u z ρu x u z ρu y u z ρu 2 z + p ( E + p ) u z α 2 u z           . The source term R ( Q ) has all elements equal to zero except the last one R ( Q ) 7 = α 2 ( ∇ · u ) + K ∇ · u , whic h app ears due to rewriting ( 5 ) in conserv ative form [ 33 ] and incorporating the presen t source term. F or the system ( 8 ), initial condition Q ( x , t = 0) = Q 0 ( x ) ov er the entire domain x ∈ D as w ell as b oundary conditions at x ∈ ∂ D for all times t ≥ 0 need to be provided to complete the full sp eciﬁcation of the m ulti-phase ﬂow problem. The metho d of lines is applied to obtain a semi-discrete representation of ( 8 ), where space contin uous op erators are approximated using the FVM for uniform s tructured grids. The approac h yields a system of ordinary diﬀeren tial equations d V ( t ) dt = L ( V ( t )) , (9) where V is a vector of cell a v erage v alues and L ( · ) is a discrete op erator that ap- pro ximates the conv ective ﬂuxes and the given sources in the gov erning system. The temp oral discretization of ( 9 ) is obtained by an explicit third-order lo w-storage Runge- Kutta sc heme [ 76 ]. The computation of the numerical ﬂuxes is based on a Godunov- t yp e scheme using the appro ximate HLLC Riemann solver originally in tro duced for single-phase ﬂow problems in [ 73 ]. The Riemann initial states are determined b y a sho c k capturing third- or ﬁfth-order accurate WENO reconstruction (see [ 32 ]). F ol- lo wing [ 33 ], the reconstruction is carried out using primitiv e v ariables, and the HLLC Riemann solv er is adapted to ( 5 ) to preven t oscillations at interface. The solution is adv anced with a time–step size that satisﬁes the Courant-F riedrichs-Lewy (CFL) condition. F or the co eﬃcient w eights in the Runge-Kutta stages, the v alues suggested in [ 25 ] are used, resulting in a total v ariation diminishing sc heme. 6 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS 2.3. Cubism-MPCF. The FVM used for the spatio-temp oral numerical dis- cretization of non-linear system of conserv ation laws in ( 8 ) is implemented in the op en source softw are Cubism-MPCF [ 64 , 27 , 28 , 61 ]. The applied scheme entails three computational kernels: computation of CFL-limited time-step size ∆ t based on a global reduction, ev aluation of the approximate Riemann problem corresponding to the ev aluation of the righ t-hand side L in ( 9 ) for each time step, and appropri- ate Runge-Kutta up date steps. The solver is parallelized with a h ybrid paradigm using the MPI and Op enMP programming mo dels. The soft ware is split into three abstraction la y ers: cluster, no de, and core [ 27 ]. The realization of the Cartesian do- main decomp osition and the inter-rank MPI communication is accomplished on the cluster lay er. On the no de la y er, the thread lev el parallelism exploits the Op enMP standard using dynamic w ork scheduling. Spatial blo cking tec hniques are used to in- crease lo cality of the data, with in tra-rank ghost cells obtained from loading fractions of the surrounding blocks, and in ter-rank ghost cells obtained from a global buﬀer pro vided by the cluster lay er. On the core lay er, the actual computational k ernels are executed, exploiting data lev el parallelism and instruction lev el parallelism, whic h are enabled b y the con version b etw een the arra y–of–structures and structure–of–arrays la yout. F or the simulations rep orted here, main parts of the computations were ex- ecuted in mixed precision arithmetic. More details on softw are design regarding the parallelization and optimization strategies used in Cubism-MPCF can b e found in [ 27 , 28 , 64 , 61 ]. Cubism-MPCF has demonstrated state–of–the–art p erformance in terms of ﬂoating p oint op erations, memory traﬃc and storage, exhibiting almost p erfect ov erlap of communication and computation [ 27 , 61 ]. The softw are has b een optimized to take adv an tage of the IBM BlueGene/Q (BGQ) and Cra y X C30 plat- forms to sim ulate cavitation collapse dynamics using up to 13 T rillion computational cells with v ery eﬃcient strong and weak scaling up to the full size of MIRA (Argonne National Laboratory) and Piz Dain t (Swiss Sup ercomputing Cen ter) supercomputers [ 64 , 28 , 61 ]. 2.4. Multi-Level Mon te Carlo (MLMC) metho d. In this section, w e intro- duce the MLMC framew ork for UQ. W e also present a nov el and improv ed numerical sampling metho d for approximating the statistical QoI. This study is grounded on the theoretical probabilistic framew ork for non-linear systems of conserv ation la ws introduced in [ 43 , 68 ]. Uncertain ties in the system of conserv ation laws ( 8 ), suc h as uncertain initial data at time t = 0 for the v ector of conserv ed quantities Q , are mo deled as r andom ﬁelds [ 43 , 68 ]. They dep end on the spatial and temp oral v ariables x and t , as well as on the sto chastic parameter ω ∈ Ω, whic h represen ts v ariability in the cloud conﬁguration. F or instance, for uncertain initial data, we assume (10) Q 0 ( x , ω ) = Q ( x , 0 , ω ) ∈ R 7 , x ∈ D, ω ∈ Ω , i.e., at ev ery spatial coordinate x , initial data Q 0 ( x , ω ) is a random v ector con taining 7 v alues, one for each equation in ( 8 ). W e further assume that Q 0 is suc h that at ev ery spatial p oint x the statistical moments such as exp ectation and v ariance exist and are deﬁned by (11) E [ Q 0 ]( x ) = Z Ω Q ( x , ω ) dω and (12) V [ Q 0 ]( x ) = Z Ω  Q 0 ( x , ω ) − E [ Q 0 ]( x )  2 dω . OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 7 Suc h uncertainties, for instance in initial data Q 0 , are propagated according to the dynamics go v erned by ( 8 ). Hence, the resulting evolv ed solution Q ( x , t, ω ) for t > 0 is also a random ﬁeld, called r andom entr opy solution ; see [ 43 , 68 , 42 ] for precise form ulation and details. 2.4.1. The classical MLMC. The statistical moments of the QoIs, suc h as exp ectation E [ q ], are obtained through sampling b y the MLMC metho dology . Multi- lev el metho ds employ a hierarch y of spatial discretizations of the computational do- main D , or, equiv alently , a hierarch y of numerical deterministic solv ers as describ ed in subsection 2.2 , ordered by increasing “level of accuracy” ` = 0 , . . . , L . On eac h suc h discretization lev el ` , and for a giv en statistical realization (a “sample”) ω ∈ Ω, a numerical appro ximation of the QoI q ( ω ) using the applied FVM will b e denoted b y q ` ( ω ). The classical MLMC estimator [ 21 ] provides accurate and eﬃcient estimates for statistical momen ts of q in terms of the telescoping sum of numerical approximations q ` ( ω ) o ver all lev els. In particular, an appro ximation E MLMC [ q L ] of E [ q ] is constructed from the approximate telescoping sum (13) E [ q ] ≈ E [ q L ] = E [ q 0 ] + L X ` =1  E [ q ` ] − E [ q ` − 1 ]  b y additionally appro ximating all exact mathematical expectations using Mon te Carlo sampling with level-dep endent num b er M ` of samples for QoIs q i ` , leading to (14) E [ q L ] ≈ E MLMC [ q L ] = 1 M 0 M 0 X i =1 q i 0 + L X ` =1 1 M ` M ` X i =1  q i ` − q i ` − 1  . W e note that eac h sample q i ` is obtained by solving the gov erning system ( 8 ) using the FVM metho d from subsection 2.2 with discretization (n um b er of mesh cells and time steps) corresp onding to level ` . Assuming that samples for q i 0 and diﬀerences q i ` − q i ` − 1 are drawn independently (ho wev er, sample pairs q i ` and q i ` − 1 ar e statistically dep enden t), the statistical mean square error of the (here assumed to b e unbiased) standard MLMC estimator is given by the sum of sample-reduced v ariances of diﬀer- ences b etw een every t wo consecutive lev els, (15) ε 2 = E h E MLMC [ q L ] − E [ q L ]  2 i = V [ q 0 ] M 0 + L X ` =1 V [ q ` − q ` − 1 ] M ` . Assuming that V [ q ` ] ≈ V [ q ], the MLMC sampling error in ( 15 ) can b e approximated in terms of correlation co eﬃcients of every t wo consecutiv e lev els, i.e., (16) ε 2 ≈ V [ q ] 1 M 0 + 2 L X ` =1 1 − Cor[ q ` , q ` − 1 ] M ` ! . Note, that strongly correlated QoIs on t wo consecutiv e lev els lead to signiﬁcan t reduc- tion in the required num b er of samples on lev els ` > 0. Optimal n um b er of samples M ` for eac h level ` = 0 , . . . L can be obtained using empirical or appro ximate asymptotic estimates on V [ q 0 ] and V [ q ` − q ` − 1 ] by minimizing the amount of total computational w ork M 0 W 0 + · · · + M L W L for a prescrib ed error tolerance τ suc h that ε ≤ τ in ( 15 ), as describ ed in [ 21 ]. Here, W 0 denotes the amount of computational w ork needed 8 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS to compute a single sample (statistical realization) on a given resolution lev el 0. F or lev els ` > 0, W ` denotes the amount of computational w ork needed to compute a pair of such samples on resolution levels ` and ` − 1. Number of samples M ` w as shown to decrease exponentially with increasing level ` , and hence suc h reduction directly translates into large computational sa vings ov er single-lev el MC sampling methods, as rep orted in [ 42 , 43 , 44 , 70 , 67 ]. 2.4.2. Optimal ﬁdelity MLMC: con trol v ariates for t w o-lev el Mon te Carlo. W e presen t a nov el metho d for reducing the v ariance and further increasing the eﬃciency of the classical MLMC metho d. The backbone of MLMC is the hierar- c hical v ariance reduction tec hnique. Assuming only t wo lev els, a coarser lev el ` − 1 and a ﬁner lev el ` , statistical momen ts at level ` use sim ulations from coarser discretization lev el q ` − 1 as a c ontr ol variate with “known” E [ q ` − 1 ] and the pr e determine d coupling co eﬃcien t α ` = 1. The coupled statistical QoI q ∗ ` is given b y: (17) q ∗ ` = q ` + α  E [ q ` − 1 ] − q ` − 1  . The v ariance reduction that is achiev ed in ( 17 ) by replacing q ` with q ∗ ` dep ends on the correlation b etw een q ` and q ` − 1 , (18) V [ q ∗ ` ] = V h q ` + α  E [ q ` − 1 ] − q ` − 1 i = V [ q ` ] + α 2 V [ q ` − 1 ] − 2 α Cov[ q ` , q ` − 1 ] . The expectation E [ q ` − 1 ] is considered “kno wn” since it can be appro ximated b y sam- pling estimator E M ` − 1 [ q ` − 1 ] that is computationally aﬀordable due to the low er reso- lution of the solv er on coarser lev el ` − 1. Statistical estimators using M ` samples are applied to q ∗ ` instead of q ` , leading to the building blo c k of the MLMC estimator: (19) E [ q ` ] = E [ q ∗ ` ] ≈ E M ` ,M ` − 1 [ q ∗ ` ] = E M ` [ q ` ] + α  E M ` − 1 [ q ` − 1 ] − E M ` [ q ` − 1 ]  . In standard MLMC, the co eﬃcient α in ( 17 )–( 19 ) is set to unity [ 21 , 22 ]. This constrain t limits the eﬃciency of the v ariance reduction tec hnique. In particular, assuming that v ariances at b oth levels are comparable, i.e. V [ q ` ] ≈ V [ q ` − 1 ], the standard MLMC estimator fails to reduce the v ariance if the correlation co eﬃcien t Cor[ q ` , q ` − 1 ] drops b elow threshold of 1 / 2, (20) V [ q ∗ ` ] ≈ 2 V [ q ` ] − 2 Cor[ q ` , q ` − 1 ] V [ q ` ] ≥ 2 V [ q ` ] − 2 1 2 V [ q ` ] = V [ q ` ] . Hence, in standard MLMC, suc h moderate lev el correlations Cor[ q ` , q ` − 1 ] ≤ 1 / 2 w ould not pro vide an y v ariance reduction; on the con trary , v ariance w ould b e increased. T o a void this, the optimal α minimizing the v ariance of q ∗ ` as in ( 18 ) can b e used: (21) α = Co v[ q ` , q ` − 1 ] V [ q ` − 1 ] ≈ Cor[ q ` , q ` − 1 ] . A consequence of ( 21 ) is that the pr e determine d c hoice of α = 1 in ( 17 ) is optimal only under very restrictive conditions: p erfectly correlated levels with correlation co eﬃcien t Cor[ q ` , q ` − 1 ] = 1 and assumption that coarser lev el estimates E [ q ` − 1 ] are already av ailable (hence no computation is needed, i.e., W ` − 1 = 0). Note, that for optimal α as in ( 21 ), v ariance is always reduced in ( 18 ), ev en for 0 < Cor[ q ` , q ` − 1 ] < 1 / 2, (22) V [ q ∗ ` ] =  1 − Cor[ q ` , q ` − 1 ] 2  V [ q ` ] < V [ q ` ] . OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 9 F or W ` − 1 6 = 0, it is necessary to obtain an estimate for E [ q ` − 1 ] in ( 17 ) as w ell. In suc h case, using the independence of estimators E M ` and E M ` − 1 and the cen tral limit theorem, the v ariance of the t wo-lev el estimator E M ` ,M ` − 1 [ q ∗ ` ] as in ( 19 ) is giv en by (23) V  E M ` ,M ` − 1 [ q ∗ ` ]  = V [ q ` − αq ` − 1 ] M ` + V [ αq ` − 1 ] M ` − 1 . Giv en a constrain t on the a v ailable computational budget B , the num b er of a v ailable samples on each lev el depends on the required computation cost at that level, i.e., (24) M ` = B W ` . Here, w e ha ve not y et considered the non-uniform distribution of computational bud- get B across levels ` and ` − 1. Such distribution will b e addressed in terms of opti- mal num b er of samples M ` for each level ` in the forth coming subsubsection 2.4.3 . Plugging ( 24 ) to ( 23 ) and multiplying by the computation budget B , the total com- putational cost for v ariance reduction in E M ` ,M ` − 1 [ q ∗ ` ] is given b y (25) C [ q ∗ ` ] = V  E M ` ,M ` − 1 [ q ∗ ` ]  B = V [ αq ` − 1 ] W ` − 1 + V [ q ` − αq ` − 1 ] W ` , where v ariances of αq ` − 1 and q ` − αq ` − 1 are weigh ted by the corresponding computa- tional costs W ` − 1 and W ` , resp ectively . In order to ﬁnd optimal α , the computational v ariance reduction cost C [ q ∗ ` ] in ( 25 ) is minimized instead of V [ q ∗ ` ] in ( 18 ). The re- sulting optimal co eﬃcient α is given b y (26) α = W ` W ` + W ` − 1 Co v[ q ` , q ` − 1 ] V [ q ` − 1 ] ≈ W ` W ` + W ` − 1 Cor[ q ` , q ` − 1 ] . W e note, that ( 26 ) reduces to standard control v ariate co eﬃcient ( 21 ) for W ` − 1 = 0. 2.4.3. Optimal ﬁdelit y MLMC: control v ariates for multi-lev el Monte Carlo. W e generalize the concept of optimal control v ariate co eﬃcients from subsub- section 2.4.2 to arbitrary n umber of levels in MLMC estimator. In particular, con trol v ariate contributions from m ultiple lev els can b e used, (27) q ∗ L =  q L − α L − 1 q L − 1  +  α L − 1 q L − 1 − α L − 2 q L − 2  + · · · +  α 1 q 1 − α 0 q 0  + α 0 q 0 , allo wing for generalization of the telescoping sum in ( 13 ) to the con trol v ariate setting in tro duced in ( 17 ), where the co eﬃcient for the ﬁnest lev el ` = L is ﬁxed to α L = 1: (28) E [ q L ] = E [ q ∗ L ] = α 0 E [ q 0 ] + L X ` =1  α ` E [ q ` ] − α ` − 1 E [ q ` − 1 ]  . Using the m utual statistical indep endence of α 0 q 0 and diﬀerences α ` q ` − α ` − 1 q ` − 1 together with the central limit theorem analogously as in subsubsection 2.4.2 , the total computational v ariance reduction cost, generalizing ( 25 ) for q ∗ L o ver all lev els as in ( 27 ), is given b y (29) C [ q ∗ L ] = V [ α 0 q 0 ] W 0 + L X ` =1 V [ α ` q ` − α ` − 1 q ` − 1 ] W ` = α 2 0 V [ q 0 ] W 0 + L X ` =1  α 2 ` V [ q ` ] + α 2 ` − 1 V [ q ` − 1 ] − 2 α ` α ` − 1 Co v[ q ` , q ` − 1 ]  W ` . 10 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS Minimization of C [ q ∗ L ] p ertains to solving a linear system of equations, (30) ∂ ∂ α ` C [ q ∗ L ] = 0 , ` = 0 , . . . , L − 1 . Denoting σ 2 ` = V [ q ` ], W `,` − 1 = W ` + W ` − 1 , and σ 2 `,` − 1 = Cov[ q ` , q ` − 1 ], system ( 30 ) can b e written in a form of a tridiagonal matrix (31)       σ 2 0 W 1 , 0 − σ 2 1 , 0 W 1 − σ 2 1 , 0 W 1 . . . . . . . . . . . . − σ 2 L − 1 ,L − 2 W L − 1 − σ 2 L − 1 ,L − 2 W L − 1 σ 2 L − 1 W L,L − 1              α 0 α 1 . . . α L − 2 α L − 1        =        0 0 . . . 0 σ 2 L,L − 1 W L        . F or a special case of tw o levels with L = 1, the solution of ( 31 ) reduces to ( 26 ). Gen- eralizing ( 14 ), the unbiased Optimal Fidelity Multi-L evel Monte Carlo (OF-MLMC) estimator is then given b y (32) E [ q L ] ≈ E OF-MLMC [ q L ] = α 0 1 M 0 M 0 X i =1 q i 0 + L X ` =1 1 M ` M ` X i =1  α ` q i ` − α ` − 1 q i ` − 1  , with optimal con trol v ariate co eﬃcients α ` giv en b y ( 31 ). W e assume that samples for α 0 q i 0 and diﬀerences α ` q i ` − α ` − 1 q i ` − 1 are drawn indep endently (how ev er, individual q i ` and q i ` − 1 ar e statistically dep endent). The statistical mean square error ε 2 OF of the OF-MLMC estimator is given by the sum of sample-reduced v ariances of α ` -w eighted diﬀerences b etw een every t wo consecutive lev els, (33) ε 2 OF = E h E OF-MLMC [ q L ] − E [ q L ]  2 i = V [ α 0 q 0 ] M 0 + L X ` =1 V [ α ` q ` − α ` − 1 q ` − 1 ] M ` =: ˜ σ 2 0 M 0 + L X ` =1 ˜ σ 2 ` M ` . Giv en computational costs W ` of ev aluating a single approximation q i ` for eac h level ` = 0 , . . . , L and a desired tolerance τ > 0, the total computational cost of OF-MLMC can b e minimized under constraint of ε 2 OF ≤ τ by choosing optimal n umber of MC samples M ` on each lev el according to [ 48 ], yielding (34) M ` =     1 τ 2 s ˜ σ 2 ` W ` L X k =0 q ˜ σ 2 k W k     . Alternativ ely , given av ailable total computational budget B instead of a desired tol- erance τ , the OF-MLMC error ε 2 OF is minimized by c ho osing optimal M ` according to (35) M ` =     B s ˜ σ 2 ` W ` . L X k =0 q ˜ σ 2 k W k     . OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 11 2.4.4. Adaptive optimal ﬁdelity MLMC algorithm. The OF-MLMC algo- rithm proceeds iterativ ely , with eac h iteration impro ving the accuracy of the estimated statistics for the quantities of interest, such as the exp ectation E [ q ]. Eac h iteration also improv es the accuracy of auxiliary parameters, such as σ 2 ` , σ 2 `,` − 1 , W ` , ˜ σ 2 ` , and the optimality of the n um b er of samples M ` for each lev el ` = 0 , . . . , L . A single iteration of the algorithm consists of the following 8 steps. (1) W arm-up: Begin with level-dep endent num b er of warm-up samples, (36) M ` =  W L W ` 1 2 ( L − ` )  , ` = L, L − 1 , . . . , 0 . The choice of M ` as in ( 36 ) preven ts eﬃciency deterioration of the ﬁnal OF- MLCM estimator b y ensuring that the total computational budget for the w arm-up iteration do es not exceed 2 W L ; at the same time, it allows to pre- scrib e suﬃciently many samples on the coarser levels, where ˜ σ 2 ` is exp ected to b e large. In our particular application, computational work for each level scales as W ` = O (2 4 ` ), and hence the amount of required warm-up samples is given by 1 , 8 , . . . , 2 3 L for ` = L, L − 1 , . . . , 0. W e note, that constan t (lev el- indep enden t) num b er of w arm-up samples can b e very ineﬃcient [ 48 , 15 ]. (2) Solver: Ev aluate approximations q i ` for eac h level ` = 0 , . . . , L and sample i = 1 , . . . , M ` . (3) Indicators: Using q i ` , estimate σ 2 ` , σ 2 `,` − 1 , and W ` for ` = 0 , . . . , L . Op- tionally , empirical estimates of σ 2 `,` − 1 could b e used within Ba yesian infer- ence framework to ﬁt an exp onential decay mo del for σ 2 `,` − 1 w.r.t. levels ` = 1 , . . . , L − 1. Assuming Gaussian errors, this reduces to a least-squares line ﬁt to the natural logarithm of indicators σ 2 `,` − 1 . (4) Co eﬃcients: Compute con trol v ariate co eﬃcients α ` from estimated σ 2 ` and σ 2 `,` − 1 using ( 31 ). (5) Errors: Using q i ` and α ` , estimate α ` -w eighted cov ariances ˜ σ 2 ` and total sampling error ˆ ε 2 OF ≈ ε 2 OF as in ( 33 ). (6) Estimator: If the required tolerance is reached, i.e., ˆ ε OF ≤ τ , or if the pre- scrib ed computational budget B is sp ent, then ev aluate OF-MLMC estimator ( 32 ) and terminate the algorithm. Otherwise, pro ceed to the optimization step. (7) Optimization: Compute optimal required n umber of samples ˆ M ` from ˜ σ ` and W ` using either ( 34 ) or ( 35 ), resp ectiv ely . R emark. If we obtain ˆ M ` < M ` for some level ` , w e ke ep the already computed samples, i.e., we set ˆ M ` = M ` . In order to adjust for such a constraint in the optimization problem, w e subtract the corresp onding sampling error ˜ σ 2 ` / M ` from the required tolerance τ 2 , or subtract the corresp onding amoun t of computational budget M ` W ` from B , resp ectively . Afterwards, the n umber of samples ˆ M ` for the r emaining levels (where ˆ M ` = M ` w as not enforced) are re-optimized according to either ( 34 ) or ( 35 ), resp ectively . W e rep eat this pro cedure until ˆ M ` ≥ M ` is satisﬁed for all levels ` = 0 , . . . , L . (8) Iterations: Go bac k to step (2) and con tinue the algorithm with the up dated n umber of samples ˆ M ` . If the empirical estimates in steps (3)–(5) of the ab ov e adaptiv e OF-MLMC algorithm are suﬃciently accurate, it will terminate after tw o iterations - the initial w arm-up it- eration and one additional iteration with already optimal α ` and M ` . A more detailed discussion of the sp ecial cases within OF-MLMC algorithm, alternative approaches, 12 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS related work and possible extensions is provided in Appendix A . 2.5. PyMLMC. The OF-MLMC algorithm is distributed as an op en source li- brary PyMLMC [ 69 ]. A diagram of the soft ware comp onents is sho wn in Figure 3 . PyMLMC provides a mo dular framew ork for nativ e execution of deterministic solvers in their respective “sandb ox” directories. This allo ws maxim um ﬂexibilit y for pro- gramming languages, distributed and/or shared memory architectures, and supp ort for many-core accelerators. Due to the lac k of comm unication among suc h sandb oxed executions for each sample, the load balancing across samples can b e relay ed to the submission system (e.g. Slurm, LSF, LoadLev eler, Cobalt) of the compute cluster. Nested (across and within samples) parallelism is used, where few large parallel jobs are submitted for ﬁne levels, and, in parallel, many small (p ossibly even serial) jobs are submitted for coarse levels. T o increase the eﬃciency and reduce the stress on submission systems, job batching (submitting a single job that computes multiple samples subsequently) and job merging (submitting a single job that computes m ulti- ple samples concurren tly) or a com bination of b oth is implemented. Once all samples (or at least some of them) are computed, statistical estimators are constructed as a p ost-pro cessing step using the NumPy and SciPy libraries. The “sandboxing” frame- w ork enables an y additional required statistical estimates or QoIs to be ev aluated at a later stage without the need to re-execute an y of the computationally exp ensive sam- pling runs. The amount of required disk space in m ulti-lev el metho ds scales linearly w.r.t. the amount of required computational budget. In particular, the total required disk space for all samples on all lev els is of the same order as a single statistical esti- mator on the full three-dimensional domain. Hence, it is highly adv an tageous to keep all computed samples for the aforemen tioned p ost-pro cessing ﬂexibility purposes. solver MPCF wrapper another solver samples scheduler statistics indicators errors Cubism- MPCF results ultimate goal static TORC conﬁdence intervals mean histogram other wrapper variances runtimes optimization median standard deviations correlations coefﬁcients KDE ??? Fig. 3 . Scheme of the OF-MLMC implementation libr ary PyMLMC. In the present work (see section 3 ), w e veriﬁed the eﬃciency of the (nested) parallelization of the OF-MLMC coupled with the Cubism MPCF solv er on the entire MIRA sup ercomputer (Argonne National Lab oratory) consisting of half a million cores. W e note that large (exa-)scale simulations on massively parallel computing platforms are sub ject to pro cessor failures at run-time [ 11 ]. Exploiting the natural fault tolerance in OF-MLMC-FVM due to indep endent sampling, a F ault T oleran t MLMC (FT-MLMC) metho d was implemented in [ 53 ] and was shown to p erform in agreemen t with theoretical analysis in the presence of simulated, comp ound P oisson OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 13 distributed, random hard failures of compute cores. Such FT-mechanisms are also a v ailable in PyMLMC, and hav e successfully excluded one corrupted sample on the coarsest level in the sim ulations rep orted in section 3 . 3. Numerical simulations and analysis. The initialization of the ca vities emplo ys exp erimental ﬁndings indicating log-normal distribution for their radii [ 47 ], whereas the p osition vectors are generated according to a uniform distribution as there is no prior knowledge. 3.1. Spherical cloud of 500 gas cavities with log-normally distributed radii. F or a cubic domain D = [0 mm, 100 mm ] 3 , we consider a cloud of 500 bub- bles lo cated at the cen ter (50 mm, 50 mm, 50 mm ) > of the domain with a radius of R cloud = 20 mm . The log-normal distribution for the radii of the cavit y is clipp ed so as to contain bubbles only within the range of r min = 0 . 8 mm to r max = 1 . 2 mm . The resulting cloud gas volume conten t (w.r.t. to the v olume of the sphere with ra- dius R cloud ) is approximately 5% and the resulting cloud interaction parameter β is appro ximately 3, where β = α ( R cloud R avg ) 2 with cloud gas v olume fraction α , cloud radius R cloud and a verage cavit y radius R avg (refer to [ 9 ] for a deriv ation). W e note that b oth of these quantities dep end on a statistical realization of the random cloud. An illustration of the cloud geometry is shown in Figure 4 . The cloud is initially at pressure equilibrium with the surrounding water phase at p 2 , 0 = 0 . 5 M P a . Throughout the en tire domain, the density of the gas phase is set to ρ 2 , 0 = 5 kg /m 3 and the densit y of the liquid is set to ρ 1 , 0 = 1000 kg /m 3 . Starting 1 mm aw ay from the surface of the cloud, there is a smo oth pressure increase to wards the prescrib ed am bient pressure of p ∞ = 10 M P a , follo wing the setup prop osed in [ 71 ]. The resulting pressure ratio is p ∞ /p 2 , 0 = 20. A t the b oundaries, non-reﬂecting, c haracteristic-based b oundary conditions are applied, together with a p e nalt y term for the prescribed far-ﬁeld pressure of p ∞ = 10 M P a [ 59 ]. A single statistical realization (sample) of the pressure ﬁeld and cavit y in terfaces across the tw o-dimensional slice at z = 50 mm computed using Cubism-MPCF with resolution of 1024 3 mesh cells is depicted in Figure 4 . Ca vities are observed to collapse in wards, emitting pressure w av es that fo cus near the center of the cloud at time t = 46 . 4 µs . 30 40 50 60 70 x [ m m ] 30 40 50 60 70 y [ m m ] v a p o r f r a c t i o n s l i c e | 0 . 0 0 µ s | L 2 T 0 S 0 0 20 40 60 80 100 120 140 160 180 200 30 40 50 60 70 x [ m m ] 30 40 50 60 70 y [ m m ] v a p o r f r a c t i o n s l i c e | 4 6 . 4 0 µ s | L 2 T 0 S 0 0 50 100 150 200 250 300 350 400 450 500 Fig. 4 . Single statistic al r ealization (sample) of the pressur e ﬁeld and c avity interfac es in the two-dimensional slic e at z = 50 mm at the initial (left, t = 0 µs ) and at the collapse (right, t = 46 . 4 µs ) stages. Cavities c ol lapse inwards, emitting pr essur e waves that fo cus ne ar the c enter of the cloud. Ther e is a two or ders of magnitude diﬀer enc e in pr essur es at initial and c ol lapse stages. 14 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS W e consider four levels of spatial and temp oral resolutions. The coarsest mesh consists of 512 3 cells with t wo intermediate meshes of 1024 3 and 2048 3 resolutions, and the ﬁnest mesh with 4096 3 cells. The time-step size decreases according to a prescrib ed CFL condition with CFL num b er set to 0 . 3, resulting in approximately 2 0 250 and 20 0 000 time steps for the coarsest and ﬁnest levels, respectively . W e note that the n um b er of uncertain ty sources in this simulation is very large: for each realization of a cloud, random three-dimensional spatial co ordinates together with a random positive radius for all 500 ca vities are needed, comprising in total 2 0 000 indep enden t sources of uncertaint y . F or eac h statistical sample of a collapsing cloud conﬁguration and on eac h res- olution level, simulations were p erformed for approximately 70 µs in ph ysical time. Dep ending on the random conﬁguration of the cloud, the main collapse o ccurred at appro ximately 40 − 50 µs , follow ed by reb ound and relaxation stages after 50 µs . The obtained results are discussed in the following sections. 3.2. Performance of the OF-MLMC. W e quantify the computational gains of the OF-MLMC method b y comparing it to standard MLMC and plain MC sampling metho ds. The chosen quan tity of interest q for ( 28 ) is the pressure as sampled by a sensor p c placed at the center of the cloud and emulated as (37) p c ( t ) = 1 |B r ( c ) | Z B r ( c ) p ( x , t ) d x , where pressure p is a veraged o ver a sphere B r ( c ) around the cen ter of the cloud located at c = (50 mm, 50 mm, 50 mm ) > with radius r = 0 . 5 mm , (38) B r ( c ) = { x ∈ D : k x − c k 2 ≤ r } . F or this particular c hoice of the QoI q = p c , estimated correlations b et ween lev els, implicitly used in ( 31 ), and the resulting optimal con trol v ariate co eﬃcients from ( 31 ) are depicted in Figure 5 . Due to relatively high correlations b etw een resolution levels, optimal con trol v ariate coeﬃcients exhibit only moderate deviations from unity , with the largest b eing at level 1 with deviation of 30%. 0 1 2 3 mesh level 0.0 0.2 0.4 0.6 0.8 1.0 c o r r e l a t i o n o f Q ` a n d Q ` − 1 0 1 2 3 mesh level 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 coefficient Fig. 5 . Estimate d c orr elations b etwe en levels (left) and the resulting optimal c ontr ol variate c o eﬃcients (right) as in ( 31 ) . In left plot, cir cles depict me asur ements of corr elations with the asso ciate d err or b ars, and solid line depicts the r esulting inferr e d values. Relatively high c orr elations and mo der ate deviations (up to 30%) in optimal c ontrol variate c o eﬃcients ar e observe d. Estimated v ariances of lev el diﬀerences, required in ( 31 ), and sampling errors for eac h level, computed in ( 33 ), are depicted in Figure 6 . V ariances of diﬀerences are OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 15 decreasing for ﬁner levels of resolution, whic h requires a smaller num b er of MC sam- ples M ` in order to reduce the statistical error on ﬁner resolution levels, where sam- pling entails computationally v ery exp ensive n umerical simulations. Measuremen ts of v ariances of diﬀerences σ 2 `,` − 1 are plotted as circles, with the asso ciated error bars, estimated from the v ariance of the estimator and the num b er of samples used in the w arm-up stage. These measurements, together with assumed Gaussian error mo del from the error bars, are used within Ba yesian inference framework to ﬁt an exp onen- tial decay mo del for σ 2 `,` − 1 w.r.t. “diﬀerence” levels ` = 1 , 2 , 3. Solid line depicts the resulting inferred v alues, whic h are later used in ( 31 ), whereas dashed line depicts the adjusted v alues ˜ σ 2 `,` − 1 from ( 33 ), where optimal con trol v ariate co eﬃcients are ap- plied. The resulting statistical errors at b oth MLMC iterations (warm-up and ﬁnal) are decreasing w.r.t. the increasing resolution level. At the ﬁnal MLMC iteration, the errors are signiﬁcantly decreased on all levels when compared to the warm-up iteration, resulting in the total statistical error estimate of approximately 1 . 3 · 10 − 2 . 0 1 2 3 mesh level 1 0 - 2 1 0 - 1 1 0 0 1 0 1 | r e l . m e a n o f α ` Q ` − α ` − 1 Q ` − 1 | 0 1 2 3 mesh level 1 0 - 4 1 0 - 3 1 0 - 2 1 0 - 1 1 0 0 r e l . s t d . d e v . o f α ` Q ` − α ` − 1 Q ` − 1 0 1 2 3 mesh level 1 0 - 4 1 0 - 3 1 0 - 2 1 0 - 1 q V a r ( α ` Q ` − α ` − 1 Q ` − 1 ) / M ` Fig. 6 . Estimate d variances of diﬀerenc es σ 2 `,` − 1 (left) re quir e d in ( 31 ) and the r esulting es- timate d sampling errors ε ` (right) deﬁne d in ( 33 ) on e ach level. In left plot, cir cles depict me a- sur ements σ 2 `,` − 1 with the asso ciated error b ars, solid line depicts the resulting inferr e d values, and dashe d line depicts the adjuste d values ˜ σ 2 `,` − 1 fr om ( 33 ) , wher e optimal c ontr ol variate co eﬃ- cients ar e applie d. In right plot, br own line indic ates estimate d err ors on e ach level after the initial MLMC warmup iter ation, wher e as or ange line and the c orr esp onding r e gion indic ate the estimate d err or impr ovements after the ﬁnal MLMC iter ation. The ﬁnal total statistic al err or is estimated to b e appr oximately 1 . 3 · 10 − 2 . F or the optimization of samples num b ers M ` on each lev el, a prescrib ed budget of 16 million core hours was used and the optimal n umber of samples w as determined b y ( 35 ). Estimated num b er of samples M ` for the w arm-up and ﬁnal iterations of the OF-MLMC algorithm are depicted in Figure 7 , together with the resulting estimated computational budget on eac h lev el. W e observe that a signiﬁcantly larger num b er of samples is required on the coarser levels of resolution owing to a strong reduction in lev el diﬀerence v ariances ˜ σ 2 `,` − 1 , which are also highest at the coarsest resolution lev els. How ever, the required computational budget is comparable across all lev els; suc h m ulti-level approac h achiev es a signiﬁcan t (more than t w o orders of magnitude) reduction in statistical error (i.e., in the v ariance of the statistical estimators), while at the same time k eeping the deterministic error (bias) small, which is controlled b y the resolution of the ﬁnest lev el. In T able 1 , a comparison of MC, MLMC, and OF-MLMC metho ds and estimated computational sp eedups is provided for a target statistical error of 1 . 3 · 10 − 2 . Note, that the n umber of samples on the same resolution as the ﬁnest lev el in (OF-)MLMC 16 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS 0 1 2 3 mesh level 1 0 0 1 0 1 1 0 2 1 0 3 1 0 4 number of samples 0 1 2 3 mesh level 0 1 2 3 4 5 6 budget [million core hours] Fig. 7 . Estimate d numb er of samples M ` for the warm-up and ﬁnal iter ations of the OF-MLMC algorithm (left) and the r esulting estimate d c omputational budget (right) on e ach level. Br own lines and the c orr esp onding r e gions indic ate numb er of samples an computational budget during the initial MLMC warmup iteration, wher eas gr e en lines indic ate the total ac counts of b oth quantities at the ﬁnal MLMC iter ation; gr e en r e gions indic ating the diﬀer ences b etwe en these two iterations. We observe that despite thousand of samples r e quir e d on c oarser levels of r esolution, the r equir e d c orr esp onding c omputational budget is c omp arable among al l levels. and resulting computational budget required for MC simulations are estimated b y (39) M MC =  σ L ε OF  , W MC = M MC W L . OF-MLMC is estimated to be more than t wo orders of magnitude faster than the plain MC metho d, and even more than three times faster than standard MLMC method without optimized control v ariate co eﬃcients. The ov erall computational budget of OF-MLMC was only appro ximately 8 times larger than a single FVM solv e on the ﬁnest resolution level. { M ` } ` =0 ,...,L budget B error ε sp eedup OF-MLMC 6400, 384, 40, 2 16.6 M CPU hours 1 . 3 · 10 − 2 176.8 MLMC 4352, 258, 32, 3 ∼ 50 M CPU hours 1 . 3 · 10 − 2 50.6 MC ∼ 2 million ∼ 2 B CPU hours 1 . 3 · 10 − 2 - T a ble 1 Comp arison of MC, MLMC, and OF-MLMC metho ds and estimate d computational sp e e dups over standar d MC. OF-MLMC is estimated to b e mor e than two orders of magnitude faster than plain MC, and even mor e than thr e e times faster than standar d MLMC. 3.3. Statistics for temporal quan tities of interest. Multi-level Mon te Carlo statistical estimates are depicted in Figure 8 , Figure 9 and Figure 10 . The statistical spread of the maximum (in physical domain) pressure is esp ecially wide at its p eak (in time) v alue, implying a large uncertaint y in the damage p otential of the ca vita- tion collapse. T o the b est of our kno wledge, such full Probability Density F unctions (PDFs) are rep orted here for the ﬁrst time when using the MLMC metho dology for non-linear systems of conserv ations laws. T o obtain such estimates, level-dependent k ernel density estimators were used, with maximum bandwidth determined using a w ell-known “solve-the-equation” metho d [ 66 , 63 ]. Such empirical PDFs are signiﬁ- can tly more v aluable in engineering, compared to less informative mean/median and deviations/p ercen tiles estimators, especially for bi-modal probabilit y distributions of- ten encoun tered in suc h non-linear systems of conserv ations la ws due to the presence of sho cks [ 43 ]. OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 17 0 10 20 30 40 50 60 70 t i m e [ µ s ] 0.00 0.01 0.02 0.03 0.04 0.05 v a p o r f r a c t i o n s e n s o r 5 [ − ] 0 10 20 30 40 50 60 70 t i m e [ µ s ] 0 1000 2000 3000 4000 5000 6000 7000 m a x p r e s s u r e [ M P a ] 1 0 - 3 1 0 - 2 1 0 - 1 1 0 0 probability density Fig. 8 . Unc ertainties in the cloud gas volume (mean values with 50% and 90% c onﬁdenc e intervals, left) and glob al maximum pr essur e (scale d empirical histo gram, right) within the cloud during the c ol lapse. Sinc e al l initial cloud conﬁgur ations c ontain the same numb er of e qual ly-size d c avities, very low unc ertainties ar e observed in the evolution of the total gas volume. However, the statistic al spr e ad of the p e ak pr essur e is esp e cial ly wide at its maximum value, implying a large unc ertainty in the damage p otential of the c avitation collapse. In Figure 8 , uncertain ties in the cloud gas v olume (represen ted by gas fraction sensor #5, lo cated at the center c = (50 mm, 50 mm, 50 mm ) > with 20 mm radius (hence containing the en tire cloud) and global maximum pressure within the cloud are measured during the entire collapse of 70 µs . As all initial cloud conﬁgurations con tain the same n um b er of equally-sized cavities, v ery low uncertain ties are observed in the evolution of the total cloud gas volume. How ever, the statistical spread of the p eak pressure is esp ecially wide at its maximum v alue, indicating a strong necessity for uncertaint y quantiﬁcation in suc h complex multi-phase ﬂo ws. 0 10 20 30 40 50 60 70 t i m e [ µ s ] 0.0 0.2 0.4 0.6 0.8 1.0 v a p o r f r a c t i o n s e n s o r 4 [ − ] 0 10 20 30 40 50 60 70 t i m e [ µ s ] 0 100 200 300 400 500 600 700 p r e s s u r e s e n s o r 4 [ M P a ] Fig. 9 . Unc ertainties (me an values with 50% and 90% c onﬁdenc e intervals) in the gas volume (left) and pr essur e (right) sensor at the c enter of the cloud during the collapse. Notic e, that the statistic al spr ead of the p e ak sensor pr essur e is esp ecially wide at its maximum value and the p ost- c ol lapse incr e ase in the gas fraction during the r ebound stage. In Figure 9 , uncertain ties are measured in the gas v olume fraction sensor α c and pressure sensor p c at the center of the cloud as deﬁned in ( 37 ). Similarly as for the observ ations of p eak pressure b ehavior, the statistical spread of the p eak sensor pressure is esp ecially wide at its maxim um v alue. The post-collapse increase in the gas fraction indicates the presence of a reb ound stage. During this stage, the post-collapse gas fraction consistently (in terms of conﬁdence in terv als) exceeds pre-collapse v alues, indicating the presence of induced outgoing rarefaction wa ves. In Figure 10 , uncertainties in the maxim um sp eed of sound and the p eak pressure lo cation distance from the cen ter of the cloud are measured during the entire collapse 18 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS 0 10 20 30 40 50 60 70 t i m e [ µ s ] 0 1 2 3 4 5 6 m a x s p e e d o f s o u n d [ m m / µ s ] 1 0 - 3 1 0 - 2 1 0 - 1 1 0 0 1 0 1 probability density 0 10 20 30 40 50 60 70 t i m e [ µ s ] 0 5 10 15 20 d i s t . o f m a x p r e s s u r e [ m m ] 1 0 - 3 1 0 - 2 1 0 - 1 1 0 0 1 0 1 probability density Fig. 10 . Unc ertainties (pr ob ability density functions) in the maximum sp e e d of sound (left) and p eak pr essur e lo c ation (right) during the c ol lapse. Dashe d line indic ates the initial lo c ation of the cloud surfac e. of 70 µs . The uncertain ty in the maxim um sp eed of sound is a direct consequence of large uncertainties in the global peak pressure. How ever, on the contrary , the uncertain ty in the distance of the p eak pressure lo cation from the cloud center is m uch smaller, i.e., the temp oral-spatial proﬁle of the pressure wa ve evolution as it tra vels from the surface of the cloud to wards the cen ter has a m uch low er uncertain ty (when compared to the large observed uncertain ties in the global maxim um pressure estimates). 3.4. Statistics for spatial quantities of in terest. In this section, w e plot the statistical estimates of QoIs extracted along one-dimensional lines that are cut as a straigh t line through the center of the cloud in the three-dimensional computational domain. W e note that radial symmetry is assumed in the statistical distribution of ca vities within the cloud, and hence suc h one-dimensional statistical estimates through the center of the cloud are suﬃcient to represen t the global uncertain ty in the entire three-dimensional domain. The ob jective of extracted one-dimensional line plots is to provide a b etter insight into uncertaint y structures at the center of the cloud by plotting all statistical estimates in a single plot. The line is cut at a sp eciﬁc physical sim ulation time, when the peak pressure is observ ed, and hence is slightly diﬀeren t for ev ery sample. In order to reduce v olatility in global maxim um pressure measuremen ts and hence the choice of the collapse time, we smo othen the observed p eak pressure measuremen ts with a Gaussian kernel of width 0 . 5 mm by the means of fast F ourier transform. Statistical estimates for such extracted one-dimensional lines for pressure at diﬀerent stages of the collapse pro cess are provided in Figure 11 . The uncertain ties are estimated using the MLMC metho dology in the extracted pressure along the line in x -direction (with co ordinates y = 50 mm and z = 50 mm ﬁxed) at the pre-collapse time t = 43 µs and at the time of largest p eak pressure, whic h o ccurs approximately at 46 µs . The time of largest p eak pressure dep ends on the initial c ould conﬁguration and hence is a random v ariable, v arying in eac h statis- tical realization. W e observe that the resulting uncertaint y in encountered pressures increases signiﬁcan tly at the ﬁnal collapse stage, and largest spreads are observed near the epicenter of the cloud ca vitation collapse, where the damage potential is the highest. 3.5. Analysis of linear and non-linear dep endencies. Statistical estimates rep orted in the previous sections indicate that even though the initial cloud setup is very similar for all realizations, i.e., equal coun t of cavit y bubbles, identical cloud OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 19 0 20 40 60 80 100 x [ m m ] 0 100 200 300 400 500 p r e s s u r e [ M P a ] 1 0 - 3 1 0 - 2 1 0 - 1 probability density 0 20 40 60 80 100 x [ m m ] 0 100 200 300 400 500 p r e s s u r e [ M P a ] 1 0 - 3 1 0 - 2 1 0 - 1 probability density Fig. 11 . MLMC-estimated unc ertainties (PDFs) of the extr acte d pr essur e line at the pre- c ol lapse time ( t = 43 µs , left) and at the time of lar gest p e ak pr essure ( ∼ 46 µs , right). The r esulting unc ertainty in enc ountere d pr essur es incr eases signiﬁc antly at the ﬁnal collapse stage, with lar gest spr e ads observe d ne ar the epic enter of the cloud c avitation c ol lapse. radius and ca vity radii ranges (whic h resulted in v ery small uncertain ties for the cloud v olume reported in Figure 8 ), and equal initial gas and liquid pressures, the resulting p eak pressure uncertaint y is v ery large, as seen in Figure 9 and Figure 11 . Hence, it is only the actual conﬁguration (or distribution) of the cavities inside the cloud that can ha ve suc h amplifying (or attenuating) eﬀect on the peak pressures. The main scop e of this section is to inv estigate v arious quan tities of in terest that could p otentially explain the cause of such non-linear eﬀects. The set of selected candidate metrics for the cloud conﬁguration includes skewness (asymmetry) of the initial spatial ca vity distribution inside the cloud, cloud in teraction parameter β , and distance (from the cen ter of the cloud) of the central ca vity (i.e. the cavit y closest to the cen ter). Cloud skewness is a measure of asymmetry of the cloud and is estimate d b y a statistical centered third moment of the distribution of cavit y lo cations along eac h of the three spatial axes. All quantities from this set of candidate metrics are tested for linear statistical correlations with the observed v alues of p eak pressure, p eak pressure distance from the cen ter of the cloud, p eak pressure at the sensor at the cen ter of the domain, and collapse time (when largest p eak pressure o ccurs). In addition, w e hav e also tested for statistical correlations among QoIs themselv es, such as p eak pressure lo cation and the location of the cen ter-most cavit y in the cloud. The results are provided as a Hin ton diagram in Figure 12 . W e observe multiple signiﬁcant direct and inv erse linear statistical correlations b et ween candidate metrics and QoIs: • mild inverse correlation b etw een cloud skewness and pr essur e sensor read- ings, mainly a consequence of the cen tral sensor placement within the cloud; • strong correlation b etw een the lo c ation of c entr al c avity and lo c ation of p e ak pr essur e (w.r.t. cloud center), conﬁrming prior observ ations in [ 62 ] that p eak pressures in the cloud are observed within cavities that are near the center of the cloud; • strong inverse correlations b etw een p e ak pr essur e lo c ation and p e ak pr essur e magnitude , indicating that highest p eak pressures are most probable near the cen ter of a cloud; • mo derate correlation b et w een β and c ol lapse time , since large β v alues can b e a consequence of large gas fraction in the cloud. Despite numerous correlations explaining the statistical spread of observ ed pres- sures, the inﬂuence of cloud interaction parameter β remains undetermined. T o this 20 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS [1] [2] [3] [4] [5] [6] [7] cloud skewness [1] beta [2] central cavity distance [3] peak preassure [4] pressure at sensor [5] peak pressure distance [6] collapse time [7] [1] [2] [3] [4] [5] [6] [7] cloud skewness [1] beta [2] central cavity distance [3] peak preassure [4] pressure at sensor [5] peak pressure distance [6] collapse time [7] Fig. 12 . Hinton diagr am indic ating pair-wise line ar statistical c orr elations b etween c andidate metrics and sele cte d quantities of inter est. White squar es indic ate p ositive and black squar es indic ate ne gative c orr elations. The size of the squar e indicates the magnitude of the corr elation in the interval [0 , 1] . Cavities within the entire cloud and only in the cor e of the cloud ar e c onsider e d in the left and right plots, r esp e ctively. end, we consider cloud skewness and β only for the c or e of the cloud. W e hav e iden tiﬁed the core of the cloud to b e a region around the cen ter of the cloud where uncertain ties in p eak pressure are largest, resulting in a spherical core with radius of 10 mm , based on results in Figure 11 . In this case, correlations inv olving resp ective metrics such as cloud sk ewness and β for the core of the cloud, are observ ed to be more signiﬁcant: • mild direct correlation b etw een β and pr essur e sensor readings, indicating stronger collapse for clouds with higher cloud interaction parameter β due to stronger pressure ampliﬁcation; • mild inv erse correlation b etw een cloud sk ewness and β . Ov erall, such insigh t into statistical correlations provides a very informativ e descrip- tion of inter-dependencies b etw een cloud prop erties and observed collapse pressures, suggesting direct causal relations w.r.t. cloud non-skewness, interaction parameter β , and proximit y of the cen tral cavit y to the cloud cen ter. Due to a non-linear nature of the multi-phase equations and the collapse process, non-linear dep endencies among candidate metrics and QoIs could b e present, whic h migh t b e inaccurately captured by the estimates of linear statistical correlations in Figure 12 . W e tested this h yp othesis b y estimating the full joint probability distribu- tion for the pairs of signiﬁcantly correlated candidate metrics and pressure b ehavior observ ations. In Figure 13 and Figure 14 , w e pro vide results for a selected subset of tested correlation pairs, where strongest and most relev ant correlations w ere observ ed. Join t probability distributions are consistent with linear correlation estimates obtained in Figure 12 , and additionally pro vide v aluable insight into non-linear de- p endencies among QoIs. Obtained results provide a go o d global o verview of causal links betw een cloud structure and collapse pressures and motiv ate further analysis to determine the complex mechanisms gov erning the dynamics of suc h large and com- plex cloud ca vitation phenomena. W e w ould also like to refer to an ongoing extensive (deterministic) parameter study [ 62 ] which inv estigates such causal links for a wider range of cloud sizes, cavit y counts and cloud in teraction parameter β v alues. OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 21 4 5 6 7 8 9 b e t a [ − ] 43 44 45 46 47 48 49 t i m e [ µ s ] 1 0 - 2 1 0 - 1 probability density 4 5 6 7 8 9 b e t a [ − ] 0 100 200 300 400 500 600 p r e s s u r e s e n s o r 4 [ M P a ] 1 0 - 4 1 0 - 3 probability density Fig. 13 . Joint PDFs of cloud inter action p ar ameter β with c ol lapse time (left) and the r esulting pr essur e sensor re ading (right). Higher values of cloud interaction p ar ameter β are mor e likely to c ause lar ger and mor e delaye d collapse pr essur es. 0 1 2 3 4 5 d i s t . o f c a v i t y [ m m ] 0 2 4 6 8 10 d i s t . o f m a x p r e s s u r e [ m m ] 1 0 - 3 1 0 - 2 1 0 - 1 probability density 0 2 4 6 8 10 d i s t . o f m a x p r e s s u r e [ m m ] 0 1000 2000 3000 4000 5000 m a x p r e s s u r e [ M P a ] 1 0 - 5 1 0 - 4 probability density Fig. 14 . Joint PDFs of loc ation of c entr al cavity and lo c ation of p eak pr essur e (left), and loc ation of p e ak pr essur e and p e ak pr essur e magnitude (right). The lo c ation of the c entr al c avity bubble c orr elates str ongly with p e ak pr essure lo c ation, which itself exhibits str ong inv erse c orr elations with the magnitude of the c ol lapse pr essur es, explaining a wide c onﬁdenc e intervals observe d in Figur e 8 and Figur e 11 . 4. Summary and conclusions. W e hav e presented a non-intrusiv e multi-lev el Mon te Carlo methodology for uncertaint y quantiﬁcation in multi-phase cloud cavi- tation collapse ﬂo ws, together with no vel optimal control v ariate coeﬃcients which main tain the eﬃciency of the algorithm even for w eak correlations among resolution lev els and delivers signiﬁcant v ariance reduction improv ements. W e hav e conducted n umerical exp erimen ts for ca vitating clouds containing 500 cloud cavities, which are randomly (uniformly) distributed within the sp eciﬁed 20 mm radius, and the radii of the cavit y follow a log-normal distribution. The results of these numerical exp er- imen ts hav e revealed signiﬁcant uncertainties in the magnitude of the p eak pressure pulse, emphasizing the relev ance of uncertaint y quan tiﬁcation in cavitating ﬂows. F urthermore, statistical correlation and joint probability densit y function estimates ha ve rev ealed potential underlying causes of this phenomenon. In particular, spatial arrangemen t characteristics of the cloud and its core, such as skewness, cloud inter- action parameter β , and the p osition of the central cavit y hav e been observ ed to ha v e a signiﬁcan t inﬂuence on the resulting pressure ampliﬁcation intensities during the collapse pro cess. All numerical exp eriments were p erformed b y coupling an open source PyMLMC framew ork with Cubism-MPCF, a high p erformance p eta-scale ﬁnite volume solver. Ev olution of collapsing clouds has been simulated b y explicit time stepping sub ject to a CFL stabilit y constrain t on a hierarch y of uniform, structured spatial meshes. Ef- 22 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS ﬁcien t m ulti-lev el Mon te Carlo sampling has b een observed to exhibit more than t w o orders of magnitude in estimated computational speedup when compared to standard Mon te Carlo metho ds, with an additional factor 3 estimated sp eedup due to optimal con trol v ariate co eﬃcien ts. In the present w ork, we ha ve observ ed the eﬃcien t scaling of the prop osed hybrid OF-MLMC-FVM metho d up to the entire MIRA sup ercom- puter consisting of half a million cores. Considering that fault-tolerance mitigation mec hanisms are implemen ted in PyMLMC and ha ve been successfully used, w e expect it to scale linearly and b e suitable for the exa-scale era of numerical computing. The proposed OF-MLMC-FVM can deal with a v ery large n um b er of sources of uncertain ty . In the problems presented here, 2’000 sources of uncertaint y are needed to fully describ e the random initial conﬁguration of the collapsing cloud. T o the b est of our knowledge, curren tly no other metho ds (particularly deterministic metho ds such as quasi Monte Carlo, sto chastic Galerkin, sto chastic collo cation, PGD, ANOV A, or sto c hastic FVM) are able to handle this many sources of uncertain ty , in particular for non-linear problems with solutions which exhibit path-wise lo w regularity and p ossibly non-smo oth dep endence on random input ﬁelds. F urthermore, the prop osed metho dology is well suited for multi-lev el extensions of Marko v Chain Monte Carlo metho ds for Bay esian inv ersion [ 30 , 36 ]. The present m ulti-level metho dology is a p ow erful general purp ose technique for quan tifying uncertain ty in complex ﬂo ws go verned by h yp erbolic systems of non-linear conserv ation laws suc h as cloud cavitation ﬂo w problems. App endix A. Discussion. Here, w e wish to commen t on p ossible shortcomings in the progress of the adaptiv e OF-MLMC algorithm describ ed in subsubsection 2.4.4 . These remarks are meant to assist the application of OF-MLMC and are not directly relev ant for the n umerical results presented in section 3 . A.1. Number of “w arm-up” samples on the ﬁnest resolution lev el. No- tice, that in order to ha v e empiric al estimates of σ 2 ` , σ 2 `,` − 1 , and ˜ σ 2 ` , at least two samples would b e required on eac h level ` = 0 , . . . , L . Enforcing M L ≥ 2 migh t b e v ery ineﬃcien t in the cases when only one sample is actually needed, since in presen tly considered applications the most computationally exp ensive samples are actually at the ﬁnest lev el ` = L . T o a void this, initially only one “w arm-up” sample on the ﬁnest mesh level ` = L could b e computed, i.e., M L = 1. F or subsequen t optimization steps of each M ` , v ariance of level diﬀerence σ 2 `,` − 1 for level ` = L is inferred from av ailable measuremen ts on lo w er resolution lev els using Ba yesian inference, as describ ed in step (3) of the OF-MLMC algorithm. If more than one sample is actually required on the ﬁnest resolution lev el, optimization step (7) of the adaptive OF-MLMC algorithm ab o ve will adjust M L to the correct optimal v alue and additional empirical estimate w ould b e av ailable for ev en more accurate inference. A.2. Iterative control of indicator errors. Final OF-MLMC error ε OF could b e underestimated by ˆ ε OF and b e actually lar ger than the prescrib ed tolerance τ , since we only ensure the estimate d total error ˆ ε OF ≤ τ . Since ˆ ε OF is based on the r andomize d statistical estimators, it is also r andom and has a spread around its mean E [ ˆ ε OF ] = ε OF . W e note, that probability P [ ε OF > ˆ ε OF ] of the resulting MLMC error ε OF b eing lar ger than the estimated error ˆ ε OF (whic h is forced to b e smaller than the prescrib ed tolerance τ ) can be reduced b y suﬃciently incr e asing the n umber of samples ˆ M ` according to the estimates of fourth c enter e d moment (kurtosis) κ 2 ` of the level diﬀerences. Such estimates would provide the varianc e estimates ˆ κ 2 ` ≈ OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 23 E [( ˆ σ 2 ` ) 2 − E [ ˆ σ 2 ` ] 2 ] of the empiric al varianc e estimators ˆ σ 2 ` . Then, increasing ˆ M ` b y sev eral standard deviations ˆ κ ` of ˆ σ 2 ` , the required p ercentile lev el α of conﬁdence in terv al [ α/ 2 , 1 − α/ 2] of ˆ ε OF can b e reduced b elow the prescrib ed tolerance τ , this w ay ensuring the required high probabilit y P [ ε ≤ τ ] ≥ 1 − α of true MLMC error ε not exc e e ding the prescribed tolerance τ . A contin uation MLMC metho d incorporating similar techniques for estimation and control of error conﬁdence in terv als w as already prop osed in [ 15 ]. Ho wev er, up dating estimates after e ach computed sample could be v ery ineﬃcient for large HPC application, since such incremental tec hniques require hea vy synchroni zation and would make eﬃcien t load balancing on distributed man y- core systems very c hallenging. A.3. Sample “recycling”. Alternative optimal co eﬃcien ts for each level in MLMC estimator was suggested in [ 54 , 55 ], where a multi-ﬁdelit y Mon te Carlo method is describ ed. There, some statistical realizations (samples) are r e-use d , i.e., the same result is used in b oth estimates E M ` [ α ` q ` ] and − E M ` +1 [ α ` q ` ], eac h con tributing to a separate diﬀerence in the telescoping MLMC estimator ( 32 ). Suc h “recycling” strategy requires less sampling, how ever, error analysis complexity is highly increased due to additional statistical dep endencies, absent in OF-MLMC method. On the other hand, for “recycled” sampling, the resulting error minimization problem is separable in terms of co eﬃcien ts α ` and n umber of samples M ` , and hence no linear system ( 31 ) needs to b e solved [ 54 ]. The linear system ( 31 ) is, ho w ever, very small, sparse, and needs to b e solved only a few times, hence is not a b ottlenec k of this algorithm. A.4. T ruly optimal num b er of MC samples on eac h resolution level. A discrete optimization problem could b e formulated, a voiding round-oﬀ op erations in ( 34 ) or ( 35 ), and pro viding truly optimal inte ger M ` , as suggested in [ 52 ]. W e note, ho wev er, that suc h round-oﬀs do not inﬂuence the eﬃciency of the method on coarser lev els, where the num b er of samples is large. F urthermore, round-oﬀ ineﬃciencies are most probably ov er-shadow ed by the used approximate estimators for ˜ σ 2 ` . Ac kno wledgmen ts. Authors ac knowledge supp ort from the following organi- zations: Innov ativ e and Nov el Computational Impact on Theory and Exp eriment (INCITE) program, for aw arding computer time under the pro ject CloudPredict; Ar- gonne Leadership Computing F acilit y , whic h is a DOE Oﬃce of Science User F acilit y supp orted under Con tract DE-A C02-06CH11357, for providing access and support for MIRA, CETUS, and COOLEY systems. Partnership for Adv anced Computing in Eu- rop e (PRA CE) pro jects PRA091 and Pra09 2376, together with J ¨ ulic h and CINECA Sup ercomputing Centers; Swiss National Sup ercomputing Cen ter (CSCS) for compu- tational resources gran t under pro ject ID s500. J ˇ S w ould lik e to thank Stephen W u for his contribution and fruitful discussions on the optimal v ariance reduction tec hniques in the m ulti-level Mon te Carlo estimators, and Timothy J. Barth for hosting him as a scien tiﬁc visitor at NASA Ames Researc h Center (California, United States) and for collab orations on robust kernel densit y estimators. REFERENCES [1] N. A. Adams and S. J. Schmidt , Sho cks in c avitating ﬂows , in Bubble Dynamics and Sho ck W aves, C. F. Delale, ed., vol. 8 of Sho ck W av e Science and T echnology Reference Library , Springer Berlin Heidelb erg, 2013, pp. 235–256. [2] M. R. Baer and J. W. Nunzia to , A two-phase mixtur e theory for the deﬂagration-to- detonation tr ansition (DDT) in r e active gr anular materials , International Journal of Mul- tiphase Flow, 12 (1986), pp. 861–889. 24 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS [3] A. Bar th, C. Schw ab, and N. Zollinger , Multi-level MC metho d for el liptic PDEs with sto chastic c o eﬃcients , Numerische Mathematik, 119 (2011), pp. 123–161. [4] J. R. Blake, M. C. Hooton, P. B. R obinson, and R. P. Tong , Col lapsing c avities, tor oidal bubbles and jet impact , Phil osophical T ransactions of the Roy al Society of London. Series A: Mathematical, Physical and Engineering Sciences, 355 (1997), pp. 537–550. [5] L. Bonfiglio, P. Perdikaris, S. Brizzolara, and G. E. Karniad akis , A multi-ﬁdelity fr ame- work for investigating the p erformanc e of sup er-cavitating hydr ofoils under uncertain ﬂow c onditions , AIAA Non-Deterministic Approaches, (2017). [6] N. Bremond, M. Arora, C.-D. Ohl, and D. Lohse , Controlled multibubble surfac e cavitation , Physical Review Letters, 96 (2006). [7] C. E. Brennen , Cavitation and bubble dynamics , Oxford Universit y Press, USA, 1995. [8] C. E. Brennen , Cavitation in medicine , Interface F o cus, 5 (2015), p. 20150022. [9] C. E. Brennen, G. Reisman, and Y.-C. W ang , Shock waves in cloud c avitation . Twen ty-First Symposium on Nav al Hydro dynamics, 1997. [10] E. A. Brujan, T. Ikeda, and Y. Ma tsumoto , Sho ck wave emission fr om a cloud of bubbles , Soft Matter, 8 (2012), p. 5777. [11] F. Cappello , F ault toleranc e in petasc ale/exasc ale systems: Curr ent know le dge, chal lenges and r ese ar ch opp ortunities , International Journal of High P erformance Computer Applica- tions, 23 (2009), pp. 212–226. [12] G. L. Chahine , Pr essures gener ated by a bubble cloud c ol lapse , Chemical Engineering Com- munications, 28 (1984), pp. 355–367. [13] Q. Y. Chen, D. Gottlieb, and J. S. Hestha ven , Unc ertainty analysis for steady ﬂow in a dual thr o at nozzle , Journal of Computational Physics, 204 (2005), pp. 378–398. [14] K. A. Cliffe, M. B. Giles, R. Scheichl, and A. L. Teckentrup , Multilevel Monte Carlo metho ds and applic ations to el liptic PDEs with random c oeﬃcients , Computing and Visu- alization in Science, 14 (2011), pp. 3–15. [15] N. Collier, A.-L. Haji-Ali, F. Nobile, E. von Schwerin, and R. Tempone , A c ontinuation multilevel Monte Carlo algorithm , BIT Numerical Mathematics, 55 (2015), pp. 399–432. [16] P. M. Congedo, E. Goncal ves, and M. G. Rodio , Ab out the unc ertainty quantiﬁc ation of turbulenc e and c avitation models in c avitating ﬂows simulations , European Journal of Mechanics - B/Fluids, 53 (2015), pp. 190 – 204. [17] C. C. Coussios and R. A. R oy , Applic ations of ac oustics and c avitation to noninvasive ther apy and drug delivery , Annual Review of Fluid Mec hanics, 40 (2008), pp. 395–420. [18] D. P. Schmidt et al. , Assessment of the pr e diction c ap abilities of a homo gene ous c avitation mo del for the c ol lapse char acteristics of a vapour-bubble cloud , in WIMRC 3rd Interna- tional Cavitation F orum, Cov entry , U.K., 2011. [19] T. J. Dodwell, C. Ketelsen, R. Scheichl, and A. L. Teckentrup , A hier archic al multilevel Markov Chain Monte Carlo algorithm with applic ations to unc ertainty quantiﬁcation in subsurfac e ﬂow , SIAM/ASA Journal on Uncertain ty Quan tiﬁcation, 3 (2015), pp. 1075– 1108. [20] D. Elsak out, M. Christie, and G. Lord , Multilevel Markov Chain Monte Carlo (MLMCMC) for unc ertainty quantiﬁc ation , in So ciety of P etroleum Engineers - SPE North Africa T ech- nical Conference and Exhibition 2015, NA TC 2015, 2015. [21] M. B. Giles , Multi-level Monte Carlo p ath simulation , Oper. Res., 56 (2008), pp. 607–617. [22] M. B. Giles , Multilevel Monte Carlo metho ds , Acta Numerica, (2015). [23] M. B. Giles and C. Reisinger , Sto chastic ﬁnite diﬀer ences and multilevel Monte Carlo for a class of SPDEs in ﬁnanc e , Preprint, Oxford-Man Institute of Quantitativ e Finance and Mathematical Institute, Universit y of Oxford, (2011). [24] D. Gottlieb and D. Xiu , Galerkin metho d for wave e quations with uncertain c o eﬃcients , Commun. Comput. Ph ys., 3 (2008), pp. 505–518. [25] S. Gottlieb and C.-W. Shu , T otal variation diminishing Runge-Kutta schemes , Mathematics of Computation of the American Mathematical So ciety , 67 (1998), pp. 73–85. [26] P. E. Hadjidoukas, P. Angelikopoulos, C. P ap adimitriou, and P. Koumoutsak os , Π 4U: A high performanc e c omputing framework for Bayesian unc ertainty quantiﬁc ation of c omplex mo dels , Journal of Computational Physics, 284 (2015), pp. 1–21. [27] P. E. Hadjidoukas, D. Rossinelli, B. Hejazialhosseini, and P. K oumoutsakos , F r om 11 to 14.4 PFLOPs: Performanc e optimization for ﬁnite volume ﬂow solver , in Proceedings of the 3rd International Conference on Exascale Applications and Softw are, 2015. [28] P. E. Hadjidoukas, D. R ossinelli, F. Wermelinger, J. Sukys, U. Rasthofer, C. Conti, B. Hejazialhosseini, and P. K oumoutsakos , High thr oughput simulations of two-phase ﬂows on Blue Gene/Q , in Parallel Computing: On the Road to Exascale, Proceedings of the International Conference on Parallel Computing, P arCo 2015, 1-4 Septem b er 2015, OF-MLMC FOR UNCER T AINTY QUANTIFICA TION IN CLOUD CA VIT A TION 25 Edinburgh, Scotland, UK, 2015, pp. 767–776. [29] S. Heinrich , Multilevel Monte Carlo metho ds , Large-scale scientiﬁc computing, Third in terna- tional conference LSSC 2001, Sozop ol, Bulgaria, 2001, Lecture Notes in Computer Science, Springer V erlag, 2170 (2001), pp. 58–67. [30] V. H. Hoang, C. Schw ab, and A. M. Stuar t , Sp arse MCMC gp c ﬁnite element metho ds for b ayesian inverse problems , Inverse Problems, (2014). [31] X. Y. Hu, B. C. Khoo, N. A. Adams, and F. L. Huang , A c onservative interfac e metho d for c ompr essible ﬂows , Journal of Computational Physics, 219 (2006), pp. 553–578. [32] G.-S. Jiang and C.-W. Shu , Eﬃcient implementation of weighte d ENO schemes , Journal of computational physics, 126 (1996), pp. 202–228. [33] E. Johnsen and T. Colonius , Implementation of WENO schemes in c ompr essible multic om- p onent ﬂow problems , Journal of Computational Physics, 219 (2006), pp. 715 – 732. [34] E. Johnsen and T. Colonius , Numeric al simulations of non-spheric al bubble c ol lapse , Journal of Fluid Mechanics, 629 (2009), pp. 231–262. [35] A. K. Kapila, R. Menik off, J. B. Bdzil, S. F. Son, and D. S. Stew ar t , Two-phase mo deling of deﬂagr ation-to-detonation transition in granular materials: R e duc e d equations , Physics of Fluids, 13 (2001), pp. 3002–3024. [36] C. Ketelsen, R. Scheichl, and A. L. Teckentrup , A hier ar chic al multilevel Markov Chain Monte Carlo algorithm with applic ations to unc ertainty quantiﬁc ation in subsurfac e ﬂow , arXiv:1303.7343, (2013). [37] E. Lauer, X. Y. Hu, S. Hickel, and N. A. Adams , Numeric al investigation of c ol lapsing c avity arr ays , Physics of Fluids, 24 (2012), p. 052104. [38] S. Li , Tiny bubbles chal lenge giant turbines: Thr ee gor ges puzzle , Interface F o cus, 5 (2015), p. 20150020. [39] G. Lin, C. H. Su, and G. E. Karniadakis , The sto chastic piston pr oblem , Pro c. Natl. Acad. Sci. USA, 101 (2004), pp. 15840–15845. [40] X. Ma and N. Zabaras , An adaptive hier ar chical sp arse grid c ol loc ation algorithm for the solution of sto chastic diﬀer ential e quations , J. Comp. Phys., 228 (2009), pp. 3084–3113. [41] R. Menikoff and B. J. Plohr , The Riemann pr oblem for ﬂuid ﬂow of re al materials , Reviews of Mo dern Physics, 61 (1989), pp. 75–130. [42] S. Mishra and C. Schw ab , Sp arse tensor multi-level Monte Carlo ﬁnite volume metho ds for hyp erb olic c onservation laws with r andom initial data , Math. Comp., 280 (2012), pp. 1979– 2018. [43] S. Mishra, C. Schw ab, and J. ˇ Sukys , Multi-level Monte Carlo ﬁnite volume metho ds for nonline ar systems of c onservation laws in multi-dimensions , Journal of Computational Physics, 231 (2012), pp. 3365–3388. [44] S. Mishra, C. Schw ab, and J. ˇ Sukys , Multi-level Monte Carlo ﬁnite volume metho ds for shal- low water e quations with unc ertain top o gr aphy in multi-dimensions , SIAM J. Sci. Comput., 24 (2012), pp. B761–B784. [45] S. Mishra, C. Schw ab, and J. ˇ Sukys , Multi-level Monte Carlo ﬁnite volume metho ds for unc ertainty quantiﬁcation of ac oustic wave prop agation in r andom heter ogene ous layer e d me dium , Journal of Computational Physics, 312 (2016), pp. 192–217. [46] K. A. Mørch , On the c ol lapse of c avity clusters in ﬂow c avitation , in Ca vitation and Inhomo- geneities in Underwater Acoustics, W. Lauterb orn, ed., Springer, 1980, pp. 95–100. [47] K. A. Mørch , Energy consider ations on the c ol lapse of c avity clusters , Applied Scien tiﬁc Research, 38 (1982), pp. 313–321. [48] F. M ¨ uller , Sto chastic metho ds for unc ertainty quantiﬁc ation in subsurfac e ﬂow and tr ansp ort pr oblems , Dissertation ETH No. 21724, 2014. [49] F. M ¨ uller, P. Jenny, and D. W. Meyer , Multilevel Monte Carlo for two phase ﬂow and BuckleyL ever ett tr ansport in r andom heter ogene ous por ous media , Journal of Computa- tional Physics, 250 (2013), pp. 685–702. [50] A. Murrone and H. Guillard , A ﬁve e quation r e duc e d mo del for c ompr essible two phase ﬂow pr oblems , Journal of Computational Physics, 202 (2005), pp. 664–698. [51] L. P arussini, D. Venturi, P. Perdikaris, and G. E. Karniadakis , Multi-ﬁdelity Gaussian pr o c ess r egr ession for pre diction of random ﬁelds , Journal of Computational Ph ysics, 336 (2017), pp. 36–50. [52] S. P a uli , Optimal numb er of multilevel Monte Carlo levels and their fault toler ant applic ation , Dissertation ETH, 2014. [53] S. P auli, M. Kohler, , and P. Arbenz , A fault tolerant implementation of multi-level monte c arlo metho ds , Parallel Computing – Accelerating Computational Science and Engineering (CSE), Adv ances in Parallel Computing, IOS Press, 23 (2014), pp. 471–480. [54] B. Peherstorfer, K. Willco x, and M. Gunzbur ger , Optimal mo del management for mul- 26 ˇ SUKYS, RASTHOFER, WERMELINGER, HADJIDOUKAS, AND KOUMOUTSAK OS tiﬁdelity Monte Carlo estimation , ACDL T echnical Report TR15-2, (2015). [55] B. Peherstorfer, K. Willco x, and M. Gunzburger , Survey of multiﬁdelity methods in unc ertainty pr op agation , infer enc e , and optimization , Preprin t, (2016), pp. 1–57. [56] P. Perdikaris, D. Venturi, J. O. Ro yset, and G. E. Karniadakis , Multi-ﬁdelity mo del ling via r e cursive c o-kriging and Gaussian Markov r andom ﬁelds , Proceedings of the Roy al Society of London A: Mathematical, Physical and Engineering Sciences, 471 (2015). [57] G. Perigaud and R. Saurel , A c ompr essible ﬂow mo del with c apil lary eﬀects , Journal of Computational Physics, 209 (2005), pp. 139–178. [58] G. Poette, B. Despr ´ es, and D. Lucor , Unc ertainty quantiﬁc ation for systems of c onservation laws , J. Comput. Phys., 228 (2009), pp. 2443–2467. [59] T. J. Poinsot and S. K. Lele , Boundary c onditions for dire ct simulations of c ompressible visc ous ﬂows , Journal of Computational Physics, 101 (1992), pp. 104–129. [60] M. Raissi, P. Perdikaris, and G. E. Karniadakis , Inferring solutions of diﬀer ential e quations using noisy multi-ﬁdelity data , Journal of Computational Physics, 335 (2017), pp. 736–746. [61] U. Rasthofer, F. Wermelinger, P. Hadijdoukas, and P. Koumoutsak os , L arge sc ale simulation of cloud cavitation c ol lapse , In pro ceedings of The In ternational Conference on Computational Science (accepted), (2017). [62] U. Rasthofer, F. Wermelinger, J. ˇ Sukys, P. Hadijdoukas, and P. Koumoutsakos , Pa- r ameter study in cloud c avitation collapse , Manuscript in progress, (2017). [63] V. C. Ra ykar and R. Duraisw ami , F ast optimal b andwidth sele ction for kernel density es- timation , In Pro ceedings of the sixth SIAM International Conference on Data Mining, Bethesda, (2006), pp. 524–528. [64] D. R ossinelli, B. Hejazialhosseini, P. E. Hadjidoukas, C. Bekas, A. Curioni, A. Ber tsch, S. Futral, S. J. Schmidt, N. A. Adams, and P. Koumoutsakos , 11 PFLOP/s simula- tions of cloud c avitation c ol lapse , in Pro ceedings of the International Conference on High Performance Computing, Netw orking, Storage and Analysis, SC ’13, New Y ork, NY, USA, 2013, ACM, pp. 3:1–3:13. [65] D. P. Schmidt and M. L. Corradini , The internal ﬂow of diesel fuel injector nozzles: A r eview , International Journal of Engine Research, 2 (2001), pp. 1–22. [66] S. J. Shea ther and M. C. Jones , A r eliable data-b ase d bandwidth sele ction metho d for kernel density estimation , Journal of the Royal Statistical So ciety . Series B (Metho dological), 53 (1991), pp. 683–690. [67] J. ˇ Sukys , A daptive load b alancing for massively p arallel multi-level Monte Carlo solvers , PP AM 2013, Part I, LNCS. Springer, Berlin Heidelb erg, 8384 (2014), pp. 47–56. [68] J. ˇ Sukys , Robust multi-level Monte Carlo Finite V olume metho ds for systems of hyp erbolic c onservation laws with r andom input data , Dissertation ETH No. 21990, 2014. [69] J. ˇ Sukys , PyMLMC , Av ailable from: https://github.com/cselab/PyMLMC , (2017). [70] J. ˇ Sukys, S. Mishra, and C. Schw ab , Static lo ad b alancing for multi-level Monte Carlo ﬁnite volume solvers , PP AM 2011, Part I, LNCS. Springer, Heidelberg, 7203 (2012), pp. 245–254. [71] A. Tiw ari, J. B. Freund, and C. P ant ano , A diﬀuse interfac e mo del with immiscibility pr eservation , Journal of Computational Physics, 252 (2013), pp. 290–309. [72] A. Tiw ari, C. P ant ano, and J. B. Freund , Gr owth-and-collapse dynamics of smal l bubble clusters ne ar a wal l , Journal of Fluid Mechanics, 775 (2015), pp. 1–23. [73] E. F. Toro, M. Spruce, and W. Speares , R estor ation of the c ontact surfac e in the HLL- Riemann solver , Sho ck W av es, 4 (1994), pp. 25–34. [74] J. Tr yoen, O. L. Maitre, M. Ndjinga, and A. Ern , Intrusive proje ction metho ds with upwinding for unc ertain non-line ar hyp erb olic systems , J. Comput. Ph ys., 229 (2010), pp. 6485–6511. [75] X. W an and G. E. Karniad akis , L ong-term behaviour of p olynomial chaos in sto chastic ﬂow simulations , Comput. Meth. Appl. Mech. Engg., 195 (2006), pp. 5582–5596. [76] J. Williamson , L ow-stor age Runge-Kutta schemes , Journal of Computational Physics, 35 (1980), pp. 48–56. [77] J. A. S. Witteveen, A. Loeven, and H. Bijl , An adaptive stochastic ﬁnite element appro ach b ase d on Newton-Cotes quadr atur e in simplex elements , Computer & Fluids, 38 (2009), pp. 1270–1288. [78] D. Xiu and J. S. Hestha ven , High-or der c ol loc ation metho ds for diﬀer ential e quations with r andom inputs , SIAM J. Sci. Comput., 27 (2005), pp. 1118–1139. [79] K. Y amamoto , Investigation of bubble clouds in a cavitating jet , in Mathematical Fluid Dy- namics, Present and F uture, Springer Nature, 2016, pp. 349–373.

Optimal fidelity multi-level Monte Carlo for quantification of uncertainty in simulations of cloud cavitation collapse

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment