Monte Carlo Type Method Research Articles

Statistical Methods for the Study of Computer Experiment Failures: Application to a Fuel-Coolant-Interaction Simulation Code

In the framework of risk assessment in nuclear accidents, simulation tools are widely used to understand and model physical phenomena. These simulation tools take into account a large number of uncertain input parameters. We often use Monte Carlo–type methods to explore their range of variation: The input space is randomly sampled, and a code run is performed on each sampled point. However, some of these code runs may fail to converge. Analyzing these code failures to understand which of the inputs have the most influence on them leads to a better understanding of how the code works. It also intends to improve the robustness of the simulation software and code computations. For this purpose, we propose two complementary approaches performing a statistical analysis of the code failures. The first approach is based on goodness-of-fit tests and compares conditional probability distributions according to code failures to a reference one. A second approach, based on a dependence measure named the Hilbert-Schmidt Independence Criterion, provides another way to measure the global dependence between the inputs and the code failures. The development of this methodology is carried out in the context of severe nuclear accidents. More especially, the presented methods are applied for the study of the simulation code MC3D, which simulates the fuel-coolant interaction in a severe nuclear accident context.

Statistically equivalent surrogate material models: Impact of random imperfections on the elasto-plastic response

Manufactured materials usually contain random imperfections due to the fabrication process, e.g., the 3D-printing, casting, etc. These imperfections affect significantly the effective material properties and result in uncertainties in the mechanical response. Numerical analysis of the effects of the imperfections and the uncertainty quantification (UQ) can be often done by use of digital stochastic surrogate material models. In this work, we present a new flexible class of surrogate models depending on a small number of parameters and a calibration strategy ensuring that the constructed model fits to the available observation data, with special focus on two-phase materials. The surrogate models are constructed as the level-set of a linear combination of an intensity field representing the topological shape and a Gaussian perturbation representing the imperfections, allowing for fast sampling strategies. The mathematical design parameters of the model are related to physical ones and thus easy to interpret. The calibration of the model parameters is performed using progressive batching sub-sampled quasi-Newton minimization, using a designed distance measure between the synthetic samples and the data. Then, employing a fast sampling algorithm, an arbitrary number of synthetic samples can be generated to use in Monte Carlo type methods for prediction of effective material properties. In particular, we illustrate the method in application to UQ of the elasto-plastic response of an imperfect octet-truss lattice which plays an important role in additive manufacturing. To this end, we study the effective material properties of the lattice unit cell under elasto-plastic deformations and investigate the sensitivity of the effective Young’s modulus to the imperfections.

On a Monte Carlo scheme for some linear stochastic partial differential equations

Abstract The aim of this paper is to study the simulation of the expectation for the solution of linear stochastic partial differential equation driven by the space-time white noise with the bounded measurable coefficient and different boundary conditions. We first propose a Monte Carlo type method for the expectation of the solution of a linear stochastic partial differential equation and prove an upper bound for its weak rate error. In addition, we prove the central limit theorem for the proposed method in order to obtain confidence intervals for it. As an application, the Monte Carlo scheme applies to the stochastic heat equation with various boundary conditions, and we provide the result of numerical experiments which confirm the theoretical results in this paper.

A Bayesian Monte Carlo method for computing the Shapley value: Application to weighted voting and bin packing games

The Shapley value is one of the most important solution concepts in cooperative game theory, it satisfies uniqueness and fairness properties while its main drawback lies in its high computational complexity. A number of approximation methods have been provided to overcome such intractability. Monte Carlo type methods are the most general and practical approaches approximating the Shapley value thanks to coalitions sampling. Their efficiency depends on the associated sampling and approximation steps. Most of the contributions focus on the sampling step that suffers from a loss of information and the difficulty to deliver a value satisfying certain expected properties such as efficiency and monotonicity. In this paper, we propose an improvement of the approximation step. It associates a Bayesian approach to that of Monte Carlo, to derive a Shapley value approximation, preserving certain expected properties essential for handling real-world applications. We develop our approach for a class of games with binary marginal contributions. We apply it for the calculation of the Shapley value for the weighted voting and bin packing games.

Dinucleotide evolutionary dynamics in influenza A virus.

