Accuracy and Convergence of Coupled Finite-Volume / Monte-Carlo Codes for Plasma Edge Simulations

Present computational techniques for coupled ﬁnite-volume/Monte-Carlo codes for plasma edge modeling un- der ITER or DEMO conditions face serious challenges with respect to computational time and accuracy. In this paper, scaling laws for different error contributions are assessed and practical procedures for error estimation are proposed. First results on a 1D and a 2D test case are discussed.


Introduction
Plasma edge simulations for plasma-neutral transport in nuclear fusion reactors are often performed using coupled finite-volume/Monte-Carlo (FV/MC) codes [1][2][3]. High recycling and detached regimes as targeted for ITER and DEMO heavily rely on increased plasma-neutral interactions to reduce the heat load and erosion of plasma-facing components to acceptable levels. As a result, plasma-neutral interactions are that strong that achieving accurate simulations within an affordable computational time is a challenging task. Moreover, to cope with the inherent statistical nature of Monte Carlo codes, different coupling techniques can be used. To assess the computational performance of different FV/MC coupling techniques, determination of simulation accuracy is indispensable. Unravelling numerical error contributions may facilitate guidance for computational effort allocation via an appropriate choice of numerical parameters.
In this paper, scaling laws for error contributions are deduced from a simple 1D plasma edge model. FV/FV results provide a reference solution to assess error contibutions from statistical noise, discretisation, finite sampling and incomplete convergence. The ultimate goal is to apply this error framework to B2-EIRENE simulations for ITER cases. To this end, practical procedures based on the derived scaling laws are proposed and tested on a 1D test case with characteristic ITER divertor conditions and for the discretisation error on a simple B2-EIRENE test case.

1D plasma edge test case
The plasma edge equations are solved on a numerical domain representing a flux tube in the poloidal cross section of a tokamak. It starts at an upstream position (x = 0) and ends at the target (x = 0.1 m). The plasma is described with the following equations for mass and parallel momentum balance to determine the ion density n i and parallel velocity u respectively: with b x is the magnetic field pitch, m i the ion mass, p = n(T i + T e ) the plasma pressure, η the parallel viscosity and S n and S m sources induced by the interaction from neutral particles with the plasma. At the upstream boundary, the density is set at n i = 10 20 m −3 , while a zero-gradient is imposed for the velocity. At the target, a sonic condition is imposed. Ion and electron temperatures for this simple test case are fixed at 40 eV throughout the domain. The neutrals are described with the following kinetic equation: In this equation, Φ n = v n f n (x, ω) represents the neutral transport flux. ω denotes the unit vector in the direction of the neutral particle velocity. It is assumed in this model that neutrals move mono-energetically in two directions. The total cross section Σ t for ionization and charge exchange is set constant according to the value at T = 40 eV. The neutral source is situated at the target with R = 0.99. Neutrals are absorped upstream and perfectly reflected at the target. For further details, the reader is referred to [6].
The fluid equations are discretized on an equidistant staggered grid with X = 100 grid cells with a finite volume approach. The kinetic equations are implemented with both a MC code and a FV code. Using the FV implementation for the kinetic equation, reference solutions without influence of MC noise are obtained. The exact, analytical solution of the system of equations is indicated with Φ EX . This solution is approximated numerically with a converged FV/FV solution on a very fine grid (9900 cells). A second solution Φ F V that includes the discretization error on the plasma grid with 100 cells but retains the fine grid for the neutral particle transport is employed for assessment of the additional error contributions in FV/MC simulations. Three coupling procedures are investigated: Correlated Sampling (CS), Random Noise (RN) and Robbins Monro (RM). To examine the errors on these solutions, they will be compared to both solutions (Φ EX and Φ F V ) obtained with the FV/FV code. In Correlated Sampling [7], the MC code uses the same seed in each iteration. Because the same sequences of random numbers are used, the particle trajectories slightly differ. The statistical noise is frozen. The simulation is deterministic once the seed is fixed and the residuals of the FV code can converge to machine accuracy (see left graph in Fig. 1). The result of one run is indicated withΦ P , where Φ can be any variable (n i , u p , n n , . . . ) and P is the number of MC particles in each MC run. Several runs are performed; each with a different seed. The expected value of all results is taken as the solution of the FV/MC system: With Random Noise [1], the MC code is run with different seeds each iteration. Due to the statistical noise, the residuals of the FV code cannot decrease to machine accuracy (see middle graph in Fig. 1). This coupling procedure is most often used. Usually, only the result of the final iteration is taken as the solution. In this paper, however, the solution is obtained from the average of many www.cpp-journal.org , with Φ P i the result of iteration i and I the total number of iterations. This way, expected values can be obtained (see section 3.3). For the Robbins Monro procedure [4], where averaged statistical values are used during the iterations, we observed that the required computational time is much higher than with Random Noise and Correlated Sampling. Therefore, it will not be discussed further.
A solution Φ obtained with the FV/MC system differs from the exact solution Φ EX . Several factors contribute to the total error. Because of the finite amount of grid cells X, a discretization error d is present. A convergence error c exists if the residuals of the plasma edge code are non-zero. The error due to the statistical noise of the finite amount of MC particles P is divided into two contributions: a statistical error s and a deterministic error or so-called finite sampling bias b . Therefore, a result Φ can be written as

