Monte Carlo Particles Research Articles

A dynamic compartment model for spatially heterogeneous reactors: Scalar and Monte-Carlo particle mixing

The coupling of multidimensional population balance models with transport equations is often necessary for the simulation of spatially heterogeneous reactors. Because of numerical intractability, a compromise on both the hydrodynamics and the population description has to be made. In this work, the Monte-Carlo method is used to calculate the local integral source term due to the dispersed phase (isodense particles) and a compartment model is used to describe the fluid motion in the reactor. The number density function n(x,t,ξ) is approximated by the number of Monte-Carlo particles per compartment and their movement in the space of compartment uses a stochastic transport algorithm adapted from existing works. The numerical tool is thoroughly and successfully validated in the context of mixing and first order reaction in a 70 L stirred tank equipped with a Rushton turbine. The dynamic method makes use of time resolved CFD simulations to built time dependent flow matrices reflecting the periodic nature of the flow field. The dynamic approach is finally compared to other usual compartment model approaches in terms of complexity, accuracy and rapidity which reveals the benefit of considering the flow unsteadiness within a compartiment model approach for stirred tanks.

Efficient simulation of stochastic interactions among representative Monte Carlo particles

Context. Interaction processes between discrete particles are often modelled with stochastic methods such as the Representative Particle Monte Carlo (RPMC) method which simulate mutual interactions (e.g. chemical reactions, collisions, gravitational stirring) only for a representative subset of n particles instead of all N particles in the system. However, in the traditionally employed computational scheme the memory requirements and the simulation runtime scale quadratically with the number of representative particles. Aims. We want to develop a computational scheme that has significantly lower memory requirements and computational costs than the traditional scheme, so that highly resolved simulations with stochastic processes such as the RPMC method become feasible. Results. In this paper we propose the bucketing scheme, a hybrid sampling scheme that groups similar particles together and combines rejection sampling with a coarsened variant of the traditional discrete inverse transform sampling. For a v-partite bucket grouping, the storage requirements scale with n and v2, and the computational cost per fixed time increment scales with n ⋅ v, both thus being much less sensitive to the number of representative particles n. Extensive performance testing demonstrates the higher efficiency and the favourable scaling characteristics of the bucketing scheme compared to the traditional approach, while being statistically equivalent and not introducing any new requirements or approximations. With this improvement, the RPMC method can be efficiently applied not only with very high resolution but also in scenarios where the number of representative particles increases over time, and the simulation of high-frequency interactions (such as gravitational stirring) as a Monte Carlo process becomes viable.

Enhancing discharge estimation from SWOT satellite data in a tropical tidal river environment

The Surface Water and Ocean Topography (SWOT) mission aims to provide essential data on river width, height and slope in order to estimate worldwide river discharge accurately. This mission offers a powerful tool for monitoring river discharge in dynamic coastal areas, like the Saigon-Dongnai estuary in Southern Vietnam. However, estimating discharge of tidally-influenced rivers using SWOT measurements can be challenging when hydraulic variables have the same order of magnitude as SWOT measurement errors. In this paper we present a methodology to enhance discharge estimation accuracy from SWOT measurements based on simulated SWOT products at the 200 meter node resolution and varying river reach size. We assess measurement error variability and its impact on discharge estimation by employing a Monte Carlo analysis. Our approach significantly improved discharge estimation in the Saigon tidal river, reducing RMSE from 1400 m3/s to 180 m3/s and increasing R² from 0.31 to 0.95. Notably, the percentage of Monte Carlo particles meeting the 30% rRMSE threshold rose from 0% to 79%. This study underscores the feasibility of obtaining reliable discharge estimates from SWOT data in complex coastal areas where hydraulic variables are of the same order of magnitude as SWOT errors. Additionally, the proposed methodology to improve discharge estimation from SWOT measurements is widely adaptable as it can be applied to similar regions and can be combined with any discharge estimation method.

Accurate Fourier-space statistics for line intensity mapping: Cartesian grid sampling without aliased power

