Mohd Pierre Jolivet
email: pierre.jolivet@enseeiht.fr

Afeef Badri

Y Favennec

M A Badri

Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning

Keywords: radiative transfer, matrix-free solver, algebraic multigrid

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

Introduction

Radiative transfer is the phenomenon of energy transfer in the form of electromagnetic radiation, often characterized by its radiative intensity I. This intensity within any participating medium is affected by absorption, emission, reflection, and scattering processes. The radiative transfer equation (RTE) models such processes at the continuous level. The RTE is commonly used in the fields of thermal radiation [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF][START_REF] Ghattassi | Galerkin method for solving combined radiative and conductive heat transfer[END_REF], neutronics [START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF][START_REF] Lewis | Computational methods of neutron transport[END_REF], astrophysics [START_REF] Richling | Radiative transfer with finite elements -I. Basic method and tests[END_REF][START_REF] Meinköhn | Radiative transfer with finite elements -II. Lyα line transfer in moving media[END_REF], and optical imaging [START_REF] Abdoulaev | Three-dimensional optical tomography with the equation of radiative transfer[END_REF][START_REF] Tarvainen | Gauss-Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation[END_REF], to name a few. This equation, in its steady-state monochromatic form, is stated as the following integro-differential equation:

( s • ∇ + β( x)) I( x, s) = σ( x) S I( x, s ′ )φ( s, s ′ )d s ′ + κ( x)I b ( x), (1)
where x = (x, y, z) ∈ Ω ⊂ R 3 and d s is the differential solid angle around the direction s = [cos ϕ sin θ, sin ϕ sin θ, cos θ] T , (ϕ, θ) ∈ [0, 2π] × [0, π] being the azimuthal and the zenith angles, respectively. β is the extinction coefficient computed as the sum κ + σ, where κ (resp. σ) is the absorption (resp. scattering) coefficient. These three coefficients are piecewise constant in Ω and their values depend on the various participating media in the domain. φ is a scattering phase function that quantifies the probability of a photon traveling in s ′ being scattered towards s. Finally, I b is the radiation from a black body specified by Planck's law.

Efficient solution strategy for the monochromatic steady-state RTE is still an active area of research for cases that involve scattering and reflections in complex geometries [START_REF] Kanschat | A Robust Multigrid Preconditioner for S N DG Approximation of Monochromatic, Isotropic Radiation Transport Problems[END_REF][START_REF] Ruh | Radiative transport in large arteries[END_REF]. The main challenge, for methods presented in this article, is the presence of in-scattering from the surface integral in eq. [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF], which leads to coupled systems that are difficult to solve due to memory and convergence issues. Associated with the multigroup approximation [START_REF] Till | Application of linear multifrequency-grey acceleration to preconditioned Krylov iterations for thermal radiation transport[END_REF][START_REF] Haut | A new multigroup method for cross-sections that vary rapidly in energy[END_REF] or spectral models [START_REF] Clements | Evaluation of FSK models for radiative heat transfer under oxyfuel conditions[END_REF][START_REF] Kez | A comprehensive evaluation of different radiation models in a gas turbine combustor under conditions of oxy-fuel combustion with dry recycle[END_REF], monochromatic RTE solvers are a stepping stone for designing multi-frequency RTE solvers in different scientific communities, e.g., combustion [START_REF] Viskanta | Radiation heat transfer in combustion systems[END_REF][START_REF] Poitou | Analysis of radiation modeling for turbulent combustion: development of a methodology to couple turbulent combustion and radiative heat transfer in LES[END_REF].

In most of this paper, the boundary ∂Ω is assumed to be nonreflective. If n is the outward normal vector at a given point of ∂Ω, the following inflow boundary condition is used to provide closure to eq. ( 1): One needs two schemes for discretizing eq. ( 1) both in space x and angles s. In this work, we couple the finite element method (FEM) for the spatial discretization with the discrete ordinates method (DOM) based on the sphere triangulation for the angular discretization [START_REF] Kanschat | Solution of radiative transfer problems with finite elements[END_REF][START_REF] Camminady | Highly uniform quadrature sets for the discrete ordinates method[END_REF]. The terminology DOM is used throughout the paper.

I = I inflow ∀( x, s) ∈ Γ -= {( x, s) : x ∈ ∂Ω, s • n < 0}. (2)
Angular discretization. The DOM, here performed according to the unit sphere triangulation, either based on Thurgood's method [START_REF] Thurgood | The T N quadrature set for the discrete ordinates method[END_REF] or icosahedron refinements [START_REF] Richling | Radiative transfer with finite elements -I. Basic method and tests[END_REF], cf. fig. 1b, provides a set of directions and weights

{( s k , ω k )} N d k=1
. Such a scheme is second-order accurate when combined with piecewise linear finite elements in space [START_REF] Badri | Vectorial finite elements for solving the radiative transfer equation[END_REF]. The continuous intensity is now approximated by a set of discrete intensities {I k } N d k=1 . The surface integral in eq. ( 1) is thus evaluated as, for all discrete directions k ∈ 1; N d :

S I( x, s )φ( s k , s )d s ≈ N d j=1 ω j I j ( x)φ k,j .
The integro-differential eq. ( 1) may now be stated as a set of N d continuous partial differential equations:

∀k ∈ 1; N d , ( s k • ∇ + β( x))I k ( x) = σ( x) N d j=1 ω j I j ( x)φ k,j + κ( x)I b ( x). (3)
In a similar fashion the inflow boundary condition from eq. ( 2) has to be rewritten in discrete angular form.

In this paper, collimated beams of intensity I inflow impinge the boundary along prescribed entrant directions in L ⊂ 1; N d . Hence, ∀k ∈ 1; N d :

I k = I inflow ∀ x ∈ ∂Ω : s k • n < 0, if k ∈ L, 0 otherwise. ( 4
)
Spatial discretization. With the FEM, given a spatial mesh Ω h , cf. fig. 1a, and a set of appropriate functional

spaces V h k N d
k=1 , the following variational equations may be defined using the Galerkin method. The finite element trial functions for the discrete radiative intensities {I h k ∈ V h k } N d k=1 have to satisfy for all test functions

{J h k ∈ V h k } N d k=1 : ∀k ∈ 1; N d , Ω h ( s k • ∇ + β( x))I h k ( x) -σ( x) N d j=1 ω j I h j ( x)φ k,j J h k ( x)d x = Ω h κ( x)I b ( x)J h k ( x)d x,
with the boundary conditions defined in eq. ( 4). In this paper, it is assumed that for all directions, the same functional space V h is used to search for all the discrete radiative intensities, i.e., for all k in 1;

N d , V h k = V h .
In particular, V h is the space of piecewise linear functions throughout the paper. Because of the convectiondominant term in eq. ( 1), an artificial diffusivity is introduced, and the following test functions derived from the SUPG method [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] are used: J h k + γ s k • ∇J h k . An effective piecewise constant value of γ was determined in prior works [START_REF] Avila | Spatial approximation of the radiation transport equation using a subgrid-scale finite element method[END_REF]:

γ( x) = 4 h 2 + β( x) 2 + 4σ( x) hN d -1/2
, where h is the size of each spatial element. By integrating by parts, the system of variational formulations becomes: ∀k ∈ 1; N d ,

Ω h β( x)I h k ( x) -σ( x) N d j=1 ω j I h j ( x)φ k,j J h k ( x) + γ( x) s k • ∇J h k ( x) -s k • ∇J h k ( x) I h k ( x) + γ( x) s k • ∇I h k ( x) s k • ∇J h k ( x) d x = Ω h κ( x)I b ( x) J h k ( x) + γ( x) s k • ∇J h k ( x) d x - Γ - s k • n J h k ( x)I inflow ( x)d x - ∂Ω h \Γ - s k • n I h k ( x)J h k ( x)d x. (5)
Now that the RTE has been discretized both in space and angles, it is possible to highlight the major contributions of this paper. Indeed, we present:

❼ an efficient way to assemble linear systems arising from the discretization of eq. ( 5), ❼ a scalable approach for solving such linear systems using a mixed matrix-ready/matrix-free approach, ❼ how to handle physics involving specular reflections.

The paper is organized as follows. In section 2, related work and limitations of current implementations of radiative transfer equation solvers are presented. In section 3, the flexibility of modern finite element domain specific languages is leveraged to efficiently discretize eq. ( 5) numerically. In section 4, methods to iteratively solve the large systems generated in the previous section on thousands of processes are presented. Validation and scaling analysis of our solver is presented in section 5. The capabilities of our solver when dealing with specular reflections are shown in section 6.

Related work

Most RTE solvers are either physics-based (stochastic) or deterministic. Concerning the former, Monte

Carlo [START_REF] Howell | Monte Carlo solution of thermal transfer through radiant media between gray walls[END_REF] or ray-tracing methods [START_REF] Liu | Discrete curved ray-tracing method for radiative transfer in an absorbing-emitting semitransparent slab with variable spatial refractive index[END_REF] are well-suited to parallel computing [START_REF] Badal | A real-time radiation dose monitoring system for patients and staff during interventional fluoroscopy using a GPU-accelerated Monte Carlo simulator and an automatic 3D localization system based on a depth camera[END_REF]. However, computations involving optically thick media are numerically challenging for such methods [START_REF] Min | Radiative transfer in very optically thick circumstellar disks[END_REF]. Moreover, coupling these with other physics, e.g., fluid mechanics, is also demanding [START_REF] Fournier | Radiative, conductive and convective heattransfers in a single Monte Carlo algorithm[END_REF].

For deterministic solvers, historically, the surface integral from eq. ( 1) was dealt with by the introduction of the DOM [START_REF] Chandrasekhar | Radiative transfer[END_REF]. Though new approaches for discretizing in angles exist [START_REF] Fujiwara | An Accurate Quadrature Rule on the Sphere for Fast Computation of the Radiative Transport Equation[END_REF], it is still nowadays the most used technique for semi-discretizing the integro-differential RTE into a set of coupled PDEs on which standard spatial discretization schemes are applicable. In the field of radiative transfer, the DOM is often combined with the finite volume method (FVM) [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF].

As an alternative, one may also use the FEM. Many recent works in the radiative transfer community support its use together with the DOM to improve the efficiency of radiative transport equation solvers [START_REF] Richling | Radiative transfer with finite elements -I. Basic method and tests[END_REF][START_REF] Lu | A parallel adaptive finite element method for the simulation of photon migration with the radiative transfer-based model[END_REF][START_REF] Widmer | An efficient sparse finite element solver for the radiative transfer equation[END_REF][START_REF] Kanschat | A Robust Multigrid Preconditioner for S N DG Approximation of Monochromatic, Isotropic Radiation Transport Problems[END_REF], or to deal with complex geometries [START_REF] Hardy | 3D numerical modelling of the propagation of radiative intensity through a X-ray tomographied ligament[END_REF][START_REF] Kang | Finite element formulation of the first-and second-order discrete ordinates equations for radiative heat transfer calculation in three-dimensional participating media[END_REF] or coupled physics [START_REF] Ghattassi | Galerkin method for solving combined radiative and conductive heat transfer[END_REF].

