J Duriez

S Bonelli

Precision and computational costs of Level Set-Discrete Element Method (LS-DEM) with respect to DEM

Keywords: computational cost, particle shape, Level Set-Discrete Element Method (LS-DEM)

The Level Set-Discrete Element Method (LS-DEM) extends DEM towards arbitrary grain shapes by storing distance-to-surface values on a grid for each Discrete Element (DE), together with considering boundary nodes located onto the DE's surface. Both these ingredients are shown to affect the precision and computational costs of LS-DEM, considering various numerical simulations at the contact-and packing-scales for ideal spherical and superellipsoid shapes. In the case of a triaxial compression for spherical particles, approaching with a reasonable precision the reference result obtained in classical DEM requires the grid spacing to be smaller than one tenth of particle size, as well as using a couple thousands of boundary nodes. Computational costs in terms of memory (RAM) or evaluation time then increase in LS-DEM by two or three orders of magnitude. Simple OpenMP parallel simulations nevertheless significantly reduce the increase in time cost, possibly dividing the latter by 20.

Introduction

At the micro-scale considered by Discrete Element Methods (DEM), granular soils reveal diverse grain's shapes, that constitute one ingredient of their discrete nature. This shape enters soil classification and is directly used in geotechnical ˚Corresponding author

Email address: jerome.duriez@inrae.fr (J. Duriez) engineering for the ballast foundations of railtracks, which rely over angular, not spherical, particles. Outside of this pratical example, particle shape has been recognised as influencing the mechanical behavior of granular materials since several studies often adopting DEM approaches. In an early 2D study on rotating cylinders and heap configurations [START_REF] Pöschel | Static friction phenomena in granular materials: Coulomb law versus particle geometry[END_REF], a non-spherical shape was shown to contribute even more to macro-behavior than contact friction in the sense non spherical particles in frictionless interaction revealed a higher slope stability than spherical particles in frictional interaction.

For a given frictional interaction, a higher shear strength of non-spherical particles has also been found for biaxial configurations in other 2D studies [START_REF] Szarf | Influence of the grains shape on the mechanical behavior of granular materials[END_REF][START_REF] Jerves | Effects of grain morphology on critical state: a computational analysis[END_REF], together with a shape influence onto the critical state line [START_REF] Jerves | Effects of grain morphology on critical state: a computational analysis[END_REF].

Investigating the mechanical influence of shape in real 3D conditions remains however technically challenging. While experimental studies require a proper particle-scale characterization of the complex shapes exhibited in nature [START_REF] Vlahinić | Towards a more accurate characterization of granular media: extracting quantitative descriptors from tomographic images[END_REF][START_REF] Wang | Superellipsoid-based study on reproducing 3D particle geometry from 2D projections[END_REF], those same real shapes have to be correctly introduced in the numerical world for DEM approaches. This induces a much more complex contact treatment in the DEM workflow, as opposed to the use of spherical particles which entails straightforward definitions of contact normals and relative displacements from the branch vector and the radii of contacting spheres. These complex contact treatments may obey several strategies which are partially listed in the following. First, rigid clusters of spheres [START_REF] Pöschel | Static friction phenomena in granular materials: Coulomb law versus particle geometry[END_REF][START_REF] Szarf | Influence of the grains shape on the mechanical behavior of granular materials[END_REF][START_REF] Garcia | A clustered overlapping sphere algorithm to represent real particles in discrete element modelling[END_REF] enable the DEM practitioner to get much closer to real shapes, making these rigid clusters probably the second most-commonly used shape for Discrete Elements, just after spheres. These clusters nevertheless inherently include some unrealistic local roundness that may affect the mechanical description [START_REF] Cho | Particle shape effects on packing density, stiffness, and strength: Natural and crushed sands[END_REF].

Convex polyedra [START_REF] Eliáš | Simulation of railway ballast using crushable polyhedral particles[END_REF][START_REF] Gladkyy | DEM simulation of polyhedral particle cracking using a combined Mohr-Coulomb-Weibull failure criterion[END_REF] now constitute another quite classical shape enhancement since [START_REF] Cundall | Formulation of a three-dimensional distinct element model-Part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF], thanks to a variety of algorithms such as searching for surface points with a common normal and/or minimizing interparticle distance [START_REF] Dubois | Numerical modeling of granular media composed of polyhedral particles[END_REF]. As described by [START_REF] Zhao | A poly-superellipsoid-based approach on particle morphology for DEM modeling of granular media[END_REF], some of those algorithms can also be adapted to superellipsoids and quite general convex shapes without any edges. A last DEM variant to be mentioned is the Level Set-DEM (LS-DEM) proposed in 3D by [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF]. LS-DEM appears as promising in terms of versatility, since it does not include any inherent requirement for convexity and may apply directly to X-ray tomography images of soil samples [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF]. Level Set concepts were initially proposed to study time evolutions of surfaces [START_REF] Sethian | Level set methods and fast marching methods[END_REF], and applied in this sense to geotechnics by [START_REF] Golay | Interfacial erosion: A three-dimensional numerical model[END_REF][START_REF] Golay | Numerical modelling of interfacial soil erosion with viscous incompressible flows[END_REF] for flow-induced interfacial soil erosion. In the sense of LS-DEM, those Level Set concepts are used for defining in space distance fields to particles' surfaces, that are at the heart of contact treatment.