ABSTRACT Estimators for n-point clustering statistics in Fourier-space demand that modern surveys of large-scale structure be transformed to Cartesian coordinates to perform Fast Fourier Transforms (FFTs). In this work, we explore this transformation in the context of pixelized line intensity maps (LIM), highlighting potential biasing effects on power-spectrum measurements. Current analyses often avoid a complete resampling of the data by approximating survey geometry as rectangular in Cartesian space, an increasingly inaccurate assumption for modern wide-sky surveys. Our simulations of a $20\, {\times }\, 20\, \text{deg}^2$ 21 cm LIM survey at $0.34\, {\lt }\, z\, {\lt }\, 0.54$ show this assumption biases power-spectrum measurements by ${\gt }\, 20~{{\ \rm per\ cent}}$ across all scales. We therefore present a more robust framework for regridding the voxel intensities on to a 3D FFT field by coordinate transforming large numbers of Monte-Carlo sampling particles. Whilst this unbiases power-spectrum measurements on large scales, smaller scale discrepancies remain, caused by structure smoothing and aliasing from separations unresolved by the grid. To correct these effects, we introduce modelling techniques, higher order particle assignments, and interlaced FFT grids to suppress the aliased power. Using a piecewise cubic spline (PCS) particle assignment and an interlaced FFT field, we achieve sub-per cent accuracy up to 80 per cent of the Nyquist frequency for our 21 cm LIM simulations. We find a more subtle hierarchical improvement in results for higher order assignment schemes, relative to the gains made for galaxy surveys, which we attribute to the extra complexity in LIM from additional discretizing steps. python code accompanying this paper is available at github.com/stevecunnington/gridimp.

Unified gas-kinetic particle method for frequency-dependent radiation transport

This paper proposes a unified gas-kinetic particle (UGKP) method for the frequency-dependent photon transport process. The photon transport is a typical multiscale process governed by the nonlinear radiative transfer equations (RTE). The flow regime of photon transport varies from the ballistic regime to the diffusive regime with respect to optical depth and photon frequency. The UGKP method is an asymptotic preserving (AP) scheme and a regime adaptive scheme. The teleportation error is significantly reduced, and the computational efficiency is remarkably improved in the diffusive regime. Distinguished from the standard multigroup treatment, the proposed UGKP method solves the frequency space in a non-discretized way. Therefore the Rosseland diffusion system can be precisely preserved in the optically thick regime. Taking advantage of the local integral solution of RTE, the distribution of the emitted photon can be constructed from its macroscopic moments. The Monte Carlo particles in the UGKP method need only be tracked before their first collision events, and a re-sampling process is performed to close the photon distribution for each time step. The high computational cost of calculating excessive scattering events can be saved, especially in the optically thick regime. The proposed UGKP method is implicit and removes the light speed constraint on the time step. The particle tracking approach combining the implicit formulation makes the proposed UGKP method an efficient solution algorithm for frequency-dependent radiative transfer problems. We demonstrate with numerical examples the capability of the proposed multi-frequency UGKP method.

Using the Graphics Processing Unit to Evaluate American-Style Derivatives

In this article, the authors apply graphics processing unit (GPU) computation to an American option pricing problem via Monte Carlo (MC) simulations and particle swarm optimization (PSO). Given that computations in both MC and PSO can be vectorized and made independent, the valuation can be readily performed on GPUs. As a result, we can increase the accuracy of the valuation by increasing MC paths and particles without spending more time. For example, with a large number of particles (but allocated to GPUs), convergence can be reached in very few steps. The method introduced in this article can be extended to a wide variety of exotic derivatives or a large portfolio of diverse derivatives (known as an eigen portfolio). This is helpful in both trading and risk management.

Numerical error analysis of SOLPS-ITER simulations of EAST

Plasma edge simulations with codes like SOLPS-ITER are widely employed to interpret fusion experiments. However, numerical errors appearing in such simulations are rarely investigated, despite their potential large impact on simulation results. These errors consist of the statistical error and the bias, both resulting from the finite number of employed EIRENE Monte Carlo particles and incomplete convergence, and the discretization error due to the finite resolution of the computational grids. In this contribution, the resulting numerical errors on simulations of pure deuterium and neon seeded H-mode EAST discharges are examined. The statistical error can be kept small compared to other numerical error contributions by averaging the plasma profiles. This allows investigating the bias and discretization errors using Richardson extrapolation. It is shown that grid refinement and the number of employed Monte Carlo particles have the largest influence on the result, in agreement with similar studies of an ITER deuterium case. For the first time, numerical error bars on the entire simulated target profiles are determined showing that the largest numerical error is 17.9%, mainly due to the plasma grid discretization. On top, also numerical errors on simulated neutral pressures are investigated in detail, for which the statistical error is dominant. The analysis demonstrates which setup is needed to keep numerical errors limited: the SOLPS-ITER averaging procedure should be employed including enough EIRENE particles, and the involved grids should be sufficiently refined to reduce discretization errors.

Preliminary investigation of neutron shielding compounds for ARC-class tokamaks

