Particle conservation in numerical models of the tokamak plasma edge

TL;DR

This paper discusses numerical diagnostics and algorithms to ensure particle conservation in Monte-Carlo models of tokamak edge plasma, addressing errors caused by statistical methods and residuals, and proposing solutions for efficient, accurate simulations.

Contribution

It introduces diagnostics and algorithms to identify and reduce particle balance errors in Monte-Carlo tokamak edge models, enabling larger time-steps without sacrificing accuracy.

Findings

Diagnostics can unambiguously identify residual-induced errors.

Algorithms allow large time-steps while maintaining particle conservation.

Techniques improve the accuracy and efficiency of plasma edge modeling.

Abstract

The test particle Monte-Carlo models for neutral particles are often used in the tokamak edge modelling codes. The drawback of this approach is that the self-consistent solution suffers from random error introduced by the statistical method. A particular case where the onset of nonphysical solutions can be clearly identified is violation of the global particle balance due to non-converged residuals. There are techniques which can reduce the residuals - such as internal iterations in the code B2-EIRENE - but they may pose severe restrictions on the time-step and slow down the computations. Numerical diagnostics described in the paper can be used to unambiguously identify when the too large error in the global particle balance is due to finite-volume residuals, and their reduction is absolutely necessary. Algorithms which reduce the error while allowing large time-step are also discussed.

Tables

Table 1. Table 1: Relative errors of the particle and power balances in the model runs of Section  IV

case	$\Delta P,\%$	$\Delta {\Gamma }^D,\%$	$\Delta {\Gamma }^{He},\%$	$\Delta {\Gamma }^C,\%$
$m$=20, $\Delta t$=3e-7	0.74	1.39	6.77	0.023
$m$=1, $\Delta t$=1e-4	0.32	91.5	99.3	4.6
$m$=1+99, $\Delta t$=1e-4∗	1.8	16.1	14.0	0.53
∗this case is discussed in Section V

Table 2. Table 2: Individual terms of the error in particle balance (in %), Equation ( 10 ), model runs of Section  IV

case	${\Delta }_R^D$	${\Delta }_S^D$	${\Delta }_T^D$	${\Delta }_R^{He}$	${\Delta }_S^{He}$	${\Delta }_T^{He}$
$m$=20, $\Delta t$=3e-7	0.02	0.14	-1.54	3.10	0.67	-10.63
$m$=1, $\Delta t$=1e-4	91.6	1.06	-1.98	100.7	-0.11	0.15
$m$=1+99, $\Delta t$=1e-4∗	0.17	16.3	0.19	-0.10	4.52	8.92
∗this case is discussed in Section V

Equations

$$\mathrm{F}\left(\varphi \right)=S\left(\varphi \right),\varphi =\left\{n_{\alpha },u_{\alpha },T_e,T_i\right\}$$

$$\mathrm{F}\left({\varphi }_k\right)=\tilde{S}\left({\varphi }_{k-1}\right)+D\left({\varphi }_k,{\varphi }_{k-1}\right)$$

$$R=\tilde{S}\left({\varphi }_k^m|{\varphi }_{k-1}\right)+D\left({\varphi }_k^m,{\varphi }_{k-1}\right)-\mathrm{F}\left({\varphi }_k^m\right)$$

$$\mathrm{F}\left({\varphi }_k^m\right)=\tilde{S}\left({\varphi }_k^m|{\varphi }_{k-1}\right)+D\left({\varphi }_k^m,{\varphi }_{k-1}\right)+R$$

$${\tilde{S}}_{\alpha }^s\left({\varphi }_k^j|{\varphi }_{k-1}\right)=\frac{Q_{\beta }^s\left({\varphi }_k^j\right)}{Q_{\beta }^s\left({\varphi }_k^0\right)}{\tilde{S}}_{\alpha }^s\left({\varphi }_{k-1}\right)$$

$$\mathrm{F}\left({\varphi }_k^m\right)=\tilde{S}\left({\varphi }_k^m\right)+$$

