E Agullo

L Giraud

V Joncquières

G Marait

L Poirel

O Vermorel

W Villafana

A note on the strong parallel scalability of numerically scalable Poisson linear solvers

Keywords: Strong parallel scalability, numerical scalability, Poisson solvers, PIC model, fluid model, plasma simulation Scalabilité parallèle forte, scalabilité numérique, solveur de Poisson, modème PIC, modèle fluide, simulation de plasma

In the context of a parallel plasma physics simulation code, we perform a qualitative performance study between two natural candidates for the parallel solution of 3D Poisson problems that are multigrid and domain decomposition. We selected one representative of each of these numerical techniques implemented in state of the art parallel packages and show that depending on the regime used in terms of number of unknowns per computing cores the best alternative in terms of time to solution varies. Those results show the interest of having both types of numerical solvers integrated in a simulation code that can be used in very different configurations in terms of selected modelisations, problem sizes and parallel computing platforms.

Une note sur le passage à l'échelle parallèle des solveurs linéaires de Poisson numériquement scalables

Résumé : Dans le contexte d'un code de simulation parallèle de la physique des plasmas, nous effectuons une étude qualitative des performances entre deux candidats naturels pour la résolution parallèle de problèmes de Poisson en 3D que sontles méthodes multigrilles et les méthodes de décomposition de domaine. Nous avons sélectionné un représentant de chacune de ces techniques numériques implémentées dans des progiciels parallèles de pointe. Nous montrons que selon le régime utilisé en termes de nombre d'inconnues par coeurs de calcul la meilleure alternative en termes de temps de résolution varie. Ces résultats montrent l'intérêt de disposer de ces deux types de solveurs numériques intégrés pour de la simulation intensive qui peut être utilisée dans des configurations très différentes en termes de modélisations choisies, taille des problèmes et plates-formes de calcul haute performance.

Introduction

Numerical simulation has become a major tool in the understanding of complex phenomena as it allows for detailed studies at fine time and spatial scales once a suited model has been designed. In many engineering applications in industry or academia, such simulations involve at some stage the solution of sparse linear systems which size increases when the discretization step decreases or the bounded domain of calculation increases. When the problem is well-posed, the matrix associated with the discretization is non-singular and the solution of the linear system is not a theoretical issue. However, it often turns out to be a challenge in practice when the model has to be implemented in a parallel numerical software where the time to solution becomes a key criterion of performance. In the literature, the parallel scalability of the linear algebra solvers is often studied in details by the designers of the numerical schemes in a stand-alone setting. However, in practice, the solvers are only one component of a more complex software environment for large simulation codes; consequently they are used in various regimes that depend on the context of the simulations.

In this work, we investigate the parallel performance of parallel iterative schemes for the solution of a 3D Poisson problem in a plasma physics simulation. In such a simulation, most of the parallel computation kernels do scale very well and the linear solver, that is not the most computing intensive part, possibly becomes the bottleneck of the parallel performance. In that context, we perform parallel experiments with two natural candidates for solving 3D Poisson problems efficiently in parallel, namely a multigrid solver and a 2-level domain decomposition solver. On a realistic test example, we show that depending on the regime the solvers are used, in terms of number of degrees of freedom per computing cores, the time to solution ranking of the two solvers varies while the domain decomposition solver is always the most memory consuming.

The paper is organized as follows. In Section 2 the plasma simulation code used in this study is presented. The set of equations used is briefly described in order to highlight the challenges raised by their numerical resolution. Section 3 reviews the different numerical techniques available for the resolution of linear systems arising from the discretization of 3D elliptic Partial Differential Equations (PDE). The multigrid solver and the 2-level domain decomposition solver used for this parallel performance study are introduced in Section 4. The main results obtained with the two solvers on a typical plasma thruster configuration are discussed. Finally, some concluding remarks are provided in Section 5.

The need of a scalable Poisson solver for space propulsion numerical simulation

Electric propulsion is an increasingly popular technology to power a satellite for station-keeping or orbit raising missions. For electric propulsion companies, the design and qualification phases of electric thrusters are long (several years) and expensive and some physical phenomena occurring there are still not understood (loss of efficiency, appearance of spurious instabilities, anomalous erosion of the discharge chamber wall, ...). The AVIP code has been developed in order to help the understanding of plasma physics phenomena inside space plasma thrusters, in particular Hall effect thrusters [START_REF] Joncquières | Modélisation et simulation numérique des moteurs à effet Hall[END_REF]. This solver uses unstructured meshes in a high performance computing framework based on the AVBP code [START_REF] Gourdain | High performance parallel computing of flows in complex geometries: I[END_REF], state-of-the-art CFD code for reactive flows. AVIP solves the equations of charged particles (ions and electrons) and neutrals inside a plasma thruster chamber using two formalisms:

• an Eulerian fluid formalism describing the evolution of particle macroscopic quantities with a fluid equation system [START_REF] Joncquières | Modélisation et simulation numérique des moteurs à effet Hall[END_REF][START_REF] Joncquières | A 10-moment fluid numerical solver of plasma with sheaths in a Hall effect thruster[END_REF], which can be written as :

Inria ∂ W α ∂t = F convection ( W α ) + F Lorentz,α ( W α ) + S collisions ( W α ) (1)
with W α the vector of conservative variables for each species α (n for neutral, i for ion and e for electron) inside the plasma :

W α =   ρ α = m α n α ρ u α ρE α   (2)
m α , ρ α , n α , u α and E α represent respectively the mass, mass density, number density, velocity vector and total energy of species α.

In Eq. 1, convections forces are denoted by F convection and collision processes between different populations of particles are represented through the term S collisions . The Lorentz force F Lorentz,α represents the dynamics of charged particles (ions or electrons) in the electromagnetic field and is expressed by:

F Lorentz,α = q( E + u α × B) (3)
with q the charge of the considered species, E and B the electric and magnetic fields. In practice q is 0 for neutrals, -e the elementary charge for electrons and a positive multiple of the elementary charge e for ions (singly ionized or more).

• a Lagrangian PIC (Particle-In-Cell) formalism in which the velocity space is discretized by considering a large number of particles belonging to species α. Each individual particle is tracked down by applying and integrating the Second Newton's Law:

d u p,α dt = F Lorentz,p,α (4)
d x p,α dt = u p,α (5)
with x p,α and u p,α the position and velocity of each individual particle. Here the Lorentz force F Lorentz,p,α has the same expression has in Eq. 3 but the macroscopic speed u α is replaced by the individual particle speed u p,α .

The system of fluid equations Eq. 1 is recognized as more complex to formulate than the system of PIC equations Eqs. 4 and 5 because it usually requires the introduction of many models to close the problem. For example, assumptions about the particle velocity distribution must be used to derive the Euler-type fluid system, which can potentially bias the simulation results. In this sense, the PIC approach can be considered more precise than the fluid approach and is generally the preferred approach for reproducing physics with a higher level of detail. However the number of particles required in a PIC approach to describe the overall plasma physics can become very large, which makes a PIC formulation more computationally challenging than a fluid formulation. Indeed the average number density of charged particles inside a thruster chamber is around n e,i = 10 18 m -3 particles. Considering a thruster chamber of volume V c ≈ 60 cm -3 , the number of particles to track is therefore around 10 13 -10 14 which greatly exceeds the current calculation capabilities. Physical particles are then grouped in numerical macro-particles to follow their trajectories all at once. Typical 3D PIC simulations use ∼ 10 9 numerical macroparticles. Despite this numerical trick, the computational time to treat all these particles is very large and high performances on large clusters are required to run realistic simulations of Hall effect thrusters. In practice, the fluid and PIC formalisms complement each other since PIC simulations can be used to assess the fluid model which aims to compute greater time ranges.

The Lorentz force F Lorentz appearing in both formalisms (Eqs. 4 and 1) is a function of the electric field E and magnetic field B. These electromagnetic fields E and B are related to each other via the well-known Maxwell equations. However in several space thrusters, the magnetic field B can be considered constant and the set of Maxwell equations reduces to a Poisson equation for the electric field E. The Poisson equation describes the behavior of the electric potential φ generated by the difference of electric charges inside the domain (i.e., the difference between ion and electron densities) :

∆φ = - e 0 (n i -n e ) (6)
with 0 the dielectric permittivity.

Once the Poisson equation has been solved, the electric field can then be calculated with the relation:

E = -∇φ.
The electric field E and the density of charged particles n e and n i are strongly coupled in this set of equations. Any error on the electric field can lead to a misprediction of the electron velocity that will directly affect the difference of charges. A high precision on the Lorentz force prediction and the electric field is therefore required to correctly reproduce the behavior of the plasma. Concerning the Poisson equation, great precision on the gradient of the unknown is required.

In general the meshes used in a 3D configuration of a thruster are relatively homogeneous. Boundary layers, commonly known as electric sheaths, exist nears walls and are characterized by high electric fields and large gradients of particle densities. Their discretization requires a local refinement of the mesh [START_REF] Joncquières | A 10-moment fluid numerical solver of plasma with sheaths in a Hall effect thruster[END_REF] but this difference in mesh size induces only a moderate heterogeneity in the entries of the discretization matrix.

In the present case, both AVIP-fluid and AVIP-PIC use a node-centered formalism based on a finite volume method. Therefore, the Poisson equation in Eq. 6 is discretized at the nodes as shown in Fig. 1. The details of the computation over the nodal volume is left to the reader [START_REF] Joncquières | Modélisation et simulation numérique des moteurs à effet Hall[END_REF].

Eventually the discretization of Eq. 6 leads to a linear system with a symmetric positive matrix A which solution φ h is the potential φ at the discretization nodes:

φ h = - e 0 (n i -n e ) = f h . (7)
A short time to solution is always a crucial feature for a numerical solver, but it is even more true for AVIP which aims to deal with real-size realistic configurations. In practice, the computational time of an AVIP simulation (regardless of the formalism used, fluid or PIC) can be divided into three main physical operators: particle convection, collisions between particles (Monte-Carlo algorithm for AVIP-PIC, implicit solver for AVIP-Fluid) and Poisson equation.

The accuracy and the cost of a fluid approach is directly governed by the size of the mesh. The Lagrangian approach is based on a different data structure and its cost responds to a different logic. Accuracy is mainly driven by the number of numerical particles tracked and not by the eulerian mesh resolution. In other words, it means that increased precision is usually achieved by increasing the number of particles while keeping the mesh unchanged, in order to improve the statistical representation of the dispersed phase. For all Lagrangian operators, increasing the number of cores is an obvious solution to reduce the time to solution considering their very good scalability. The challenge then is to solve the Poisson equation on an even smaller number of degrees of freedom (dof) per core. In practice, the objective is to solve the Poisson equation for numbers of dof per core of the order of a few hundred, numbers for which standard linear solvers are not necessarily designed and efficient. Thus, a Poisson solver able to work under different regimes in terms of number of unknowns per computing cores is absolutely essential for anyone who wishes to use both fluid and PIC approaches in a same code.

Numerical linear solvers

The solution of sparse linear systems of equations has been studied for a long time and still remains an active field of research. Many numerical packages have been developed during the last decades to efficiently cope with the features of the ever evolving computer architectures. Among the oldest and general purpose approaches are the ones referred to as direct methods. They are based on Gaussian elimination and rely on either multifrontal or supernodal schemes to address sparsity and parallelism. State of the art parallel packages that implement these techniques are: MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF], PARDISO [START_REF] Schenk | PARDISO: A High-Performance Serial and Parallel Sparse Linear Solver in Semiconductor Device Simulation[END_REF], PaStIX [START_REF] Hénon | PASTIX: A high-performance parallel direct solver for sparse symmetric positive definite systems[END_REF] or SuperLU [START_REF] Li | SuperLU DIST: A Scalable Distributed-memory Sparse Direct Solver for Unsymmetric Linear Systems[END_REF] to name a few. The direct methods are numerically robust, but their robustness comes at the price of large memory consumption and floating point operation (flop) count for the solution of linear systems arising from 3D PDE. For instance, the solution of a PDE defined on a 3D Cartesian grid of size n = N 3 requires a memory consumption of O(n 2/3 ) and a flop count of O(n 2 ). In order to reduce these costs, recent developments attempt to exploit data sparsity as introduced in the H-matrix community [START_REF] Grasedyck | Construction and Arithmetics of H-Matrices[END_REF]. Alternatively, iterative methods have a low memory footprint but need to be combined with acceleration techniques to control the number of iterations. Modern iterative approaches are based on Krylov subspace and the acceleration techniques are referred to as preconditioning. For the solution of linear systems arising from the discretization of elliptic PDE, numerically scalable preconditionners have been designed, that can be split into two main families.