Results
The statistical error can be quantified by the standard deviation σ Φ . The sample standard deviation with N the number of iterations (Random Noise) or the number of runs (Correlated Sampling) is expected to be a good measure for σ Φ . In Fig. 2 (left) the standard deviation is depicted as a function of the number of MC particles P for Φ being the L2-norm over the grid for the ion density. It is observed that σ Φ is inversely proportional to the square root of the number of MC particles per iteration σ ∝ 1/ √ P . When the mean of R runs is taken (CS), the statistical error has a standard deviation σ ∝ 1/ √ P R, according to the central limit theorem. The same proportionality applies with RN when the mean of I iterations is taken: σ ∝ 1/ √ P I. By averaging over many runs (CS) or iterations (RN), the statistical error can be made negligible compared to the deterministic errors. In order to examine the deterministic errors, this will be done for all solutions employed in the remainder of this section.
x ) 2 taken over the grid, with X the number of grid cells. Different solutions, each with a certain number P of MC particles are presented. For P > 20 the error saturates at a constant level. This is the discretization error d , which is found to be equal for all coupling techniques in the limit of infinite number of particles.
To be able to analyse the other error contributions, Φ F V is taken for comparison as it is a solution that includes the discretization error. Thus differences between Φ F V and FV/MC solutions are governed by the use of a MC method for the neutral part. It should be noted that for this reference solution the plasma variables are computed on the original plasma FV grid, whereas a much finer grid (9900 cells) is used in the neutral FV code to mimic the fact that there is no spatial discretization error in the MC code. To return to the plasma simulations, the source terms are integrated from the fine grid to the coarse. The blue lines in Fig. 2 (right) show the error for solutions obtained with different number of MC particles P . For CS, this error is the finite sampling bias originating from the non-linearities in the modeling equations. It is inversely proportional to the number of MC particles P : b ∝ 1/P , . With RN, the solutions over which is averaged are not converged, therefore, the deterministic error consists of a finite sampling bias and a convergence error. It is observed that the sum of these error is also inversely proportional to the number of MC particles P : b + c ∝ 1/P . For CS, these scaling laws indicate that the smallest error for a given amount of computational time is obtained with one run with as many MC particles as affordable, rather than multiple runs with fewer MC particles. Indeed, assuming that computational time scales with ∝ RP , the same statistical error is obtained (∝ 1/ √ RP ). However, with high P , the finite sampling bias (∝ 1/P ) will be lower than with low P . For RN, it can be seen that averaging the results over iterations is recommended. Then the statistical error can be made negligible and only a much lower deterministic error remains. Taking the result from the final iteration only, results in a much higher statistical error.