$$+\left[R+\tilde{S}\left({\varphi }_k^m|{\varphi }_{k-1}\right)-\tilde{S}\left({\varphi }_k^m\right)+D\left({\varphi }_k^m,{\varphi }_{k-1}\right)\right]$$

$$\sum _i\sum _{{\alpha }^{\prime }}\left[{\tilde{S}}_i^{{\alpha }^{\prime }}\left({\varphi }_k^m\right)-{\mathrm{F}}_i^{{\alpha }^{\prime }}\left({\varphi }_k^m\right)\right]=\sum _i\sum _{{\alpha }^{\prime }}\left[-R_i^{{\alpha }^{\prime }}+{\tilde{S}}_i^{{\alpha }^{\prime }}\left({\varphi }_k^m\right)-{\tilde{S}}_i^{{\alpha }^{\prime }}\left({\varphi }_k^m|{\varphi }_{k-1}\right)-D_i^{{\alpha }^{\prime }}\left({\varphi }_k^m,{\varphi }_{k-1}\right)\right]$$

$$\Delta {\Gamma }^{\beta }=\frac{{\Gamma }_{puff}^{\beta }+{\Gamma }_{core}^{\beta }+{\Gamma }_{spt}^{\beta }-{\Gamma }_{pump}^{\beta }-{\Gamma }_{leak}^{\beta }}{{\Gamma }_{puff}^{\beta }+{\Gamma }_{core}^{\beta }+{\Gamma }_{spt}^{\beta }}$$

$${\Gamma }_{in}^{\beta }={\Gamma }_{puff}^{\beta }+{\Gamma }_{core}^{\beta }+{\Gamma }_{spt}^{\beta }$$

$${\Delta }^{\beta }={\Delta }_R^{\beta }+{\Delta }_S^{\beta }+{\Delta }_T^{\beta }$$

$${\Delta }_R^{\beta }=\frac{-{\sum }_i{\sum }_{{\alpha }^{\prime }}{R}_i^{{\alpha }^{\prime }}}{{\Gamma }_{in}^{\beta }}$$

$${\Delta }_S^{\beta }=\frac{{\sum }_i{\sum }_{{\alpha }^{\prime }}\left[{\tilde{S}}_i^{{\alpha }^{\prime }}\left({\varphi }_k^m\right)-{\tilde{S}}_i^{{\alpha }^{\prime }}\left({\varphi }_k^m|{\varphi }_{k-1}\right)\right]}{{\Gamma }_{in}^{\beta }}$$

$${\Delta }_T^{\beta }=\frac{1}{{\Gamma }_{in}^{\beta }}\sum _i\sum _{{\alpha }^{\prime }}\frac{n_{k-1}^{{\alpha }^{\prime },i}-n_k^{{\alpha }^{\prime },i}}{\Delta t}$$

$$\sum _i\sum _{{\alpha }^{\prime }}{\mathrm{F}}_i^{{\alpha }^{\prime }}\left({\varphi }_k^m\right)={\Gamma }_{out}^{\beta +}-{\Gamma }_{core}^{\beta }$$

$$\sum _i\sum _{{\alpha }^{\prime }}{\tilde{S}}_i^{{\alpha }^{\prime }}\left({\varphi }_k^m\right)={\Gamma }_{out}^{\beta +}+{\Gamma }_{puff}^{\beta }+{\Gamma }_{spt}^{\beta }-{\Gamma }_{pump}^{\beta }-{\Gamma }_{leak}^{\beta }$$

$$\ln X={\tau }_X^{-1}t+C\Rightarrow \frac{1}{X}\frac{dX}{dt}=\frac{1}{{\tau }_X}$$

$$N_{\beta }=\int \sum _{{\alpha }^{\prime }}{n}_{{\alpha }^{\prime }}dV,E_e=\frac{3}{2}\int n_e{T}_{e}dV$$

$$E_i=\int \left(\frac{3}{2}\sum _{\alpha }{n}_{\alpha }{T}_i+\frac{1}{2}\sum _{\alpha }{m}_{\alpha }{n}_{\alpha }{v}_{\alpha }^2\right)dV$$

