Giovanni Samaey

Frederic Legoll, Tony Lelievre, Giovanni Samaey

We introduce a micro-macro parareal algorithm for the time-parallel integration of multiscale-in-time systems. The algorithm first computes a cheap, but inaccurate, solution using a coarse propagator (simulating an approximate slow macroscopic model), which is iteratively corrected using a fine-scale propagator (accurately simulating the full microscopic dynamics). This correction is done in parallel over many subintervals, thereby reducing the wall-clock time needed to obtain the solution, compared to the integration of the full microscopic model. We provide a numerical analysis of the algorithm for a prototypical example of a micro-macro model, namely singularly perturbed ordinary differential equations. We show that the computed solution converges to the full microscopic solution (when the parareal iterations proceed) only if special care is taken during the coupling of the microscopic and macroscopic levels of description. The convergence rate depends on the modeling error of the approximate macroscopic model. We illustrate these results with numerical experiments.

Mathias Rousset, Giovanni Samaey

We discuss variance reduced simulations for an individual-based model of chemotaxis of bacteria with internal dynamics. The variance reduction is achieved via a coupling of this model with a simpler process in which the internal dynamics has been replaced by a direct gradient sensing of the chemoattractants concentrations. In the companion paper \cite{limits}, we have rigorously shown, using a pathwise probabilistic technique, that both processes converge towards the same advection-diffusion process in the diffusive asymptotics. In this work, a direct coupling is achieved between paths of individual bacteria simulated by both models, by using the same sets of random numbers in both simulations. This coupling is used to construct a hybrid scheme with reduced variance. We first compute a deterministic solution of the kinetic density description of the direct gradient sensing model; the deviations due to the presence of internal dynamics are then evaluated via the coupled individual-based simulations. We show that the resulting variance reduction is \emph{asymptotic}, in the sense that, in the diffusive asymptotics, the difference between the two processes has a variance which vanishes according to the small parameter.

Mathias Rousset, Giovanni Samaey

We discuss velocity-jump models for chemotaxis of bacteria with an internal state that allows the velocity jump rate to depend on the memory of the chemoattractant concentration along their path of motion. Using probabilistic techniques, we provide a pathwise result that shows that the considered process converges to an advection-diffusion process in the (long-time) diffusion limit. We also (re-)prove using the same approach that the same limiting equation arises for a related, simpler process with direct sensing of the chemoattractant gradient. Additionally, we propose a time discretization technique that retains these diffusion limits exactly, i.e., without error that depends on the time discretization. In the companion paper \cite{variance}, these results are used to construct a coupling technique that allows numerical simulation of the process with internal state with asymptotic variance reduction, in the sense that the variance vanishes in the diffusion limit.

Tony Lelièvre, Giovanni Samaey, Przemysław Zieliński

We analyse convergence of a micro-macro acceleration method for the Monte Carlo simulation of stochastic differential equations with time-scale separation between the (fast) evolution of individual trajectories and the (slow) evolution of the macroscopic function of interest. We consider a class of methods, presented in [Debrabant, K., Samaey, G., Zieliński, P. A micro-macro acceleration method for the Monte Carlo simulation of stochastic differential equations. SINUM, 55 (2017) no. 6, 2745-2786], that performs short bursts of path simulations, combined with the extrapolation of a few macroscopic state variables forward in time. After extrapolation, a new microscopic state is then constructed, consistent with the extrapolated variable and minimising the perturbation caused by the extrapolation. In the present paper, we study a specific method in which this perturbation is minimised in a relative entropy sense. We discuss why relative entropy is a useful metric, both from a theoretical and practical point of view, and rigorously study local errors and numerical stability of the resulting method as a function of the extrapolation time step and the number of macroscopic state variables. Using these results, we discuss convergence to the full microscopic dynamics, in the limit when the extrapolation time step tends to zero and the number of macroscopic state variables tends to infinity.

Ward Melis, Thomas Rey, Giovanni Samaey