Regardless of the spatial discretization, the fully discrete set of linear equations is commonly solved by a direction-by-direction iterative scheme known as source iteration [START_REF] Fiveland | Discrete ordinates solutions of the radiative transport equation for rectangular enclosures[END_REF][START_REF] Adams | Fast iterative methods for discrete ordinates particle transport calculations[END_REF][START_REF] Mihalas | Stellar atmospheres[END_REF]. It provides a memory-efficient strategy needed to handle the high-dimensionality of RTE in its discrete form eq. ( 5). Then, an optimal sweeping sequence is performed to compute solution of linear systems. Source iteration is widely popular in parallel deterministic solvers: Uintah [START_REF] Kumar | Scalable data management of the Uintah simulation framework for next-generation engineering problems with radiation[END_REF] from University of Utah, Denvo [START_REF] Evans | Denovo: a new three-dimensional parallel discrete ordinates code in SCALE[END_REF] from ORNL, Ardra [START_REF] Hannebutte | Ardra: Scalable parallel code system to perform neutron-and radiation-transport calculations[END_REF] from LLNL, PDT [START_REF] Adams | Provably optimal parallel transport sweeps on regular grids[END_REF] from TAMU. However, they are most often optimized for regular geometries that have to be (semi-)structurally meshed to facilitate optimal marching (sweeping sequence) [START_REF] Pautz | Parallel Deterministic Transport Sweeps of Structured and Unstructured Meshes with Overloaded Mesh Decompositions[END_REF], which is also sensitive to boundary conditions, unless, e.g., orthogonality conditions between reflective surfaces hold [START_REF] Adams | Provably optimal parallel transport sweeps on semi-structured grids[END_REF]. These schemes may be extended to unstructured grids [START_REF] Plimpton | Parallel S N sweeps on unstructured grids: algorithms for prioritization, grid partitioning, and cycle detection[END_REF][START_REF] Asllanaj | Solution of radiative heat transfer in 2D geometries by a modified finite volume method based on a cell vertex scheme using unstructured triangular meshes[END_REF][START_REF] Wang | Diffusion synthetic acceleration for high-order discontinuous finite element S N transport schemes and application to locally refined unstructured meshes[END_REF][START_REF] Haut | An Efficient Sweep-Based Solver for the S N Equations on High-Order Meshes[END_REF][START_REF] Vermaak | Massively parallel transport sweeps on meshes with cyclic dependencies[END_REF] but with understandably lower parallel efficiencies.

Krylov methods have been often used with source iteration. They are particularly suited to situations where there is a tight coupling between the discrete directions, e.g., strongly scattering or reflecting media [START_REF] Seaıd | Efficient numerical methods for radiation in gas turbines[END_REF].

Currently, there are different ways to leverage Krylov methods for solving the discrete RTE, e.g.:

❼ direct application on the full linear system yielded by the spatial and DOM discretizations [START_REF] Godoy | On the use of flux limiters in the discrete ordinates method for 3D radiation calculations in absorbing and scattering media[END_REF];

❼ in conjunction with sweep-based iterative schemes recasted as preconditioned operators [START_REF] Patton | Application of preconditioned GMRES to the numerical solution of the neutron transport equation[END_REF].

The latter choice is often less memory-demanding since the workload can be implemented direction-by-direction.

Contrary to direction-by-direction schemes, all the discrete radiative intensities are considered simultaneously in the present work. Consequently, RTE solvers based on direct application of Krylov methods are most often limited by the high dimensionality of eq. ( 5) and thus require a large amount of memory to store the according linear system. As such, these solvers usually tackle small-scale problems [START_REF] Krishnaprakas | Use of GCG methods for the efficient solution of matrix problems arising from the FVM formulation of radiative transfer[END_REF][START_REF] Godoy | On the use of flux limiters in the discrete ordinates method for 3D radiation calculations in absorbing and scattering media[END_REF] (up to 512 MPI processes on a regularly-decomposed Cartesian grid of 262 thousand elements with up to 120 directions) [START_REF] Lorca | Multilevel Schwarz Preconditioners for Singularly Perturbed Symmetric Reaction-Diffusion Systems[END_REF] (up to four directions in 2D) or have to rely on tools such as adaptive mesh refinement in space [START_REF] Charest | Solution of the equation of radiative transfer using a Newton-Krylov approach and adaptive mesh refinement[END_REF] or in angles [START_REF] Adigun | A Haar wavelet method for angularly discretising the Boltzmann transport equation[END_REF], which may complexify a distributed solver due to the need of proper load balance. To deal with larger problems, at least two modern ways to solve the RTE have emerged. First, the use of reduced order models [START_REF] Soucasse | An angular reduced order model for radiative transfer in non grey media[END_REF] or sparse grids may be used to deal efficiently with the "curse of dimensionality" [START_REF] Widmer | Sparse adaptive finite elements for radiative transfer[END_REF][START_REF] Grella | Sparse tensor spherical harmonics approximation in radiative transfer[END_REF] (up to 3 • 10 6 unknowns solved sequentially). Second, Jacobian-free Newton-Krylov methods [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF] have been extended to the RTE [START_REF] Godoy | Parallel Jacobian-free Newton-Krylov solution of the discrete ordinates method with flux limiters for 3D radiative transfer[END_REF] (up to 512 MPI processes and 8 threads per MPI process on a Cartesian grid of 65 thousand elements with up to 120 directions). These two methods are based on simplifications of the RTE by approximating either the interactions between the unknowns in space and in angles for the latter, or the global operator for the former.

Efficient assembly

Nowadays, the finite element method is available in a large number of professional or academic domain specific (embedded) languages: deal.II [START_REF] Bangerth | II: a general-purpose object-oriented finite element library[END_REF], FEniCS [START_REF] Logg | Automated solution of differential equations by the finite element method: the FEniCS book[END_REF], FreeFEM [START_REF] Hecht | New development in FreeFem++[END_REF], etc. The goal of the next section is to show that these software typically cannot cope with the dimension complexity of the RTE if the corresponding linear system is assembled explicitly. However, we will use one of them, FreeFEM, to discretize auxiliary problems and show how this relates to the original PDE.

Vectorial finite elements

Naively, N d linear systems need to be assembled corresponding to the N d variational formulations of eq. ( 5).

However, by introducing the vectorial finite element space V h defined as the Cartesian product of N d finite element spaces V h , eq. ( 5) may be formulated in terms of a single vectorial unknown I h and a single vectorial

test function J h . Let I h = [I h 1 • • • I h N d ] T , J h = [J h 1 • • • J h N d ] T , S = [ s 1 • • • s N d ] T , and {Φ i,j = ω j φ i,j } N d i,j=1
. If we denote by ⊗ (resp. ⊙) the entrywise product (resp. scalar product) of two vectorial functions, and by ⋄ the product of a matrix by a vectorial function, then, the bilinear part of eq. ( 5) becomes, for all J h ∈ V h :

a(I h , J h ) = Ω h β( x)I h ( x) -σ( x)Φ ⋄ I h ( x) ⊗ J h ( x) + γ( x)S ⊙ ∇J h ( x) (6)
+ γ( x) S ⊙ ∇I h ( x) ⊗ S ⊙ ∇J h ( x) -S ⊙ ∇J h ( x) ⊗ I h ( x) d x + ∂Ω h \Γ - S ⊙ n ⊗ J h ( x) ⊗ I h ( x)d x.
The same goes for the linear part of eq. ( 5) such that one may write for all J h ∈ V h : a(I h , J h ) = l(J h ).

Matrix-ready assembly

Let the linear system yielded by the discretization of the previous equation be written as:

AX = B. (7)
If the number of degrees of freedom in V h is n, since the number of components in the vectorial finite element space V h is equal to the number of directions N d in the angular mesh, A is of order andB, then:

N d • n. Furthermore, if it is assumed that each component of V h is interleaved in A, X,
❼ the sparse matrix A is made of dense N d × N d blocks, ❼ B = [b 1,1 , . . . , b N d ,1 , b 1,2 , . . . , b N d ,2 , . . . , b 1,n , . . . , b N d ,n ] T .
Increasing the number of directions for the DOM means that the number of nonzero coefficients in A will increase as well. In practice, the efficiency of two traditional finite element languages, FEniCS and FreeFEM, is assessed in fig. 2. For this test case, the RTE is discretized in a unit two-dimensional square using a fixed unstructured spatial mesh made of 6,768 triangles. The resulting finite element space V h has 3,487 degrees of freedom. The time spent to assemble the corresponding linear systems, while increasing the number of directions N d , is displayed. The results have been averaged over five consecutive runs, performed after cache warming.

With FEniCS, prior to any computations, a first run to JIT compile the solver is also performed. This is done in sequential, on a single core of an Intel Core i7-6567U. Clearly, these timings are not satisfactory since both the spatial and angular meshes are tiny. They could not be used to resolve any nontrivial physics.

Semi-matrix-free assembly

The two previous general-purpose languages do not succeed at exploiting the structure of eq. ( 6). Instead, they assemble the matrix from eq. ( 7) coefficient by coefficient, which has thus a cost proportional to the order of the linear system N d • n. If it is assumed that the sparse matrix A does not have to be assembled explicitly, this cost may be greatly reduced. Indeed, if one looks for example at the following term from the bilinear form eq. ( 6):

Ω h β( x)I h ( x) ⊗ J h ( x)d
x, assembled as A 1 , one can instead discretize this bilinear form

Ω h β( x)I h ( x)J h ( x)d x, assembled as A 1
, where both I h and J h are in the finite element space V h (of dimension n), and not in V h (of dimension N d • n). The action of A 1 on a vector X of order N d • n, is nothing else than the action of A 1 on all N d components of order n of X:

∀(k, i) ∈ 1; N d × 1; n , (A 1 X) k+N d •(i-1) = A 1 [x k,1 • • • x k,n ] T i .
Now, looking at:

a 2 (I h , J h ) = Ω h β( x)I h ( x) ⊗ γ( x)S ⊙ ∇J h ( x) -S ⊙ ∇J h ( x) ⊗ I h ( x) d x = Ω h (γ( x)β( x) -1)I h ( x) ⊗ S ⊙ ∇J h ( x) d x,
assembled as A 2 , it is instead possible to assemble the operators A 2a , A 2b , and A 2c , corresponding to the discretization of the following bilinear forms:

a 2a (I h , J h ) = Ω h (γ( x)β( x) -1)I h ( x) ∂J h ( x) ∂x d x a 2b (I h , J h ) = Ω h (γ( x)β( x) -1)I h ( x) ∂J h ( x) ∂y d x a 2c (I h , J h ) = Ω h (γ( x)β( x) -1)I h ( x) ∂J h ( x) ∂z d x.
The action of A 2 on a vector X of order N d • n may be computed as:

∀(k, i) ∈ 1; N d × 1; n , (A 2 X) k+N d •(i-1) = s k,1 A 2a + s k,2 A 2b + s k,3 A 2c [x k,1 • • • x k,n ] T i ,
where s k,j is the jth coordinate of the discrete direction s k , for j ∈ 1; 3 . Another part of the bilinear form is:

a 3 (I h , J h ) = Ω h γ( x) S ⊙ ∇I h ( x) ⊗ S ⊙ ∇J h ( x) d x,
assembled as A 3 . It is instead possible to assemble the operators A 3a , A 3b , A 3c , A 3d , A 3e , and A 3f , corresponding to the discretization of the following bilinear forms:

a 3a (I h , J h ) = Ω h γ( x) ∂I h ( x) ∂x ∂J h ( x) ∂x d x a 3b (I h , J h ) = Ω h γ( x) ∂I h ( x) ∂y ∂J h ( x) ∂y d x a 3c (I h , J h ) = Ω h γ( x) ∂I h ( x) ∂z ∂J h ( x) ∂z d x a 3d (I h , J h ) = Ω h γ( x) ∂I h ( x) ∂x ∂J h ( x) ∂y + ∂I h ∂y ∂J h ∂x d x a 3e (I h , J h ) = Ω h γ( x) ∂I h ( x) ∂x ∂J h ( x) ∂z + ∂I h ∂z ∂J h ∂x d x a 3f (I h , J h ) = Ω h γ( x) ∂I h ( x) ∂y ∂J h ( x) ∂z + ∂I h ∂z ∂J h ∂y d x.
The action of A 3 on a vector X of order N d • n may be computed as: 7) for a two-dimensional domain with P 1 finite elements either explicitly or using our semi-matrix-free approach.

∀(k, i) ∈ 1; N d × 1; n , (A 3 X) k+N d •(i-1) = s 2 k,1 A 3a + s 2 k,2 A 3b + s 2 k,3 A 3c + s k,1 s k,2 A 3d + s k,1 s k,3 A 3e + s k,2 s k,3 A 3f [x k,1 • • • x k,n ] T i . 2
Eventually, the last part of eq. ( 6) with a volume integral assembled as A 4 to deal with is:

a 4 (I h , J h ) = Ω h (σ( x)Φ ⋄ I h ( x)) ⊗ (J h ( x) + γ( x)S ⊙ ∇J h ( x))d x,
for which we instead assemble the operators A 4a , A 4b , A 4c , and A 4d , corresponding to the discretization of the following bilinear forms:

a 4a (I h , J h ) = Ω h σ( x)I h ( x)J h ( x)d x a 4b (I h , J h ) = Ω h σ( x)γ( x)I h ( x) ∂J h ( x) ∂x d x a 4c (I h , J h ) = Ω h σ( x)γ( x)I h ( x) ∂J h ( x) ∂y d x a 4d (I h , J h ) = Ω h σ( x)γ( x)I h ( x) ∂J h ( x) ∂z d x.
Now, one can formally view the column vector X of dimension N d • n as a column-major matrix X of dimensions

N d × n. Let Y = ΦX be of dimensions N d × n as well.
The action of A 4 on a vector X of order N d • n may be computed as:

∀(k, i) ∈ 1; N d × 1; n , (A 4 X) k+N d •(i-1) = s k,1 A 4b + s k,2 A 4c + s k,3 A 4d + A 4a [Y k,1 • • • Y k,n ] T i .
Additional work is needed to take into account the surface integral from eq. ( 6). It cannot be easily decomposed, neither with FreeFEM nor with FEniCS, in a tensor form, like for the other parts of the variational formulation, because of the Boolean expression in the definition of Γ -in eq. ( 2). Thus, for the last line of eq. ( 6), the operators {S k } N d k=1 associated with the following bilinear forms have to be evaluated:

∀k ∈ 1; N d , ζ k (I h , J h ) = ∂Ω h \Γ - s k • n I h k ( x)J h k ( x)d x. (10)
Once again, this would require a number of finite element assemblies that grows linearly with N d . Hence, it is not tractable for large number of directions. To keep a low number of finite element assemblies, we color the

N be boundary elements of ∂Ω h \ Γ -into subsets {F k } χ k=1
such that two boundary elements, i.e., faces in 3D or edges in 2D, of F k , ∀k ∈ 1; χ , do not share a single vertex. This coloring is computed using a simple greedy algorithm. The following bilinear forms are then assembled into matrices

{S ′ k } χ k=1 : ∀k ∈ 1; χ , ζ ′ k (I h , J h ) = F k I h ( x)J h ( x)d x. (11)
Furthermore, given a direction k ∈ 1; N d , the integrand s k • n is constant per boundary element, so these scalar products can be precomputed and stored in ψ ∈ R N d ×N be . For a spatial degree of freedom i ∈ 1; n , π(i)

returns the subset of boundary elements to which i is associated. For a boundary element j ∈ 1; N be , π ′ (j) returns the color of the given face. The assembled operators {S k } N d k=1 which are the discretized counterparts of eq. ( 10) are thus computed as follows: Keep in mind that the matrices {S ′ k } χ k=1 are very sparse since they correspond to the discretization of a surface variational formulation. As a consequence, the matrices {S k } N d k=1 are also sparse, even though eq. ( 12) uses dense notations for the row and column indices.

∀(k, i, j) ∈ 1; N d × 1; n 2 , S k i,j = N d k=1 m∈π(i)∩π(j) ψ k,m S ′ π ′ (m) i,j . (12)
For realistic geometries, it is expected that the number of unknowns associated with degrees of freedom on the surface ∂Ω h \ Γ -will be lower than those associated with degrees of freedom in the volume Ω h . The cost of evaluating these surface integrals is thus expected to be much lower than the cost of evaluating the fixed number of reference matrices with volume integrals. Using the same test case as for matrix-ready assemblies in section 3.2, FreeFEM is now used to assemble the reference matrices A 1 , A 2a,2b , A 3a,3b,3c,3d , and A 4a,4b,4c , as well as {S ′ k } χ k=1 . Note that for two-dimensional test cases, the matrices A 2c , A 3c,3e,3f , and A 4d do not exist (only nine reference matrices). As displayed in fig. 2, the semi-matrix-free approach is clearly faster than matrix-ready assemblies by multiple orders of magnitude. Moreover, the cost of this part of the solver is now only loosely dependent on N d because of the surface integrals. If the coefficient matrix of order N d • n from eq. ( 7) is not needed explicitly, we have derived a way to compute implicitly its action on a vector by assembling only reference matrices of order n.

Iterative scheme/preconditioning

In the previous paragraph, it was shown how one may assemble auxiliary matrices of a much lower dimension than the original coefficient matrix from eq. ( 7), to discretize the RTE using the DOM coupled with the FEM.

The resulting system, however, still needs to be solved. This is the focus of this section.

Matrix-vector products

Because the coefficient matrix from eq. ( 7) is not known explicitly, it seems natural to use for the solution phase methods that only require the action of the matrix on vectors. Krylov methods [START_REF] Van Der | Iterative Krylov methods for large linear systems[END_REF] which are widely used in conjunction with the FEM seem like good candidates. Because they iterate by computing matrix-vector products, it is of paramount importance to design an efficient routine for our semi-matrix-free approach. To do so, a key observation is that, thanks to the compact support of finite element basis functions, the sparsity pattern of all fourteen reference matrices is the same. Thus, no matter the format used to store these sparse matrices, their pattern may be traversed only once. In the case of FreeFEM (and FEniCS), sparse matrices may be stored using the compressed sparse row (CSR) format, which is also widely used in common linear algebra backends such as PETSc [START_REF] Balay | PETSc web page[END_REF]. Hence, if all reference matrices have m nonzero coefficients, we will assume that two arrays I (of size n + 1) and J (of size m) are used to traverse the matrices. We also denote by A ⋆ (i) the ith nonzero coefficient of the reference matrix ⋆. The same strategy and notations are used for the matrices associated with eq. ( 10) (m S , I S , J S , and S ⋆ (i)). All the previous operations from section 3.3 may now be implemented as in fig. 3. If A from eq. ( 7) would have been explicitly assembled, its action on a vector would require 2 • m • N 2 d double-precision floating-point operations. Our semi-matrix-free matrix-vector product performs the following number of FLOP:

N d • n • (N d -1) line 2 of fig. 3 + 28 • m • N d lines 4-11 + 2 • m S • N d lines 12-16 . ( 13
)
Its efficiency is assessed in table 1. This time, a cube meshed with 235 thousand tetrahedra is used. The corresponding finite element space has 44 thousand degrees of freedom. It is assumed that the matrices A ⋆ and S ⋆ have already been assembled. The following results have been averaged over five consecutive runs, performed after cache warming. The time spent in the dense matrix-matrix product (line 2) as well as in the CSR traversal (lines 3-17) are reported in the second and third column respectively, with increasing number of directions N d , reported in the first column. The number of FLOP performed in our product is reported in the fourth column.

The numbers in parentheses are the actual throughputs in GFLOP/s of the matrix-vector product. In the fifth

1: function y = MatMultAdd(x, Φ, I, J , A ⋆ , I S , J S , S ⋆ , s) 2: gemm("N", "N", N d , n, N d , 1, Φ, N d , x, N d , 0, Y , N d ) 3: do i = 0, n ⊲ row indices 4: do j = I(i), I(i + 1)
⊲ column indices of the volume reference matrix 5:

do k = 0, N d ⊲ direction indices 6: b 1 = (A 2a (j) + s k,1 A 3a (j) + s k,2 A 3d (j) + s k,3 A 3e (j)) s k,1 7: b 2 = (A 2b (j) + s k,2 A 3b (j) + s k,3 A 3f (j)) s k,2 8: b 3 = (A 2c (j) + s k,3 A 3c (j)) s k,3 9
:

y(i • N d + k) = y(i • N d + k) + A 1 (j) + b 1 + b 2 + b 3 x(J (j) • N d + k) + s k,1 A 4b (k) + s k,2 A 4c (k) + s k,3 A 4d (k) + A 4a (k) Y (J (j) • N d + k)
y(i • N d + k) = y(i • N d + k) + S k (j)x(J S (j) • N d + k) 15:
end do 16:

end do 17: end do 18: end function ⊲ return y = y + Ax column, knowing that the number of nonzero entries m in each reference matrix is 6.29 • 10 5 , the number of double-precision floating-point operations that would have required a matrix-ready matrix-vector product is reported. For small number of directions (N d 320), the performance of our algorithm is constrained by the single traversal of the sparsity pattern of the reference matrices, lines 3-17. With the CSR format, this is known to be an inefficient operation for discretizations on unstructured meshes [START_REF] Goumas | Performance evaluation of the sparse matrix-vector multiplication on modern architectures[END_REF]. When the number of directions becomes greater, the performance is driven by the one of the matrix-matrix product, line 2. Moreover, the vectorization of the inner loops lines 5 and 13 is most efficient for large values of N d . In general, one may notice that our approach is one order of magnitude more efficient, needing fewer operations, for a wide range of directions. Note that for the case with N d = 5,120, if A were to be assembled explicitly using double-precision floating-point numbers, it would require 6.29 • 10 5 • 5,120 2 • 8 ≈ 132 TB of memory (not accounting for integers used to represent the sparsity pattern of A). With the semi-matrix-free approach, this operator is assembled on an ordinary desktop computer. Though the Intel MKL (with AVX-512) is used to perform the matrix-matrix product, we are focused here once again on the single-thread performance of our algorithm on an Intel Xeon Platinum 8168.

Semi-matrix-free preconditioning

Krylov methods usually require a preconditioner so that they converge in a timely manner to an accurate solution. Because the coefficient matrix from eq. ( 7) is not assembled explicitly, standard preconditioners such as the restricted additive Schwarz method [START_REF] Cai | A restricted additive Schwarz preconditioner for general sparse linear systems[END_REF], or algebraic multigrid methods [START_REF] Stüben | A review of algebraic multigrid[END_REF], are not usable. Matrix-free solvers typically rely on either: a) basic preconditioners such as the Jacobi method, or b) more sophisticated tools like geometric multigrid methods [START_REF] May | A scalable, matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow[END_REF] or element-wise preconditioning [START_REF] Axelsson | Preconditioning of boundary value problems using elementwise Schur complements[END_REF]. In the case of the RTE, the former certainly makes a Krylov method such as the GMRES [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] converge, but with an insufficient convergence rate. When using different high-level languages, the latter is not easily interoperable since it requires a tight link with the language data structures, e.g., hierarchical meshes. In the next paragraph, an auxiliary operator will be explicitly assembled so that it can be used to precondition eq. ( 7) using standard techniques.

The diffusion synthetic acceleration (DSA) is a common tool in the radiative transfer community [START_REF] Larsen | Diffusion-synthetic acceleration methods for discrete-ordinates problems[END_REF][START_REF] Haut | Diffusion synthetic acceleration preconditioning for discontinuous galerkin discretizations of s n transport on high-order curved meshes[END_REF] to improve the convergence of iterative schemes when dealing with scattering media with large optical depths. It is based on a simplification of eq. ( 1) using the diffusion approximation. To build an efficient preconditioner, we use another approximation by assuming that the phase function matrix φ in eq. ( 3) is now diagonal:

∀k ∈ 1; N d , ( s k • ∇ + β( x))I k ( x) = σ( x)ω k I k ( x)φ k,k .
The following physical phenomena are encompassed in these pruned PDEs: transport 11:

s k •∇I k ( x)
I B (N d • i + k) = j shifted + k • m row ⊲ number
J B (idx) = J (j) • N d + k ⊲ column index of the pruned preconditioner 12: b 1 = (A 2a (j) + s k,1 A 3a (j) + s k,2 A 3d (j) + s k,3 A 3e (j)) s k,1 13
:

b 2 = (A 2b (j) + s k,2 A 3b (j) + s k,3 A 3f (j)) s k,2
14:

b 3 = (A 2c (j) + s k,3 A 3c (j)) s k,3
15:

B(idx) = B(idx) + A 1 (j) + b 1 + b 2 + b 3 + s k,1 A 4b (k) + s k,2 A 4c (k) + s k,3 A 3d (k) + A 4a (k) Φ k,k 16:
end do 17:

end do 18: end do components of I h in eq. ( 6) are only coupled through Φ. The assembly cost is thus divided by N d . Zero offdiagonal coefficients from the dense blocks of B are obviously not stored. Thus, there is no block structure in B and it is assumed that it is stored in a standard sparse matrix using the CSR format. Its arrays are noted I B , J B , and B. In fig. 4, the procedure used to allocate and assemble the pruned preconditioner is displayed.

19: I B (N d • n) = N d • m ⊲ total
For the same test case as in table 1 and on the same hardware, the time spent in this routine is shown in the second column of table 2, averaged over two runs. In the third (resp. fourth) column, the order of the system (resp. number of nonzero entries) is displayed. Choosing an efficient preconditioner to approximate the action of B -1 is not trivial. In our experiments, we used hypre algebraic multigrid solver BoomerAMG [START_REF] Henson | BoomerAMG: A parallel algebraic multigrid solver and preconditioner[END_REF] with the following parameters: ext+i interpolation, HMIS coarsening, two sweeps of l 1 Gauss-Seidel smoothing per level, and aggressive coarsening on all levels. Similar parameters are used in the neutron transport community [START_REF] Hanuš | A study of various thermal upscattering acceleration schemes for massively parallel transport sweeps[END_REF].

In the fifth column of table 2, the time spent in BoomerAMG for constructing the preconditioner is reported.

Obviously, with a fixed spatial mesh but a refined angular mesh, a steady increase in setup time is observed.

In the following section, a strategy to mitigate this issue will be presented, based on a hybrid multigrid in space-domain decomposition in angles preconditioner.

Large-scale numerical results

Both the matrix-vector product and the preconditioner have been defined in the previous paragraphs. In this section, the performance of the complete solver will be assessed for large-scale simulations. One application of the solver is its use to improve the design of heterogeneous participating media for material science. Thus, for a given geometry, multiple problems have to be solved, for example because of changing physical properties.

In [START_REF] Badri | Efficient finite element strategies for solving the radiative transfer equation[END_REF], this solver is used to study the response of open-cell foams made of SiC or Al 2 O 3 to radiations, which represent two different classes of materials (opaque or semi-transparent). We will thus focus on the strong scaling of the solver. However, since the transport community is heavily focused on weak scaling, there is also such an analysis in section 6.4. Before presenting any scaling tests of the solver, it is first validated in the next paragraph.

Validation

Problems 3i and 3ii from Kobayashi benchmarks [START_REF] Kobayashi | 3D radiation transport benchmark problems and results for simple geometries with void region[END_REF] are used to validate the solver. These test cases correspond to monochromatic emitting source radiation problems within heterogeneous materials with purely reflecting walls. They have already been used to validate radiation transport solvers in the literature [START_REF] Brown | Spherical harmonic solutions to the 3D Kobayashi benchmark suite[END_REF][START_REF] Zmijarevic | Deterministic solutions for 3D Kobayashi benchmarks[END_REF][START_REF] Azmy | TORT solutions for the 3D radiation transport benchmarks for simple geometries with void region[END_REF]. The material geometry for both problems 3i and 3ii is a 60 × 100 × 60 cm 3 rectangular cuboid which consists of three regions: source, void, and shield. They are respectively highlighted in blue, red, and transparent in fig. 5a. The radiative exchange in both of these problems is driven by the source region, a 10 × 10 × 10 cm 3 cube located at the origin of the coordinate system. It emits a black-body radiation of unit strength. The xz-plane for y = 0, xy-plane for z = 0, and yz-plane for x = 0, i.e., front, bottom, and left planes, are specularly reflecting, while the three other planes are vacuum walls. Problem 3i involves absorption, emission, and reflection but neglects scattering. The coefficients (κ, σ) are constant per region: (0.1 cm -1 , 0 cm -1 ) in the source, (10 -4 cm -1 , 0 cm -1 ) in the void, and (0.1 cm -1 , 0 cm -1 ) in the shield region. Problem 3ii involves all modes of radiation: absorption, emission, scattering, and reflection. The coefficients are now: (0.05 cm -1 , 0.05 cm -1 ),

(5 • 10 -5 cm -1 , 5 • 10 -5 cm -1 ), and (0.05 cm -1 , 0.05 cm -1 ), respectively.

Problems 3i and 3ii are prone to ray effects and not trivial to solve with the DOM [START_REF] Zmijarevic | Deterministic solutions for 3D Kobayashi benchmarks[END_REF]. As such, these benchmarks require a fine spatial and angular meshes, cf. fig. 1. This mitigates errors due to ray effects. Two sets of meshes were used with our solver to further quantify ray effects and validate our results. A coarse set made of a spatial mesh with 1.41 • 10 5 nodes and 8.15 • 10 5 elements and an angular mesh with 320 directions.

And a fine set made of a spatial mesh with 7.88 • 10 5 nodes and 4.71 • 10 6 elements and an angular mesh with 5,120 directions. As suggested in Kobayashi benchmark document [START_REF] Kobayashi | 3D radiation transport benchmark problems and results for simple geometries with void region[END_REF], densities are extracted at certain point locations within the materials. Table 3 provides these specific densities obtained with four different approaches:

1. numerical evaluation of an exact expression only available for problem 3i [START_REF] Kobayashi | 3D radiation transport benchmark problems and results for simple geometries with void region[END_REF], column "reference", 2. Monte Carlo solver [START_REF] Nagaya | MVP/GMVP 2: general purpose Monte Carlo codes for neutron and photon transport calculations based on continuous energy and multigroup methods[END_REF], column "MC", 3. spherical harmonics P 21 solver [77], column "P 21 ", 4. our solver with both sets of meshes, column "DOM c " and "DOM f " for the coarse and the fine set, respectively.

The data provided in the table suggests that the density computed by our solver when using the fine set of meshes are in close agreement with the ones provided by other approaches. Moreover for both problems 3i and 3ii, the accuracy of the results improves significantly when switching from the coarse to the fine set of meshes.

Contours of the density in logscale are presented in fig. 5b (resp. fig. 5c) for problem 3i (resp. 3ii).

Parallelism

Workload distribution. The finite element method is most often parallelized on distributed-memory systems using a coarse-grained approach such as domain decomposition. One of the many linear algebra backends that does support this workload distribution is PETSc, that will be used in all the following tests. Considering that an input mesh is split into P domains using a mesh partitioner such as Metis [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF], FreeFEM is in charge of generating a global unknown numbering of V h . This numbering is such that all P processes will own a contiguous chunk of rows of all reference matrices A ⋆ (one-dimensional row-wise distribution).

Table 3: density (W cm -2 ) of Kobayashi problems 3i and 3ii computed with: an exact expression [START_REF] Kobayashi | 3D radiation transport benchmark problems and results for simple geometries with void region[END_REF], GMVP [START_REF] Nagaya | MVP/GMVP 2: general purpose Monte Carlo codes for neutron and photon transport calculations based on continuous energy and multigroup methods[END_REF], Ardra [START_REF] Brown | Spherical harmonic solutions to the 3D Kobayashi benchmark suite[END_REF] or our solver, for x = z = 5 cm. Matrix-vector product. As it is the case with standard matrix types from PETSc (and other backends like hypre [START_REF] Falgout | hypre: a library of high performance preconditioners[END_REF]), it is furthermore assumed that the local representation of each matrix A ⋆ is split into two matrices

