RoSSBi3D: a 3D and bi-fluid code for protoplanetary discs

Our code was developed assuming the gas is completely inviscid to conserve the equilibrium solution, to treat compressibility at best, and to respect the supersonic character of the Keplerian flow. These conditions correspond to the ones expected in an ideal dead zone of PPDs. Thus, in this outline, the equations governing gas dynamics in a 3D PPD are Euler’s equation coupled with the conservation of mass and energy; they read:

where $\mathcal{E}=\frac{\rho\@vec{V}^{2}}{2}+\frac{P}{\gamma-1}$ is the gas energy and $\@vec{f}_{V}$ are all volume forces applied to the gas. $\mathcal{S}_{\rho}$ refers to source or sink terms in the continuity equation while $\mathcal{S}_{\mathcal{E}}$ stands for heat transfer by conduction, radiation, or cooling. The variable $S$ is the surface enclosing the volume $\mathcal{V}$. In current version of RoSSBi3D, the volume forces applied to the gas are the central’s object gravitation, $\@vec{f}_{star}=-\frac{GM_{\odot}\,\rho}{\left(r^{2}+z^{2}\right)^{\frac{3}{2}}}\left(r\@vec{e}_{r}+z\@vec{e}_{z}\right)$, aerodynamic forces with solid particles, self-gravity (SG)²²2SG is not yet supported in the 3D version and the indirect potential due to a shift of the barycenter of mass, since the frame is centered on the star’s position. The equations governing dust dynamics are slightly different, and it was preferred to devote the entire Sect. 2.9 to solid particles. The conservative form provided by Eqs. 1 is attractive if the solution is discontinuous, which occurs when the flow is inviscid. This is one of the primary reasons that led us to use the FVM that we describe in the next section.

We summarise the set of Eqs. 1 under the compact form:

where $E_{p}$, $F_{p}$, and $S_{p}$ are the $p^{th}$ component of the Euler variables, flux, and source vectors. The expression of these quantities in cylindrical coordinates is provided in Appendix B. We use a cell-centred formulation that is variables are defined by their mean value over finite volumes. This implies that triple integrals are simply approximated on elementary volumes as follows: $\iiint\limits_{\mathcal{V}_{i,j,k}}f(E_{p})\,d\tau\approx\overline{f}(E_{p})\,\mathcal{V}_{i,j,k}$ where $\overline{f}(E_{p})$ is the mean value of function $f$ in the elementary cell of respective volume $\mathcal{V}_{i,j,k}$. For conciseness, we will use the notation $f(E_{p})=\overline{f}(E_{p})$ throughout the rest of this paper. On the other hand, the incoming fluxes are decomposed as elementary fluxes in each cell surface element. Considering the above statements, the flux and source terms of Eq. 2 can be rewritten as:

where $\@vec{n}_{s}$ is the outward unit normal vector to surface $A_{s}$. Accounting for the above approximations, the conservative form of the evolution Eqs. 2 can be simply written:

The 2D and 3D versions of RoSSBi3D solve the above equations in a polar and cylindrical coordinate system, respectively.

For the sake of conciseness, in next sections we will only treat the 1D case, which is easily generalised to three dimensions.

In the current version, we use a uniform structured grid in the azimuth and vertical directions while logarithmic in the radial direction. The choice of a logarithmic mesh in the radial direction stems from the 2D code. Indeed, a logarithmic mesh is necessary for computing SG thanks to Fast Fourier methods (Baruteau & Masset, 2008; Surville, 2013; Rendon Restrepo & Barge, 2023). Even if the logarithmic cutting is unnecessary in the 3D version, we kept it in order to save the effort during the coding stage.

In our code, we connected the numerical domain to the stationary flow in isothermal equilibrium (see Eq. 25) thanks to a set of 2 ghost cells in each direction. A sudden change in the flow at the domain boundaries could lead to waves reflection, which could affect the computation domain, and even instabilities can arise from an improper choice of BC. Therefore, in the radial direction, we apply wake-killing and sheared BC in order to avoid (i) spurious reflections which can perturb the computation domain and (ii) the triggering of Papaloizou-Pringle like instabilities in the disc inner edge (Papaloizou & Pringle, 1984; Goldreich & Narayan, 1985). We naturally use periodic BC in the azimuthal direction to close the physical domain. Finally, motivated by the work of Meheut et al. (2012b); Lin (2013), we implemented two different kinds of BC for the vertical direction: outflow BC or sheared damped BC.