We present high-order, fully explicit projective integration schemes for nonlinear collisional kinetic equations such as the BGK and Boltzmann equation. The methods first take a few small (inner) steps with a simple, explicit method (such as direct forward Euler) to damp out the stiff components of the solution. Then, the time derivative is estimated and used in an (outer) Runge-Kutta method of arbitrary order. The procedure can be recursively repeated on a hierarchy of projective levels to construct telescopic projective integration methods. Based on the spectrum of the linearized collision operator, we deduce that the computational cost of the method is essentially independent of the stiffness of the problem: with an appropriate choice of inner step size, the time step restriction on the outer time step, as well as the number of inner time steps, is independent of the stiffness of the (collisional) source term. In some cases, the number of levels in the telescopic hierarchy depends logarithmically on the stiffness. We illustrate the method with numerical results in one and two spatial dimensions.

Ward Melis, Thomas Rey, Giovanni Samaey

We present a high-order, fully explicit, asymptotic-preserving projective integration scheme for the nonlinear BGK equation. The method first takes a few small (inner) steps with a simple, explicit method (such as direct forward Euler) to damp out the stiff components of the solution. Then, the time derivative is estimated and used in an (outer) Runge-Kutta method of arbitrary order. Based on the spectrum of the linearized BGK operator, we deduce that, with an appropriate choice of inner step size, the time step restriction on the outer time step as well as the number of inner time steps is independent of the stiffness of the BGK source term. We illustrate the method with numerical results in one and two spatial dimensions.

Pieter Vanmechelen, Geert Lombaert, Giovanni Samaey

We compare three random field discretization strategies for probabilistic identification of spatially varying material parameters in high-resolution finite element models. These strategies are (i) the Karhunen-Loève expansion, (ii) a wavelet expansion, and (iii) local average subdivision. The methods are assessed in the context of multilevel Markov chain Monte Carlo applied to plane stress elasticity with high-resolution displacement observations. Emphasis is placed on numerical efficiency, initialization cost, Markov chain mixing, and cost-to-error behaviour as the discretization resolution increases. While all approaches yield comparable posterior estimates, significant differences are observed in multilevel variance reduction and sampling efficiency. In particular, local average subdivision exhibits improved mixing and lower cost-to-error ratios at fine resolutions, despite its higher nominal parameter dimension. The results provide practical guidance for selecting stochastic field representations in uncertainty quantification in finite element simulations of heterogeneous materials.

Hannes Vandecasteele, Przemysław Zieliński, Giovanni Samaey

We discuss through multiple numerical examples the accuracy and efficiency of a micro-macro acceleration method for stiff stochastic differential equations (SDEs) with a time-scale separation between the fast microscopic dynamics and the evolution of some slow macroscopic state variables. The algorithm interleaves a short simulation of the stiff SDE with extrapolation of the macroscopic state variables over a longer time interval. After extrapolation, we obtain the reconstructed microscopic state via a matching procedure: we compute the probability distribution that is consistent with the extrapolated state variables, while minimally altering the microscopic distribution that was available just before the extrapolation. In this work, we numerically study the accuracy and efficiency of micro-macro acceleration as a function of the extrapolation time step and as a function of the chosen macroscopic state variables. Additionally, we compare the effect of different hierarchies of macroscopic state variables. We illustrate that the method can take significantly larger time steps than the inner microscopic integrator, while simultaneously being more accurate than approximate macroscopic models.

Afroja Parvin, Giovanni Samaey, Julian Koellermeier

We derive the first-moment Shallow Water Exner Moment model with sediment entrainment and deposition (SWEMED1) and show that the full source term has a fully-settled water-at-rest equilibrium manifold. We prove that the model is only weakly hyperbolic at this equilibrium, which prevents the use of Yong's structural stability framework. However, a linear spectral analysis and numerical results do not indicate instability. Based on numerical results, we introduce a fast-slow scaling of the source term, and for the fast limit, we derive a new suspended water-at-rest equilibrium manifold, which has a different structure but is still only weakly hyperbolic. Our results show that the remaining obstruction is linked to the transport closure of the SWEMED1, and we give a constructive direction for the derivation of new closures leading to models with more desirable analytical properties.