$$\Delta P=\frac{P_{SOL}-P_{PFC}^{+}-P_{PFC}^n-P_{rad}-P_{core}^n}{P_{SOL}}$$

Full text

I Introduction

A combination of a 2D finite-volume plasma transport code with a kinetic Monte-Carlo model for neutral particles is typically applied for numerical modelling of the tokamak edge and divertor plasmas. A well known example of such modelling tool is the code package B2-EIRENE ReiterJNM1992 ; Reiter_EIRENE2005 (SOLPS) widely used in the field. The Monte-Carlo method allows physically accurate description of atomic and molecular kinetics in complex geometries, but has a disadvantage of random error - statistical noise in the calculated quantity. There were always concerns that this statistical noise can have detrimental impact on the coupled solution Maddison1994 .

In the present paper one specific noise related issue which can lead to pathological solutions is addressed - violation of the global particle balance. It is shown that the error in the steady-state particle balance can be presented as a sum of three terms. Those are the operator splitting error, residual of the fluid solver, and the time-derivative. Whereas the first term can be effectively reduced by the source re-scaling, the reduction of residuals may require iterative solution of the discretized fluid equations after each call of the Monte-Carlo model. This can, in turn, pose severe restrictions on the time-step and lead to a very long overall run-time. E.g. in the ITER modelling studies KukushaFED2011 one model run could take several months of wall-clock time. Special diagnostics for monitoring of the particle balance allow to clearly identify the cases when reduction of residuals is absolutely necessary, and the corresponding measures must be taken.

This paper presents in condensed form the most important findings from the dedicated studiy of the SOLPS code KotovJul2014 . Prototypes of the numerical diagnostics were implemented and tested in the code SOLPS4.3 which is the legacy version of B2-EIRENE used in the past for the ITER design modelling KukushaFED2011 . The approach itself is thought to be applicable to any finite-volume edge code. The numerical convergence is analyzed here only in terms of the global balances and criteria of the (quasi-)steady-state. It is not attempted to use the stricter methods of analysis proposed recently for the combination of fluid and Monte-Carlo models in GhoosJCP2016 . Only steady-state solutions are considered.

The rest of the paper is organized as follows. In the next section a finite-volume fluid code with source terms calculated by Monte-Carlo is described in general terms. In Section III the diagnostics for monitoring of the particle balance are introduced. An example of calculations with different error (residual) reduction techniques is discussed in Section IV. Further methods which can be used to reduce the residuals and the associated error in the particle balance are outlined in Section V. Last section summaries the conclusions.

II Coupling of a finite-volume and a Monte-Carlo models

Here only minimal information about numerical procedure of the code B2-EIRENE is given which is required for the subsequent discussion. The plasma transport code B2 BraamsPhDThesis ; BraamsNET solves a set of 2D (axi-symmetric) equations for particle conservation, parallel momentum balance, electron and ion energy. The full set of equations can be found in Ref. BraamsNET, , Chapter 2. The computational domain comprises the scrape-off-layer (SOL) region outside of the 1st magnetic separatrix, and the edge of the core plasma inside the separatrix.

Finite-volume discretization of the differential equations  PatankarBook1980 ; BraamsNET leads to a set of algebraic equations which can be symbolically written as:

[TABLE]

Here $\phi$ is the solution vector: $n_{\alpha}$ is the number density and $u_{\alpha}$ is the parallel velocity of the ion fluid $\alpha$, $T_{e}$ and $T_{i}$ are the electron and ion temperatures. The discrete variables are defined in the cell centers or on the cell faces of the grid. ${\rm F}\left(\phi\right)$ is the non-linear vector function, $S\left(\phi\right)$ are the source terms calculated by the test particle Monte-Carlo method in each grid cell.

To find the solution of Equations (1) a discrete time-derivative $D$ is added to the equations, and iterations over time are performed. On each time-iteration $k$ the solution $\phi_{k}$ of the following set of equation has to be found:

[TABLE]