In order to match high accuracy and reasonable computation effort, we use an explicit, second order Runge-Kutta Method (RKM2) time-stepping scheme for time integration. Crudely, this RKM2 is based on the first prediction of a differential equation solution which permits performing a correction and thus increases the accuracy of the final solution. The most general form of the explicit RKM2 applied to our problem is written:

where $\Delta t$ is the time step and $t_{n}$ is the time at step $n$. For our calculations, we used $a=0$, $b=1$, $\alpha=\frac{1}{2}$ and $\beta=\frac{1}{2}$.

In order to ensure the scheme stability during the time integration, a CFL-like condition must be respected at each time step (Courant et al., 1928; Toro, 2009, chapter 6.3). In the absence of magnetic fields and viscosity, this condition is fulfilled as long as the time step is less than:

where $c_{s}$ is the sound speed in the center of the cell and $N_{CFL}$ is a constant parameter, which is equaled to ${1}/{2}$ in order to improve stability. Physically, the quantity $\max\left\{||\@vec{V}||+c_{s}\right\}$ is the highest wave speed in the whole computational domain which allows interpreting the CFL requirement as a condition for avoiding multiple wave overlap within cell $i$ when the Riemann problem is solved. Indeed, only two solutions of the Riemann problem can affect cell $i$: the solution coming from left cell $i-1$ and the solution coming from the right cell $i+1$ (Toro, 2009, chapter 6.1).

From the intrinsic definition of FVM, a numerical solution possesses discontinuities at inter-cells positions, which forces the fluxes at cell boundaries to be computed. Indeed, flux terms are step functions, which are not continuous from one cell to the other. At each time step, $t_{n}$, a correct evaluation of these fluxes is done thanks to the solution of the following initial value problem (IVP):

This corresponds to a local Riemann problem, the solution of which is executed numerically thanks to an exact Riemann solver which has been taken from Toro (2009, chapter 4). Even if $\mathscr{G}_{p}(t,E_{q})$ contains source terms, which does not pose calculation problems, the Riemann solver is more sensitive to the scheme permitting to estimate fluxes. Indeed, a difficulty is that Euler variables are only defined on cell’s centers, but flux terms need to be computed on cell’s frontiers. A workaround is to improve the FVM with a more advanced spatial scheme, which could permit better estimation of fluxes between cells and, thus, obtaining accurate solutions even in the presence of shocks. Therefore, for this version of RoSSBi3D, we use (1) a three-point parabolic interpolation of Euler variables, which allows accurate computation of flux terms at cell’s boundaries, combined with (2) a MinMod flux limiter in order to avoid undesired oscillations.

All steps of this section are summarised below (for cell $i$). First, we perform the (1) Euler variables interpolation at the cell frontiers. For instance, the left and right Euler variables at frontier $i+1/2$, $E_{p}^{l}\left(i+1/2\right)$, and $E_{p}^{r}\left(i+1/2\right)$, respectively, are estimated thanks to the parabolic interpolation described in Sect. 2.7. Then, (2) fluxes are evaluated in the right and left side of each frontier. Still for frontier $i+1/2$, we get:

Finally, we proceed to the (3) Riemann problem resolution. When fluxes are estimated for the first time at frontiers $i-1/2$ and $i+1/2$, the exact Riemann solver allows computing exactly the Euler variables at frontiers and then reevaluating fluxes:

where Riemann stands for the solver. We highlight that the choice of the exact Riemann solver is justified by the absence of viscosity in the simulation box. Indeed, this solver was specially designed for the treatment of Euler’s equations.

The equilibrium solution presented in Eq. 25 could be used to improve the numerical scheme described in Sect. 2.6. First, we rewrite the Euler variables under the form:

where $E_{p}^{0}$ are the stationary solution Euler variables. $E_{p}^{\prime}$ stands for deviations to the steady axisymmetric flow: for the stationary solution, we get $E_{p}^{\prime}=0$. Thanks to this formulation, the Euler variables at cell $i$ frontiers are:

Such formulation has two advantages. First, the stability of the stationary solution is improved since we use the exact value of the Euler variables at the cell frontiers and subsequently the exact flux terms. Second, the accuracy of flux terms is improved thanks to a more accurate interpolation of Euler’s variables described in 2.6. Within this maneuver a parabolic interpolation on $E_{p}^{\prime}$ permits obtaining the left and right Euler’s variables at cell’s frontiers:

where:

In order to avoid undesired oscillations, we use a MinMod flux limiter which, when necessary, reduces the curvature of the interpolation function (Surville, 2013, Appendix).

One sensitive issue in PPD is that the stationary solution is not in a numerical equilibrium. Indeed, the scheme induces residual terms, which affect the stationary solution during the time evolution steps. For instance, the evolution of the radial velocity, written in terms of Euler variables is:

where the right hand side stems from the computation of flux and source terms (see Eq. 4 and Appendix B). For an ideal infinite mesh resolution, the r.h.s term of above equation vanishes, but in practise it is not the case. Indeed, this term depends on the accuracy of $\left[rP\right]^{(i+\frac{1}{2},k)}_{(i-\frac{1}{2},k)}$ computation, which introduces a small residual term affecting the stationary solution after few tens of orbits.

At 2D, Surville (2013, Chapter III.2.1) corrected this problem using an interpolation of Euler variables and performing the exact computation of the stationary solution source terms. However, for the current 3D code, this is a more difficult task since the same method would imply the analytic integration of source terms in the radial and vertical directions. For above reason, we adopted another technique for the 3D code proposed by Hervé Guillard, which consists on performing an empty run which allows to compute the r.h.s residual term for the equilibrium solution. This residual term is then subtracted at each time step from the right hand side of Eq. 15. Such strategy was also adopted by Meheut et al. (2010, Sect. 4.4). Other sophisticated techniques exist, but they do not show significant improvements while the present one has the advantage to be easy and rapid to implement. In Fig. 1, we report the error in different physical quantities for the equilibrium solution explored in Sect. 4.2, depending on whether the numerical scheme is balanced or not. The effectiveness of the exposed method is confirmed: the errors of the unbalanced scheme are 5 orders of magnitude higher than the balanced scheme at least.

In PPDs, particles interact with the gas mainly through friction forces that look like a headwind depending on the size of the particles and the gas free mean path. From the kinetic theory of gases, the mean-free-path of a gas molecule, $\lambda$, is defined as the length travelled between two successive collisions:

where $m_{u}$, $\mu$, and $\sigma_{eff}$ are the atomic mass, the number of atoms of the gas, and the molecule cross-section, respectively. The two regimes describing drag forces are the Epstein and Stokes regime. The former describes particles with sizes smaller than the gas mean-free-path interacting with the gas through molecular collisions and is equivalent to a ram pressure force. The latter deals with particles with a size being greater than the gas mean-free-path and experiencing a viscous fluid friction. For ideal spherical particles, the laws of friction expressions are well known (Weidenschilling, 1977; Youdin, 2010); the norm of the volume drag force is $||\@vec{f}_{drag}||=\rho_{d}\,{||\Delta\@vec{v}||}/{t}_{stop}$ with:

where $r_{p}$, $\rho_{s}$, and $\Delta\@vec{v}$ are the radius, the internal density of particles, and particles velocity relative to the gas, respectively.

Solid material is considered as a pressureless fluid interacting with gas through the aerodynamic forces introduced in the above paragraph. Therefore the equations governing dust dynamics are similar to the ones of gas except that the pressure term vanishes, which implies that there is no energy equation to be solved for this phase. The absence of pressure implies that the Riemann solver should be adapted to solid particles, but essentially it is almost the same as for gas.