Maarten Volkaerts, Marie Cloet, Hans Dierckx et al.

We present a Bayesian approach to estimate the parameters of mathematical models of cardiac electrophysiology with quantified uncertainty. Such models capture the dynamics of the electrical signal that coordinates the muscle cell contraction in the heart wall and can support cardiac arrhythmia treatment. We consider an illustrative case motivated by a cardiac arrhythmia, namely, by intra-atrial reentrant tachycardia. We estimate a low-dimensional geometrical parameter that describes the boundary of an electrically non-conducting region in the heart tissue from synthetic electrical measurements outside of the tissue. Instead of relying on a deterministic fit for this region, we estimate a posterior distribution on the geometrical parameter using Bayesian inference that captures the uncertainty due to measurement errors. We propose a likelihood based on a set of quantities that characterize the data for improved accuracy. To efficiently approximate the posterior distribution, we propose a compressed likelihood function and an adapted Metropolis-Hastings (MH) algorithm. We obtain an algorithm that strongly decreases the number of samples by using an adaptive proposal strategy. Our algorithm also gives attention to the impact of discretization errors on inference outcomes, as these introduce artificial discontinuities in the posterior if not properly addressed. We account for discretization errors in the likelihood and in the accept-reject step of our adapted MH algorithm to improve the robustness of our estimates and to further increase the sampling efficiency. All of these elements combined give us a method that efficiently estimates the non-conducting parameters with uncertainty. We perform several experiments with different amounts of measurement noise and illustrate how this translates into the posterior distributions.

Oskar Lappi, Emil Løvbak, Thijs Steel et al.

Particle-based kinetic Monte Carlo simulations of neutral particles is one of the major computational bottlenecks in tokamak scrape-off layer simulations. This computational cost comes from the need to resolve individual collision events in high-collisional regimes. However, in such regimes, one can approximate the high-collisional kinetic dynamics with computationally cheaper diffusion. Asymptotic-preserving schemes make use of this limit to perform simulations in these regimes, without a blow-up in computational cost as incurred by standard kinetic approaches. One such scheme is Kinetic-diffusion Monte Carlo. In this paper, we present a first extension of this scheme to the two-dimensional setting and its implementation in the Eiron particle code. We then demonstrate that this implementation produces a significant speedup over kinetic simulations in high-collisional cases.

Afroja Parvin, Giovanni Samaey, Julian Koellermeier

In this paper, we develop the Hyperbolic Shallow Water Exner Moment model with Erosion and Deposition (HSWEMED), extending the shallow water moment framework to capture coupled morphodynamics with erosion and deposition. HSWEMED introduces a suspended-sediment concentration equation, couples concentration-dependent mixture density with the momentum and higher-order moment equations, and includes source terms due to erosion and deposition. Starting from the incompressible Navier-Stokes equations for a water-sediment mixture, we derive a coupled system consisting of the shallow water equations, moment equations for polynomial velocity coefficients, a depth-averaged suspended-sediment equation, and an Exner equation for bedload transport with erosion-deposition coupling. Although the transported scalar is depth-averaged, we reconstruct a low-order vertical concentration profile consistent with the moment representation of velocity, providing the near-bed concentration needed in the closure. We prove hyperbolicity through hyperbolic regularization and derive dissipative energy balance relations for lower-order models. Numerical results are obtained with a path-conservative finite-volume scheme based on a Lax-Friedrichs-type flux. Several dam-break tests, including wet/dry front cases, are validated against laboratory experiments, showing improved accuracy over existing shallow water moment models. The proposed HSWEMED provides a mathematically well-posed and computationally efficient framework for morphodynamic simulations.

Kristian Debrabant, Giovanni Samaey

