Micro-macro acceleration and multilevel Monte Carlo methods for kinetic equations

Since the 1960's people have made regular claims to the effect of: "A working nuclear fusion reactor is 20 years away!" This persistent underestimaton of the challenges in the way of building a working fusion reactor can be explained by the following dichotomy between theory and practice: On the one hand, the physical reaction is, in principle, quite simple, i.e., heat hydrogen gas to the point where it ionizes and the ions collide with sufficient kinetic energy to fuse. On the other hand, many unexpected phenomena appear in practice, hampering the maintenance of a stable plasma for sufficiently long reactor operation to produce net energy. A major challenge in tokamak reactor design is the shape of the divertor, which eliminates particles that have escaped the plasma's magnetic confinement. Due to the high cost of testing prototypes, in-silico design is the most cost-effective strategy for such components.

Mathematically, we consider an optimal design problem, constrained by a system of partial differential equations, modeling the physics in the reactor. The main component of these physics is the ion plasma, where the fusion reaction takes place. However, there is a significant population of neutral particles that interact with the plasma. Whereas a fluid approximation can often be applied to the plasma model,  the neutral particles are typically modeled with a kinetic equation, due to their relatively low density. Kinetic equations model the particle distribution over a position-velocity phase-space with twice as many dimensions as the problem at hand, making deterministic grid-based methods infeasible for three dimensional simulations. The plasma source terms, based on the neutral particle distribution are however only spatially dependent and do not suffer from this high dimensionality.

Monte Carlo simulations avoid the construction of such high-dimensional grids, by tracing stochastic particle trajectories through the phase-space. Each trajectory produces a sample from a distribution with, as expected value, an approximation of the source term being estimated. One can then average over a large number of such trajectories to compute a source term estimate. Despite avoiding the construction of a high-dimensional grid, Monte Carlo simulations are still expensive as one must often simulate many trajectories to limit the noise on the computed results. High collision rates between the neutral particles and the background plasma also increase the individual cost per trajectory. As such, these simulations are the main bottleneck in codes for simulating the tokamak plasma edge.

In this thesis we introduce a class of asymptotic-preserving multilevel Monte Carlo methods to reduce this computational cost. These methods first compute an initial estimate with large simulation time steps. We make use of an asymptotic-preserving scheme to stably perform these large time steps using diffusion as an approximation for many aggregated collisions. As these simulations are cheap, many trajectories can be simulated. The initial estimate thus has low cost, but a large bias due to the approximate nature of the trajectories. Next, we generate pairs of trajectories, one with the same coarse time steps and one with finer time steps with more kinetic behavior. The difference between the results computed by these trajectories, corresponds with an approximate sample for the bias of the initial coarse estimate. By correlating these particle pairs, the variance of their differences is small, meaning fewer fine trajectories are needed to get a low variance difference estimator. Evaluating further difference estimators with finer pairs of particles produces a hierarchy of such estimators. Combined, they produce a multilevel Monte Carlo estimate with the accuracy of the finest discretization, at significantly reduced computational cost.

After having developed efficient computational tools to evaluate designs, we wish to perform optimal design. This requires computing the gradient of an objective function, which has been evaluated, in part, using a Monte Carlo simulation. This gradient can be estimated accurately and efficiently using the discrete adjoint approach. This approach simulates the adjoint equations to the discretized PDE, which map perturbations on the residual to perturbations in the plasma-neutral state-space. In the Monte Carlo setting, the discrete adjoint Monte Carlo solver simulates the same particle trajectories as the original Monte Carlo simulation, but backwards in time. These time-reversed trajectories often have infeasible memory requirements for all but toy-problems. As such, a checkpointing approach, where snapshots of the simulation are stored at intermediate moments in time, is often applied in practice. The relevant portions of the trajectories are then recomputed from their corresponding checkpoint during the adjoint simulation, replacing the explosion in memory usage with additional computation due to recomputing paths and checkpointing overhead. In this thesis, we additionally present a novel, checkpointing free, approach for the adjoint simulation, where we generate the particle trajectories in a fully reversed manner, by reversing the underlying random number generator.