Test case description
Eventually, the aim is to estimate the errors for ITER cases with B2-EIRENE. To advance in this direction, an error assessment based on the scaling laws is performed for a 1D case with ITER characteristics. To test the feasibility of these methods for B2-EIRENE, the discretization error is also assessed for a simple B2-EIRENE case. The 1D test case uses data from a B2-EIRENE simulation of a partially detached ITER-F12 case (Fig. 3 (left)).
The figure in the middle shows the flux line studied on the original grid (red curve). It starts at a poloidal distance of 1.069 m upstream from the right divertor target. Continuity and momentum equations are solved, while profiles from B2-EIRENE are imposed (Fig. 3 (right)). Errors for the target flux are examined. For CS, one simulation Φ CS is run with P = 10 6 MC particles per MC run (approximately 20 MC runs were performed). For RN, a solution Φ RN is obtained with P = 2 × 10 3 MC particles per iteration (10 4 iterations were performed). Both solutions require approximately the same computational time. In the next subsections the numerical errors of these solutions are assessed.
To test the procedure for the discretization error in B2-EIRENE a simple toroidal slab geometry (100×20 cells) with top and bottom targets under low recycling coefficients is used (separatrix at n i = 3 × 10 19 m −3 , T i = T e = 400 eV and wall conditions at n i = 10 17 m −3 , T i = T e = 2 eV and R = 0.3). For this case the density profile between the targets in the middle of the domain as well as the volume averaged density over the domain are studied. Reference cases are run with 140,000 MC particles.

Discretization errors
In FV codes, the discretization error can be estimated by comparing solutions computed with different grid spacings [5] a coarse grid n 2 h, a finer grid nh and a very fine grid h. The discretization error on the finest grid is then estimated from Richardson's extrapolation, assuming the error d scales with (Δx) p . This leads to . Given the use of upwind schemes for the discretization of convection terms in the plasma equations p is expected to be close to one. Therefore, the value for p is taken as a criterion to www.cpp-journal.org check grid convergence. This method provides good results for the simple 1D edge case solved with the FV/FV code. However, it needs to be adapted when using a FV/MC code in order to cope with the noise introduced by MC. A first approach is to compare solutions of three different grid spacings (n = 3 is used). For each grid, the solution Φ is calculated as the mean of R solutions obtained using correlated sampling. The individual runs are converged sufficiently long to achieve negligible convergence errors. To reduce the influence of the finite number of MC particles, a high number of Monte Carlo particles (P = 10 4 for the simple 1D case of section 2) is used. This method works well for the 1D case and is in principle possible with any FV-MC code. However, it easily becomes impractical: many runs with many MC particles are needed for an accurate estimate.
To determine the discretization error in more elaborate (2D) cases, we propose to use simulations that only converge the FV plasma code further from the FV/MC results under investigation. Sources are taken from the last iteration of the FV/MC run on a coarse grid and interpolated for finer grids using a piecewise linear approach while keeping the integrated sources constant. The discretization error is then computed with Richardson's equation. For the 1D ITER case, this results in a discretization error on the target flux of d ≈ 2.7 × 10 22 m −2 s −1 (relative error of 1.3%). Figure 4 shows discretization errors for density profiles for the 1D ITER case (left) and the B2-EIRENE slab case (right). For the 1D ITER case a relative error between 3 and 15% is obtained. For the 2D case relative errors for n i are approximately 4%. Locations with p < 0.8 are indicated with a cross. It can be concluded that further grid refinement near the target is needed to achieve valuable estimates at these positions.