We present and analyze a micro/macro acceleration technique for the Monte Carlo simulation of stochastic differential equations (SDEs) in which there is a separation between the (fast) time-scale on which individual trajectories of the SDE need to be simulated and the (slow) time-scale on which we want to observe the (macroscopic) function of interest. The method performs short bursts of microscopic simulation using an ensemble of SDE realizations, after which the ensemble is restricted to a number of macroscopic state variables. The resulting macroscopic state is then extrapolated forward in time and the ensemble is projected onto the extrapolated macroscopic state. We provide a first analysis of its convergence in terms of extrapolation time step and number of macroscopic state variables. The effects of the different approximations on the resulting error are illustrated via numerical experiments.

Zhirui Tang, Niels Horsten, Giovanni Samaey

Neutrals in the plasma edge are commonly modeled by kinetic equations, with quantities of interest given by macroscopic quantities such as density, velocity, and temperature. In reactor-relevant regimes, fully kinetic descriptions solved by Monte Carlo (MC) methods, although accurate, become computationally expensive, whereas fluid-limit approximations are computationally more efficient but may lose accuracy due to boundary effects or low-collisional regimes. Hybrid fluid-kinetic approaches aim to combine the strengths of both descriptions. However, existing simulation methods face challenges, including interface handling in domain decomposition, unphysical assumptions, and iterative coupling in distribution decomposition. In this work, we propose a distribution-decomposition hybrid model constructed at the particle level based on the kinetic-diffusion Monte Carlo (KDMC) method. The model inherits key properties of KDMC: it is asymptotic-preserving and does not require iterative coupling between the fluid and kinetic components. To improve the accuracy of the fluid-part quantities estimation, a Navier-Stokes-type fluid system is derived via Hilbert-Chapman-Enskog expansions, tailored for KDMC. In the considered one-dimensional tests, the resulting fluid system has comparable accuracy to the AFN model used in SOLPS-ITER while requiring substantially fewer nonlinear iterations. Additionally, a tunable reflective boundary condition is introduced that allows balancing accuracy and efficiency. The model exhibits at least 500 times speedup over the kinetic MC, while maintaining relative L2 errors around 10% in a charge exchange (CX)-dominant test case. In non-CX-dominant regimes, the accuracy becomes increasingly sensitive to boundary treatment due to the inherent limitations of the fluid approximation near the boundary, motivating further refinement of the KDMC boundary conditions.

Klaas Willems, Erik Arlemark, Giovanni Samaey et al.

Polyatomic kinetic models are essential for accurately capturing the thermodynamic behavior of real gases, as internal energy modes significantly influence transport coefficients, relaxation processes, and non-equilibrium effects that cannot be represented by monoatomic models. The polyatomic ESBGK model describes molecular collisions as a relaxation towards a generalized Gaussian distribution with an anisotropic covariance matrix and an exponentially decaying internal energy distribution. We present a new semi-Lagrangian scheme for the polyatomic Ellipsoidal Statistical BGK (ESBGK) model of the Boltzmann equation. The semi-Lagrangian framework, being deterministic and grid-based, removes the time-step restriction associated with the linear transport term by following the method of characteristics. The potentially stiff relaxation term is treated using an implicit A-stable linear multistep method which, owing to the structure of the BGK operator, can be reformulated into a cheap time-stepping scheme. This yields a highly efficient and numerically stable method. The numerical method is asymptotic preserving and stiffly accurate, meaning the scheme asymptotically converges to a scheme for the Euler equations in the vanishing Knudsen limit. In addition, we prove that the first-order scheme, asymptotically converges to the compressible Navier-Stokes equation with correct transport coefficients. Finally, we propose inflow and outflow boundary conditions suitable for BGK-type kinetic equations. We perform simulations of the Fourier and Couette test case to compare the BGK model with Direct Simulation Monte Carlo (DSMC). To conclude, we demonstrate the method on a challenging orifice flow test case with moving boundaries.

William De Deyn, Michael Herty, Giovanni Samaey