The current version of RoSSBi3D supports SG for 2D simulations and accounts for the correction and generalisation of the smoothing length method to bi-fluids proposed by Rendon Restrepo & Barge (2023). On the other hand, the 3D version lacks SG calculation. Indeed, this is a long-term task requiring care in the theoretical and programming stages. Yet, we already thought about it and we provide a discussion on different possibilities for future versions in Sect. 5.3.

The presence of a disc asymmetry can lead to an offset of the mass barycenter (Zhu & Baruteau, 2016; Regály & Vorobyov, 2017). We have taken this effect into account through the gradient of the indirect gravity potential, which can be read in cylindrical coordinates:

where $A_{disc}=\int\limits_{disc}\frac{\tilde{r}^{\prime}\cos{(\theta^{\prime})}}{\left(\tilde{r}^{\prime 2}+\tilde{z}^{\prime 2}\right)^{3/2}}\,d\tilde{m}^{\prime}$ and $B_{disc}=\int\limits_{disc}\frac{\tilde{r}^{\prime}\sin{(\theta^{\prime})}}{\left(\tilde{r}^{\prime 2}+\tilde{z}^{\prime 2}\right)^{3/2}}\,d\tilde{m}^{\prime}$. We assumed that the density vertical stratification is symmetric with respect to the $z=0$ plane that permits the vertical component of the indirect gravity to be cancelled.

We presented the new bi-fluid RoSSBi3D code, which is the natural extension of RoSSBi to 3D. This 3D version is based on a FVM scheme, is third order in space, second order in time, stable, balanced, and accounts for inherent boundary instabilities, which could arise in PPDs simulations. The current version of RoSSBi3D supports the options gathered in Table 1. In the next section, we show how we succeeded in making a fast and scalable 2D and 3D code well-suitable for HPC.

Domain decomposition is one key to reach scalability and high-performance at high resolutions. As an instructive example, let us consider a 2D square computational box containing $N$ cells in each direction where the computational effort (workload) is shared among $P$ processors. In a crude approach, the workload computed by a processor is proportional to the volume of data the assigned cells contain ($\mathcal{V}_{comp}$) while communication (the data to be transferred among processors) is proportional to their surface ($\mathcal{S}_{comm}$). One may use the communication-to-computation ratio, $CC=\mathcal{S}_{comm}/\mathcal{V}_{comp}$³³3In numerical code metrics, we often use granularity defined by $G=T_{comp}/T_{comm}$ where $T_{comp}$ and $T_{comm}$ are the measured computation and communication times, respectively., which helps to roughly estimate the bandwidth requirements. Since the ratio between the available bandwidth and the computing power of the modern parallel HPC systems is typically low, it is also important to keep the $CC$ ratio low during the algorithm design. For distributing the $P$ processors in this regular grid, we consider 1D and 2D domain decompositions the features of which are gathered in Table 2. The $CC$ ratio is lower in the 2D decomposition compared to the 1D decomposition since it grows more slowly with $P$. Therefore, communication costs are smaller in the 2D case and more suitable for architectures with low bandwidth. For the reasons above, it was decided to (1) change the 1D grid decomposition to a 2D grid decomposition in the 2D code and (2) use a 3D grid decomposition in the 3D code. Fig. 2 illustrates how the computational grid could be decomposed and distributed among 128 processors.

Many bottlenecks intrinsically related to the structure of the code may suppress the potential of parallelisation. Therefore, we took advantage of the offer by Performance Optimisation and Productivity CoE (POP) that helps users and developers to optimise parallel software and understand performance issues. The initial performance assessment of the RoSSBi3D code was requested to identify areas of improvement in large resolution simulations. The assessment was conducted using Barbora cluster at IT4Innovations, National Supercomputing Center at VSB – Technical University of Ostrava. The cluster consists of 192 powerful x86-64 compute nodes, each equipped with 2x 18-core Intel Cascade Lake 6240 at 2.6 GHz processors and 192 GB of DDR4 memory, interconnected by high-speed InfiniBand HDR100 network. As the vast majority of the current supercomputers, Barbora is operated by a Linux OS, specifically Red Hat Enterprise Linux 7.9. To build and install the RoSSBi3D code, the GCC 9.3.0 toolchain together with OpenMPI 4.0.3 was used. In addition, the BSC performance toolset including Extrae 3.8.3, Paraver 4.9.0, Dimemas 5.4.2, and Basic Analysis 0.3.9 was employed to collect and analyse performance data.