The “time derivative” is defined such that $D\left(\phi_{k},\phi_{k-1}=\phi_{k}\right)=0$. E.g. for the particle continuity $D\left(n_{k},n_{k-1}\right)=\left(n_{k}-n_{k-1}\right)/\Delta t$, where $\Delta t$ is the time-step. The notation with tilde $\tilde{S}\left(\phi_{k-1}\right)$ underlines that this source term is calculated by Monte-Carlo and contains random error, as opposite to the “exact” value $S\left(\phi\right)$ which would be obtained with the infinite number of test particles.

In the code B2 the set of non-linear algebraic equations (2) is solved by simple iterations and block Gauss-Seidel algorithm (splitting by equations). The so called “internal iterations” of B2 are described in detail in Ref. BraamsNET, , Chapter 3, one may also refer to Ref. KotovJul2014, , Chapter 1.2. Approximate solution $\phi^{m}_{k}$ obtained at the end of internal iteration $m$ can be inserted back into Equation (2) to find the residual:

[TABLE]

That is, the found $\phi^{m}_{k}$ fulfills the equation:

[TABLE]

By comparing with Equation (1) one can see that the difference between $S\left(\phi\right)$ and the right hand side of Equation (4) can be seen as generalization of the common residual $R$.

In the simplest procedure the source terms are calculated at the beginning of internal iterations and are fixed afterward. That is, they stay as $\tilde{S}\left(\phi_{k-1}\right)$. However, certain modifications of the sources can be made in the iterative solver to adjust them with the changed plasma solution $\phi^{m}_{k}$. This modification is reflected in the notation as $\tilde{S}\left(\phi^{m}_{k}|\phi_{k-1}\right)$.

II.0.1 Measures to ensure particle conservation

Critical importance of very high accuracy in the global particle balances for the reactor-scale edge modelling was recognized back at the early stages of the ITER analysis KukushaFED2011 ; KukushaPPCF2002 . To reach this high accuracy the Monte-Carlo neutral transport code must ensure perfect particle conservation in its solution. The internal balance in the neutral solver is usually achieved by re-scaling of the volumetric ion sources estimated by the statistical procedure to make them entirely consistent with the primary sources of neutral particles. To increase accuracy the particles originating from the different primary sources $s$ are sampled independently from each other - the source is split into independent “strata”. The primary sources of neutrals are: i) recombination of ions on the solid surfaces - “recycling”; ii) volumetric recombination in plasma; iii) gas puff; iv) erosion. The strength of recycling sources is proportional to the ion fluxes.

If the volumetric ion sources $\tilde{S}\left(\phi_{k-1}\right)$ stay fixed, but the fluxes of neutralized (recycled) ions change in the course of internal iterations, then an imbalance in the sinks and sources occurs. To compensate for this inconsistency the sources of ions $\alpha$ coming from recycling strata $s$: $\tilde{S}^{s}_{\alpha}\left(\phi_{k-1}\right)$, must be re-scaled as follows:

[TABLE]

Here $j$ is the index of internal iteration, $\phi^{0}_{k}=\phi_{k-1}$, $Q$ is the total flux of neutralized ions to which the source $S^{s}_{\alpha}$ is proportional. E.g. if $\alpha$ is He+ then $Q^{s}_{\beta}$ is the sum of the fluxes of He+ and He*++*.

III Monitoring of the particle balance

Numerical diagnostic for monitoring of the steady-state global particle balance can be derived from Equation (4) by transforming it into the form:

[TABLE]

Error (inconsistency) in the particle balance is defined separately for each ion species $\beta$. “Ion species” here is the chemical element as opposite to “ion fluids” which are charged states of an element. E.g. species Carbon includes 6 ion fluids from C+ to C6+.

Equation (6) is applied to the discretized continuity equation for each ion fluid $\alpha$ in each cell $i$. Then the sum is calculated:

[TABLE]

Here $\sum_{i}$ is the sum over all grid cells, $\sum_{\alpha^{\prime}}$ is the sum over all ion fluids which belong to ion species $\beta$. It is readily seen that zero left hand side of Equation (7) means perfect balance between volumetric sources and fluxes, and the right hand side is the error in the global particle balance of species $\beta$.