One can finally think about introducing more complex contact laws as an indirect description of particle's shape [START_REF] Wensrich | Rolling friction as a technique for modelling particle shape in DEM[END_REF]Aboul Hosn et al., 2017). However, this strategy obviously induces additional model parameters and increased calibration efforts that diminish the appealing mechanical simplicity of DEM.

Advocating therefore for a direct description of particle's shape through e.g.

LS-DEM, the present manuscript then aims to discuss associated technical aspects in terms of obtained precision and increased computational costs, in the case of an implementation based on the YADE code ( Šmilauer et al., 2015). Detailed information in these technical aspects seem lacking until now, even though one can await significant costs from the mentions of gigabytes RAM footprint in [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF] or superprocessors with 480 cores in [START_REF] Kawamoto | All you need is shape: Predicting shear banding in sand with LS-DEM[END_REF].

Section 2 presents the YADE implementation of LS-DEM based on the principles given by [START_REF] Jerves | Effects of grain morphology on critical state: a computational analysis[END_REF]; [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF]. Section 3 discusses the variable precision of LS-DEM in describing contact-or packing-scale configurations adopting spherical or superellipsoid shapes: ideal spherical shapes are in particular considered for the precision analysis to ground on reference results obtained using DEM. LS-DEM precision is then connected with computational costs in Section 4, before that parallel scalability is examined in Section 5 in order to alleviate time costs.

Outline of LS-DEM

Shape description

Describing shape, i.e. particle morphology, in LS-DEM relies on the signed distance function φp xq that returns, for any point x in space, the shortest distance from x to the surface at hand, with the convention of negative distances when x lies inside the surface. The surface of a Discrete Element (DE) then corresponds to the zero level set of the function φ, while the exterior (resp.

inner) to the surface obeys φ ą 0 (resp. φ ă 0).

In this sense, LS-DEM is similar to the potential particles approach proposed by [START_REF] Houlsby | Potential particles: a method for modelling non-circular particles in DEM[END_REF]; [START_REF] Boon | A new contact detection algorithm for three-dimensional non-spherical particles[END_REF] where the sign of a potential function f defines the position of any point with respect to particle's surface, with f " 0 along the surface. Potential particles however require convex shapes and polynomial equations for the potential f , unlike LS-DEM.

In LS-DEM, the signed distance function φ is actually defined in a discrete fashion, storing φ-values on a cartesian body-centered grid, for each DE (Figure 1). This minor requirement of a discrete distance field, instead of an analytical equation, confers LS-DEM a great versatility to mimic real shapes, as exemplified by [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF][START_REF] Kawamoto | All you need is shape: Predicting shear banding in sand with LS-DEM[END_REF]. From the knowledge of φ-values at each node of the grid, φp xq is also defined for any point x within the grid extents from trilinear interpolation of φ-values at the eight surrounding grid nodes. In addition to defining particle's surface, and serving for contact treatment as described in the following section 2.2, this distance fied also enables one to define inertial quantities for DE summing mass and inertia contributions of all grid voxels that are considered inside a particle.