The scaling of the RoSSBi code was evaluated on 1 – 16 compute nodes, i.e 36 – 576 cores, respectively. During the initial stage of the analysis, it turned out that full-length simulations comprising millions of time steps produce an unmanageable amount of collected performance data. Using high-level techniques for detecting the computation structure, it was proved that the performance characteristics do not evolve in time significantly. Therefore, further analysis was done with the simulation duration being reduced to the order of low hundreds of time steps. The detailed computation structure analysis showed that the code had issues with computation granularity that were indicated by the very short computation bursts – the sequences of user-code instructions outside the MPI library runtime. Most of the bursts were shorter than 500 $\mu$s. The scalability tests showed rapid speedup degradation on more than 2 nodes as depicted in Fig. 3 (original version).

The following investigation revealed significant growth of the communication phase that consisted of a sequence of loops with many short, mostly collective, MPI calls. Moreover, the loops overhead also caused the increased number of instructions issued. Each MPI call or user instruction causes an additional latency that increases the total runtime and effectively degrades the scalability. To narrow this limiting factor, the code was refactored in the way of fusing the loops and collapsing the thousands of the small MPI messages transferred in each time step into only tens of larger messages. This enhancement significantly improved efficiency of the inter-node communication through the network and also reduced the instruction overhead. Another important optimisation was improving the insufficient parallelisation of 3D domains in all directions (see Sect. 3.1). Indeed, the original version of RoSSBi3D was only parallelised in the azimuthal direction while the updated version used a 3D domain decomposition.

The follow-on performance assessment demonstrated that the applied optimisations significantly improved the scaling as shown in Fig. 3 (updated version). During the analysis, we found that also the MPI process configuration and domain decomposition, as explained in Sect. 3.1, has a non-negligible impact on performance. For example, we obtained the lower speedup value (13.62) of the updated version in Fig. 3 with the 8x36x2 configuration of processes, while the higher value (15.13) using the 8x18x4 configuration. Table 3 reports up to 4.9x speedup on the same largest test case. Also increased average CPU frequency was observed. For future scaling to even larger test cases, some optimisation opportunities in the communication pattern were also identified. The single-core performance might be enhanced by analysing and optimising the vectorisation opportunities.

We also compared the performances of RoSSBi3D when compiled with $\{$OpenMPI library + GCC compiler$\}$ and with $\{$IntelMPI library + ifort compiler $\}$ for a standard test on 64 processors in CeSAM. The test reports a 1.86 speedup when compiled with Intel compilers that we recommend for Intel architectures. This is particularly interesting since Intel^® compilers are free to use now.

The explosion test is a 2D analogue of the sod shock-tube test which permits both the accuracy of a numerical scheme and its ability to keep the cylindrical symmetry of the numerical solution to be assessed (Toro, 2009, chap. 17). This test consists of a 2D Riemann problem in which, as an initial state, a high pressure region and a low pressure region are separated by a diaphragm:

For this test, we take the ratio of specific heats $\gamma=1.4$, and the number of cells in the radial and azimuthal direction is $(N_{r},N_{\theta})=(1600,3200)$. In Fig. 4, we show this test results at $t=2.5$, which are in good agreement with the reference solution of Toro (2009, Fig. 17.4).

We use as initial state the equilibrium solution exhibited in Eq. 25 and checked that this solution remained stationary during 300 orbits. The error evolution of all physical variables, shown in Fig. 1, remained smaller than $10^{-9}$, which is satisfactory.

PPDs are subject to instabilities, such as the RWI (Lovelace et al., 1999b; Li et al., 2000). Numerical and theoretical work was already considered at 3D for this instability by few authors and in this section, we aim to retrieve Meheut et al. (2010, 2012a) results. As an initial state, we introduced a radial Gaussian overdensity bump in the stationary solution provided by Eq. 25