Statistical error
The statistical error can be estimated via the sample standard deviation s. When using CS simulations, only one run with a high amount of MC particles P h is recommended. However, while 1 run with a high amount of MC particles minimizes the variance, it does not allow this variance to be accurately estimated. Using the scaling laws and at least 30 runs with a lower amount of MC particles P l , the statistical error of the CS simulation with high P becomes: s ≈ σ P h ≈ s P l P l P h , with s P l the sample standard deviation of the runs with a low amount of MC particles P l . For 1D ITER case, 200 additional runs with 1 MC particle have been executed, which results in an estimation of the statistical error for Φ CS of s = 2.8 × 10 21 m −2 s −1 .
When only the final RN solution is taken as the solution, it is expected that the statistical error is dominant. The latter can be estimated as the sample standard deviation s of the iterations. For the target flux in the 1D ITER case this results in s = 5.8 × 10 21 m −2 s −1 . For the B2-EIRENE case the densities have a statistical error ranging from 0.1 to 0.02%. However, it is recommended to average the iterations until the statistical error s becomes negligible. This can be determined by monitoring the residual of the average (see right graph in Fig. 1): when s is dominant, this residual decreases. When s is negligible, this residual remains constant. A similar result is obtained for the 2D case.

Finite sampling bias and convergence error
To estimate the remaining deterministic errors b + c (with c ≈ 0 for CS), we elaborate a deviation similar to Richardson's extrapolation. Starting from the scaling law as obtained in section 2, i.e. b + c ∝ 1/P , 3 different solutions Φ P , Φ P f and Φ P f 2 with a decreasing amount of MC particles (P , fP and f 2 P ) are used. This results in: Also in this case the statistical noise hampers the evaluation of equation (3). However, this cannot be circumvented as this variance is inherently associated with the bias error b . In a first approach again a sufficiently low statistical error s can be aimed at. This is checked for the 1D case using additional runs with P = 1, 2 and 4. For RN, this method works well as statistical errors can be decreased by increasing the number of iterations. In the case of CS, it requires many runs to have a sufficiently low statistical error estimate. A more practical procedure obtained by choosing a very low reduction factor f 1 and statistically testing for p = 1. When p = 1 is acceptable, In addition to the solution Φ P only a few other solutions with an intermediate amount of MC particles Φ P l have to be acquired for testing the acceptance of p = 1. Finally, it should be checked that b + c > s for this additional solution. It can be shown that for f < 0.1 this condition will be either reached or the error on b + c will be statistically small enough to show that s is the dominant error. For the 1D test case with CS, the fast estimate using again additional 200 runs with P = 1 the resulting bias on the target flux amounts to b = 7.6 · 10 17 m −2 s −1 . For Φ RN the estimates with equation (3) and with the second approach give simular results and lead for the target flux to b + c = 1.5 × 10 20 m −2 s −1 for an iteration averaged solution. The error due to statistical noise is approximately 40 times lower than in case the final iteration of the RN simulation is taken. Further, it it is noticed that for this 1D ITER case the averaging RN simulation ( = 1.5 × 10 20 m −2 s −1 ) is approximately 20 times more accurate than Φ CS ( = 2.8 × 10 21 m −2 s −1 ). Finally it should be noted that this procedure is also implemented in B2-EIRENE, leading to similar results.

Conclusion
A procedure to estimate the numerical error contributions is elaborated based on the scaling laws. With Random Noise, a significant reduction of the error induced by statistical noise can be achieved by averaging over iterations instead of taking only the final solution. Another option is to run with Correlated Sampling. In this case, the optimal solution is to execute one run with as many MC particles as you can afford instead of averaging over multiple runs with a lower amount of MC particles. It is shown that the discretization error can be estimated with grid refinement in the plasma equations only. A similar refinement procedure, using simulations with a decreasing amount of MC particles, is elaborated for the finite sampling bias and the convergence error. Finally, the statistical error can be assessed with the standard deviation. The assessments were elaborated based on a simple test case and next used to assess error contributions for a 1D test case under ITER divertor conditions at one hand and for a simple B2-EIRENE test case on the other. Future work will focus on its application for B2-EIRENE ITER simulations and possible speed-up while keeping the required accuracy.