Problem 3i

{D i } P i=1 and {O i } P i=1
, where the unknowns associated with the {D i } P i=1 are purely local to each domain, whereas those associated with the {O i } P i=1 involve interprocess communication for computing matrix-vector products, cf. [START_REF] Bienz | Analyzing the performance of a sparse matrix vector multiply for extreme scale computers[END_REF]. By using a matrix of type "shell", one can wrap the matrix-vector product from fig. 3 so that PETSc will use it in its Krylov solvers. If the reference matrices are split using the previous decomposition, it is trivial to extend the algorithm so that it overlaps communication due to interprocess unknowns shared with neighboring domains with local computation.

Pruned preconditioner. In the previous paragraphs, matrix-free system eq. ( 7) is solved by using an assembled preconditioner B defined in section 4. 2, setting up BoomerAMG on the complete pruned operator B becomes costly for large numbers of directions. To decrease the setup cost, we propose to modify the preconditioner by first using a domain decomposition-like approach in angles, and then using multigrid in space. Instead of using a single instance of BoomerAMG with the single operator B, each structured diagonal block of dimension N d will be decomposed into q subblocks, where q is a common divisor of the number of directions N d and the number of processing elements P . This may be interpreted as nonoverlapping domain decomposition according to the angles. Then, q instances of BoomerAMG are created for each subsystem of dimension N d q • n. This is standard multigrid in space. Suitable restriction and interpolation operators will be defined later on to go from the full space to a single subsystem and vice versa. This is somehow similar to the anglesets strategy from [START_REF] Adams | Provably optimal parallel transport sweeps on regular grids[END_REF]. From an algebraic point of view, the modified pruned preconditioner is now written as:

B -1 = q k=1 R T k (R k BR T k ) -1 R k ,
where q is a common divisor of the number of directions N d and the number of processing elements P , and R k is the restriction on the kth angleset. Formally, for each spatial unknown, R k selects from all N d directions only the M q = N d /q ones from the appropriate angleset. In matrix form, R k is a full-rank diagonal Boolean matrix written as:

∀k ∈ 1; q , R k =    J k . . . J k    ∈ {0, 1} n•Mq×n•N d , with J k = ✵ Mq×(k-1)•Mq ✶ Mq ✵ Mq×(q-k)•Mq ∈ {0, 1} Mq×N d .
Within the PETSc framework, this construc- tion of B -1 could be performed using an additive fieldsplit preconditioner (with as many fields as anglesets), but this would have two main bottlenecks:

❼ the operator B would have to be assembled explicitly first, ❼ q would need to be low for the preconditioner assembly to be efficient.

Instead, as it is standard when doing domain decomposition in space, the respective R k BR T k q k=1 are directly assembled on q different MPI communicators of size P/q. This can be done by looping on the appropriate directions from a given angleset line 7 from fig. 4, instead of looping on all N d directions. For better memory efficiency, we use PETSc matrices of type "hypre" to store these operators, so that there is no copy or conversion from the default PETSc type "aij" to "hypre". Once the problem has been decomposed in angles, q instances of BoomerAMG are called concurrently to approximate the action of each (R k BR T k ) -1 q k=1 (multigrid in space), using the same parameters as in section 3.3.

Test case and scaling

Hardware and software settings. Results were obtained on a partition of Irène composed of 1,656 nodes with two 24-core Intel Skylake clocked at 2.7 GHz. The interconnect is an InfiniBand EDR full fat tree and the MPI implementation exploited was IntelMPI version 2019.0.5.281. All binaries and shared libraries were compiled with Intel compilers and Math Kernel Library support (for dense linear algebra computations) version 19.0.5.281.

PETSc version 3.12.4 with hypre version 2.18.2 were used.

Geometry and physics. For the spatial domain, a three-dimensional centimetric rectangular cuboid of dimension

[0; 1] × [0; 1] × [0; 3] is used.
In the rectangular cuboid lies a collection of multiple small spheres with varying physical parameters, cf. fig. 6. The resulting mesh generated by Gmsh [START_REF] Geuzaine | Gmsh: a 3-D finite element mesh generator with built-in pre-and postprocessing facilities[END_REF] is made of 1.9 million tetrahedra, which are uniformly refined once, for a total of 8 × 1.9 = 15 million tetrahedra. Before the refinement step, the unstructured mesh is partitioned in parallel using ParMETIS [START_REF] Karypis | Multilevel k-way partitioning scheme for irregular graphs[END_REF]. The background coefficients used in the rectangular cuboid are: κ = 10 -6 cm -1 and σ = 10 -6 cm -1 , i.e., a nearly transparent medium. In the spheres, the values of this pair of coefficients vary between: {10 -6 , 10}, {10, 10 -6 }, {0.1, 0.1}, {10, 10} . These values represent a mix of highly scattering, highly absorbing, and semi-transparent materials. They are considered homogenized with φ modeled using the Henyey-Greenstein phase function [START_REF] Henyey | Diffuse radiation in the galaxy[END_REF], with the anisotropy factor g = 0.5. As inflow boundary conditions, the rectangular cuboid is impinged with a collimated external radiative source I inflow on its left face:

I inflow = I 0 ✶ ( x, s)∈{0}×[0.25;0.75]×[0.25;0.75]×{[1,0,0] T } ,
where I 0 = 1 W cm -2 sr -1 . ✶ represents the vectorial characteristic function of a set. Since emission is neglected and the condition on angles reads s ∈ {[1, 0, 0] T }, the right-hand side B of eq. ( 7) is nonzero for a single component I h associated with this direction. In total, N d = 5,120 directions are used to discretize the unit sphere, cf. fig. 1b. This is an extremely fine angular mesh. Such a resolution will be necessary in a following test when dealing with specular reflections on spheres in order to solve the physics accurately [START_REF] Hardy | 3D numerical modelling of the propagation of radiative intensity through a X-ray tomographied ligament[END_REF]. These directions are grouped in q = 32 anglesets of dimension N d /q = 160. To the best of our knowledge, radiative transfer equation solvers using Krylov methods have never been able to deal with such an angular resolution, mostly due to memory limitations. After a simple post-processing step, it is easy to compute the spatial density profile D, i.e., at the continuous level, the integral over the complete solid angle of the radiative intensity, by computing for each spatial degree of freedom: ∀j ∈ 1; n , D j ≈ N d i=1 ω i x i,j , assuming x is the solution vector from eq. ( 7).

For this test, the density is shown in fig. 7. Notice the jumps in the density due to the presence of spheres with nonzero extinction coefficients along the collimated radiation.

# of MPI processes Time to solution (s) Performance. In the following numerical experiments and as standard in the radiative transfer community [START_REF] Godoy | On the use of flux limiters in the discrete ordinates method for 3D radiation calculations in absorbing and scattering media[END_REF][START_REF] Grella | Sparse tensor spherical harmonics approximation in radiative transfer[END_REF][START_REF] Krishnamoorthy | Parallel computations of radiative heat transfer using the discrete ordinates method[END_REF], convergence is assumed to be reached when the relative residual norm is lower than 10 -6 . Such a stopping criterion provided accurate final residuals 1 . Since the global operator is ill-conditioned, the long-recurrence Krylov method GMRES(30) 2 using double-precision floating-point numbers is preferred over short-recurrence methods such as BiCGSTAB, cf. [START_REF] Kanschat | Solution of radiative transfer problems with finite elements[END_REF] for a comparison with radiative transfer operators.

The performance of the complete solver is assessed in the strong scaling regime on 1,536 up to 24,576 MPI processes in fig. 8. The following steps of our method are timed, and averaged over two consecutive runs (same node distribution): assembly of the reference matrices A ⋆ and of the surface integrals (section 3.3), of the preconditioner B (section 4.2), and preconditioner setup (all three steps are summed and referenced to as the "setup" phase), and solution phase. In fig. 8a, each color represents the fraction of the total time spent in these two phases. For example, with 1,536 MPI processes, the height of the red bar is 8b). These two colors are used to differentiate at a glance both phases, in order to estimate the setup time to solve time ratio. This may be useful to infer how to amortize setup costs when doing successive linear solves. For the exact timings, interested readers are referred to the right-hand side table. Overall, even with a rather high number of anglesets and aggressive multigrid parameters, the preconditioner remains very robust in terms of iterations on a wide range of MPI process counts. The test case was chosen so that it saturates the node memory when running on 1,536 MPI processes.

For moderate number of MPI processes, computations (and thus cache effects) prevail over communications, which explains ideal speedups. At scale, the total time starts to be dominated by the solution phase. Indeed, on 24,576 processes, the average number of tetrahedra per process is a mere 610, so that the communication cost in the matrix-vector product cannot be fully overlapped and is no more negligible compared to the highly efficient computation, cf. table 1. For this largest configuration, the average grid and operator complexities of (R k BR T k ) -1 32 k=1 , as returned by hypre, are 1.09 and 1.17 respectively. These low numbers imply that most of the communications take place in the matrix-vector products. Note that their volume of communication is proportional to N d , while their computational cost is proportional to N 2 d , cf. eq. ( 13).

Specular reflections

In the previous sections, a solver for the radiative transfer equation was derived using a mixed matrixready/matrix-free approach. The goal of the following paragraphs is to show how this solver may be extended to deal with specular reflections. Such phenomena appear when radiations reach surfaces with different reflective indices. The radiative intensity is thus modified since a portion of the incident radiations is reflected back into the domain, according to Fresnel equations.

Assembly

To take into account these reflections, the bilinear form from eq. ( 6) has to be modified. Indeed, it is now assumed that the inflow boundary condition for all {I h k } N d k=1 is given, for all ( x, s k ) ∈ Γ -, by:

I h k = I k,inflow + j∈ 1;N d sj • n>0 ρ n,k,j I h j . ( 14
)
1 checked with the PETSc option -ksp view final residual 2 the default restart value of PETSc

13

Describing precisely how the coefficients ρ are computed is beyond the scope of this paper. There are various ways to do so, for example using the 1-DP (used in this paper), 3-DP, or the partition methods [START_REF] Hardy | Specular reflection treatment for the 3D radiative transfer equation solved with the discrete ordinates method[END_REF]. Their values depend on the reflective index of the surface material, the direction of the normal n at the point of reflection, as well as the incident and reflected directions s k and s j . It is assumed that the reflective surface Γ - may be decomposed into a union of surfaces with p different normals, i.e., Γ -= ∪ p r=1 Γ - r . For example, if Ω h is a hexahedron, p is equal to six. With the new inflow boundary condition, for all components of I h , the bilinear form a of eq. ( 6) has to be augmented by the following surface integral, for each component:

a(I h , J h ) += p r=1 Γ - r (P r ⋄ I h ( x)) ⊗ J h ( x)d x, (15)
with the matrices {P r } p r=1 being defined as:

∀r ∈ 1; p , P r =    ρ n,1,1 • • • ρ n,1,N d . . . . . . . . . ρ n,N d ,1 • • • ρ n,N d ,N d    .
For all ( x, s k ) ∈ Γ -, s k • n < 0, cf. the definition of Γ -given eq. ( 2). To take into account the conditional sum from eq. ( 14), s j • n > 0, the diagonals of {P r } p r=1 are set to zero. Moreover, since a large number of directions are used, these matrices are relatively sparse. In practice, they are stored using the coordinate list format. As in section 3.3, instead of assembling the vectorial bilinear form (15), the following variational formulations are assembled:

∀r ∈ 1; p , q r (I h , J h ) = Γ - r I h ( x)J h ( x)d x.
Let n Γ -(resp. {n Γ - r } p r=1 ) be the number of unknowns from V h associated with elements of Γ -(resp. {Γ - r } p r=1 ).

The assembly cost of eq. ( 15) is proportional to N d • n Γ -, while the total assembly cost of all reference matrices from the previous equation is proportional to p r=1 n Γ - r .

Matrix-vector products

The initial algorithm from fig. 3 has to be extended to take into account the reference matrices from the previous paragraph. It is assumed that the bilinear forms {q r } p r=1 are assembled using the CSR format {I qr , J qr , Q r } p r=1 . This time, the sparsity pattern of the reference matrices corresponding to the surface integral may be different. For a surface r ∈ 1; p , and for all unknown indices i ∈ 1; n Γ - r , we assume that n Γ - r (i) returns the unknown index in the finite element space V h . The extended matrix-vector algorithm is presented in fig. 9. For each reflective surface, it involves first gathering the relevant values from the input vector of size N d • n, then computing a sparse matrix-dense matrix product, and eventually, scattering back the result into the output vector. When dealing with meshes with a large amount of different reflective surfaces (large p), this algorithm may not scale well. However, it is important to keep in mind that for realistic three-dimensional geometries, the number of unknowns associated with surface elements n Γ -is expected to be much lower than the global number of unknowns n. Thus, the most time-consuming part of the matrix-vector procedure should still be the algorithm from fig. 3.

Because of the definition of {P r } p r=1 , eq. ( 15) would not introduce any diagonal values in the dense N d × N d blocks of A from eq. ( 7) if it were to be assembled explicitly. Thus, the definition of the preconditioning matrix B from section 4.2 will remain the same for the numerical experiments.

Test case and scaling

The same geometry with the same background coefficients from section 5.3 will be used. The main difference now is that all the spheres inside the rectangular cuboid model reflecting materials with no transmissivity. As such, their insides are not meshed anymore, and their reflective indices are picked from the following set: {1.8, 2.0, 2.2, 2.5}. The number of directions is still set to N d = 5,120, and with a mesh of 14.5 million tetrahedra, the discrete system is of order 14 × 10 9 . The outer surface of the rectangular cuboid is nonreflective (reflective index set to one). The same preconditioner as before is used: hybrid multigrid in space-domain decomposition in angles with the same BoomerAMG parameters and q = 32 anglesets. In the density now displayed in fig. 10, one can see that at the surface of the spheres in the rectangular cuboid, part of the radiative intensity is reflected back into the spatial domain. The impact of reflections that cause more variations in the density profile compared to fig. 7 may also be observed. Notice that for this test case, the total number of reflective surfaces p is equal to 407. Depending on how the global domain is decomposed, each process has a different number of local reflective surfaces {p i } P i=1 . Accordingly, they only load the corresponding matrices {P r } pi r=1 that are generated in a pre-processing step. As these numbers decrease, so does the local cost of the matrix-vector product (outer loop line 2 from fig. 9). For this configuration, the performance of the solver, cf. fig. 11, delivers similar speedups as in the previous section (without reflective surface).

do i = 0, n Γ - r 6: do k = 0, N d 7: Y (N d • i + k) = x(N d • n Γ - r (i) + k) 8:
end do do j = I qr (i), I qr (i + 1)

13:

do k = 0, N d 14: idx = N d • n Γ - r (i) + k 15: y(idx) = y(idx) + tmp(N d • i + k)Q r (j) 16:
end do (b) Timings of the solver. The last column provides the number of anglesets, equal to the number of hypre instances working concurrently.

Figure 12: weak scaling analysis of the complete solver for a system of 1.7 million unknowns per MPI process (with six reflective surfaces).

Spatial mesh refinement and weak scaling

Results were obtained on a partition of Irène composed of 2,292 nodes with two 64-core AMD Rome clocked at 2.6 GHz. The interconnect is an InfiniBand HDR full fat tree. A simple three-dimensional unstructured cube being impinged with a collimated radiative beam is now considered. It is again generated by Gmsh, with approximately 750k tetrahedra. The reflective index on the six faces of the cube is 1.5. The number of directions is now set to N d = 1,280. This saturates the memory available on one node of Irène, so it will be used next to perform a weak scaling analysis. Inside the cube, the following coefficients are used: (κ, σ) = (0.085, 20.10), as suggested in the literature [START_REF] Yang | Radiative Properties of Ceramic Al 2 O 3 , AlN and Si 3 N 4 -II: Modeling[END_REF]. After using ParMETIS to decompose this initial mesh in 128 × i 3 for i ∈ 1, 6 , and assigning one subdomain per MPI rank, uniform spatial mesh refinement is performed concurrently on all processes. Each tetrahedron is split into i 3 elements, by splitting all edges of the initial elements in i.

No refinement is performed for i = 1. The smallest problem generated for i = 1 has a total of 216 million unknowns, while the largest problem, running on 216 compute nodes for i = 6, has 45.7 billion unknowns.

Conclusion

The radiative transfer equation is a high-dimensional integro-differential equation. In this paper, a discretization in both space and angles was used to solve it numerically. Our solver is able to discretize the RTE efficiently by using a high-level finite element language, FreeFEM. By using such a language, most of the burden inherent of the finite element method, e.g., mesh data structure, is hidden from the user. Unfortunately, these languages do not exploit the structure of the RTE and cannot be used out of the box to directly discretize it.

The mixed matrix-ready/matrix-free approach introduced in this paper has the two following important properties. First, it relies on a low number of reference matrices of a much lower dimension than the original system. These may be assembled fast using almost any domain specific language at a relatively low memory and computational cost. Second, the matrix-vector product uses a dense matrix-matrix product, which may be highly optimized, and it is possible to assemble explicitly a preconditioner to reduce the number of iterations needed by Krylov methods to converge.

By using a hybrid multigrid in space-domain decomposition in angles preconditioner, we are able to discretize and solve a 45.7 billion unknown system in less than three minutes on 27,648 MPI processes. Moreover, we show how our solver may be extended to deal with physics involving specular reflections on both exterior or interior surfaces with no orthogonality condition. Thanks to the flexibility of modern domain specific languages and of PETSc, it would be easy to couple our solver with other fields, e.g., through a "fieldsplit" preconditioner 3 .

Overall, we exhibit large-scale deterministic results on fine unstructured and arbitrarily-decomposed spatial meshes and large numbers of directions, with a mix of heterogeneous scattering, absorbing, and reflecting materials.

  Angular mesh (N d = 320).

Figure 1 :

 1 Figure 1: DOM and FEM coupling strategy, each spatial unknown is associated with a set of N d angular unknowns.

Figure 2 :

 2 Figure 2: time spent assembling A from eq. (7) for a two-dimensional domain with P 1 finite elements either explicitly or using our semi-matrix-free approach.

  do j = I S (i), I S (i + 1)⊲ column indices of the surface reference matrix 13:do k = 0, N d ⊲ direction indices 14:

Figure 3 :

 3 Figure3: matrix-vector product using the semi-matrix-free discretization.

7 :

 7 , absorption κ( x)I k ( x), out-scattering σ( x)I k ( x), and unidirectional in-scattering σ( x)ω k I k ( x)φ k,k . This construction makes the preconditioner able to handle a wide range of radiative transfer problems with absorbing, scattering, transparent, or semi-transparent media, as displayed in the following numerical experiments. Assembling explicitly the pruned coefficient matrix B from the approximated RTE becomes much cheaper than assembling A, since the dense N d × N d blocks now become diagonal N d × N d blocks. Indeed, the various 1: function M = PrunedPC(x, Φ, I, J , A ⋆ , s) 2: allocate I B to store N d • n + 1 integers 3: allocate J B to store N d • m integers 4: allocate B to store N d • m scalars 5: do i = 0, n ⊲ row indices 6: m row = I(i + 1) -I(i); j shifted = N d • I(i) do k = 0, N d ⊲ direction indices 8:

of nonzeros of this row 9 :

 9 do j = I(i), I(i + 1)⊲ column indices of the volume reference matrix 10:idx = j shifted + k • m row + (j -I(i))

  Density profile of problem 3i. (c) Density profile of problem 3ii.

Figure 5 :

 5 Figure 5: Kobayashi benchmark.

2 .

 2 The strategy used to assemble B was to prune all off-diagonal nonzeros in the N d × N d dense blocks of A. The pruned operator B of dimension N d • n thus still remains structured, with blocks of dimension N d , but they are now diagonal instead of dense. As displayed in the previous section in table

Figure 6 :

 6 Figure 6: Gmsh geometry for the large-scale experiments.

2 D (W cm - 2 )Figure 7 :

 227 Figure 7: density profile D. The rectangular cuboid is impinged on the front face and the radiation travels up to the back face, while being partially absorbed in the spheres.

( a )

 a Time to solution. Each color represents the fraction of the total time spent in the respective tasks. In parenthesis, number of iterations.