Alternative way of writing the particle balance uses formulation via fluxes KukushaFED2011 ; KukushaPPCF2002 :

[TABLE]

Here $\Gamma^{\beta}_{puff}$ is the strength of external particle source - gas puff, $\Gamma^{\beta}_{core}$ is the ion flux through the core grid boundary, $\Gamma^{\beta}_{spt}$ is the flux sputtered (eroded) from the solid surfaces, $\Gamma^{\beta}_{pump}$ is the flux (of both ions and neutrals) absorbed on solid surfaces - pumped flux, $\Gamma^{\beta}_{leak}$ is the flux of atoms which leak to the core.

The final steady-state solution has to self-adjust in such way that the rate with which the particles are removed from the system $\Gamma^{\beta}_{pump}+\Gamma^{\beta}_{leak}$ becomes equal to the particle input:

[TABLE]

That is, $\Gamma^{\beta}_{in}$ serves as a scale with which the particle balance error has to be compared. The numerical solution can be considered as physically meaningful only if this error $<<\Gamma^{\beta}_{in}$.

Coming back to Equation (7), its right hand side yields the following expression for the relative error:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

First term $\Delta^{\beta}_{R}$ contains residuals calculated with Equation (3) after the end of internal iterations. This is the error in the solution of the set of nonlinear finite-volume equations on each time-iteration. The term $\Delta^{\beta}_{S}$ is due to inconsistency of the neutral-related sources calculated on the “old” and “new” plasma. It can be called an operator splitting error. This term can become large if, e.g., the re-scaling procedure, Equation (3), is not implemented. Last term $\Delta^{\beta}_{T}$ is the time derivative which is considered as error when a stationary solution is looked for.

If the plasma fluxes in Equation (8) are taken from the solution $\phi^{m}_{k}$, and the neutral fluxes are calculated on the same plasma, then it is easy to show that Equations (8) and (10) must yield exactly same result when one extra condition is fulfilled. This condition is the discrete analogue of the divergence theorem:

[TABLE]

Here $\Gamma^{\beta+}_{out}$ is the total flux of ions of species $\beta$ to the grid boundaries. The total ion source is calculated as the total source of neutral particles minus their pumped and leaked fluxes:

[TABLE]

Volume recombination does not appear in Equation (12) because atoms originating from recombination which re-ionize back in plasma do not contribute to the net source, and particles which are removed from the system are already included in $\Gamma^{\beta}_{pump}$ and $\Gamma^{\beta}_{leak}$. Subtracting Equation (11) from Equation (12) yields the nominator of Equation (8).

In practice it makes sense to use both diagnostics in parallel. Incorrect particle balance in the solution for neutrals or a mistake in the transfer of ion fluxes to the Monte-Carlo code manifests itself as non-physical particle sinks or sources. The diagnostic of Equation (10) may not be able to detect them because it does not distinguish between “legitimate” and “illegitimate” sources and sinks of neutrals. This distinction is made in Equation (8). The two diagnostics are complimentary to each other and enable an additional consistency check.

IV An example of case study