Radiation shielding design is a crucial step in the development of a nuclear fusion reactor. Superconducting magnets are indeed a highly engineered component whose properties are particularly sensitive to irradiation damage. In addition, novel magnet technologies, such as High Temperature Superconductors (HTS), and the attractiveness of compact tokamak designs (like Affordable Robust Compact - ARC - reactor), require an even higher precision in the shielding design. In this work most promising radiation shielding solutions for compact tokamaks are explored. The study applies a Monte Carlo particles transport technique assisted by the OpenMC code. A D-shaped model is built and ARC reactor core is considered as main case study. A detailed analysis of the radiation environment is carried out. The study then assesses the most promising compounds on the basis of most suitable cross sections, material shielding properties and space constraints. In addition, the work addresses the problem of energy deposition on the shield providing an analytical model for active cooling necessities. Results highlighted the necessity of heavy elements and high density compounds as a considerable photon flux is generated in neutron-material interactions. Tungsten borides (WxBy) ceramics seem to be particularly promising for shielding radiation in a compact configuration.

A cell-based population control of Monte Carlo particles for the global variance reduction for transport equations

We present a population control method with sampling and regulation steps for Monte Carlo particles involved in the numerical simulation of a transport equation. We recall in the first section the difficulties related to the variance reduction methods in the general framework of transport equations; we continue with a brief presentation of the mathematical tools invoked when solving the radiative transport equations and we focus on the importance of the emission and control of existing Monte Carlo particles.The next part discusses several novel methods based on the cell-based population control method proposed in [27]. To this end, we analyze theoretically two types of splitting: one is conservative in energy (at the particle level) and the other is not. Thanks to these results, a new algorithm is introduced that uses cell-based population control and a spatial distribution. A numerical comparison of the different types of splitting is proposed in a simplified framework, then the various algorithms presented are compared against two benchmarks: the propagation of a Marshak wave and the propagation of two waves having different intensity and speed scales. To carry out these last tests, we use the multi-physics code FCI2 [13].

New variance reduction techniques for Monte-Carlo calculation based on first-collision source method

ABSTRACT A new variance-reduction method, called the first-collision source (FCS) method, is proposed to optimize Monte-Carlo source sampling to reduce variance in deep-penetration problems. In the FCS method, Monte-Carlo samples source particles from the first-collision source distribution rather than the original source distribution. The first-collision source is generated from the original source through a transport process and a scatter process, thus it has wider phase space distribution compared with the original source. Sampling from the first-collision source is capable of obtaining source particles near the aim phase space (specific space area, energy segment, and angle segment which are most related to the aim response) to increase particle popularity in aim phase space. Moreover, a new method called the FCS-CADIS method has also been proposed to further improve the sampling efficiency and performance of variance reduction. The FCS-CADIS method combines the FCS method and the Consistent Adjoint Driven Importance Sampling (CADIS) method. Both of the FCS method and the FCS-CADIS method have been implemented in the hybrid Monte-Carlo-Deterministic particle-transport code NECP-MCX. Three problems are used to assess the performance of the FCS method and the FCS-CADIS method. The results show satisfactory improvement of Monte-Carlo calculation efficiency.

Neutrino Transport with the Monte Carlo Method. II. Quantum Kinetic Equations

Abstract Neutrinos have a unique quantum feature as flavor conversions. Recent studies suggested that collective neutrino oscillations play important roles in high-energy astrophysical phenomena. The quantum kinetic equation (QKE) is capable of describing the neutrino flavor conversion, transport, and matter collision self-consistently. However, we have experienced many technical difficulties in their numerical implementation. In this paper, we present a new QKE solver based on a Monte Carlo (MC) approach. This is an upgraded version of our classical MC neutrino transport solver; in essence, a flavor degree of freedom including mixing state is added into each MC particle. This extension requires updating numerical treatments of collision terms, in particular for scattering processes. We deal with the technical problem by generating a new MC particle at each scattering event. To reduce statistical noise inherent in MC methods, we develop the effective mean free path method. This suppresses a sudden change of flavor state due to collisions without increasing the number of MC particles. We present a suite of code tests to validate these new modules with comparison to the results reported in previous studies. Our QKE-MC solver is developed with fundamentally different philosophy and design from other deterministic and mesh methods, suggesting that it will be complementary to others and potentially provide new insights into physical processes of neutrino dynamics.

Adaptive global weight window generator based on particles density uniformity for Monte Carlo particles transport simulation

Variance reduction techniques are necessary for the Monte Carlo (MC) calculations in which obtaining a detailed flux distribution for a large and complex model is required. A new method for generating mesh weight window (WW) parameters is presented, named global weight window generator (GWWG). In the method, MC calculation is performed to generate the importance of each mesh voxel in the user-set region, and then the corresponding WW parameters are obtained. The importance is related to whether the simulated particles can be uniformly transported to each mesh voxel, which is called particles density uniformity. The GWWG method does not rely on user’s design experience and can significantly improve the calculation efficiency. The tests for this method have been conducted on ITER reference neutronics models. In the calculation of ITER C-Lite model, the FoM was improved comparing to the analog FoM by the factor of 637.4.