Figure 8 :

 8 Figure8: strong scaling analysis of the complete solver for a system of thirteen billion unknowns (without reflection).

  .1 = 64% of the height of the complete bar (figures taken from table

Figure 9 : 5 D

 95 Figure 9: matrix-vector product for taking reflections into account.

Figure 10 :

 10 Figure 10: density profile D when dealing with reflective surfaces. The rectangular cuboid is impinged on the front face and the radiation travels up to the back face, while being reflected on the surface of the spheres.

Figure 11 :

 11 Figure11: strong scaling analysis of the complete solver for a system of fourteen billion unknowns (with reflections).

1

 1

  While the number of iterations to solve the global problem slightly increases, from 31 to 39, see third column in fig.12b, one can notice that the overall efficiency remains above 80%, see fig.12a. The same steps as in the previous strong scaling experiments are timed: assembly of the reference matrices A ⋆ and of the surface integrals (section 3.3), of the preconditioner B (section 4.2), preconditioner setup, and solution phase. To keep the preconditioner setup tractable, the number of anglesets q is increased with the global number of processes as shown in the rightmost column of fig.12b. On this range of process numbers, the global preconditioner remains numerically efficient, and it converges in fewer iterates than in our previous strong scaling experiments with spherical inclusions and internal reflections, see figs. 8b and 11b for a comparison.

Table 1 :

 1 single-thread performance of the matrix-vector product with an increasing number of directions, for a fixed spatial domain. n = 44,204, m = 6.29 • 10 5 , and m S = 89,483.

	N d	gemm (ms)	CSR (ms)	GFLOP MF (/s)	GFLOP MR
	80	26	115	2 (12)	8
	128	48	125	3 (17)	21
	320	184	262	10 (22)	129
	512	383	395	21 (26)	330
	1,280	1,960	1,000	95 (32)	2,060
	2,048	4,579	1,574	222 (36)	5,274
	5,120	27,281	4,116	1,250 (39)	32,964

Table 2 :

 2 number of nonzeros 20: end function ⊲ return B = CSR(I B , J B , B) time spent assembling the pruned preconditioner. n B and m B represent the order and the number of nonzero entries respectively in B.

	Figure 4: explicit assembly of the pruned preconditioner.
	N d	Assembly (s)	n B	m B	Setup (s)
	80	1.9	4 • 10 6	5 • 10 7	9.1
	128	2.9	6 • 10 6	8 • 10 7	15.4
	320	7.5	1 • 10 7	2 • 10 8	42.5
	512	12.1	2 • 10 7	3 • 10 8	72.7
	1,280	31.0	6 • 10 7	8 • 10 8	198.2
	2,048	53.0	9 • 10 7	1 • 10 9	317.5
	5,120	102.1	2 • 10 8	3 • 10 9	1,181.7

  {P r , I qr , J qr , Q r } p r=1 ) 2: do r = 1, p

	3: 4:	allocate Y to store N d • n Γ -r scalars allocate tmp to store N d • n Γ -r scalars
	5:

1: function y = MatMultAdd(x,

  a) Efficiency and number of unknowns with an increasing number of uniform mesh refinements. In parenthesis, number of iterations.

	Efficiency relative to 128 processes	0% 20% 40% 60% 80% 100%	2 8 (31)	1 , 0 2 4 (32)	3 , 4 5 6 8 , 1 9 2 1 6 , 0 0 0 2 7 , 6 4 8 (34) (34) (37) (39) 0.2 45.8	# of unknowns (billions)	P 128 1,024 3,456 8,192 16,000 27,648	Setup + solve (s) 133.0 140.4 145.5 143.8 149.6 157.1	# of iterations 31 32 34 34 37 39	q 2 8 16 32 64 128
				# of MPI processes