We study Consensus-Based Optimization (CBO) for two-layer neural network training. We compare the performance of CBO against Adam on two test cases and demonstrate how a hybrid approach, combining CBO with Adam, provides faster convergence than CBO. Additionally, in the context of multi-task learning, we recast CBO into a formulation that offers less memory overhead. The CBO method allows for a mean-field limit formulation, which we couple with the mean-field limit of the neural network. To this end, we first reformulate CBO within the optimal transport framework. In the limit of infinitely many particles, we define the corresponding dynamics on the Wasserstein-over-Wasserstein space and show that the variance decreases monotonically.

G. Dimarco, L. Pareschi, G. Samaey

We develop a new Monte Carlo method that solves hyperbolic transport equations with stiff terms, characterized by a (small) scaling parameter. In particular, we focus on systems which lead to a reduced problem of parabolic type in the limit when the scaling parameter tends to zero. Classical Monte Carlo methods suffer of severe time step limitations in these situations, due to the fact that the characteristic speeds go to infinity in the diffusion limit. This makes the problem a real challenge, since the scaling parameter may differ by several orders of magnitude in the domain. To circumvent these time step limitations, we construct a new, asymptotic-preserving Monte Carlo method that is stable independently of the scaling parameter and degenerates to a standard probabilistic approach for solving the limiting equation in the diffusion limit. The method uses an implicit time discretization to formulate a modified equation in which the characteristic speeds do not grow indefinitely when the scaling factor tends to zero. The resulting modified equation can readily be discretized by a Monte Carlo scheme, in which the particles combine a finite propagation speed with a time-step dependent diffusion term. We show the performance of the method by comparing it with standard (deterministic) approaches in the literature.

Ward Melis, Giovanni Samaey

We study a variance reduction strategy based on control variables for simulating the averaged macroscopic behavior of a stochastic slow-fast system. We assume that this averaged behavior can be written in terms of a few slow degrees of freedom, and that the fast dynamics is ergodic for every fixed value of the slow variable. The time derivative for the averaged dynamics can then be approximated by a Markov chain Monte Carlo method. The variance-reduced scheme that is introduced here uses the previous time instant as a control variable. We analyze the variance and bias of the proposed estimator and illustrate its performance when applied to a linear and nonlinear model problem.

Pauline Lafitte, Ward Melis, Giovanni Samaey

We present a general, high-order, fully explicit relaxation scheme which can be applied to any system of nonlinear hyperbolic conservation laws in multiple dimensions. The scheme consists of two steps. In a first (relaxation) step, the nonlinear hyperbolic conservation law is approximated by a kinetic equation with stiff BGK source term. Then, this kinetic equation is integrated in time using a projective integration method. After taking a few small (inner) steps with a simple, explicit method (such as direct forward Euler) to damp out the stiff components of the solution, the time derivative is estimated and used in an (outer) Runge-Kutta method of arbitrary order. We show that, with an appropriate choice of inner step size, the time step restriction on the outer time step is similar to the CFL condition for the hyperbolic conservation law. Moreover, the number of inner time steps is also independent of the stiffness of the BGK source term. We discuss stability and consistency, and illustrate with numerical results (linear advection, Burgers' equation and the shallow water and Euler equations) in one and two spatial dimensions.

Ward Melis, Giovanni Samaey

We study a general, high-order, fully explicit numerical method for simulating kinetic equations with a BGK-type collision model with multiple relaxation times. In that case, the problem is stiff and its spectrum consists of multiple separated eigenvalue clusters. Projective integration methods are explicit integration schemes that first take a few small (inner) steps with a simple, explicit method, after which the solution is extrapolated forward in time over a large (outer) time step. These are very efficient schemes, provided there are only two clusters of eigenvalues. Telescopic projective integration methods generalize the idea of projective integration methods by constructing a hierarchy of projective levels. Here, we show how telescopic projective integration methods can be used to efficiently integrate kinetic equations with multiple relaxation times. We show that the required number of projective levels depends on the number of clusters, which in turn depends on the stiffness of the BGK source term. The size of the outer level time step only depends on the slowest time scale present in the model and is independent of the stiffness of the problem. We discuss stability and illustrate the approach with simulations in one and two spatial dimensions.