An example discussed here is based on a SOLPS4.3 run from the data-base of ITER simulations KukushaNF2009 (case #1568vk4, see Ref. KotovJul2014, , Chapter 4.2). The model plasma consists of all charged states of D, He and C. Power entering the computational domain from the core is equal to $P_{SOL}$=80 MW, 47 % of $P_{SOL}$ is radiated, mainly by C ions. The D particle content is controlled by the gas puff $\Gamma^{D}_{puff}$=1.17e22 D-at$\cdot$s*-1* and ion flux from the core $\Gamma^{D}_{core}$=0.91e22 s*-1*. Influx of He ions from the core is set to $\Gamma^{He}_{core}$=2.1e20 s*-1*. All plasma facing components in the model are covered by carbon. The pump is modelled by an absorbing surface in divertor beneath the dome. The solution represents a relatively hot attached plasma in front of divertor targets, with insignificant parallel momentum losses and volume recombination.

In the ITER modelling studies KukushaFED2011 ; KukushaPPCF2002 ; KukushaNF2009 the B2-EIRENE code was always applied with internal iterations in the fluid solver. In the model run in question $m$=20 internal iterations are used, the time-step is $\Delta t$=3e-7 sec. Significant increase of the time-step is not possible: with $\Delta t>$1e-6 sec a numerical instability develops and no stationary solution can be found. It turns out that $\Delta t$ can be increased by orders of magnitude if no internal iterations are applied, that is $m=1$. In this case no visible instability develops even with $\Delta t$=1e-4 sec. However, solutions obtained with and without internal iterations - they are shown in Figure 1 - strongly deviate from each other.

Strictly-speaking, in the presence of Monte-Carlo noise in the source terms the solutions never reach true steady-state. One can only speak about quasi-steady-state solution which randomly oscillates around some average. As applied to the B2-EIRENE runs, the “quasi-steady-state” is defined through characteristic decay times of selected parameters derived from their time-traces, see Appendix. In practice the run is regarded as converged if the condition of quasi-steady-state is fulfilled, and if errors in the global power and particle balances are small.

Errors in the balances for the two model runs considered here are given in Table 1. The power balance error $\Delta P$ is defined by Equation (13), $\Delta\Gamma$ are calculated using Equation (8). One can see that $\Delta P$ is small in both cases. The situation is completely different for the particle balance. Whereas in the simulation made with $m=20$ both \Delta\Gamma^{D,He}$$<10 %, in the “fast” run the error approaches 100 %. That is, in the solution obtained with $m=1$ the pumped fluxes are negligible compared to $\Gamma^{D,He}_{in}$.

Individual terms of the error are shown in Table 2 for both recycling species D and He. There is a good agreement between $\Delta\Gamma$ and $\Delta_{R}+\Delta_{S}+\Delta_{T}$ calculated independently by two diagnostics. From Table 2 the reason of the large error in $m=1$ case immediately becomes clear. While $\Delta^{D}_{S}$ and $\Delta^{D}_{T}$ always remain relatively small, $\Delta_{R}$ becomes very large if the code is operated without internal iterations after each Monte-Carlo call.

Particle balance is much more difficult to converge than the power balance because of different relation between the controlling flux and internal sources and sinks in the system. For the power the sources and sinks in plasma are smaller than $P_{SOL}$. In contrast, the particles are “recycled” between the plasma and solid surfaces and the total volumetric ion sources by far exceed $\Gamma_{in}$. In the present example $\sum_{i}\tilde{S}^{D}$=4.3e24 s*-1* and $\sum_{i}\tilde{S}^{He}$=8.1e22 s*-1*. Those numbers are by more than two orders of magnitude larger than $\Gamma_{in}$ of those species. This problem does not appear for C because in the present model all incident C particles are absorbed on the surfaces - this species does not recycle.

In the B2-EIRENE model run above extra measures for reduction of residuals on each time-iteration were absolutely necessary. Only in this case a solution can be obtained which is correct in terms of the global balances. The techniques such as internal iterations in B2 can impose severe limitation on the time-step, and it is not attractive from the run-time point of view to operate the code in this mode. Experience has shown that the use of B2-EIRENE with $m$=1 does not always lead to deviations as dramatic as that shown in Figure 1. E.g. in ITER cases with single fluid (D only) plasma $\Delta\Gamma^{D}$ was found to be sufficiently small both with and without internal iterations, and the obtained solutions are close to each other, see example in Ref. KotovJul2014, , Chapter 4.1.

The multi-fluid simulation analyzed here clearly demonstrates that this must not always be the case. This example emphasizes that in each simulation the particle balance has to be carefully monitored with the special diagnostics. Too large error detected by the diagnostic is an unequivocal indication that the residual reduction techniques must be applied irrespective of the run-time penalty which they impose.

V Reduction of residuals

A series of studies was undertaken with the B2-EIRENE code to find algorithms which would deliver sufficiently good accuracy without penalizing the run-time. Their outcome may be of general interest for developers and users of other edge modelling codes as well. Main results are briefly summarized in this section.

As a simplest remedy to the particle balance problem a “0D correction” was first tried, see Ref. KotovJul2014, , Chapter 5.5. The ion density in the whole computational domain is multiplied by a constant factor calculated in such way that with the corrected ion fluxes $\Delta_{R}$ automatically becomes zero. It was found that this method cannot be used because it always produces solutions oscillating in time, and no stationary solutions.

Much more success was achieved with a correction based on iterative relaxation of the finite-volume continuity equations. Technical details of the implementation in B2 can be found in Ref. KotovJul2014, , Chapter 5.2-5.4. This algorithm works as follows. The whole set of equations for particle, momentum and energy balances is relaxed only on the first internal iteration. On subsequent iterations only equations for particle continuity are relaxed. To be precise, in the code B2 those are pressure correction equations where both the density and velocity fields are modified. (B2 uses compressible version of the Patankar’s SIMPLE algorithm, see Ref. PatankarBook1980, , Chapter 6.7 and Ref. BraamsNET, , Chapter 3.) Nevertheless, correction of the particle balance via relaxation of the pressure correction equations was found to be very reliable. Tests performed for the same ITER model as in Section IV showed that such iterations robustly converge with time-steps up to $\Delta t$=1e-4 sec.

Results obtained with this algorithm can be found in the last row in Tables 1 and 2. The run was performed with 99 iterations for continuity equations after one full internal iteration, which is reflected in the designation $m$=1+99. Despite increased $\Delta_{S}$ the method leads to significant reduction of the total error $\Delta$ due to reduction of $\Delta_{R}$. As expected, the main disadvantage of this procedure is that it increases residuals of other equations. Closer investigations (Ref. KotovJul2014, , Chapter 6.1) showed that especially the parallel momentum balance suffers. However, comparison of the solutions obtained with full internal iterations and with the reduced scheme demonstrates that they are close to each other: the $m$=1+99 case is shown by dashed lines in Figure 1. Moreover, the tests demonstrated that this result holds even for the ITER model with high density detached divertor, see Ref. KotovJul2014, , Chapter 6.4. Hence, the method can be suggested for use in a two stage approach for fast finding of the initial approximation to the solution which is then refined on the second “slow” stage by more accurate techniques.

As a next step a scheme was proposed where coupled continuity and parallel momentum balance equations are iterated - without equations for temperatures. This kind of “incomplete internal iterations” was implemented in B2-EIRENE and tested as well, but the results were found to be unsatisfactory: Ref. KotovApr2015, , Chapter 2.3. Tests showed that similar to the full internal iterations the “incomplete iterations” are prone to numerical instabilities with large time-steps, and therefore bring no advantages. The SIMPLE pressure correction which introduces extra non-linearity is a possible reason of this behavior. The scheme could be improved if monolithic coupling of the continuity and momentum equation would be applied instead. That is, when corrections for both the density and the velocity fields are calculated simultaneously in a one set of linear equations.

A fairly simple technique which increases accuracy and can be easily implemented in any code is time-averaging of source terms, Ref. KotovJul2014, , Chapter 3. Although this algorithm can be helpful in many cases, it was found to be not always efficient enought in reducing $\Delta\Gamma$, in particular with impurities, see example in Ref. KotovJul2014, , Chapter 4.2. In KawashimaPFR2006 a more advanced “piling method” is described which do not reset the whole history as the calculation of the new average starts.

Finally, the brute force method can always be applied to decrease both the statistical error in the source terms and the residuals - massive increase of the number of test particles. Applicability of this solution strongly depends on the available computing hardware. The test particle Monte-Carlo algorithm is easy to parallelize, and the increased number of particles does not necessarily mean the increased wall-clock run time. Experience KotovJul2014 has indicated that the pure “brute force compensation” of the particle balance issue described in Section IV is likely to require $>>$100 processors to be practical.

VI Conclusions

Use of the test particle Monte-Carlo for neutrals in the tokamak edge modelling codes has an unpleasant side effect of random error in the source terms. If no special measures are taken, then this persistent statistical noise leads to residuals of the discretized fluid equations which do not converge, but saturate at a certain level. In the present paper one particular well identified issue caused by the saturated residuals has been described. It has been shown that too large finite-volume residuals can cause crude violation of the global particle balance. In turn, for the system in question - the tokamak edge and divertor plasma - violation of the particle conservation may have a very strong (“zero order”) non-local impact on the whole numerical solution.

There are computational techniques which can effectively reduce the residuals. E.g. in the code B2 which uses splitting by equations an extra loop of simple iterations on each time-iteration is applied. However, severe restriction imposed by those internal iterations on the time-step leads to a very long overall model run-time when this option is used. With numerical diagnostics proposed in this work it can be unambiguously identified when the too large error in the particle balance is caused by the saturated residuals, and the residual reduction techniques must be applied to obtain the physically meaningful solution. The diagnostics can be implemented in any finite-volume edge code.

The problem describe here would become less of an issue if solving the set of non-linear equations on each time-iteration would not require reduced time-step. If such solvers are not feasible, then the accuracy and run-time drawbacks may even outweight the very advantage of using the kinetic test particle Monte-Carlo in the self-consistent models. The drawbacks can be partly compensated by reducing the statistical error which is, in principle, only a matter of available computing resources. Emerging heterogeneous CPU-booster architectures DEEP could be particularly well suited for the combination of a fluid and a Monte-Carlo code. While the serial finite-volume part runs on CPU, the Monte-Carlo part can make use of massive parallelization on hundreds of processing units on the accelerator.

Acknowledgements.

This work was performed under EFDA Work Programme 2013 “Assessment Studies for SOLPS Optimisation” (WP13-SOL).

Appendix A Practical convergence criteria applied to the tokamak edge modelling code B2-EIRENE

Characteristic time-scale $\tau_{X}$ of the parameter $X$ is calculated from its time-trace $X(t_{k})$ by fitting it with a linear function:

[TABLE]

In the present paper the number of last time-iterations used for the fit was equal to $\max{\left(2000,N_{p}^{(5\;\mu s)}\right)}$, where $N_{p}^{(5\;\mu s)}$ is the number of points which cover last 5 $\mu s$ of physical time. Least-square method is applied to find the parameters $\tau_{X}$ and $C$. Same data-points were used to calculate average $\Delta\Gamma$ and $\Delta P$ in Table 1 and $\Delta_{R,S,T}$ in Table 2.

The control parameters for which $\tau_{X}$ are calculated are the total amount of ions $N_{\beta}$ of species $\beta$, total diamagnetic energy in electrons $E_{e}$ and ions $E_{i}$:

[TABLE]

[TABLE]

as well as plasma parameters averaged along the magnetic separatrix: $<n_{e}>^{sep}$, $<T_{e}>^{sep}$, $<T_{i}>^{sep}$. Here the integration is performed over the whole computational grid, $V$ is geometrical volume, $\sum_{\alpha}$ is the sum over all ion fluids, $n_{e}$ is the electron density, $m_{\alpha}$ is the atomic mass of ions, $v_{\alpha}$ is their average macroscopic velocity. The B2-EIRENE solutions analyzed in this paper were regarded as stationary when $\tau_{X}>$3 sec for all the parameters listed above. For $N_{D}$ and $N_{He}$ $\tau_{X}>$15 sec.

Besides this condition of steady-state the errors in the global particle and power balances are checked. The error in particle balance is expressed by Equation (8). Relative error in the power balance is defined as follows:

[TABLE]

Here $P_{SOL}$ is the power influx into the computational domain from the core plasma, $P^{+}_{PFC}$ is the power deposited by charged particles to the Plasma Facing Components (PFC), $P^{n}_{PFC}$ is the power deposited to PFC by neutrals, $P_{rad}$ is the power radiated by both charged and neutral particles, $P^{n}_{core}$ is the power transferred by neutrals back to the core.