To facilitate the trigger of the unstable modes we introduced a small perturbation in the radial velocity field:

where $\epsilon=10^{-4}$ and $m$ is the perturbation mode. We ran 5 simulation tests for $m=[2,3,4,5,6]$ the outline of which is the same as (Meheut et al., 2012b). In Fig. 5 we depicted the exponential growth rate obtained for each mode. Fig. 6 shows the perturbed density evolution in the $z\leq 0$ region for the $m=4$ mode up to instability saturation and in Fig. 7, we illustrated in detail the vortex generated by the $m=4$ mode and its vertical distribution. We found that the growth rate values and the curve trend are in agreement with Meheut et al. (2012b, Fig. 7).

Vortices are gaseous structures whose main interest is their strong ability to capture solid particles with nearly optimal Stokes numbers (Barge & Sommeria, 1995; Tanga et al., 1996). Here we propose to find this trapping mechanism in a 3D vortex using RoSSBi3D.

For this test, we use the vortex obtained after RWI saturation in the $m=4$ mode of Sect. 4.3. Fig. 7 shows this 3D gaseous vortex used as an input for this test. Then we introduce uniformly dust in the radial and azimutal direction but with a vertical Gaussian distribution:

where $H_{d,0}/H_{g}=0.01$ and $Z_{ini}=10^{-4}$. We chose particles interacting with the gas only through aerodynamic forces (friction and dust back-reaction) and with a nearly optimal Stokes number ($St=0.5$). In Fig. 8, the dust density is exhibited at different times. We do not show the gaseous vortex evolution since it is nearly unaffected by the dust back-reaction.

In the early stages, we observe that dust settles in the midplane and starts to be captured by the vortex. After $\sim 10$ orbits almost all dust lies exactly in the midplane: the settling time is consistent with the theoretical prediction $\tau_{sed}/t_{0}\sim\left(St+\frac{1}{St}\right)=2.5$. Because of the absence of viscosity the dust scale vanishes (at the limit of the vertical resolution). The gaseous vortex collects dust during the whole simulation and after 110 orbits the dust density is increased by a factor 2000 in the vortex core which is in agreement with the theoretical expectations.

According to Gheller et al. (2015) and Wei et al. (2020), GPU technology is promising for solving Euler’s equations. However, the current code architecture and the observed workload on processors in tests conducted do not permit confirming yet if RoSSBi3D performance would benefit notably from GPU acceleration. Another audit or proof of concept analysis would be necessary to elucidate this question, but this is not a current priority. Furthermore, the parallelisation philosophy in the CUDA language is different from MPI and this would imply a significant rewriting of the code. Alternatively, the OpenACC language could be employed for prototyping GPU acceleration since it is way less intrusive than CUDA and can potentially reach similar performances. Finally, the metrics and scalability of the code illustrated in Table 3 and Fig. 9 are satisfactory, and for the current use there is no real need to migrate to GPU.

Currently the calculation of 2D SG is based on convolution products computed thanks to fast Fourier transforms. In Fig. 10, we compare the scalability of the 2D version of RoSSBi3D with and without SG with respect to the number of processors. We found that with SG, the scalability starts to show limitations from 16 nodes (576 processors), whereas such behaviour is not observed without SG. This suggests that SG calculation is responsible for the loss of performances for a large number of nodes. Possible improvements that permit alleviating this issue consists in performing in-place transforms, which is particularly interesting if each processor has to treat large amounts of data. Another idea is to carrying out a single transform on a unique multidimensional array instead of doing successive transforms on multiple arrays. This could be beneficial for bi-fluid simulations where 6 Fourier transforms need to be performed at each time step, compared to 2 for a single fluid simulation (Rendon Restrepo & Barge, 2023).

SG is a rather scarce option in 3D hydrodynamic codes. There are only a few codes with this option, e.g.: the PENCIL Code (Pencil Code Collaboration et al., 2021), RAMSES (Teyssier, 2002), NIRVANA (Ziegler, 2004) or FLASH (Fryxell et al., 2000). The scarcity mainly comes from computational costs and trade offs between accuracy and computational efficiency. In the next paragraph, we summarise few of the main methods for computing SG and identify those that could help us.