The first ones are multigrid techniques initially introduced on regular grids and later extended to unstructured matrices with the naming Algebraic Multi-Grid (AMG). In a two-level setting, the general idea is to smooth the high frequencies of the error via simple relaxation techniques and solve for the low frequency modes of the error on a coarser grid where those modes can be captured in a lower dimensional space. Applying recursively the two-grid idea to solve for the low frequencies leads to the multigrid method. In that context a sequence of linear systems of decreasing size are considered, on each of them a few relaxation steps, called pre-smoothing, are performed before restricting the residual on the next coarser space, where the error equation is handled. Once the coarsest level is reached, the error equation is solved and the solution is prolongated to the preceding finer level to correct the current iterate on this level; few more relaxation steps, called post-smoothing, might be applied before prolonging again the resulting approximated error to the next finer level. This procedure is applied until the finest level is reached and a new iterate for the solution of the linear system is produced. For elliptic problems, the multigrid algorithms have unique computational and numerical features such as a convergence that is independent of the mesh size and an asymptotic computational cost (memory, flop) that is linear with respect to the problem size. Their parallel implementation is based on a matrix partitioning, all the steps (smoothing, restriction, prolongation) can be implemented with a few neighbor-to-neighbor communications as long as the coarse problems are fine enough to keep all the computing cores busy. Intrinsically, when going down the multigrid hierarchy the computational granularity reduces. This is a very brief and incomplete overview of multigrid principles, we refer the reader to the rich bibliography and textbooks such as [START_REF] Briggs | A Multigrid Tutorial (2Nd Ed[END_REF][START_REF] Hackbusch | Multi-Grid Methods and Applications[END_REF][START_REF] Trottenberg | Multigrid[END_REF].

Another class of techniques for the solution of PDEs is known as domain decomposition methods. The main motivation is to naturally express some parallelism while overwhelming the unaffordable calculation costs of sparse direct methods if they were applied to the whole problem. In these approaches, the robustness of the direct methods is still exploited locally for the solution of the sub-problems resulting from the decomposition and an iterative process is designed to "glue" the local solutions so that each of them eventually converges towards the global solution restricted to the individual partition. The decomposition of the domain can have overlapping sub-domains, leading to a class of techniques originally called Schwarz methods. This refers to the alternating method introduced by H. A. Schwarz as a constructive procedure to prove the existence of the solution to a Poisson problem in a complex domain [START_REF] Schwarz | Über Einen Grenzübergang Durch Alternirendes Verfahren[END_REF]; he built the solution as the limit of a series of solutions defined on a set of overlapping domains where each has a canonical geometry and whose union forms the complex domain. This pioneering idea has inspired and still inspires research activities and his name is associated with a large class of techniques. Alternatively, sub-domains without overlap can be considered leading to techniques often referred to as sub-structuring approaches where local solutions are composed to form a reduced linear system, only defined on the interface between the sub-domains. Within this approach, the reduced system, also called the Schur complement system, is solved iteratively. The natural parallelism introduced by the splitting of the domain has a defect due to the associated induced decoupling, the number of iterations increases with the number of sub-domains. A numerical remedy to design numerically scalable preconditioners, i.e., with a convergence that does not depend on the number of sub-domains, consists in introducing a coarse space to recouple contributions from all the sub-domains represented in a low dimensional space. In the parallel design of the solver a special attention should be paid in the implementation of the coarse space correction that might become a bottleneck for the parallel scalability. This illustrates the Inria fragile trade-off to be found between numerical and parallel scalability of the linear solvers. For a more detailed presentation of the domain decomposition techniques and associated theory we refer to the following list of books and reference therein [START_REF] Victorita Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation[END_REF][START_REF] Tarek | Domain Decomposition Methods for the Numerical Solution of Partial Differential Equations[END_REF][START_REF] Quarteroni | Domain Decomposition Methods for Partial Differential Equations[END_REF].

Numerical experiments 4.1 Experimental framework

The parallel performance study has been conducted with one representative of each of the two scalable numerical linear algebra solvers. Namely, for the representative of the multigrid we chose the GAMG algebraic multigrid solver available in the PETSc library [START_REF] Balay | PETSc Web page[END_REF] named AMG in the paper; it implements a smoothed aggregation algebraic multigrid method. As representative of domain decomposition solver, we selected MaPHyS [START_REF] Agullo | Robust preconditioners via generalized eigenproblems for hybrid sparse linear solvers[END_REF][START_REF] Poirel | Algebraic domain decomposition methods for hybrid (iterative/direct) solvers[END_REF] that implements an additive Schwarz preconditioner for the Schur complement with a coarse space correction scheme that follows GENEO [START_REF] Spillane | Achieving robustness through coarse space enrichment in the two level Schwarz framework[END_REF][START_REF] Spillane | Méthodes de décomposition de domaine robustes pour les problèmes symétriques définis positifs[END_REF] ideas. This specific version of Maphys [START_REF] Poirel | Algebraic domain decomposition methods for hybrid (iterative/direct) solvers[END_REF] will be referred to as DDM in the sequel. For both solvers, we tuned the control parameters following the advice of the package developers through a trial and correct iterative loop. The selected parameters that depart from the default are listed in Appendix A, while we refer the reader to the user's guide of the two packages for more details.

The same stopping criterion is used in both linear solvers, which relies on a standard scaled residual norm:

||f h -Aφ h || L2 ≤ f h L2 . (8)
Solver performances are evaluated on a standard AVIP configuration of a 3D Hall effect thruster cavity. The geometry represents a 5 degree azimuthal sector of the ionization chamber and the plume, with dimensions of a standard SPT-100 Hall thruster (15mm high and 34 mm long) [START_REF] Morozov | Fundamentals of stationary plasma thruster theory[END_REF]. It corresponds to an azimuthal extrusion of the 2D case simulated in [START_REF] Joncquières | A fluid formalism for low-temperature plasma flows inside hall effect thrusters in an HPC unstructured solver[END_REF]. The geometry and the boundary conditions applied for ions, electrons and the electric potential are shown in Fig. 2. Only periodic and Dirichlet boundary conditions are used for the Poisson equation. Both fluid and PIC simulations are performed on the same mesh made of 27 million tetrahedral elements (around 5 million dof). The cell size is relatively homogeneous with a minimum of 50µm near the ceramic walls (blue hatches in Fig. 2) and around 80µm elsewhere.

Figure 3 shows some typical plasma quantities inside the computational domain for the fluid simulation at steady state. In this configuration, the propellant (neutral xenon atoms) is injected at the anode (see Fig. 2) while electrons are injected at the cathode. As the neutral xenon atoms diffuse into the channel of the thruster (Fig. 3a), they are ionized by collisions with circulating electrons (Fig. 3b). The difference of electric charges creates an electric field in the axial (z) direction that accelerates the ions in the thruster exit plane in order to provide the thrust (Fig. 3d).

RR n°9423

Parallel performance and discussion

All simulations have been performed on the Occigen supercomputer hosted by CINES in Montpellier (France) [START_REF] Occigen | [END_REF]. In order to get reliable and avoid possible outliers, each run was repeated three times and subsequent performances were averaged out. All the parallel elapsed times reported in this section have been measured on the first 50 time steps of the AVIP solver using either the fluid or the PIC model. A typical simulation requires around two million time steps, but 50 are enough to get representative values of the parallel performance of all the numerical kernels. For all the simulations the number of unknowns involved in the Poisson problem is about 5 millions.

The fluid solver

The computational complexities of the computational kernels involved in the AVIP-Fluid solver are fully governed by the mesh size. In Fig. 4, we display the time spent in each of these numerical kernels. Regarding the performance of the two parallel linear solvers, it can be observed that when the number of unknowns per core is larger than 5 000 (i.e., using less than 1000 cores on our examples), the multigrid solver is faster than the domain decomposition one; after this crossover point the tendency is reversed and multigrid can be up-to 2.5 slower than the domain decomposition solver while it was more than 4 times faster on 360 cores. The time spent in the linear system solution with the multigrid solver reaches a minimum value for a computing core number around 800, while the minimum for the domain decomposition approach is slightly lower and observed for a number of cores around 1500. Eventually, the time for the two solvers increases with the same rate, when the solution time starts to be dominated by communications.

Regarding the other numerical kernels, the cost of the collision operator is always marginal, while the convection operator consumes most of the time when the number of computing cores is small. These two operators scale in excellent proportions to about 2000 dof per core (i.e., 2500 computing cores) but still gain speed beyond this limit.

In terms of total time to solution, one can observe that it reduces with the two solvers up-to 2500 computing cores. Beyond, the time spent in the linear solvers starts to dominate so that the solution time even increases with AMG beyond 2500 computing cores and 3500 cores for the DDM solver. Consequently, it does not make sens to use more than 2500 computing cores for AVIP-Fluid with AMG and more than 3500 with DDM for this specific case. Yet, because AVIP spends more than 80% of the CPU time in the convection operator in the high dof per core range, the initial advantage of the AMG solver is unfortunately not significative.

Inria

The PIC solver

Three PIC simulations were performed by varying the number of macro-particles per cell at initialization: 30, 60 and 120. The total number of macro-particles tracked during the computation therefore varies between approximately 2.43 10 9 and 9.72 10 9 . The computational complexity for the PIC simulations depends on the number of macro-particules involved in the simulation. In order to illustrate this effect on the parallel performance we consider two configurations with 30 and 120 macro-particules per cell at the initialisation. We report in Figs. 5 and6, the parallel strong scalability performance when the number of cores varies from 360 up-to 7200. For these simulations the mesh remains unchanged and the Poisson problem still has around 5 Mdof. In that respect, the comments made above on the parallel performance of the linear solvers remain the same.

In Fig. 5, it can be seen that the kernels related to the convection, collison and interpolation calculations scale relatively well up-to 7200 cores. The total simulation time scales reasonable well up-to 2500 cores, beyond the times to solution do still reduce with the DDM solver and stagnates with AMG. This relatively poor parallel performance is due to the lack of scalability of the two parallel linear solvers when the number of unknowns per computing core decreases.

Similar observations can be made for results presented in Fig. 6 on the parallel performance where 120 macro-particules per cell are used to get a better statistical analysis on the simulated physics. On that example, due to memory constraint, the first measurement point corresponds to a run on 1000 cores, that is the cross-point between the AMG and DDM parallel performance. Because the relative weigth of the linear solver is lower than for the fluid model, the total times to solution decrease up-to 7200 cores when either of the two solvers is used.

Inria

We report in Fig. 7 the memory consumption (mean and max value) for both linear solvers when the number of computing cores varies. The AMG solver is significantly less memory demanding with a very good scalability and a nice load balancing; futhermore the memory footprint becomes negligible when using more and more processors. The DDM approach is much more memory consuming, the mean memory scales also rather well but a significant unbalance appears. This is due to the fact that the coarse space correction that involved a linear system which size grows linearly with the number of sub-domains is solved on a fixed size MPI sub-communicator. Consequently the memory used on thus sub-communicator increases as the number of sub-domains/cores increases. Therefore, the DDM strategy can represent beyond 80 % of total allocated memory for the sub-communicator. While the parallel implementation could be improved it remains an intrinsic parallel bottleneck for a two-level numerical scheme that should be extended to multi-level techniques to possibly exhibit memory performance similar to AMG. For more detailed results on the memory consumption scalability of the two linear solvers we refer the reader to Figs. [START_REF] Schenk | PARDISO: A High-Performance Serial and Parallel Sparse Linear Solver in Semiconductor Device Simulation[END_REF]

Concluding remarks

As shown by the results of this paper, neither of the two numerically scalable linear solvers that are multigrid or domain decomposition are systematically superior to the other on parallel platforms even for the widely encountered Poisson problem. The parallel ranking between the two classes of numerically scalable linear solvers depends on the range of problem and platform sizes. In our example in plasma physics, multigrid is the most efficient with large number of unknowns per core for a moderate number of cores; it quickly suffers from a lost of computational granularity that leads to an increasing time to solution. Scalable domain decomposition methods are more expensive for large sub-domain size due to the weight of its direct solver component but significant gains can quickly be obtained thanks to the nonlinear computational cost of this direct solver ingredient.

While it might sound obvious, this study highlights that the performance agility of a complex simulation code relies on its capability to use a variety of parallel linear solvers and select then on their best performance ranges. The rule for AVIP is to use multigrid when the number of unknowns is larger than 5000 and switch to domain decomposition for lower values. One option for a possibly better scalability on a large number of computing cores would be to depart from the classical MPI SPMD programming model and would consist in buiding a sub-communicator for the solution of the Poisson problem which size will be adjusted to get the best performance of the chosen linear solver. Such an implementation in the AVIP code will be the subject of a possible future work.

A Useful parameters for reproducibility PETSc 3.11 GAMG Options Definition -gamg est ksp type cg iterative solver for the computation -gamg est ksp max it 10 max number of iterations -pc gamg process eq limit 500 max unknowns on each process -pc gamg coarse eq limit 1000 max unknowns on coarsest grid -mg levels esteig ksp type cg iterative solver for the coarse grid -mg levels esteig ksp max it 10 max unknowns for the coarse grid -mg levels pc type sor sor smoothing for the coarse grid -mg levels ksp type richardson solver used for the smoothing -ksp norm type unpreconditioned stopping criterion based on the unpreconditioned residual norm MaPHyS 0.9.

C Fluid

  of a scalable Poisson solver for space propulsion numerical simulation 3 Numerical linear solvers 4 Numerical experiments 4.1 Experimental framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Parallel performance and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 The fluid solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 The PIC solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Concluding remarks A Useful parameters for reproducibility B Benefits of the coarse grid correction C Fluid D Various performance results D.1 Elapsed times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.2 Memory consumption . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Inria

Figure 1 :

 1 Figure 1: Discretization of the Poisson equation in a 2D triangular mesh with a node-centered strategy.

Figure 2 :

 2 Figure 2: Geometry and boundary conditions for the simulation of a SPT-100 (the blue hatches represent ceramic walls) [12]. The figure does not represent the azimuthal direction θ (5 degree sector) for which periodic conditions are imposed.

Figure 3 :

 3 Figure 3: Working principle of a Hall thruster at steady state from a r -z plane view.

Figure 4 :

 4 Figure 4: Elapsed time for experiments using fluid model.

Figure 5 :

 5 Figure 5: Elapsed time for experiments with 30 macro-particles per cell.

Figure 6 :

 6 Figure 6: Elapsed time for experiments with 120 macro-particles per cell.

  -23 in Appendix D.2.

Figure 7 :

 7 Figure 7: Memory consumption in function of the number of cores for two linear solvers. The maximal memory consumption is represented with dashed lines, the mean consumption with continuous lines.

Table 1 :

 1 Parameters that depart from the default for each package.InriaB Benefits of the coarse grid correction

Figure 8 :

 8 Figure 8: Elapsed time in function of the number of cores for MaPHyS with and without the coarse grid correction.

Figure 9 :

 9 Figure 9: Elapsed time spent in each operator depending on the number of cores for a AVIP-Fluid simulation.

Figure 10 :

 10 Figure 10: Elapsed time for experiments using fluid model.

Figure 11 :

 11 Figure 11: Elapsed time for experiments using fluid model with log scale on core axis.

Figure 12 :

 12 Figure 12: Elapsed time for experiments with 30 particles per cell.

Figure 13 :

 13 Figure 13: Elapsed time for experiments with 30 particles per cell with log scale on core axis.

Figure 14 :

 14 Figure 14: Elapsed time for experiments with 60 particles per cell.

Figure 15 :

 15 Figure 15: Elapsed time for experiments with 60 particles per cell with log scale on core axis.

Inria

Figure 16 :

 16 Figure 16: Elapsed time for experiments with 120 particles per cell.

Figure 17 :

 17 Figure 17: Elapsed time for experiments with 120 particles per cell with log scale on core axis.

Figure 18 :

 18 Figure 18: Comparison of the memory consumption of the AMG and DDM solvers.

Figure 19 :

 19 Figure 19: Memory consumption with AMG and ratio of the total memory for PIC simulation with 30 particles per cell.

Inria

Figure 20 :

 20 Figure 20: Memory consumption with DDM and ratio of the total memory for PIC simulation with 30 particles per cell.

Figure 21 :

 21 Figure 21: Memory consumption with AMG and ratio of the total memory for PIC simulation with 60 particles per cell.

Figure 22 :

 22 Figure 22: Memory consumption with DDM and ratio of the total memory for PIC simulation with 60 particles per cell.

Figure 23 :

 23 Figure 23: Memory consumption with AMG and ratio of the total memory for PIC simulation with 120 particles per cell.

Figure 25 :

 25 Figure 25: Elapsed time in Poisson solution using AMG or DDM.

RR n°9423

Inria

Acknowledgements

This work is financially supported by Safran Aircraft Engines as a part of a CIFRE convention and supervized by Stephan Zurbach. This work is also financially supported by the ANR chaire industrielle POSEIDON (ANR-16-CHIN-0003-01) (LPP/CERFACS/SAFRAN). It was performed using HPC resources from GENCI (Grant A0032B10157). We are extremely thankful to the PETSc team for the great and reactive support all along this study, which allowed us to significantly improve the usage of the solver.


Pierre Bourdaud

Géraldine Lassalle

Nathalie Niquil

Marie Vagner

Importance of grey mullets on coastal food webs in a context of omega-3 deficiency: from sub-individual level to trophic dynamics

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Importance of grey mullets on coastal food webs in a context of omega-3 deficiency: What effects on the individual fish physiology? What consequences on the food web functioning?

  from sub-individual level to trophic dynamics Pierre Bourdaud, Géraldine Lassalle, Nathalie Niquil, Marie Vagner Advances in Marine Ecosystem Modelling Research (AMEMR) Symposium. 12Update of a previous Linear Inverse Model (LIM; van Oevelen et al., 2010) of MO (van der Heijden et al., and Necessity model (CaN; Planque and Moullon, 2019) of MO from LIM outputs Static mass-balanced trophic flows with uncertaintyWhat consequence on the trophic web?


Nils Krah
email: nils.krah@creatis.insa-lyon.fr

Denis Dauvergne

Jean Michel Létang

Simon Rit

Étienne Testa

Energy-adaptive calculation of the most likely path in proton CT

This note addresses an issue faced by every proton computed tomography (CT) reconstruction software: the modelling and the parametrisation of the multiple Coulomb scattering power for the estimation of the most likely path of each proton. The conventional approach uses a polynomial model parameterised as a function of depth for a given initial beam energy. This makes it cumbersome to implement a software that works for proton CT data acquired with an arbitrary beam energy or with energy modulation during acquisition.

We propose a simple way to parametrise the scattering power based on the measured proton CT list-mode data only and derive a compact expression for the most likely path (MLP) based on a conventional MLP model. Our MLP does not require any parameter. The method assumes the imaged object to be homogeneous, as most conventional MLPs, but requires no information about the material as opposed to most conventional MLP expressions which often assume water to infer energy loss. Instead, our MLP automatically adapts itself to the energy-loss which actually occurred in the object and which is one of the measurements required for proton CT reconstruction. We validate our MLP method numerically and find excellent agreement with conventional MLP methods.

Introduction

Image reconstruction in list-mode proton computed tomography (CT) is typically performed along the most likely path (MLP) that a proton has taken through the imaged object [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF][START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF][START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF][START_REF] Dedes | The role of Monte Carlo simulation in understanding the performance of proton computed tomography[END_REF]. The reason for this is that protons follow stochastic non-straight trajectories through the object, mainly because of multiple Coulomb scattering (MCS) [START_REF] Gottschalk | On the scattering power of radiotherapy protons[END_REF]. Most protons which have undergone nuclear scattering are filtered out prior to image reconstruction and MLP algorithms generally assume that only MCS bends the proton trajectory. The MLP attempts to account for MCS by providing curved integration lines to the tomographic reconstruction algorithm. It is estimated based on the position and angle measurement upstream and downstream of the object and considers the likelihood of MCS in the medium. To calculate the latter, a kinematic term essentially equivalent to the protons' scattering power at a given depth needs to be calculated which depends on the proton's kinetic energy. In most MLP formalisms, this term is integrated along depth and therefore has a significant impact on the computational speed of the reconstruction code. A practical approach to speed up the computation is to approximate the kinematic term by a polynomial as a function of depth and to perform the integration analytically [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF][START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF]. The result is again a polynomial which is computationally efficient to calculate. The polynomial coefficients need to be determined a priori for a given beam energy, e.g. based on a Monte Carlo simulation or suitable tabulated data such as NIST's PSTAR tables [START_REF] Berger | ESTAR, PSTAR, and ASTAR: Computer Programs for Calculating Stopping-Power and Range Tables for Electrons[END_REF]. While this approach is functional and has been successfully used by several authors, it has a few technical downsides: The value of the coefficients depends on the beam energy making it difficult for a reconstruction software to work with list-mode proton CT data acquired with arbitrary beam energy. Either the code would need to include a large database of such polynomial coefficients to cover an adequate range of beam energies, potentially with interpolation between energies, or the user would need to determine the coefficients and provide them as input to the reconstruction, which is arguably impractical. For image acquisition scenarios where the beam energy varies from one proton to another, e.g. when using pencil beams with varying beam energy, the approach is cumbersome because each proton requires a separate set of coefficients.

In this note, we present a method to parametrise and evaluate the kinematic term in the MLP estimation purely based on the measured data. No external parameters other than established physical constants are required. The resulting MLP expression is compact and easy to implement.

Materials and Methods

Context

This note does not propose a new formalism to estimate the MLP, but a method to parametrise existing MLP formalisms. For this purpose, we chose the MLP formalism in [START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF] because it yields a compact MLP expression.

Most MLP formalisms in proton CT use the following model for the multiple Coulomb scattering power T (u) as a function of depth in the object:

T (u) = 1 X 0 Ω 2 0 β 2 (u)p 2 (u) = 1 X 0 Ω 2 0 (E(u) + E p ) 2 c 2 (E(u) + 2E p ) 2 E 2 (u) ≡ 1 X 0 K(u), (1)
with Ω 0 = 13.6 MeV/c a constant, β = v/c the proton's velocity relative to the speed of light, p the proton's momentum, E p the proton's rest mass energy, X 0 the radiation length of the material, and E(u) the protons' kinetic energy at depth u. We defined the unit-less function K(u) to contain the energy-dependent part of the scattering power. The radiation length X 0 is depth-dependent in a heterogeneous medium, but most MLP formalisms treat objects as 2.2 Conventional parametrisation of the scattering power 2 MATERIALS AND METHODS homogeneous material and so do we. [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] also includes an empirical correction factor which we omit. Evaluating the MLP from [START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF] involves four variants of integrals over K(u),

A(u) = 1 X 0 u 0 dv K(v); C(u) = 1 X 0 u 0 dv v 0 dw K(w); B(u) = 1 X 0 u 0 dv K(v)v; D(u) = 1 X 0 u 0 dv v 0 dw K(w)w (2)
where u, v, w are variables indicating depth in the object.

The MLP, or more precisely the most likely transverse position t MLP at depth u, is calculated as

t MLP (u) = t in + θ in u + c 0 C(u) + c 1 D(u), (3)
where t in and θ in are the transverse position and angle of the proton entering the object. The factors c 0 and c 1 are independent of depth and given by

c 0 = 1 A out D out -B out C out (-B out (t out -t in -θ in u out ) + D out (θ out -θ in )) c 1 = 1 A out D out -B out C out (A out (t out -t in -θ in u out ) -C out (θ out -θ in )) , (4)
with t out and θ out the exit position and angle, respectively, and u out the exit depth. We used the short-hand notation A out = A(u out ) for brevity. Determining u out requires knowledge of the external contour of the object to be reconstructed -an assumption common to all MLP formalisms. In a proton CT reconstruction algorithm, equation 3 is evaluated independently for both perpendicular directions transverse to the beam, i.e. t MLP,x and t MLP,y (assuming the beam along z). Throughout this work, we ignore this technicality and let "transverse direction" refer to either of these two directions.

Conventional parametrisation of the scattering power

The term K(u) needs to be adequately parametrised in a proton CT reconstruction code to evaluate the integrals in equation 2. The conventional way is to obtain E(u) from a Monte Carlo simulation for a given proton beam energy and material and fit a polynomial P(u) = N n=0 a n u n to the so-derived K(u) [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF][START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF][START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF], with N = 5 suggested in [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF]. The integration can then be performed analytically. As explained in the introduction, one inconvenience of this approach is that P(u) depends on the beam energy.

New parametrisation

The new parametrisation proposed in this work starts from an observation depicted in the left panel of figure 1, namely that the reciprocal of the kinematic term in the scattering power as a function of depth in a homogeneous material is approximately linear, i.e.

1/K(u) ≈ au + b, so that K(u) ≈ 1 au + b (5)
where we choose u by convention so that protons enter the object at u in = 0. To generate figure 1, we determined the proton energy as a function of depth and calculated the scattering power as in equation 1. More specifically, we integrated u(E) = E E in 1/S w (E) dE numerically, using the stopping power S w (E) of water from NIST's PSTAR table [START_REF] Berger | ESTAR, PSTAR, and ASTAR: Computer Programs for Calculating Stopping-Power and Range Tables for Electrons[END_REF] and obtained E(u) by numerical inversion of the relation. The dots in figure 1 proton's entry and exit depths and, for this figure, we chose u out to be 2 cm less than the proton range for the given beam energy.

In practice, the coefficients a and b can be determined via equation 1 from the proton's incoming and outgoing energy,

b = 1 K(u in ) a = 1 u out 1 K(u out ) - 1 K(u in ) . (6)
The right panel in figure 1 compares the scattering power approximated according to equation 5 (dashed) with the exact scattering power from which a and b were derived (solid), as well as a polynomial fit P(u) of order N = 5 with the exact scattering power (dotted). Visually, the approximation in equation 5 is at least as good as the polynomial fit. Quantitatively, the right panel in figure 1 shows that equation 5 approximates the scattering power to within a few percent, depending on the initial energy and depth.

Using equation 5, the integrals in equation 2 can be evaluated analytically, as outlined in the Appendix, yielding:

A(u) = 1 a log a b u + 1 (7) B(u) = b a 2 a b u -log a b u + 1 (8) C(u) = b a 2 a b u + 1 log a b u + 1 - a b u (9) D(u) = b 2 a 3 1 2 a b u 2 + a b u - a b u + 1 log a b u + 1 . ( 10
)
Inserting A(u), B(u), C(u), and D(u) into equations 3 and 4 and grouping terms yields the new expression for the MLP,

t MLP (u) = t in + θ in -c 0 1 a + c 1 b a 2 u+ 1 2a c 1 u 2 + c 0 b a 2 -c 1 b 2 a 3 a b u + 1 log a b u + 1 , (11)
i.e. a quadratic function plus an additional linear-logarithmic term. This expression for the MLP does not depend on any external parametrisation of the scattering power T (u) (or the term K(u) therein) and can be fully determined via equations 1 and 6 based on measurements only, i.e. the proton's entry and exit coordinates as well as the entry and exit energy.

Implementation

The MLP in equation 11 can be implemented according to the following scheme: For each proton, project the entry and exit coordinates (position and angle) onto the external contour of the object to be reconstructed. This yields u in , u out , t in , t out , θ in and θ out . Next, determine the parameters a and b via equation 6 based on the proton's entry and exit energy. Calculate and store the parameters c 0 and c 1 (equation 4). Finally, calculate the most likely transverse position t MLP at all desired depths u using a, b, c 0 , c 1 . It can be shown that the case a = 0 (no energy loss) does not constitute a true singularity because in that case K(u) = 1/b and all integrals simplify. This case must be handled separately to avoid numerical instability.

Validation and simulation

We validated the new MLP parametrisation by comparison with conventional MLP expressions. In particular, we chose the methods from [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] and [START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF]. The former is well-established and was one of the first to appear in the literature on proton CT and the latter is the method to which the adaptive parametrisation in this work is applied.

To quantify how well the new MLP agrees with conventional MLP expressions, we performed a GATE Monte Carlo simulation [START_REF] Sarrut | A review of the use and potential of the GATE Monte Carlo simulation code for radiation therapy and dosimetry applications[END_REF] of protons propagating through water and traced their true trajectory in 1 mm depth intervals. We chose water because this is the underlying assumption of both reference MLPs. We used the emstandard_opt4 physics list to ignore nuclear interactions as these are not considered by the MLP formalisms. In practice, most nuclear events are filtered out from proton CT data, but that would have been an unnecessary complication for the purpose of this validation.

For each proton, we calculated the MLP using the three different MLP methods. We compared the MLPs among each other as well as with the true paths from Monte Carlo. Regarding the latter, we generated the distribution of the difference between MLP estimate and true position, t MLP -t MC , at a certain depth.

Results

The upper left panel in figure 2 shows a few example MLPs in a 20 cm thick box of water for protons of 200 MeV and calculated using the MLP method from this work as well as those from [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] and [START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF]. The lower left panel illustrates how much these MLPs vary among each other. The right panels depict the error of t MLP compared to the transverse position from Monte Carlo for 200 MeV (upper) and 150 MeV (lower) protons. The position was evaluated at the central depth and the water boxes were 20 cm and 12 cm thick, respectively.

The differences among MLPs are no greater than a few tens of micrometers and much smaller than the statistical uncertainty of the MLP estimate itself. The MLP uncertainty, which directly influences the spatial resolution in the reconstructed image, is essentially identical for all MLPs. We also verified these results at other beam energies and depths and came to the same conclusions.

Discussion and Conclusion

In this work, we have presented a new method to parametrise the energy dependence in the scattering power which is part of the most likely path formalism in proton CT. The method relies on the fact that the reciprocal of the scattering power as a function of depth in an object is approximately linear, which we have verified for the range of beam energies relevant for proton CT. Our new MLP parametrisation automatically adapts to the energy-loss of a proton in a given object because it uses the actually measured exit energy.

From the practical point of view, the new MLP parametrisation presented in this work greatly simplifies the implementation of a proton CT reconstruction software. Indeed, with a conventional MLP, either the software or the user need to provide the coefficients of the polynomial which approximates the scattering power as a function of depth. The new method can calculate the MLP for protons of any initial energy only relying on the proton CT list-mode data, i.e. proton coordinates and energies. This also makes reconstruction of images acquired with energy-modulated proton beams [START_REF] Dickmann | Proof of concept image artifact reduction by energy-modulated proton computed tomography (EMpCT)[END_REF] technically easier than it would be with the conventional MLP methods.

It is worth remarking that the scattering power term K(u) (equation 1) in a homogeneous object is also well approximated by a power law of the form K(u) = (R -u) δ , where R is the proton range in the given material and and δ are model parameters, similar to the model used in [START_REF] Bortfeld | An analytical approximation of depth -dose distributions for therapeutic proton beams[END_REF]. In a known material, and δ could be determined using K(u in ) and K(u out ) (which can be calculated from E in and E out ). For the example data in water shown in figure 1, we find δ ≈ -1.05, which is coherent with our inverse linear model. In practice, however, the material of an imaged object and the proton range R therein is not known. One could use the range in water as an approximation, but that would incorrectly represent the relation K(u) in materials which are on average denser or less dense than water. In mathematical terms, the power law model in an unknown material has three free parameters ( , δ, R), but only two data points (E in , E out ) are available, so the model remains underdetermined.

Interestingly, the new MLP expression in equation 11 does not make any assumption about the material of the imaged object, contrary to the conventional MLP parametrisation. For example, in an object which is less dense than water such as the lungs, protons loose less energy and K(u) does not increase as much as it would in water. Consequently, the slope parameter a in equation 6 is smaller than it would be in the same object filled with water. The conventional parameterisation, which uses a polynomial P(u) fitted to data simulated in water, would overestimate the energy loss in the less dense material. In this context, we also note that the radiation length X 0 (and other constant factors) which appear in the scattering power actually cancel out in the MLP expression, specifically in the terms c 0 C(u) and c 1 D(u) in equation 3. In other words, in a homogeneous medium, the MLP depends only on K(u) and not on the radiation length.

We point out that our MLP expression does not account for heterogeneous material composition along the particle path, and neither do the conventional MLP from [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] and similar variants. There are MLP methods available in the literature which use a priori knowledge about the object to account for inhomogeneous tissue composition [START_REF] Collins-Fekete | Extension of the Fermi-Eyges most-likely path in heterogeneous medium with prior knowledge information[END_REF][START_REF] Brooke | An inhomogeneous most likely path formalism for proton computed tomography[END_REF]. A comparison of such methods, both in terms of image quality and computation speed, was beyond the scope of this work.

The purpose of this work was to find a convenient way of parametrising the scattering power in proton CT reconstruction. We concentrated on a specific MLP formalism, namely the one in [START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF], because it leads to a compact expression for the MLP. We point out, however, that other MLP formalisms could also use equations 5 and 6 when evaluating scattering integrals, e.g. equations 7-9 in [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF].

  Figure1: Left: Inverse scattering power calculated according to equation 1 (solid). The linear approximation (dashed line) was determined as the line connecting the two dots which represent the scattering power at the object entry and exit. Centre: Reciprocal of the left plot. The dotted line was obtained as a polynomial fit of order 5 to the solid line between the two dots. Right: Relative difference of the fitted scattering powers to the reference.

  Upper left: Example of most likely paths for 200 MeV protons in water calculated with different MLP methods. Lower left: Difference between the MLPs. Right: Distributions of the difference between t MLP and true position (from Monte Carlo) at the centre of a 20 cm (upper) and 12 cm (lower) thick water box. Beam energies are indicated in the graphics. The vertical lines indicate ± one root mean square error of the distributions, with specific numbers provided in the legend.

				Occurence		Adaptive MLP, RMSE=0.685 mm Krah 2019, PMB, RMSE=0.684 mm Schulte 2008, MedPhys, RMSE=0.683 mm	200 MeV
	Difference in trans. pos. in mm	0.04 0.02 0.00 0.02	Adaptive -Schulte Krah -Schulte Krah -Adaptive	Occurence	3	2 Adaptive MLP, RMSE=1.098 mm Difference between tMLP and true path in mm 1 0 1 Krah 2019, PMB, RMSE=1.194 mm Schulte 2008, MedPhys, RMSE=1.093 mm	2 150 MeV	3
	Figure 2:

Acknowledgements

The work of Nils Krah was supported by ITMO-Cancer (CLaRyS-UFT project). This work was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). We also gratefully acknowledge the support from the CNRS/IN2P3 Computing Center (Lyon/Villeurbanne -France), for providing a significant amount of the computing resources needed for this work.

Appendix

Evaluating the integral terms A(u), B(u), C(u), and D(u) in equation 2 with the new parametrisation of the scattering power involves two similar types of integral expressions which can be found in common mathematical handbooks, namely

where we have assumed u in = 0 for simplicity and without loss of generality. A(u) is directly equivalent to equation 12. C(u) is obtained by integrating equation 12 and using the substitution v ≡ av/b + 1, so that

which is solved recalling that log(x) dx = x log(x) -x = x (log(x) -1) .

B(u) is directly given by equation 13 and D(u) follows from it by integration and substitution v ≡ av/b + 1, in analogy to C(u).