Significant biases of dinucleotide composition in many RNA viruses including influenza A virus have been reported in recent years. Previous studies have showed that a codon-usage-altered influenza mutant with elevated CpG usage is attenuated in mammalian in vitro and in vivo models. However, the relationship between dinucleotide preference and codon usage bias is not entirely clear and changes in dinucleotide usage of influenza virus during evolution at segment level are yet to be investigated. In this study, a Monte Carlo type method was applied to identify under-represented or over-represented dinucleotide motifs, among different segments and different groups, in influenza viral sequences. After excluding the potential biases caused by codon usage and amino acid sequences, CpG and UpA were found under-represented in all viral segments from all groups, whereas UpG and CpA were found over-represented. We further explored the temporal changes of usage of these dinucleotides. Our analyses revealed significant decrease of CpG frequency in Segments 1, 3, 4, and 5 in seasonal H1 virus after its re-emergence in humans in 1977. Such temporal variations were mainly contributed by the dinucleotide changes at the codon positions 3-1 and 2-3 where silent mutations played a major role. The depletions of CpG and UpA through silent mutations consequently led to over-representations of UpG and CpA. We also found that dinucleotide preference directly results in significant synonymous codon usage bias. Our study helps to provide details on understanding the evolutionary history of influenza virus and selection pressures that shape the virus genome.

A review and perspective on a convergence analysis of the direct simulation Monte Carlo and solution verification

In 1963, G. A. Bird published a research note on his investigation of a rigid sphere gas reaching translational equilibrium using a Monte Carlo type method. Since then, the method has been developed into a primary workhorse to computationally solve the Boltzmann kinetic equation. As it is increasingly applied to challenging problems in the real world, verification studies of the method have become a critical issue. In this paper, we review previous studies on this challenging subject and present a perspective on a convergence analysis of the direct simulation Monte Carlo (DSMC) method and solution verification. During this process, a verification method based on the physical laws of conservation is studied in depth. In particular, a convergence history plot on all three types of computational errors—decomposition, statistical, and round-off—is presented for two benchmark problems. Finally, future research topics to maximize the full potential of the DSMC method, pioneered by the late G. A. Bird, are suggested.

The Pleistocene contribution to the net erosion and sedimentary conditions in the outer Bear Island Trough, western Barents Sea

The Pleistocene sedimentary conditions and the glacial contribution to net erosion were determined for the outer Bear Island Trough by using a Monte Carlo-type method. The approach uses ages for glacial/interglacial periods that were based on the regional ice-sheet volume curve. The results indicate that the western Barents Sea was glaciated during four marine isotope stages: MIS 16 (635.6–624.7 ka), MIS 12 (438.7–428.0 ka), MIS 6 (138.6–134.6 ka) and MIS 2 (19.3–16.0 ka) for a total duration of 29 ka. The first glacial event (before 440 ka, MIS 16) resulted in homogeneous erosion over the study area (with an erosion rate of 24.2 ± 8.5 mm/a). After 440 ka, a change in sedimentary conditions resulted in inhomogeneous erosion rates over the study area from −12.6 ± 1.6 (i.e. net deposition) to 1.6 ± 1.8 mm/a (i.e. net erosion). The most likely values of average deposition rates during interglacial periods were modelled as 0.12 ± 0.1 mm/a. In the outer Bear Island Trough the net erosion was found to be mainly the effect of tectonic uplift and subsequent erosion prior to the glacial ages. The results show that in most parts of the study area the Pleistocene glacial contribution to the total net erosion was small: the most likely glacial contribution in the eastern part of the area reaches 100 m, which is about 9 % of the total net erosion, while in the westernmost part the glaciations did not contribute to the net erosion.

Direct discrete simulation of the kinetic Landau-Fokker-Planck equation with an alternating external electromagnetic field

A stochastic method of simulation of Coulomb interaction is considered. The main idea of the method is to approximate the nonlinear Landau kinetic collision integral by the Boltzmann integral. In its realization, the method can be attributed to a wide class of Monte Carlo-type methods. It is easily combined with the existing particle methods used to simulate collisionless plasmas. This is important for simulation of the dynamics of both laboratory and space plasmas when the mean free path of plasma particles is comparable with the plasma inhomogeneity scale length. Illustrative examples of relaxation of two-temperature plasma being subject to a high-frequency alternating electric field are given, and differences from their classical description are considered. The method satisfies the conservation laws for the number of particles, momentum, and energy and is simple and efficient in implementation.