Here, a grid voxel made of eight nodes tpi, j, kq ; i P ri 0 ; i 0 `1s, j P rj 0 ; j 0 `1s, k P rk 0 ; k 0 `1su is considered inside a particle depending on φ-value at the lowest node pi 0 , j 0 , k 0 q. A smoother description was proposed by Kawamoto et al.

(2016) but is not considered here, having in mind quasi-static simulations with no influence from the inertial quantities onto the results.

As will be discussed in more detail in section 3, the grid spacing g grid , compared with particle's characteristic length l grain obviously affects the precision of the interpolated distance field, and that of LS-DEM.

Moreover such a distance field, the contact algorithm precised below in § 2.2 introduces a second key ingredient for the method, since a LS-DEM shape also involves a set of so-called boundary nodes, being exactly located on the surface (Figure 1). These are obtained through ray tracing (e.g. [START_REF] Lin | An efficient volume-rendering algorithm with an analytic approach[END_REF]:

starting from the center of mass of a DE, as determined from the inside voxels, a half-line ray defined by its direction v is followed until crossing the DE's surface. Rays v could be chosen adopting various partitions of the pθ, ϕq space, with θ P r0; πs and ϕ P r0; 2πs being the two spherical angles. Here, boundary nodes follow a spiral path in the spirit of [START_REF] Rakhmanov | Minimal discrete energy on the sphere[END_REF], where a total number N nod of boundary nodes is located along the following spherical coordinates pθ k , ϕ k q, k P r0; N nod ´1s:

θ k " arccos ˆ´1 `1 `2k N nod ˙(1) ϕ k " πp3 ´?5qk (2)
For spheres at least, such a spiral path seeds boundary nodes more uniformly over the particle's surface, when compared with a rectangular partition of the pθ, ϕq space. As a matter of fact, it avoids an overdiscretization of the poles (θ " 0 rπs) thanks to the non-constant step in θ. For each ray direction v, and due to the trilinear description of distance within each grid voxel, the ray-surface intersection can be obtained solving the roots of a cubic polynom, giving the position of boundary nodes.

As it will be detailed in the following paragraph, no real update of the boundary nodes, nor of the distance field is needed during LS-DEM simulations:

considering rigid particles with constant shapes, both are determined once for all at the beginning of a simulation, in reference configurations of the DE.

The present shape description appears as very general and distance fields for non-convex shapes could be readily obtained through Level Set algorithms [START_REF] Sethian | Level set methods and fast marching methods[END_REF] that also apply to such cases. Ray traced boundary nodes may also follow non-convex shapes, with the only limitation being that ray tracing leads to a maximum of one boundary node per grid cell, along a given ray, due to the trilinear description of the distance field.

Kinematics of contact from Level Set shape and boundary nodes

Contact detection between two Level Set-shaped DEs first implies an approximate neighbour search that is common to all YADE simulations, following a so-called sweep and prune algorithm working on bodies' axis-aligned bounding boxes [START_REF] Dubois | Numerical modeling of granular media composed of polyhedral particles[END_REF][START_REF] Šmilauer | [END_REF]. This leads to a reduced list of potential contacts between bodies pairs. Exact determination of contact between two bodies in this list then relies on a master-slave algorithm whereby the exact determination of interparticle distance both relies on the distance field φ B to the biggest (in volume) particle, and on the boundary nodes Ý Ý Ñ ON i (with O the origin) of the smallest particle (Figure 2). For convenience, labels 1,2 will replace in the following the mention of small or big particles, with φ 2 " φ B . Contact is then obtained for at least one boundary node Ý Ý Ñ ON i showing φ 2 p Ý Ý Ñ ON i q ď 0. Boundary nodes logically need to be numerous enough to avoid bias in the LS-DEM results through missing contacts if φ 2 p Ý Ý Ñ ON i q ą 0 @N i , as it will be investigated in the following sections.

After detecting at least one boundary node of 1 touching 2, the interaction description is based on the node N c showing the greatest penetration, leading to the following interparticle overlap u n :

u n " ´minpφ 2 p Ý Ý Ñ ON i q, Ý Ý Ñ ON i P S 1 q " ´φ2 p ÝÝÑ ON c q ě 0 (3)
The current "greatest penetration" choice follows classical contact laws in DEM and corresponds to another recent LS-DEM study [START_REF] Li | Capturing the inter-particle force distribution in granular material using LS-DEM[END_REF]. On the other hand, LS-DEM was initially proposed by [START_REF] Jerves | Effects of grain morphology on critical state: a computational analysis[END_REF]; Kawamoto et al.

(2016) with a mechanical interaction at each contacting node, which used to make the model behavior directly dependent on the number of boundary nodes, in addition to the k n and k t stiffnesses discussed below. That other choice would still enable to address non-convex shapes, which is not done here.

While the overlap u n serves as the normal relative displacement, the present contact treatment does not resort to any total tangential displacement but just to an incremental one at the subsequent stage of applying the contact law, see the next § 2.3. The normal and tangential contact directions actually refer to the normal to S 1 at N c , chosen as the contact normal:

n " ∇φ 1 p ÝÝÑ ON c q (4)
For simplicity, special shapes showing pathological definitions of the normal, with tips or edges, are not considered here.

For e.g. the purpose of subsequent torque computations, a contact point x c is defined in the middle of the overlap between 1 and 2:

x c " ÝÝÑ ON c ´un 2 n (5)
Considering the rigid bodies transformations of 1 and 2, the current contact algorithm easily makes use of the initial distance field and boundary nodes, as defined in the previous § 2.1 in reference configurations.

In line with its master-slave nature, such a contact treatment is not symmetric and this could be seen as a possible source of inaccuracy in the contact model in the sense different results could have been obtained adopting other choices, using e.g. φ 2 instead of φ 1 in Eq. ( 4). It is however reasonably believed that a sufficient discretization of particle's surfaces with many boundary nodes would cancel this possible bias. One should also note that the present choice of the smallest particle for carrying the boundary nodes allows to explore distance fields (whose precision depends upon grid resolution only) with the greatest surface density in nodes.

Mechanics of contact

Once a contact is detected and kinematically described as presented in the above, classical elastic (resp. elastic-plastic) contact laws apply in the normal (resp. tangential) directions, with k n and k t the normal and tangential stiffnesses and µ the contact friction coefficient.

The repulsive normal force F n is first given by:

F n " k n u n n (6)
In the tangent plane, the frictional tangential force is incrementally computed from 0, one time step after another as per the following equation:

d F t " d ˜|| F t || F t || F t || ¸" || F t ||d ˜ F t || F t || ¸`dp|| F t ||q F t || F t || (7)
In the rhs of Eq. ( 7), the first term just accounts for a possible change in the tangential force direction (its unit vector F t {|| F t ||) while the interacting pair would move as a rigid body with possible variations in the orientation of the tangent plane. This first term is computed from the previous and current normal directions and from the angular velocities of each DE ( Šmilauer et al., 2015). On the contrary, the last term in Eq. ( 7) accounts for the force variation due to a incremental tangential relative displacement, d u t , as computed at the contact point between the two DEs. A classical elastic-plastic force-displacement relationship here applies:

dp|| F t ||q F t || F t || " k t d u t enforcing || F t || ď µ|| F n || (8)
The interaction force being determined, an associated torque is also imposed with a possible contribution of the normal force for arbitrary shapes, unlike spheres.

Equations of motion

Sustaining resultant forces and torques, each DE is classically characterized in space using xptq, the current position of its center of mass P , as well as a rotation matrix Rptq that describes its current orientation, i.e. the orientation of the local frame of eigenvectors for the inertia matrix, p e i q, i P r1; 3s, as seen in the global frame. The rotation matrix R actually transforms each vector u L of the local frame in its current counterpart in the global frame u G through classical change of basis relation u G " R u L . Newton-Euler equations for the motion of rigid bodies then rule the evolutions of v, the velocity of point P and of ω, the angular velocity of the body:

m d v dt " f (9) I d ω dt ` ω ^I ω " Γ (10)
, with f the resultant force on the DE and Γ the resultant torque computed at the center of mass P . For the purposes of deriving Eq. ( 10) Γ and ω are expressed in the local frame p e i q, where I components are constant. We recall that Eq. ( 10) would simplify to Id ω{dt " Γ for simple, isotropic, shapes with a spherical inertia matrix I " kδ (with δ the identity matrix), such as spheres or cubes.

Global damping is classically considered, modifying the resultant forces and torques in Eqs. ( 9)-( 10) in dynamic cases where those are non-zero. A damping coefficient D, taken here equal to 0.2, enters the equations such that the right hand sides of Eqs. ( 9)-( 10) actually are p1 ˘Dq f or p1 ˘Dq Γ, depending on the power of resultant forces or torques. Accelerating cases with a positive power are hindered, considering p1 ´Dq, while decelerating conditions with a negative power are amplified through the use of p1 `Dq.

Time variations of position and orientation finally follow from the above

Newton-Euler equations as per:

d x dt " v (11) dR dt " R Ω (12)
, with Ω in Eq. ( 12) being the antisymmetric matrix such that Ω x " ω ^ x, @ x. Integrating these Eqs. ( 9) to ( 12) is achieved in YADE from appropriate explicit numerical schemes and using a quaternion equivalent for the rotation matrix R ( Šmilauer et al., 2015).

Precision of LS-DEM

Materials and methods

The precision of LS-DEM in connection with boundary nodes and grid spacing is now investigated for different kinds of numerical simulation, comparing when possible LS-DEM with classical DEM serving as a numerical reference.

For comparison purposes, ideal spherical shapes are then often adopted, since they enable one to obtain such a DEM reference result. The distance fields necessary to LS-DEM are straightforward to define for spheres of given radii.

Extending towards arbitrary shapes, superellipsoids, also known as superquadrics [START_REF] Barr | Rigid physically based superquadrics[END_REF], are also considered. Generalizing ellipsoids, they constitute a convenient choice for exploring non-spherical shapes, e.g. [START_REF] Wang | Superellipsoid-based study on reproducing 3D particle geometry from 2D projections[END_REF], since they offer an analytical description through three radii r x , r y , r z distorting length along the three axes, combined with two additional exponents Shape index Half-extents (length unit) Curvature exponents (-)

r x r y r z ǫ e ǫ n 0 0.4 1 0.8 0.4 1.6 1 0.42 = 0.83 0.1 1 2 = = = 1 0.5 3 0.5 0.7 1 1.4 1.2
Table 1: Shape parameters of the four superellipsoids shown in Figure 3 ǫ e , ǫ n that modify the surface curvature. In local axes, their surface equation namely reads:

f px, y, zq " ˜ˇˇˇx r x ˇˇˇ2 ǫe `ˇˇˇy r y ˇˇˇ2 ǫe ¸ǫe ǫn `ˇˇˇz r z ˇˇˇ2 ǫn ´1 " 0 (13)
While such an analytical description is not required in LS-DEM, it aptly provides a first order approximation for the signed distance function to a superellipsoid, which is herein simply proposed as:

φ « f || ∇f || (14)
Eq. ( 14) obviously describes a zero distance, φ " 0, along the surface. It is furthermore easily verified that the Eikonal equation defining distances, || ∇φ|| " 1 [START_REF] Sethian | Level set methods and fast marching methods[END_REF], is by construction verified at the first order close to the surface. This approximation, illustrated in Figure 3, is sufficient for typical LS-DEM simulations with negligible overlaps since an accurate distance field is then necessary close to the surface only.

The Table 1 lists a chosen set of 4 shape parameters, with the corresponding 4 different superellipsoids being depicted in Figure 3. The radii r x , r y , r z shown therein will be scaled to appropriate lengths in the following.

Regardless of the shape or the modelling approach (DEM or LS-DEM) chosen thereafter, the same contact parameters and particle size distribution are used, see 1, illustrated together with their distance fields (right). Image scales are constant for each shape (on each row), and the positive range of color maps (shape's interior) is capped to 0.2 length units for convenience size distribution. Doing so, a sieve diameter is chosen for each superellipsoid as the diameter of its circumscribed sphere, i.e. twice the greatest center-boundary node distance.

k n k t {k n µ D min D max {D min (N/m) (-) (-) (cm) (-) 6ˆ10 5
0.3 0.577 6.1 3

Single contact description

The precision of LS-DEM is first analyzed for the simple case of a single contact between two spherical particles, with a possible discrepancy in size (Fig-

Isotropic reconstruction

A second examples devotes to the LS-DEM reconstruction of a dense packing of 8000 spherical particles. While the current reconstruction procedure is essentially similar to the definition a LS-DEM sample from an experimental one, e.g. through computed tomography [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF][START_REF] Kawamoto | All you need is shape: Predicting shear banding in sand with LS-DEM[END_REF], it actually here applies to DEM data describing the isotropic state of a numerical Color maps show the overlap unpx, yq of a LS-DEM interaction between a sphere 1 centered at pxc, yc, zcq and a bigger sphere 2 centered at pxc `pr `R2 q cospθq, yc `pr `R2 q sinpθq, zcq with pr, θq the polar counterparts to the cartesian px, yq. The origin of the map, x " y " 0, for instance corresponds to the center of 1 belonging the surface of 2, and to an expected overlap value equal to R 1 . White region correspond to the absence of an interaction. Each map is constructed using 401 2 colored pixels and as many relative configurations of the two spheres estimation. As a matter of fact, a 80% precision can here be obtained choosing tN nod ; D{g grid u either as t2500;50u or t1600;90u. Among the cases tested, a maximum precision of 94 % is reached for 9000 boundary nodes and a grid resolution of 90, which is another step towards validating the present LS-DEM implementation with respect to DEM and investigating the role of its technical ingredients tN nod ; D{g grid u. This is pushed further in the following section.

Triaxial compression

Another comparison between DEM and LS-DEM for spherical shapes eventually considers the triaxial compression of that same dense sample, under the confining stress σ 2 " σ 3 " 16.5 kPa and until an axial strain ε 1 " 5 %. This axial strain value is posterior to the peak in deviatoric stress q " σ 1 ´σ3 that is observed in DEM.

Again, several LS-DEM simulations are carried on, for N nod P t100;400;1600; 2500;4000u and D{g grid P t10;20;50u. Any LS-DEM simulation starts with the same sample definition than before, defining appropriate Level Set shaped bodies from the DEM data that describe the isotropic stress p ref " 16.5 kPa.

Because the same mechanical state is not directly captured within LS-DEM, confining phase is pursued further, with a servo-control of boundary walls until that reference isotropic stress p ref is re-obtained. Then, both DEM and LS-DEM simulations apply triaxial shear loading with a constant axial strain rate 9 ε 1 that corresponds to an inertial number I " 9 ε 1 D 50 a ρ{σ 3 « 10 ´4 low enough for its influence and the one of global damping to vanish. It is actually verified in DEM and LS-DEM that stresses measured along the boundary walls equal homogenized Love-Weber stresses [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF][START_REF] Weber | Recherches concernant les contraintes intergranulaires dans les milieux pulvérulents[END_REF][START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF] for static equilibrium conditions. Table 3 details relevant parameters, with a fictitious ρ " 1000 kg/m 3 density being herein adopted.

The latter could be replaced by another value provided that time step and loading rate are also modified in order to avoid divergence of the explicit scheme and maintain the same inertial number. Such changes would keep constant the total number of DEM iterations required for simulating triaxial shear until ε 1 " 5 %. On that second example, the LS-DEM precision is quantified comparing the deviator peak q max of each LS-DEM simulation with the reference DEM value q max ref « 33 kPa, through a q max {q max ref ratio that is illustrated in the Figure 7. Similar trends in precision are observed on this third example, with a joint influence of the grid resolution and the number of boundary nodes. This being said, the present DEM vs LS-DEM comparison with non-fixed DEs under deviatoric loading is more favorable than the isotropic reconstruction. Indeed, using 4000 boundary nodes and a grid resolution of 50 now enables one to reach an excellent 96% overall precision, whereas it previously led to just 85% for the isotropic example. This 85% precision would here be exceeded choosing tN nod " 400; D{g grid " 20u only. The particular case of tN nod ě 1600; D{g grid " 10u illustrates the marginal possibility for a nonmonotonous increase in precision with respect to N nod . One may think for instance to the very specific case of two spheres in contact that could be perfectly described with just one boundary node located along their branch vector.

In addition to the only consideration of peak deviatoric stress, the Figure 8 illustrates the effects of tN nod ; D{g grid u choices onto the evolutions of other average quantities according to axial strain. LS-DEM is therein also compared with DEM for what concerns the volumetric strain ε V , the anisotropy a c of the contact network, and the average contact number z c . As for the contact anisotropy a c , the latter is expressed as the difference between the axial and the lateral components of the fabric tensor F whose expression is represented in the following Eq. ( 15).

F " 1 N c ÿ c n b n (15)
For the purpose of computing F in LS-DEM, it is recalled contact normals are computed in this case from the distance gradient as per the previous Eq. ( 4). The precision in evaluating this distance gradient again depends on grid resolution.

The Figure 8 confirms that the LS-DEM evaluation of any quantity of interest tends to its DEM counterpart for tN nod ; D{g grid u reaching the order of t4000;50u. It furthemore illustrates how the dense-like behavior traits, with softening and dilation, of the present numerical sample appear as diminished when using an insufficient LS-DEM discretization in terms of boundary nodes and grid resolution. One can lastly note that LS-DEM curves are generally speaking somewhat more noisy than DEM counterparts, due to the surface discretization in boundary nodes. Such a surface discretization, when poor in particular, may indeed enhance the discontinuous i.e. sudden changes in overlap and contact forces already present in DEM due to the time discretization, possibly affecting the curves at the macro-scale.

Triaxial compression of superellipsoids

A last example devotes to a packing of 8000 superquadrics, as defined in the above § 3.1, under the same triaxial loading than the one imposed on spherical particles. After reaching the isotropic state (Figure 9) p " 16. for spheres, triaxial shear is again pursued until an axial strain ε 1 " 5 % being posterior to the deviator's peak. Among the simulation parameters, being listed in Tables 2 and3, time step is modified from the spherical case because of a possibly lower volume, hence mass, of a superellipsoid when compared to a sphere having the same circumscribed diameter. Such a LS-DEM simulation is carried on for different choices of N nod P t400; 1600;2500;4000u and 2 minpr x , r y , r z q{g grid P t10;20;50u, disregarding here the less precise case N nod " 100. Looking at the obtained peak in q, the data illustrated in the Figure 10 once again show how both the grid resolution and the boundary nodes number affect the LS-DEM results. With respect to the ideal spherical shapes considered in the above, the results also suggest that capturing more complex shapes might be more demanding in terms e.g. of boundary nodes number N nod . While using N nod ě 1600 induced fairly constant LS-DEM results

for spheres (within a 2-3 % variation, see Figure 7), the present results on superellipsoids still vary by nearly 10 % in that range, without a clear plateau.

As for the deviator strength itself, one can also note from the most precise LS-DEM simulations that the superquadrics packing exhibits a deviator strength q max « 48 kPa, which is approximately 45 % higher than the ones for spheres (where q max ref « 33 kPa) and combined with differences in initial porosity or coordination number. A greater ultimate triaxial strength at critical state is also obtained, with M " q{p « 0.76 for spheres, versus M « 1.13 for superel- lipsoids using N nod " 2500 and 2 minpr x , r y , r z q{g grid " 20 until ε 1 " 40 %.

While further discussion is left for future work, these results confirm the shape influence upon the mechanical properties.

Discussion

From the comparisons shown in the above, and with a greater focus on the more meaningful triaxial simulation with moving DEs, one could advice to use a grid resolution (l grain {g grid ) in the order of few tenths, and a couple of thousands boundary nodes at least. Even though previous LS-DEM studies [START_REF] Jerves | Effects of grain morphology on critical state: a computational analysis[END_REF][START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF][START_REF] Kawamoto | All you need is shape: Predicting shear banding in sand with LS-DEM[END_REF] did not explicitly provide such technical details, similar order of magnitudes can be inferred as follows.

Regarding the boundary nodes, the key references [START_REF] Jerves | Effects of grain morphology on critical state: a computational analysis[END_REF][START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF] formulated the same guideline in terms of node-to-node spacing, proposing therein that restricting these distances to one tenth of particle diameter would avoid bias in the results. In addition to distance considerations, a proper set of boundary nodes should obviously cover the whole direction space θ ˆϕ " r0; πs ˆr0; 2πs. Assuming this was done in [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF] with a rectangular partition, and considering that R a ∆θ 2 `∆ϕ 2 , with ∆θ, ∆ϕ the increments in the spherical angles θ, ϕ between two adjacent nodes, is an upper bound to that node-to-node distance, one can connect node-to-node spacing to the increments ∆θ, ∆ϕ, then to the total number of nodes N nod . As such, the above distance guideline quoted by [START_REF] Jerves | Effects of grain morphology on critical state: a computational analysis[END_REF]; [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF] can eventually be related to a total number of nodes N nod being in the order of 1200. The present comparisons rather confirm this order of magnitude of thousand of boundary nodes as a minimum, and they furthermore illustrate how the grid resolution articulates with N nod for what concerns the precision of the method.

As for the grid resolution itself, no exact mention of the latter seems to be found in [START_REF] Jerves | Effects of grain morphology on critical state: a computational analysis[END_REF][START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF][START_REF] Kawamoto | All you need is shape: Predicting shear banding in sand with LS-DEM[END_REF]. One can nevertheless speculate from [START_REF] Kawamoto | Level set discrete element method for three-dimensional computations with triaxial case study[END_REF] that a resolution l grain {g grid in the order of 30 or 40 was adopted therein, which also appears to be the required order of magnitude.

To conclude, LS-DEM practice certainly requires to consider grid resolution and boundary nodes as similar technical ingredients than meshes for Finite Element Methods, and eventually to check their (non-)influence onto the results.

Computational costs

The greater flexibility of LS-DEM logically comes along greater computational costs, be in terms of memory (RAM) footprint or evaluation time. These are now carefully investigated for the triaxial compression of spherical particles until ε 1 " 5% that was considered in the previous section 3. Finally, the present cost analysis also recalls the combined influence of both boundary nodes and grid resolution onto the results. It actually illustrates the possibility for different strategies of ressource managements, when seeking a given precision. Aiming to limit RAM consumption, a 95% precision could be here obtained choosing D{g grid " 20 and 4000 boundary nodes. On the other hand, choosing D{g grid " 50 and 1600 boundary nodes would show higher memory requirements, but would lead to the same precision after faster simulations.

OpenMP scalability for parallel simulations

Parallel computing is an obvious strategy to alleviate the high time costs of LS-DEM, and is available in YADE e.g. in a OpenMP shared memory framework ( Šmilauer, 2010). The OpenMP framework distributes the treatment of DEM variables among parallel threads that will collectively move forward the simulation. Typical examples include integrating motion for different DEs with different threads, or the parallel computing of interaction forces for different interactions. However, the shared memory paradim inherently requires costly safeguards to avoid conflicts between possible operations from different threads onto the same DEM variable. One can think for instance to the resultant force of one given DE contributing to different interactions, which could be modified by different threads after parallel computations of interaction forces. After performing extra-operations to avoid such pitfalls, OpenMP speedups in YADE usually do not reach the optimal value of threads number ( Šmilauer, 2010), with possible peaks in speedup around 8 threads for spherical particles [START_REF] Zhao | A poly-superellipsoid-based approach on particle morphology for DEM modeling of granular media[END_REF].

As for LS-DEM, parallel speedups are investigated hereafter for the same triaxial shear on spheres and until ε 1 " 5 % than considered in the previous sections 3.4 and 4, using 1600 boundary nodes and a grid resolution of 20 which confered LS-DEM a sufficient precison (93%). Allocating a variable number of threads, the LS-DEM simulation is executed on the server machine mentioned in the above section 4, as well as on a workstation with one 4 cores (8 threads)

Intel i7-7700, 3.60GHz processor with 8 MB of cache memory, as well as 64 GB of 2.4 GHz RAM.

Allocated threads go from 1 to 8 for the workstation, and from 1 to 100 for the server. For each thread number j (including the sequential case j " 1), simulation time t is measured repeating 3 times the simulation to account for the possible variations in time cost. Then, 9 parallel speedups can be measured for a given j, through the 9 ratios tpjq{tpj " 1q.

After averaging among these 9 measurements and quantifying error as one standard deviation, the data (Figure 12) show LS-DEM parallel simulations follow a linear speedup until 22 threads approximately. Under those conditions the workstation shows a fairly optimal speedup, while a 0.6 speedup coefficient, 40% smaller than the optimal one, is obtained on the server. Using even more threads, simulations then continue to speed up, at a lower rate, until 50 threads approximately. For that number of threads, parallel execution is more than 20 times faster than the sequential one. The simulation speed afterwards starts to decrease with the number of threads, whereby allocating more ressources eventually just increases evaluation time.

Even though the OpenMP scalability is not necessarily optimal, significant time can then be saved in a LS-DEM simulation using an appropriate number of threads between 20 and 50. Time gains are even greater in proportion than one could get for classical DEM simulations. Indeed, the maximum parallel speed-up for the DEM simulation approximates 3.5 only, which is obtained for 10 threads approximately (Figure 13). Such a scalability corresponds to the one observed for spheres by [START_REF] Zhao | A poly-superellipsoid-based approach on particle morphology for DEM modeling of granular media[END_REF]. Allocating more threads to the DEM simulation does not bring any benefit and can even be detrimental since parallel simulations using more 60 threads are eventually slower than the sequential one. This enhanced scalability of LS-DEM versus DEM relates with the former's specificity that more than 99% of a sequential simulation is spent in contact treatment, with costly loops over boundary nodes. Other parallel paradigms such as MPI, distributing memory instead of sharing it, may be even more useful and have yet to be investigated. Together with possible code and algorithmic [START_REF] Duriez | Level set representation on octree for granular material with arbitrary grain shape[END_REF] improvements, they will hopefully make geotechnical simulations with real particle's shape even more affordable.

Conclusions and perspectives

LS

Figure 1 :

 1 Figure 1: Plane view of the 3D regular grid at the roots of shape description in LS-DEM. Exact values of the signed distance function φ are known at each grid node M (the blue cross evidences just one of them). Boundary nodes N i play a role in contact treatment as described in §2.2

Figure 2 :

 2 Figure 2: Distance field (colored map) and boundary nodes (black points) serving for the LS-DEM contact algorithm, illustrated for spherical particles

Figure 3 :

 3 Figure 3: The four superellipsoids (left) defined in Table1, illustrated together with their distance fields (right). Image scales are constant for each shape (on each row), and the positive range of color maps (shape's interior) is capped to 0.2 length units for convenience

  Figure4illustrates how these three parameters affect the LS-DEM measure of an overlap between the two spherical particles.It is for instance observed in Figure4(a) that using just 100 boundary nodes (in 3D space) leads to miss interactions close to the unit circle of the map, and to an approximation between the detected overlap and the true distance to a sphere. On the other hand, the Figure4(d) confirms the true distance field can be re-obtained with a very good precision, i.e. u n " ´φ, using D{g grid " 50 and N nod " 1600, with D 2 {D 1 " 1. Thanks to the present choice of locating boundary nodes on the smallest sphere, cases with D 2 {D 1 ą 1 are described with a greater precision, see Figure4(b) vs 4(a).

Figure 4 :

 4 Figure 4: Precision of the LS-DEM contact algorithm in capturing a sphere's distance field.Color maps show the overlap unpx, yq of a LS-DEM interaction between a sphere 1 centered at pxc, yc, zcq and a bigger sphere 2 centered at pxc `pr `R2 q cospθq, yc `pr `R2 q sinpθq, zcq with pr, θq the polar counterparts to the cartesian px, yq. The origin of the map, x " y " 0, for instance corresponds to the center of 1 belonging the surface of 2, and to an expected overlap value equal to R 1 . White region correspond to the absence of an interaction. Each map is constructed using 401 2 colored pixels and as many relative configurations of the two spheres

Figure 5 :Figure 6 :

 56 Figure 5: Geometric precision of LS-DEM in terms of porosity n after reconstructing a fully determined spherical packing in isotropic state

Figure 7 :

 7 Figure 7: Precision of LS-DEM in terms of peak strength during the triaxial loading of spherical grains

Figure 8 :

 8 Figure 8: DEM vs LS-DEM comparisons during a triaxial loading of spherical grains: effects of LS-DEM discretization onto averaged quantities

Figure 9 :

 9 Figure 9: Initial (left) and sheared (right, for ε 1 " 40 %) configurations of the superellipsoids packing under triaxial loading

Figure 10 :

 10 Figure10: LS-DEM description of the peak strength for a triaxial loading imposed on superellipsoids, choosing l grain " 2 minprx, ry, rzq.

  Figure 11(b) illustrates how the present time cost is primarily affected by the number of boundary nodes, with an increasing N nod leading to longer loops for contact treatment, in the same time it globally improves precision. For a given N nod , slight variations in time cost are observed depending on the grid resolution D{g grid , which just come from the previously mentioned non-reproducible

Figure 12 :

 12 Figure 12: OpenMP speed up for the LS-DEM triaxial compression using spherical grains

Figure 13 :

 13 Figure 13: OpenMP speed up for the DEM triaxial compression on spheres

Table 2 .

 2 The distribution of particle's diameter D is uniform in number between extreme D min and D max , whose values do not necessarily correspond to any physical entity. Numerical samples made of superellipsoids include in equal proportion the 4 shapes presented in the above and conform that same particle

Table 2 :

 2 DEM and LS-DEM mechanical parameters

  sample, showing a n ref « 0.372 porosity while subjected to an hydrostatic pressure p ref " 16.5 kPa. This pressure value corresponds to a stiffness ratio κ " k n {ppD 50 q « 300 which is an intermediate value among DEM studies. One

	1.12
	can for instance mention κ-values in the order of several hundreds up to one
	1.1
	thousand in qualitative (Duriez et al., 2018) as well as quantitative (Aboul Hosn
	et al., 2017) studies. 1.08
	As such, a first DEM simulation, whose parameters were presented in Ta-1.06
	ble 2, is run to reach that mechanical state. After exporting from the DEM 1.04
	model the positions and diameters D of all spherical particles, a LS-DEM recon-
	struction is attempted using at the particle scale different numbers of boundary 1.02
	0 nodes N nod P t0;100;400;900;1200;1600;2000;2500;4000;9000u and grid resolu-20 40 60 80 100 1
	On the other hand, in terms of mean stress p{p ref data (Figure 6) illustrate
	how grid resolution and boundary nodes both contribute to the mechanical
	precision of LS-DEM. Starting from an absence of contacts and stress in the ex-
	treme case of N nod " 0, boundary nodes obviously have to be numerous enough
	for all contacts to be detected. For a given number of boundary nodes, grid
	resolution still improves precision since it contributes to more exact locations of
	these boundary nodes, closer to the true surface, as well as to a better overlap

tion D{g grid P t10;20;30;50;90u. LS-DEM spheres being so defined from known positions and radii, reconstructed porosity n can be measured and one LS-DEM iteration is finally performed in order to also reconstruct normal contact forces being responsible for the sample's mean stress p, while preventing any movements of the DE. The obtained precision in terms of porosity or mean stress can be quantified through the n{n ref or p{p ref ratios, where a value of 1 or 100% indicates a perfect LS-DEM reconstruction of the reference case.

Porosity precision is actually independent of the boundary nodes and can be seen as geometric in nature since voxellised particles volumes are fully determined from the grid resolution. As such, the Figure

5

disregards boundary nodes number N nod and evidences how spherical morphologies can be satisfactorily described with tens of grid voxels per diameter, with the error on porosity i.e. solid volumes reducing below 4% for D{g grid ě 20.

Table 3 :

 3 DEM and LS-DEM numerical parameters for the triaxial compressions

	Density	Timestep	Damping	Loading rate
	ρ (kg/m 3 )		∆t (s)	coefficient D (-)	9 ε 1 (s ´1)
		Spheres	Superquadrics
	1000	3.4 ˆ10 ´4	1.7 ˆ10 ´4	0.2	2.5 ˆ10	´3

Acknowledgements

Financial support from the French Sud region to the LS-ENROC project is warmly acknowledged, as well as fruitful discussions with Cédric Galusinski and Frédéric Golay (Université de Toulon, IMATH). Technical assistance was furthermore provided by Jérémy Verrier (INRAE) for the server acquisition.