Parareal computation of stochastic differential equations with time-scale separation: a numerical convergence study

The parareal algorithm is known to allow for a significant reduction in wall clock time for accurate numerical solutions by parallelising across the time dimension. We present and test a micro-macro version of parareal, in which the fine propagator is based on a (high-dimensional, slow-fast) stochastic microscopic model, and the coarse propagator is based on a low-dimensional approximate effective dynamics at slow time scales. At the microscopic level, we use an ensemble of Monte Carlo particles, whereas the approximate coarse propagator uses the (deterministic) Fokker-Planck equation for the slow degrees of freedom. The required coupling between microscopic and macroscopic representations of the system introduces several design options, specifically on how to generate a microscopic probability distribution consistent with a required macroscopic probability distribution and how to perform the coarse-level updating of the macroscopic probability distribution in a meaningful manner. We numerically study how these design options affect the efficiency of the algorithm in a number of situations. The choice of the coarse-level updating operator strongly impacts the result, with a superior performance if addition and subtraction of the quantile function (inverse cumulative distribution) is used. How microscopic states are generated has a less pronounced impact, provided a suitable prior microscopic state is used.

Neutrino Transport with Monte Carlo Method. I. Toward Fully Consistent Implementation of Nucleon Recoils in Core-collapse Supernova Simulations

Abstract The small energy exchange via nucleon recoils in neutrino–nucleon scattering is now supposed to be one of the important factors for successful explosions of core-collapse supernovae (CCSNe), as they can change neutrino spectra through accumulation of a large number of scatterings. In deterministic methods employed for neutrino transport in CCSN simulations, we normally cannot afford to deploy a large enough number of energy bins needed to resolve this small energy exchange, and subgrid techniques are employed one way or another. In this paper, we study quantitatively with the Monte Carlo (MC) method how well such a treatment performs. We first investigate the effects of nucleon recoils on the neutrino spectra and confirm that the average energy is reduced by ∼15% for heavy-lepton neutrinos and much smaller amounts for other types of neutrinos in a typical postbounce situation. It is also observed that the nucleon scattering dominates the electron scattering in the thermalization of neutrino spectra in all flavors. We then study possible artifacts that the coarse energy grid may produce in the deterministic methods. In order to mimic the latter calculation, we redistribute MC particles in each energy bin after a certain interval in a couple of ways and study how the results are affected, depending on the energy resolution. We also discuss the possible implications of our results for the deterministic methods.

Hadronic Light-by-Light Scattering Contribution to the Muon Anomalous Magnetic Moment from Lattice QCD.

We report the first result for the hadronic light-by-light scattering contribution to the muon anomalous magnetic moment with all errors systematically controlled. Several ensembles using 2+1 flavors of physical mass Möbius domain-wall fermions, generated by the RBC and UKQCD collaborations, are employed to take the continuum and infinite volume limits of finite volume lattice QED+QCD. We find a_{μ}^{HLbL}=7.87(3.06)_{stat}(1.77)_{sys}×10^{-10}. Our value is consistent with previous model results and leaves little room for this notoriously difficult hadronic contribution to explain the difference between the standard model and the BNL experiment.

Numerical accuracy and convergence with EMC3-EIRENE

The iterative Monte Carlo (MC) code EMC3-EIRENE is frequently used for plasma edge simulations in 3D applications. So far, a quantitative evaluation of the numerical quality of the code results remains an open issue. In this paper, we demonstrate a framework for the practical assessment of accuracy and convergence with EMC3-EIRENE. Moreover, we provide a first accuracy analysis with EMC3-EIRENE for a DIII-D divertor edge plasma case. First, we introduce post-processing averaging to efficiently reduce the variance of the statistical error. Then, we estimate the deterministic error contributions based on their theoretical reduction rates by comparing solutions with a different grid resolution, time step, or number of MC particles per iteration. Finally, using parameterized expressions for the error and the computational time, suitable numerical parameters are determined to achieve faster and/or more accurate results. We found that simulations can be more than twice as fast without losing accuracy by making use of post-processing averaging and choosing optimal parameters. In addition, we conclude that the discretization error is the dominant error contribution for the case selected in this paper, which demonstrates the importance of constructing an adequate mesh.

TURTLE: A C library for an optimistic stepping through a topography