The numerical computation of the SG potential and the associated forces are handled by two main approaches. The first one is called the difference method and is based on the resolution of the Poisson’s equation:

where $\Phi_{sg}$ is the gravitational potential. This equation discretization results in a linear system which is usually solved using iterative multigrid methods as it was done by Ziegler (2005) in NIRVANA. Nonetheless, this requires evaluating the SG potential at the numerical domain boundaries which could be done by multipole expansion. The second approach called integral method relies on the equivalent integral form:

This summation could be performed using a multipole expansion of $||\@vec{r}-\@vec{r}^{\prime}||^{-1}$, which is explained in detail by Wongwathanarat (2019). This second formulation is also compatible and appealing for some Adaptive Mesh Refinement (AMR) codes where one can take advantage of the tree structure: in a crude way, the self-gravity calculation is refined in the neighbourhood of a cell or regions of strong density gradients while long distance and slow varying density regions are averaged. This method is very efficient numerically and is used by Wünsch et al. (2018). Finally, fast Fourier methods could also be used for both formulations and different types of geometries (Müller & Chan, 2019; Moon et al., 2019; Krasnopolsky et al., 2021).

This is a development often considered to address the complex evolution of PPDs. Numerical simulations with two phase codes generally describe the solid material as a single population of spherical and solid particles. This is of course a major draw-back which biases the estimate of dust-gas friction and, consequently, the description of the disc evolution. A natural possibility to overcome this problem is to use a multi-fluid approach as, for example, the ones implemented by Hanasz et al. (2010), Benítez-Llambay et al. (2019) and Huang & Bai (2022).

However, in a view of a more realistic evolution, it is also important to describe the exchanges of mass and energy between the various fluids of the dispersed phase and between the two phases, gas and solids. This refers to the formation of planetesimals, a problem in which binary collisions and gas drag must be taken into account and result in coagulation/fragmentation or condensation/evaporation processes (Weidenschilling, 1980; Voelk et al., 1980; Weidenschilling & Cuzzi, 1993; Brauer et al., 2008; Birnstiel et al., 2011). An appropriate way to implement a multi-fluid approach could be to use a sectional method (Laurent & Massot, 2001; Laurent, 2002), which is well known in the context of the dispersed droplet evolution (Dufour & Villedieu, 2005) and can be interpreted in terms of the kinetic theory of gases (Tambour, 1985; Laurent & Massot, 2001). The advantages of this method are to ensure conservation of moment transfer between the various fluids and, possibly, to rely on the large developments of the numerical simulation of reactive flows (Drui, 2017).

In order to explain the angular momentum transport in PPDs, it is common to associate local turbulence, generated in various regions of the disc, with a pseudo-viscosity (Shakura & Sunyaev, 1973). Promising mechanisms generating turbulence include the magneto rotational instability (MRI) (Balbus & Hawley, 1991) and the VSI (Nelson et al., 2013; Stoll & Kley, 2014), among others. Such mechanisms encourage to include viscosity in the RoSSBi3D code for a realistic description of PPD. However, accounting viscous effects modifies the nature of the equations to be solved: Euler’s equations are hyperbolic whereas Navier-Stokes equations are parabolic. The presence of diffusion causes the Riemann solver to be inadequate for solving the Navier-Stokes equations, forcing an alternative numerical strategy. A possible solution could be to use an explicit Runge-Kutta-Legendre method (Meyer et al., 2014).

The confrontation of numerical models with dust thermal emission (Andrews, 2015) is key since it allows discriminating between physical processes operating on PPDs. Such comparison is made possible thanks to dust radiative transfer codes, see (Steinacker et al., 2013, and references therein). Therefore, an important step is to include radiative transfer in the post-processing pipeline from another existing code. We think that the best option is to use a free, robust, documented, and up to date existing code: the most attractive solutions seem to be RADMC-3D 2.0 or MCFOST.