Giovanni Samaey

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.

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.

Pauline Lafitte, Giovanni Samaey

We investigate a projective integration scheme for a kinetic equation in the limit of vanishing mean free path, in which the kinetic description approaches a diffusion phenomenon. The scheme first takes a few small steps with a simple, explicit method, such as a spatial centered flux/forward Euler time integration, and subsequently projects the results forward in time over a large time step on the diffusion time scale. We show that, with an appropriate choice of the inner step size, the time-step restriction on the outer time step is similar to the stability condition for the diffusion equation, whereas the required number of inner steps does not depend on the mean free path. We also provide a consistency result. The presented method is asymptotic-preserving, in the sense that the method converges to a standard finite volume scheme for the diffusion equation in the limit of vanishing mean free path. The analysis is illustrated with numerical results, and we present an application to the Su-Olson test.

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.

Kristian Debrabant, Giovanni Samaey, Przemysław Zieliński

We present and analyse a micro-macro acceleration method for the Monte Carlo simulation of stochastic differential equations with separation between the (fast) time-scale of individual trajectories and the (slow) time-scale of the macroscopic function of interest. The algorithm combines short bursts of path simulations with extrapolation of a number of macroscopic state variables forward in time. The new microscopic state, consistent with the extrapolated variables, is obtained by a matching operator that minimises the perturbation caused by the extrapolation. We provide a proof of the convergence of this method, in the absence of statistical error, and we analyse various strategies for matching, as an operator on probability measures. Finally, we present numerical experiments that illustrate the effects of the different approximations on the resulting error in macroscopic predictions.

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.

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.

Klaas Willems, Axel Klar, Giovanni Russo et al.

We present a numerical method for simulating rarefied gases that interact with moving boundaries and rigid bodies. The gas is described by the BGK equation in Lagrangian form and solved using an Arbitrary Lagrangian-Eulerian method, in which grid points move with the local mean velocity of the gas. The main advantage of the moving grid is that the algorithm can deal well with cases where the domain boundaries are time-dependent and the simulation domain contains rigid objects. Due to the irregular nature of the grid, we use a novel meshless MUSCL-like Moving Least Squares Method (MLS) for spatial discretisation coupled with a higher-order Implicit-Explicit Runge-Kutta method. To avoid spurious oscillations at discontinuities, we use the so-called Multi-dimensional Optimal Order Detection (MOOD) method with an adapted criterion to relax the discrete maximum property. Finally, we employ a new implementation of the boundary conditions that requires no iterative or extrapolation procedure. The method achieves fourth-order in 1D and second-order in 2D for simulations with moving boundaries. We demonstrate the method's effectiveness on classical test cases such as the driven square cavity, shear layer, and shock tube.

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.

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.

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.

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.

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.

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.

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.