(

Acknowledgments

This work was granted access to the HPC resources of TGCC@CEA under the allocations A0050607519, A0070607519, and A0090607519 made by GENCI.


Sophie Sanguin

Horace Roman

Arthur Foulon

Jean Gondry

Eric Havet

Dr Céline Klein
email: celinekleinfr@yahoo.fr

Deep infiltrating endometriosis with sacral plexus involvement: improving knowledge through human cadaver dissection

Keywords: human cadaver dissection, deep infiltrating endometriosis, sacral roots, sciatic nerve

A human cadaver model is a suitable means of improving a physician's level knowledge about this rare, complex type of deep infiltrating endometriosis.

Introduction

Although endometriosis is a concern for around 10% of women of childbearing age, deep infiltrating endometriosis (DIE) with sacral plexus involvement is rare [START_REF] Possover | Laparoscopic Neurolysis of the Sacral Plexus and the Sciatic Nerve for Extensive Endometriosis of the Pelvic Wall[END_REF][START_REF] Nehme-Schuster | Alcock's canal syndrome revealing endometriosis[END_REF][START_REF] Possover | Isolated infiltrative endometriosis of the sciatic nerve: a report of three patients[END_REF]. The sacral plexus involvement results from the extension of endometriotic nodules deep into the retroperitoneal spaces, with entrapment of the sacral roots S2 to S4, the hypogastric nerve, the splanchnic nerves, and the inferior hypogastric plexus. In these specific locations, women can experience urinary symptoms (dysuria), bowel symptoms (constipation, dyschesia) and neurologic pain (sciatic nerve pain, pudendal nerve pain), for 51%, 45% and 71% of them [START_REF] Roman | Postoperative Outcomes after Surgery for Deep Endometriosis of the Sacral Plexus and Sciatic Nerve: A 52-patient Consecutive Series[END_REF].

Moreover, the surgical treatment of DIE with sacral plexus involvement can also be associated with a risk of nerve damage. Indeed, surgical injury to the sacral roots can lead to the sensory or motor impairments [START_REF] Díaz-Feijoo | Nerve-Sparing Technique during Laparoscopic Radical Hysterectomy: Critical Steps[END_REF][START_REF] Possover | Laparoscopic management of neural pelvic pain in women secondary to pelvic surgery[END_REF], and injury to the splanchnic nerves can result in the loss of urinary, defecatory or sexual function, as has been reported for radical surgery [START_REF] Laterza | Bladder function after radical hysterectomy for cervical cancer[END_REF][START_REF] Jensen | Early-stage cervical carcinoma, radical hysterectomy, and sexual function. A longitudinal study[END_REF][START_REF] Sood | Anorectal dysfunction after surgical treatment for cervical cancer[END_REF].

Hence the challenge of this surgery is to avoid worse or new symptoms, like for example selfurinary catheterization in the situation of a severe bladder atonia. In addition to these functional complications, dissection of the retroperitoneal spaces (in order to gain access to sacral roots) may result in vascular injury (especially for the internal iliac vein) and severe hemorrhage. Several studies have highlighted limitations in postgraduate surgeons' knowledge of surgical anatomy [START_REF] Ripperda | Anatomic relationships of the pelvic autonomic nervous system in female cadavers: clinical applications to pelvic surgery[END_REF][START_REF] Barton | Dissection of soft-preserved cadavers in the training of gynaecological oncologists: report of the first UK workshop[END_REF][START_REF] Gordinier | The effects of a course in cadaver dissection on resident knowledge of pelvic anatomy: an experimental study[END_REF]. Gynecological surgeons are often unfamiliar with the anatomy of the pelvic nerves because the sacral roots are far from the surgical field in operations for benign gynecological diseases.

In view of (i) the rarity of DIE with sacral root involvement and (ii) the surgeons' lack of knowledge, we decided to assess these topics among residents in gynecology before and after a lecture and a cadaver dissection class.

The primary objective of the present study was to determine whether pelvic dissection of human cadavers improves a surgeon's level of knowledge about the deep pelvic anatomy and about the functional relationship with nerve entrapment in DIE or nerve injury during surgery for DIE. The secondary objective was to describe the advantages and disadvantages of pelvic dissection in cadavers.

Materials and methods

We performed a prospective, observational study of participants in an anatomy course that included human cadaver dissection. The study participants were recruited from the pool of residents in the Department of Obstetrics and Gynecology at Amiens University Medical Center (Amiens, France) and gave their consent. The study procedure had been approved by the local investigational review board. Statistical analyses of continuous variables were performed using Student's t-test. The threshold for statistical significance was set to p<0.05.

For the dissection procedure, four female cadavers (1 fresh non-preserved cadaver, and 3 embalmed cadavers; all Caucasian) were dissected in a university anatomy laboratory between January and June 2019. One cadaver had undergone a hysterectomy -probably via a vaginal approach. None of the cadavers had laparotomic scars or pelvic endometriotic lesions.

The senior surgeon and the junior surgeon on each cadaver performed measurements separately.

Dissection technique:

We dissected the cadavers by applying an open surgical technique for the treatment of DIE with sacral root involvement; a laparoscopic variant of the technique has been described previously [START_REF] Zanatta | Laparoscopic Dissection and Anatomy of Sacral Nerve Roots and Pelvic Splanchnic Nerves[END_REF][START_REF] Possover | Laparoscopic therapy for endometriosis and vascular entrapment of sacral plexus[END_REF][START_REF] Ceccaroni | Laparoscopic nerve-sparing transperitoneal approach for endometriosis infiltrating the pelvic wall and somatic nerves: anatomical considerations and surgical technique[END_REF] . Laparoscopic cadaver dissection was not technically feasible on our laboratory, so we applied an open approach with a cruciate incision. For each cadaver, only one hemipelvis was dissected. The same dissection steps were performed in all cases. All dissections were performed by a senior gynecological surgeon with expertise in DIE surgery (SS), and a resident in gynecology. Firstly, we inspected the pelvis and identified anatomic structures: the sacral promontory, the rectum, the uterus (which was then sutured attached to the anterior wall), and the uterosacral ligaments (USLs). Secondly, we opened up the retroperitoneal spaces in front of the sacral promontory, in order to identify the aortic bifurcation, the superior hypogastric plexus, the external and internal iliac vessels, the ureter, and the deep uterine vein and artery. Thirdly, we opened up the pararectal spaces: (i) Okabayashi's pararectal space between the ureter and the posterior layer of the broad ligament, and (ii) Latzko's pararectal space between the internal iliac vein and the ureter.

Through these spaces, we identified the hypogastric nerve, the pelvic splanchnic nerves emerging from the sacral nerve roots S2 to S4, and the inferior hypogastric plexus. Lastly, the internal iliac vein was dissected medially from its origin, the parietal pelvic fascia was opened up, and the piriformis muscle and the sacral nerve root S1 were identified. Dissection was then performed caudally to successively reveal the sacral nerve roots S2 to S4, and the thin nerve fibers leading to the pelvic splanchnic nerves that had been identified in the previous step (Figures 1, 2 and3).

Any difficulties related to the dissection procedures or anatomical variants were fully documented. The duration of each dissection was noted. The following distances were measured by both the senior surgeon and the resident using a flexible ruler, and then averaged: the distance between the ureter and the hypogastric nerve, and the distances between the mid-portion of the USLs and the sacral nerve roots S2, S3 and S4, respectively.

Results

The test results are summarized in Table 1. Initial levels of knowledge about pelvic innervation and the impact of nerve damage in this region were low for all 8 students (mean (range) score in the pre-test: 2.5 (0-8)) but improved significantly after the lecture (mean score: 8.5, p=0.9.10 -5 for pre-vs. post-) and after the dissection class (mean score: 14.25, p=0.6.10 -6 for pre-vs. post-). The improvement was significantly greater for the dissection class (p=0.0003).

When considering the cadaver dissection, no technical difficulties were encountered in cases #1 and #3. In case #2, pelvic adhesions complicated the dissection. In case #4, a rectal fecaloma hindered correct exposure of the Douglas pouch and the USLs; the fecaloma was removed by transrectal incision and segmental resection.

We were able to dissect and identify the above-mentioned anatomic structures rapidly and unambiguously in three of the four cadavers. In case #3, a relevant anatomical variant was present: large inferior gluteal vessels were located in front of the sacral roots between S2 and S3 (Figure 4). These vessels were not observed in the other three cases.

The distances between the anatomic structures are given in Table 2. The median (range) dissection time was 145 min (130-180). There were no inter-observer differences in the distances measured. The median (range) distance between the ureter and the hypogastric nerve was 2.25 cm (2 -2.5). The median (range) distances between the mid-portion of the USLs and the sacral roots S2, S3 and S4 were respectively 5.5 cm (5.5 -6), 5.0 cm (4.5 -6), and 5.0 cm (4.5 -6.5).

Discussion

A good understanding of pelvic anatomy is essential for gynecologic surgery and especially for difficult gynecologic surgery, such as the treatment of DIE with pelvic nerve involvement.

A combination of risk factors in this kind of surgery can lead to irremediable functional consequences and (in the event of hemorrhage) severe complications. The first risk factor is the rarity of this pathology [START_REF] Possover | Laparoscopic Neurolysis of the Sacral Plexus and the Sciatic Nerve for Extensive Endometriosis of the Pelvic Wall[END_REF][START_REF] Nehme-Schuster | Alcock's canal syndrome revealing endometriosis[END_REF][START_REF] Possover | Isolated infiltrative endometriosis of the sciatic nerve: a report of three patients[END_REF]. The second is the gynecologist's lack of knowledge of deep pelvic anatomy [START_REF] Ripperda | Anatomic relationships of the pelvic autonomic nervous system in female cadavers: clinical applications to pelvic surgery[END_REF][START_REF] Gordinier | The effects of a course in cadaver dissection on resident knowledge of pelvic anatomy: an experimental study[END_REF]. Our study confirmed this lack of knowledge in a selected population of residents in obstetrics and gynecology. Although practical anatomy can be learned in the operation room, the small number of cases and time and safety constraints during surgery prevent residents from learning properly about these topics.

Our results demonstrate that the human cadaver is a suitable model for teaching deep pelvic anatomy (and especially pelvic innervation) and the functional consequences in the event of injury. Furthermore, use of a cadaver model enables the surgeon to become familiar with vascular entrapments, anatomical variants, and the relationship between the pelvic nerves and surrounding anatomic structures.

Our results are in line with the literature data and also highlight the objectives of the French reform of the 3rd cycle of the medical studies in which surgical anatomical training is emphasized. It is well known that didactic session using cadaveric model after gaining surgical experience improves the anatomy comprehension and surgical skills [START_REF] Heisler | Importance of adequate gross anatomy education: The impact of a structured pelvic anatomy course during gynecology fellowship[END_REF]. For example, Corton et al. reported that cadaver dissection can increase the gynecologist's proficiency in surgical anatomy [START_REF] Corton | A comprehensive pelvic dissection course improves obstetrics and gynecology resident proficiency in surgical anatomy[END_REF]. In the UK, Barton et al. [START_REF] Barton | Dissection of soft-preserved cadavers in the training of gynaecological oncologists: report of the first UK workshop[END_REF] developed a cadaver-based workshop and found that it improved the trainees' knowledge of pelvic anatomy. Improving the level of anatomical knowledge is one of the prerequisites for performing surgeryespecially for difficult situations like DIE with plexus entrapment.

Three studies described a standardized laparoscopic approach for the treatment of DIE with sacral plexus entrapment [START_REF] Zanatta | Laparoscopic Dissection and Anatomy of Sacral Nerve Roots and Pelvic Splanchnic Nerves[END_REF][START_REF] Possover | Laparoscopic therapy for endometriosis and vascular entrapment of sacral plexus[END_REF][START_REF] Ceccaroni | Laparoscopic nerve-sparing transperitoneal approach for endometriosis infiltrating the pelvic wall and somatic nerves: anatomical considerations and surgical technique[END_REF]. In both cases, the researchers recommended that these cases should be managed in tertiary referral centers by skilled surgeons. Increasing the surgeon's level of skill in rare, complex procedures such as the surgical treatment of DIE requires a companionship approach, in line with Halstead's "see one, do one, teach one" principle [START_REF] Kotsis | Application of the "see one, do one, teach one" concept in surgical training[END_REF].

Since the Renaissance, human cadavers have been used as models to teach anatomy to artists and physicians [START_REF] Denton | Sciatic syndrome due to endometriosis of sciatic nerve[END_REF], and they are still used to teach medical students or residents. The main advantage of a cadaver training model for the identification of pelvic nerves and sacral roots is the high degree of reproducibility. Repetition of the multistep dissection makes it easier for the surgeon to understand anatomic features, and might improve his/her surgical skills [START_REF] Ripperda | Anatomic relationships of the pelvic autonomic nervous system in female cadavers: clinical applications to pelvic surgery[END_REF][START_REF] Sharma | A cadaveric procedural anatomy course enhances operative competence[END_REF]. Hence, cadaver dissection might reduce the operating time and increase the surgeon's level of confidence [START_REF] Sharma | A cadaveric procedural anatomy course enhances operative competence[END_REF].

Our measurements of the main anatomical landmarks were similar for the four cadavers. Indeed, we found that the distance between the ureter and the hypogastric nerve and the distances between the sacral roots S2 to S4 and the ureter or the USLs were remarkably similar from one cadaver to another. Likewise, the anatomic locations of hypogastric nerve, splanchnic nerves and inferior hypogastric plexus did not differ from one cadaver to another.

Our present results show that these landmarks can be used during the surgery of deep pelvic spaces. Anatomical variants are nevertheless possible, and may particularly concern the deep pelvic vasculature. Overlooking these situations may lead to unexpected severe hemorrhages and thus conversion to laparotomy. Awareness of anatomical variants may prepare the surgeon for unusual intraoperative findings that can complicate the surgical procedure in general and the dissection in particular.

Cadaver models nevertheless have a number of limitations. Firstly, the lack of intraoperative hemorrhages from hypogastric vessels and their branches means that young surgeons cannot learn to manage this severe complication. Secondly, female cadavers rarely present DIE, and the uterus and the adnexa tend to be atrophied in elderly cadavers. Thirdly, the open and laparoscopic approaches differ with regard to the angle of vision (influencing the view and the magnification) [START_REF] Rabischong | Nerve sparing techniques in deep endometriosis surgery to prevent urinary or digestive functional disorders: Techniques and results: CNGOF-HAS Endometriosis Guidelines[END_REF]. Fourthly, the surgical environment in the dissection of fresh or embalmed cadavers is very different from that encountered during laparoscopy in living patients. Fifthly, the dissection of only four cadavers in the present study limits the validity of our anatomic measurements. Lastly, cadaver dissection is expensive, and requires access to an anatomy laboratory and a sufficient supply of cadavers.

In order to reduce the cost of cadaver dissection (one of the main limitations of these training models), several in-house models of specific surgical techniques (such as inguinal hernia repair) have been developed [START_REF] Aljamal | Cadaveric Dissection vs Homemade Model: What is the Best Way to Teach Endoscopic Totally Extraperitoneal Inguinal Hernia Repair[END_REF]. It might well be possible to develop a similar model for the dissection of sacral roots. Recently, a French group developed a fresh cadaver perfusion/ventilation technique, which may mimic the in vivo environment [START_REF] Delpech | SimLife a new model of simulation using a pulsated revascularized and reventilated cadaver for surgical education[END_REF]. As mentioned above, cost is not the only limitation of cadaver models; the variable availability of fresh cadavers, the time required to prepare the bodies, and the need for a cadaver-dedicated laparoscope are also problematic [START_REF] Delpech | SimLife a new model of simulation using a pulsated revascularized and reventilated cadaver for surgical education[END_REF].

Lastly, the most modern approach to surgical training (three-dimensional laparoscopic computer simulators) is promising -particularly for younger surgeons [START_REF] Diesen | Effectiveness of laparoscopic computer simulator versus usage of box trainer for endoscopic surgery training of novices[END_REF][START_REF] Leblanc | A comparison of human cadaver and augmented reality simulator models for straight laparoscopic colorectal skills acquisition training[END_REF]. However, simulators do not provide tactile feedback, and are still very expensive [START_REF] Leblanc | A comparison of human cadaver and augmented reality simulator models for straight laparoscopic colorectal skills acquisition training[END_REF].

The present study presented several limitations. It will be necessary to confirm our results in a prospective study with a larger sample size of cadavers and participants. A larger cohort would allow for more robust statistical values.

Despite its limitations, a human cadaver model is useful for improving the level of knowledge of deep pelvic innervation and the functional consequences of injury to pelvic nervesespecially in cases of DIE with sacral plexus entrapment. Use of this model may also help to prepare surgeons for complex procedures, such as the surgical treatment of DIE with pelvic nerve and sacral plexus involvement. The reproducibility of cadaver dissection, the facilitated identification of anatomic landmarks, and the relatively invariant distances between anatomic landmarks can safely improve surgical skills and procedures.

  Eight junior residents (with less than four years of residency) agreed to participate in the study. Each course consisted of a one-hour lecture and a cadaver dissection class led by an expert in DIE, with assistance from a junior surgeon. The first part of the course consisted in a one-hour lecture, with a pre-test and post-test. The pre-test and post-test were respectively performed before the lecture and 1 week after the lecture. The tests assessed the candidate's knowledge of the disease and the related anatomical features. In each test, the study participant could score up to 16 points. The second part of the course involved observing a human cadaver dissection focused on the pelvic nerves. The second post-test took place one week after the cadaver dissection class.

-Figure 1 . 4 Figure 2 . 3 Figure 3 .Figure 4 .

 142334 Figure 1. A fresh cadaver

Table 1 :

 1 pre-and post-test results

	Participant	Pre-test	Post-test (after the lecture)	Post-test (after the cadaver dissection)
	1	4	12	16
	2	4	10	14
	3	8	14	16
	4	0	6	12
	5	0	4	14
	6	0	8	14
	7	4	8	16
	8	0	6	12
	Mean	2.5	8.5	14.25

Table 2 .

 2 Distances between pelvic anatomic structures, and the dissection time. (HN: hypogastric nerve, USLs: uterosacral ligaments)

		Case #1 Case #2 Case #3 Case #4 Median
	Distance (cm)
	-Ureter to HN	2	2	2.5	2.5	2.25
	-Mid-point of USLs to S2 5.5	6	5.5	-	5.67
	-Mid-point of USLs to S3 5	6	4.5	-	5.17

ACKNOWLEDGMENTS

The authors sincerely thank those who donated their bodies to science so that anatomical research could be performed. Results from such research can potentially increase mankind's overall knowledge that can then improve patient care. Therefore, these donors and their families deserve our highest gratitude.

From which sacral roots (S) are the origins of the splanchnic nerves and the sacral sympathic nerves:

What are the consequences of inferior hypogastric plexus's injury or splanchnic nerves's injury?

• Motors defect =0