TURTLE is a C library providing utilities allowing to navigate through a topography described by a Digital Elevation Model (DEM). The library has been primarily designed for the Monte Carlo transport of particles scattering over medium to long ranges, e.g. atmospheric muons. But, it can also efficiently handle ray tracing problems with very large DEMs (109 nodes or more), e.g. for neutrino simulations. The TURTLE library was built on an optimistic ray tracing algorithm, detailed in the present paper. This algorithm proceeds by trials and errors, approximating the topography within the modelling uncertainties of the DEM data. This allows to traverse a topography in constant time, i.e. independently of the number of grid nodes, and with no added memory. Detailed performance studies are provided by comparison to other ray tracing algorithms and as an application to muon transport in a Monte Carlo simulation. Program summaryProgram Title: The TURTLE libraryProgram Files doi:http://dx.doi.org/10.17632/pbnsctnbk8.1Licensing provisions: LGPL-3.0Programming language: C99Nature of problem: Transport long range Monte Carlo particles through a topography described by a high resolution Digital Elevation Model (DEM).Solution method: Mapping the ground surface on the fly using a heuristic requiring only the original DEM data.

Sparse-Lagrangian MMC modelling of the Sandia DME flame series

A series of turbulent, piloted dimethyl ether (DME)/air jet flames (Sandia DME flames D–G’), with Reynolds numbers ranging from 29,300 to 73,250, has been simulated using a sparse-Lagrangian multiple mapping conditioning (MMC) approach coupled to a large eddy simulation (LES) flow field solver. Mixing between the Monte-Carlo particles is modelled by a generalised form of MMC combined with a sparse distribution of particles leading to significant computational savings compared to what is required for conventional mixing models. This is achieved by pairwise mixing of particles that are selected dependent on their distance in an extended space comprised of a reference variable, given by the LES mixture fraction, and spatial location. The MMC-LES method successfully predicts the flame structure and composition field for the full flame series. Numerical results are compared against conditional statistics and spatially resolved experimental data acquired with Raman/Rayleigh scattering and laser-induced fluorescence measurements. They show good agreement even for flame DME-G’ where large turbulence-chemistry interactions lead to significant local extinction and large deviations from a flamelet structure. The influence of the mixing time on the predicted flame structure is investigated, and the systematic validation of the time scale models with the aid of measurements of the entire flame series has corroborated the findings of earlier DNS and single flame studies: a modified time scale model is needed to provide accurate predictions of conditional fluctuations and thus of possible deviations from a flamelet-like combustion regime.

Monte Carlo simulations of homogeneous nucleation and particle growth in the presence of background particles

The application of the Monte Carlo (MC) simulation technique for the modelling of nucleation processes with an existing background particle concentration is presented in this paper. Next to the nucleation of novel particles, the coagulation of an existing particle population as well as the condensational growth and evaporation of unstable particles (whose diameter is smaller than the critical Kelvin diameter) are included into the simulation. The usage of statistically weighted MC particles allows the description of particle size distribution (PSD), whose concentrations differ in several orders of magnitude. It is shown, that this approach allows to model the complex interplay between freshly nucleated particles and an existing background particle population. In this work, the nucleation of novel particles is modelled by three different nucleation theories discussed by [Girshick, S. L. and C.-P. Chiu (1990), The Journal of Chemical Physics 93], which comprise of (1) the classical nucleation theory, (2) a mathematical correction to (1) and (3) a self-consistency correction of (2). For the chosen simulation conditions, the resulting PSDs are independent of the used nucleation theory for longer simulation times, in which the simulations are described by the coagulation mechanism only. The time-frame is identified for which relevant discrepancies of the PSDs have to be taken into account.

Enforcing conservation at Monte Carlo level in a coupled finite volume—Monte Carlo simulation

Plasma edge simulations generally require coupling a finite‐volume (FV) fluid plasma discretization with a kinetic neutral Monte Carlo (MC) model. From the MC neutral simulation, source terms are estimated and inserted into the FV plasma simulation, transferring MC stochasticity to the FV simulation. While the FV simulation is based on the physical principle of conservation, most common MC source term estimators do not strictly conserve mass, momentum, and energy. Through the coupling, this also breaks the conservation in the FV simulation. Losing conservation can lead to numerical instability on the one hand and the inability to fulfil global mass, momentum, and energy balances on the other. Presently, the SOLPS‐ITER code rescales the neutral flux launched from the target to mitigate the instability. In this paper, we propose to instead independently rescale each of the source term estimations. We show that rescaling introduces a bias that depends on the number of MC particles and is negligible compared to the statistical error in the simulation. The resulting method restores all conservation properties in the FV/MC code, which will allow monitoring of the convergence through stabilization of the global quantities.