Quadrature Methods for the Calculation of Subgrid Microphysics Moments

Abstract Many cloud microphysical processes occur on a much smaller scale than a typical numerical grid box can resolve. In such cases, a probability density function (PDF) can act as a proxy for subgrid variability in these microphysical processes. This method is known as the assumed PDF method. By placing a density on the microphysical fields, one can use samples from this density to estimate microphysics averages. In the assumed PDF method, the calculation of such microphysical averages has primarily been done using classical Monte Carlo methods and Latin hypercube sampling. Although these techniques are fairly easy to implement and ubiquitous in the literature, they suffer from slow convergence rates as a function of the number of samples. This paper proposes using deterministic quadrature methods instead of traditional random sampling approaches to compute the microphysics statistical moments for the assumed PDF method. For smooth functions, the quadrature-based methods can achieve much greater accuracy with fewer samples by choosing tailored quadrature points and weights instead of random samples. Moreover, these techniques are fairly easy to implement and conceptually similar to Monte Carlo–type methods. As a prototypical microphysical formula, Khairoutdinov and Kogan’s autoconversion and accretion formulas are used to illustrate the benefit of using quadrature instead of Monte Carlo or Latin hypercube sampling.

Surfing the Network for Ranking by Multidamping

PageRank is one of the most commonly used techniques for ranking nodes in a network. It is a special case of a family of link-based rankings, commonly referred to as functional rankings. Functional rankings are computed as power series of a stochastic matrix derived from the adjacency matrix of the graph. This general formulation of functional rankings enables their use in diverse applications, ranging from traditional search applications to identification of spam and outliers in networks. This paper presents a novel algorithmic (re)formulation of commonly used functional rankings, such as LinearRank, TotalRank and Generalized Hyperbolic Rank. These rankings can be approximated by finite series representations. We prove that polynomials of stochastic matrices can be expressed as products of Google matrices (matrices having the form used in Google's original PageRank formulation). Individual matrices in these products are parameterized by different damping factors. For this reason, we refer to our formulation as multidamping. We demonstrate that multidamping has a number of desirable characteristics: (i) for problems such as finding the highest ranked pages, multidamping admits extremely fast approximate solutions; (ii) multidamping provides an intuitive interpretation of existing functional rankings in terms of the surfing habits of model web users; (iii) multidamping provides a natural framework based on Monte Carlo type methods that have efficient parallel and distributed implementations. It also provides the basis for constructing new link-based rankings based on inhomogeneous products of Google matrices. We present algorithms for computing damping factors for existing functional rankings analytically and numerically. We validate various benefits of multidamping on a number of real datasets.

Regularized inversion of microphysical atmospheric particle parameters: Theory and application

Retrieving the distribution of aerosols in the atmosphere via remote sensing techniques is a highly complex task that requires dealing with a wide range of different problems stemming both from Physics and Mathematics. We focus on retrieving this distribution from multi-wavelength lidar data for aerosol ensembles consisting of spherical particles via an iterative regularization technique. The optical efficiencies for spherical scatterers are examined to account for the behavior of the underlying integral equation. The ill-posedness of the problem and the conditioning of the discretized problem are analyzed. Some critical points in the model, like the assumed wavelength-independence of the refractive index and the fixed grid of investigated refractive indices, are studied with regard to their expected impact on the regularized solution. A new Monte-Carlo type method is proposed for retrieval of the refractive index. To validate the results, the developed algorithm is applied to two measurement cases of burning biomass gained from multi-wavelength Raman lidar.

A new method for the statistical simulation of the virtual values of parameters in inverse orbital dynamics problems

A Monte Carlo-type method for simulating virtual values of the parameters in inverse orbital dynamics problems for highly nonlinear cases is proposed. The method is based on imitating Fisher’s statistics employed to specify the confidence region, and is implemented by solving repeatedly nonlinear least-squares problems with various samples of simulated observations obtainable by suitable random variations.

Input uncertainty propagation methods and hazard mapping of geophysical mass flows

This paper presents several standard and new methods for characterizing the effect of input data uncertainty on model output for hazardous geophysical mass flows. Note that we do not attempt here to characterize the inherent randomness of such flow events. We focus here on the problem of characterizing uncertainty in model output due to lack of knowledge of such input for a particular event. Methods applied include classical Monte Carlo and Latin hypercube sampling and more recent stochastic collocation, polynomial chaos, spectral projection and a newly developed extension thereof named polynomial chaos quadrature. The simple and robust samplings based Monte Carlo type methods are usually computationally intractable for reasonable physical models, while the more sophisticated and computationally efficient polynomial chaos method often breaks down for complex models. The spectral projection and polynomial chaos quadrature methods discussed here produce results of quality comparable to the polynomial chaos type methods while preserving the simplicity and robustness of the Monte Carlo‐type sampling based approaches at much lower cost. The computational efficiency, however, degrades with increasing numbers of random variables. A procedure for converting the output uncertainty characterization into a map showing the probability of a hazard threshold being exceeded is also presented. The uncertainty quantification procedures are applied first in simple settings to illustrate the procedure and then subsequently applied to the 1991 block‐and‐ash flows at Colima Volcano, Mexico.

Monte Carlo simulation of DNMR spectra of coupled spin systems

A new program MC-DNMR is presented for the simulation of dynamic nuclear magnetic resonance spectra. The algorithm is a Monte Carlo type method based on the extension of single spin vector model to coupled spin systems. This extension is explained in detail and the theory is justified by examples. The main advantage of this program is the significantly smaller sizes of matrices than that in programs based on density matrix theory. So spectra of systems can be simulated that was impossible previously.

Performance of particle in cell methods on highly concurrent computational architectures

Particle in cell (PIC) methods are effective in computing Vlasov-Poisson system of equations used in simulations of magnetic fusion plasmas. PIC methods use grid based computations, for solving Poisson's equation or more generally Maxwell's equations, as well as Monte-Carlo type methods to sample the Vlasov equation. The presence of two types of discretizations, deterministic field solves and Monte-Carlo methods for the Vlasov equation, pose challenges in understanding and optimizing performance on today large scale computers which require high levels of concurrency. These challenges arises from the need to optimize two very different types of processes and the interactions between them. Modern cache based high-end computers have very deep memory hierarchies and high degrees of concurrency which must be utilized effectively to achieve good performance. The effective use of these machines requires maximizing concurrency by eliminating serial or redundant work and minimizing global communication. A related issue is minimizing the memory traffic between levels of the memory hierarchy because performance is often limited by the bandwidths and latencies of the memory system. This paper discusses some of the performance issues, particularly in regard to parallelism, of PIC methods. The gyrokinetic toroidal code (GTC) is used for these studies and a new radial grid decomposition is presented and evaluated. Scaling of the code is demonstrated on ITER sized plasmas with up to 16K Cray XT3/4 cores.

PIC-MC Code to Model Fast Electron Beam Transport Through Dense Matter

A PIC (particle-in-cell)-MC (Monte Carlo) code to model electron beam transport into dense matter is developed. The background target is treated as a cold, stationary fluid and the fast electrons as particles with the relativistic motions. The process is described by a particle-in-cell method with consideration of the influence of both the self-generated electric and magnetic fields as well as collisions between the fast electrons and the target. The collisional part of the code is solved by the Monte Carlo-type method. Furthermore by assuming that the background current balances with the fast electron current, the electric field is given by the Ohm's law and the magnetic field is calculated from the Faraday's law. Both are solved in a two-dimensional cylindrical geometry. The algorithms implemented in the code are demonstrated and the numerical experiments are performed for monoenergy homogeneous fast electron beam transport in an aluminum target when the fields, collision and angular scattering are switched on and off independently.

Sensitivity Analysis of Fatigue Behaviour of Steel Structure under In-Plane Bending

The aim of the present paper is to analyse the influence of some input quantities on the fatigue resistance of a steel structure. The steel element subjected to many times repeated in-plane bending moment was analysed. The fatigue resistance was defined as the number of cycles causing the initial crack size propagation up to the critical size. The variability influence of input random quantities on the fatigue resistance variability was studied by means of the stochastic sensitivity analysis. All input imperfections were considered to be random quantities. The Latin Hypercube Sampling (LHS) numerical simulation method (Monte Carlo type method) was used. It has been found by the stochastic sensitivity analysis that the fatigue resistance variability is most sensitive to the initial crack variability. The paper presented draws attention to the necessity of identifying the statistical characteristics of the initial crack size as exactly as possible because their random variability can largely influence the failure probability of a structure. Large diversity in applying the statistical characteristics of initial crack size by numerous specialists is illustrated by the list of international publications at the end of the present paper.

Structure of BRIJ-35 Nonionic Surfactant in Water: A Reverse Monte Carlo Study

There are some contradictions in the literature on the structure of micelles formed by the BRIJ-35 surfactant in water. One can find reasonable differences in the aggregation numbers and micellar sizes, but there is a lack of data on the intermicellar structure. In this study, we reevaluated the small-angle neutron scattering experiments performed previously on the BRIJ-35 surfactant in the concentration range of 5-200 g/dm3 at 20, 40, and 60 degrees C. The data were analyzed with a reverse Monte Carlo-type method developed recently for colloids. The micelles were modeled as spherical cores representing the hydrophobic parts and number of balls put on the cores to mimic the hydrated hydrophilic chains. The simulations provided data on the mean aggregation number and on the extent of hydration of the hydrophilic shell of the micelles. We obtained intermicellar pair-correlation functions indicating different micelle-micelle interactions from the usually assumed hard-sphere ones.

Monte-Carlo Version of Vertex Model Applied to Recrystallization Modeling

A newly developed model based on vertex concept is presented in this paper. Contrary to its standard version, which is strictly deterministic, some concepts of Monte-Carlo type method were introduced. It makes the model more flexible and allows to introduce some parameters appearing in vertex movement equations, which are not easy to express in analytical form. Initial microstructure in the model is characterized by topology, crystallographic orientations and stored energy values of the grains. The boundary energies and mobilities are anisotropic in general. Nucleation mechanism of a given type is selected at the beginning of calculations. Deformation texture, stored energy distribution and initial microstructure are input parameters of the model. The aim of the calculations is to predict the texture and microstructure modifications during recrystallization. The model was also applied to the study of the kinetics of grain growth and recrystallization. The preliminary tests of the model are presented.

Location of EPR-Active Spins Buried in Proteins from the Simulation of the Spin−Lattice Relaxation Enhancement Caused by Dy(III) Complexes

A computational method is described that simulates the EPR microwave power saturation characteristics of a protein sample containing an EPR signal from a buried free radical in a frozen solution containing paramagnetic spin-lattice relaxation-enhancement agents such as [Dy(III)HEDTA], where HEDTA is N-hydroxyethylenediaminetriacetate. The specific experiment modeled in this work is the EPR progressive power saturation experiment, where a series of spectra are recorded with increasing microwave power incident on the sample. A fit of the data series yields a characteristic power at half saturation, or P 1 / 2 , that is proportional to the spin-lattice relaxation rate of the signal observed in the sample. The addition of DyHEDTA to the protein solution increases P 1 / 2 through dipole-dipole coupling to the buried free radical by an amount strongly dependent on the distance r from the fast-relaxing DyHEDTA to the slow-relaxing radical. The measured value of P 1 / 2 reflects an ensemble of protein molecules in which each buried spin interacts with a specific distribution of DyHEDTA molecules. To simulate the effect of DyHEDTA on the P 1 / 2 of the buried radical spin, we have created a Monte Carlo-type method that models the power-saturation characteristics of the bulk sample by creating saturation profiles of individual, randomly oriented spin ensembles. These individual spin ensembles incorporate arbitrary structural data for the protein as a region where the individual molecules of DyHEDTA are spatially excluded from the calculation, as opposed to previous methods where simple geometric models were employed. The calculation method is validated mathematically and is used to simulate experimental data from the heme-nitric oxide-containing protein systems myoglobin and horseradish peroxidase to demonstrate the effectiveness of the new simulation method.