diff --git a/configure.ac b/configure.ac index 212dd5ee47bbcb6951426e2285e5b5bb00ac1aaf..80c3a8d8d4d9ce8bfb769d14a30b84ee62084e4f 100644 --- a/configure.ac +++ b/configure.ac @@ -1963,7 +1963,8 @@ AC_SUBST([SUNDIALS_INCS]) # As an example for this, see the call to AC_ARG_WITH for cooling. AC_ARG_WITH([subgrid], [AS_HELP_STRING([--with-subgrid=<subgrid>], - [Master switch for subgrid methods. Inexperienced user should start here. Options are: @<:@none, GEAR, AGORA, QLA, QLA-EAGLE, EAGLE, EAGLE-XL, SPIN_JET_EAGLE default: none@:>@] + [Master switch for subgrid methods. Inexperienced user should + start here. Options are: @<:@none, GEAR, GEAR-G3, AGORA, QLA, QLA-EAGLE, EAGLE, EAGLE-XL, SPIN_JET_EAGLE default: none@:>@] )], [with_subgrid="$withval"], [with_subgrid=none] @@ -1990,12 +1991,24 @@ case "$with_subgrid" in GEAR) with_subgrid_cooling=grackle_0 with_subgrid_chemistry=GEAR_10 + with_subgrid_pressure_floor=GEAR + with_subgrid_stars=GEAR + with_subgrid_star_formation=GEAR + with_subgrid_feedback=GEAR + with_subgrid_black_holes=none + with_subgrid_sink=GEAR + with_subgrid_extra_io=none + enable_fof=no + ;; + GEAR-G3) + with_subgrid_cooling=grackle_3 + with_subgrid_chemistry=GEAR_10 with_subgrid_pressure_floor=none with_subgrid_stars=GEAR with_subgrid_star_formation=GEAR with_subgrid_feedback=GEAR with_subgrid_black_holes=none - with_subgrid_sink=none + with_subgrid_sink=GEAR with_subgrid_extra_io=none enable_fof=no ;; diff --git a/doc/RTD/source/SubgridModels/GEAR/sinks/index.rst b/doc/RTD/source/SubgridModels/GEAR/sinks/index.rst index 20d3fffcf36e73b727ffd8446d17d116d1a8f804..699dfa53dcdb7991c914bcc4c5c91ff2382758ea 100644 --- a/doc/RTD/source/SubgridModels/GEAR/sinks/index.rst +++ b/doc/RTD/source/SubgridModels/GEAR/sinks/index.rst @@ -13,5 +13,6 @@ Sink particles and star formation in GEAR model introduction theory + sink_timesteps params output diff --git a/doc/RTD/source/SubgridModels/GEAR/sinks/output.rst b/doc/RTD/source/SubgridModels/GEAR/sinks/output.rst index 90fe867ba075db0a03136f47c7b33633f9d48315..fda4c6334cbef571c55cbb4d562830a9eef7b9c6 100644 --- a/doc/RTD/source/SubgridModels/GEAR/sinks/output.rst +++ b/doc/RTD/source/SubgridModels/GEAR/sinks/output.rst @@ -37,6 +37,14 @@ Sink particles | | | element | | | number of elements ``N`` is determined at | | | | | | | compile time by ``--with-chemistry=GEAR_N``. | +---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+ +| ``BirthScaleFactors`` | | Scale-factors when the sinks were | [-] | Only used in cosmological runs. | +| | | born | | | +| | | | | | ++---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+ +| ``BirthTimes`` | | Time when the sinks were | [U_T] | Only used in non-cosmological runs. | +| | | born | | | ++---------------------------------------+-------------------------------------+-----------+---------------------------------------------------+ + Stars diff --git a/doc/RTD/source/SubgridModels/GEAR/sinks/params.rst b/doc/RTD/source/SubgridModels/GEAR/sinks/params.rst index 61d502935b490b7af3a2c432431403c0813b1c8d..9f8409de7e95f48099c09684a618a07a6853b33e 100644 --- a/doc/RTD/source/SubgridModels/GEAR/sinks/params.rst +++ b/doc/RTD/source/SubgridModels/GEAR/sinks/params.rst @@ -44,7 +44,19 @@ The next set of parameters controls the sink formation scheme. More details are Those criteria are checked if the density and temperature criteria are successfully passed. They control the behaviour of the sink formation scheme. By default, they are all activated and set to ``1``. -The last parameter is ``disable_sink_formation`` (default: 0). It controls whether sinks are formed or not in the simulation. The main purpose is when sinks are put in initial conditions and sinks are not wanted to be added during the run. This parameter is set to ``0`` by default, i.e. sink formation is *enabled*. +The next parameter is ``disable_sink_formation`` (default: 0). It controls whether sinks are formed or not in the simulation. The main purpose is when sinks are put in initial conditions and sinks are not wanted to be added during the run. This parameter is set to ``0`` by default, i.e. sink formation is *enabled*. + +The next set of parameters deals with the sink time-steps: + +* Courant-Friedrich-Levy constant for the CFL-like time-step constraint:``CFL_condition``, +* age (in Myr) at which a sink is considered dead (no accretion) and without time-step limitations, except for 2-body encounters involving another young/old sink and gravity: ``timestep_age_threshold_unlimited_Myr`` (default: FLT_MAX), +* age (in Myr) at which sinks switch from young to old for time-stepping purposes:``timestep_age_threshold`` (default: FLT_MAX), +* maximal time-step length of young sinks (in Myr): ``max_timestep_young_Myr`` (default: FLT_MAX), +* maximal time-step length of old sinks (in Myr): ``max_timestep_old_Myr`` (default: FLT_MAX), +* number of times the IMF mass can be swallowed in a single time-step: ``n_IMF`` (default: FLT_MAX). +* tolerance parameter for SF timestep constraint: ``tolerance_SF_timestep`` (default: 0.5) + +The last parameter is ``sink_minimal_mass_Msun``. This parameter is mainly intended for low-resolution simulations with :math:`m_\text{gas} > 100 \; M_\odot`. It prevents :math:`m_\text{sink} \ll m_\text{gas}` simulations when sinks spawn stars, which can lead to gravity run away. The full section is: @@ -56,6 +68,7 @@ The full section is: temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. density_threshold_Hpcm3: 1e3 # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle. maximal_density_threshold_Hpcm3: 1e5 # If the gas density exceeds this value (in Hydrogen atoms/cm3), a sink forms regardless of temperature if all other criteria are passed. (Default: FLT_MAX) + sink_minimal_mass_Msun: 0. # (Optional) Sink minimal mass in Msun. This parameter prevents m_sink << m_gas in low resolution simulations. (Default: 0.0) stellar_particle_mass_Msun: 20 # Mass of the stellar particle representing the low mass stars (continuous IMF sampling) (in solar mass) minimal_discrete_mass_Msun: 8 # Minimal mass of stars represented by discrete particles (in solar mass) stellar_particle_mass_first_stars_Msun: 20 # Mass of the stellar particle representing the low mass stars (continuous IMF sampling) (in solar mass). First stars @@ -67,6 +80,13 @@ The full section is: sink_formation_bound_state_criterion: 1 # (Optional) Activate the bound state check for sink formation. (Default: 1) sink_formation_overlapping_sink_criterion: 1 # (Optional) Activate the overlapping sink check for sink formation. (Default: 1) disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0) + CFL_condition: 0.5 # Courant-Friedrich-Levy condition for time integration. + timestep_age_threshold_unlimited_Myr: 100. # (Optional) Age above which sinks no longer have time-step restrictions, except for 2-body encounters involving another young/old sink and gravity (in Mega-years). (Default: FLT_MAX) + timestep_age_threshold_Myr: 25. # (Optional) Age at which sink switch from young to old for time-stepping purposes (in Mega-years). (Default: FLT_MAX) + max_timestep_young_Myr: 1.0 # (Optional) Maximal time-step length of young sinks (in Mega-years). (Default: FLT_MAX) + max_timestep_old_Myr: 5.0 # (Optional) Maximal time-step length of old sinks (in Mega-years). (Default: FLT_MAX) + n_IMF: 2 # (Optional) Number of times the IMF mass can be swallowed in a single timestep. (Default: FLTM_MAX) + tolerance_SF_timestep: 0.5 # (Optional) Tolerance parameter for SF timestep constraint. (Default: 0.5) .. warning:: Some parameter choices can greatly impact the outcome of your simulations. Think twice when choosing them. diff --git a/doc/RTD/source/SubgridModels/GEAR/sinks/sink_timesteps.rst b/doc/RTD/source/SubgridModels/GEAR/sinks/sink_timesteps.rst new file mode 100644 index 0000000000000000000000000000000000000000..f31551bd1b595e40ebf4f6a5a4a64bdd7170e89f --- /dev/null +++ b/doc/RTD/source/SubgridModels/GEAR/sinks/sink_timesteps.rst @@ -0,0 +1,66 @@ +.. Sink particles in GEAR model + Darwin Roduit, 24 November 2024 + +.. _sink_GEAR_timesteps: + +Sink timesteps +~~~~~~~~~~~~~~ + +Sink particles interact with the surrounding gas through accretion. To accurately follow the local gas dynamics and assess the stability of gas hydrodynamics, we must impose timestep constraints on the sink particles. + +First, we want to ensure sinks know the *local dynamics* and thus they obey a Courant-Friedrichs-Lewy (CFL)-like timestep constraint: + +.. math:: + \Delta t_s \leq \Delta t_\text{CFL} = C_\text{CFL} \frac{r_{\text{cut, min}}}{\sqrt{c_{s,s}^2 + \| \Delta \mathbf{v}_s \|^2}} \; , + +where :math:`r_\text{cut, min} = \min(r_{\text{cut}, s}, \gamma_k \min_j(h_j))` is the minimal cut-off radius between the sink :math:`s` and the gas neighbours:math:`j`, :math:`c_{s, s}` and :math:`\Delta \mathbf{v}_s` are the gas sound-speed and the relative gas-sink velocity at the sink location. The latter two are reconstructed at the sink location with SPH interpolation. The value of :math:`C_\text{CFL}` is given in the YAML parameter file with ``GEARSink:CFL_condition``. + +Since sink particles accrete gas, they must anticipate the surrounding gas infall. We achieve this with a gas free-fall time criterion similar to `Grudic et al. (2021) <https://academic.oup.com/mnras/article/506/2/2199/6276745>`_: + +.. math:: + \Delta t_s \leq \Delta t_\text{ff} = \sqrt{ \frac{3 \pi}{32 G \rho_s} } \quad \text{with} \quad \rho_s = \frac{3 m_s}{4 \pi {r_{\text{cut, min}}}} \; , + +with :math:`m_s` the sink mass. + +These constraints ensure smooth gas accretion by removing a few particles per timestep. This is important since, in SPH, gas particles serve as interpolation points. By preventing the removal of a large number of particles, we ensure the stability of the hydrodynamic computations. + +Sink 2-body encounters +++++++++++++++++++++++ + +To accurately follow *sink mergers*, we implemented `Grudic et al. (2021) <https://academic.oup.com/mnras/article/506/2/2199/6276745>`_ two body timestep constraints between all sink particles: + +.. math:: + \Delta t_s \leq \Delta t_\text{2-body} = \frac{ t_\text{c, min} t_\text{dyn, min}}{t_\text{c, min} + t_\text{dyn, min}} \; , + +with + +.. math:: + \quad t_\text{c, min} = \min_{s \neq s'} \frac{ |\varphi^{-1}(r_{ss'}, \, H_\text{sink})| }{v_{ss'}} \quad \text{and} \quad t_\text{dyn, min} = \min_{s \neq s'} \sqrt{ \frac{ |\varphi'(r_{ss'}, \, H_\text{sink})|^{-1}} { G (m_s + m_{s'})} } \; , + +where :math:`r_{ss'}` is the sinks relative separation, :math:`v_{ss'}` the relative velocity, :math:`m_{s}` and :math:`m_{s'}` their masses and :math:`H_\text{sink}` is the sink (fixed) gravitational softening. The function :math:`\varphi(r, H)` is the potential corresponding to the Wendland C2 kernel density field (see `Schaller et al. (2024) <https://doi.org/10.1093/mnras/stae922>`_ section 4.1) and :math:`\varphi'(r, H) \equiv \frac{\mathrm{d} \varphi(r, H)}{\mathrm{d} r}` its derivative. + +Timesteps per sink's age categories ++++++++++++++++++++++++++++++++++++ + +We also implemented maximal timesteps sizes depending on the sink age; :math:`\Delta t_\text{max,s}^\text{age}`. A sink can be young, old or dead. In the first two cases, the sink's timestep is :math:`\min(\Delta t_\text{max,s}^\text{age}, \Delta t_s)`. In the last case, we impose :math:`\Delta t_\text{2-body}` only if a dead sink is involved in a two-boy encounter with an alive sink. Otherwise, the sink has no timestep constraint (apart from gravity). The parameters controlling the transition between the young and old is ``GEARSink:timestep_age_threshold``, and the one between old and dead is ``GEARSink:timestep_age_threshold_unlimited_Myr``. The maximal timesteps are given by ``GEARSink:max_timestep_young_Myr`` and ``GEARSink:max_timestep_old_Myr``. + +Notice that sink particles also satisfy a gravity timestep constraint, as do all gravitational particles in Swift. + +Star formation constraint ++++++++++++++++++++++++++ + +Although the stars' masses are sampled from the IMF, the stars' metallicities are not. If a sink accretes mass, it can create many star particles simultaneously. However, these stars will all have the same metallicity, which does not represent the actual metals' evolution during the accretion. In such a situation, the galaxies' properties are affected and do not represent the underlying physics. + +Another problem is that we can spawn many stars simultaneously, and the code may complain. Such a situation could be better. Although the constraints in the previous section will help, more is needed. Our solution is to introduce a new accretion criterion using the IMF properties. However, since our politics is that accretion should be feedback-regulated and not based on an arbitrary accretion rate, we reduce the sink time step to avoid limiting the star formation rate to an arbitrary value. + +The new accretion criterion is the following. The swallowed gas and sink mass does not exceed ``n_IMF`` times the IMF mass (see the IMF sampling section), but make sure to swallow at least one particle: :math:`M_\text{swallowed} \leq n_\text{IMF} M_\text{IMF} \text{ or } M_\text{swallowed} = 0`. + +Since we artificially restrict mass accretion, we keep track of the mass :math:`M_\text{eligible}` that would be swallowed without this criterion. Then, we compute the error :math:`\Delta M` between the restricted and unrestricted swallow. The absolute error is :math:`\Delta M = M_\text{swallowed} - M_\text{eligible}` and the relative error is :math:`| \Delta M | / M_\text{eligible}`. + +When :math:`\Delta M < 0` (i.e. :math:`M_\text{swallowed} \neq M_\text{eligible}`), we know the accretion was restricted and we can apply another time-step contraint. To compute a timestep, we convert :math:`\Delta M` to accretion rate by dividing by :math:`\Delta t_\text{s, estimated} = \min(\Delta t_\text{CFL}, \, \Delta t_\text{ff}, \Delta t_\text{2-body})`. Hence, we have the constraint: + +.. math:: + \Delta t_s \leq \Delta t_\text{SF} = \eta \cfrac{M_\text{eligible} \Delta t_\text{s, estimated}}{\Delta M} \text{ if } \Delta M < 0 \; , + +where :math:`\eta` is a tolerance parameter. This parameter corresponds to ``GEARSink:tolerance_SF_timestep`` in the code. + diff --git a/doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst b/doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst index 1414c92fb8302984f985e9f91d963a374e516aef..cbb8337a5bf3afc665f10feb76f2f9d5a26e9dd3 100644 --- a/doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst +++ b/doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst @@ -1,5 +1,5 @@ .. Sink particles in GEAR model - Darwin Roduit, 17 April 2024 + Darwin Roduit, 24 November 2024 .. _sink_GEAR_model_summary: @@ -117,6 +117,8 @@ The third criterion is mainly here to prevent two sink particles from forming at The tenth criterion prevents the formation of spurious sinks. Experiences have shown that removing gas within the accretion radius biases the hydro density estimates: the gas feels a force toward the sink. At some point, there is an equilibrium and gas particles accumulate at the edge of the accretion radius, which can then spawn sink particles that do not fall onto the primary sink and never merge. Moreover, the physical reason behind this criterion is that a sink represents a region of collapse. As a result, there is no need to have many sinks occupying the same space volume. They would compete for gas accretion without necessarily merging. This criterion is particularly meaningful in cosmological simulations to ensure proper sampling of the IMF. *This criterion can be disabled*. +Once a sink is formed, we record it birth time (or scale factor in cosmological runs). This information is used to put the sink into three categories: young, old and dead. If a sink is dead, it cannot accrete gas or sink anymore. However, a dead sink can still be swallowed by a young/old sink. Young and old sink only differ by their maximal allowed timestep. Details are provided in :ref:`sink_GEAR_timesteps`. + .. note:: However, notice that contrary to `Bate et al. <https://ui.adsabs.harvard.edu/abs/1995MNRAS.277..362B/abstract>`_, no boundary conditions for sink particles are introduced in the hydrodynamics calculations. @@ -129,12 +131,15 @@ Gas accretion Now that sink particles can populate the simulation, they need to swallow gas particles. To be accreted, gas particles need to pass a series of criteria. In the following, :math:`s` denotes a sink particle and :math:`i` is a gas particle. The criteria are the following: +#. The sink is not dead. If it is dead, it does not accrete gas. A sink is considered dead if it is older than ``timestep_age_threshold_unlimited_Myr``. #. If the gas falls within :math:`f_{\text{acc}} r_{\text{acc}}` (:math:`0 \leq f_{\text{acc}} \leq 1`), the gas is accreted without further check. #. In the region :math:`f_{\text{acc}} r_{\text{acc}} \leq |\mathbf{x}_i| \leq r_{\text{acc}}`, then, we check: #. The specific angular momentum is smaller than the one of a Keplerian orbit at :math:`r_{\text{acc}}`: :math:`|\mathbf{L}_{si}| \leq |\mathbf{L}_{\text{Kepler}}|`. #. The gas is gravitationally bound to the sink particle: :math:`E_{\text{tot}} < 0`. - #. Out of all pairs of sink-gas, the gas is the most bound to this one. This case is illustrated in the figure below. + #. The gas size is smaller or equal to the sink size: :math:`\gamma_k h_i \leq r_{\text{acc}}`. + #. Out of all pairs of sink-gas, the gas is the most bound to this one. This case is illustrated in the figure below. + #. The total swallowed mass does not exceed ``n_IMF`` times the IMF mass (see the IMF sampling section), but make sure to swallow at least one particle: :math:`M_\text{swallowed} \leq n_\text{IMF} M_\text{IMF} \text{ or } M_\text{swallowed} = 0`. The physical specific angular momenta and the total energy are given by: @@ -145,7 +150,11 @@ The physical specific angular momenta and the total energy are given by: .. note:: Here the potential is the softened potential of Swift. -Those criteria are similar to `Price et al. <https://ui.adsabs.harvard.edu/abs/2018PASA...35...31P/abstract>`_, with the addition of the internal energy. This term ensures that the gas is cold enough to be accreted. Its main purpose is to avoid gas accretion and star spawning in hot regions far from sink/star-forming regions, which can happen, e.g., if a sink leaves a galaxy. +Those criteria are similar to `Price et al. <https://ui.adsabs.harvard.edu/abs/2018PASA...35...31P/abstract>`_ and `Grudic et al. (2021) <https://academic.oup.com/mnras/article/506/2/2199/6276745>`_, with the addition of the internal energy. This term ensures that the gas is cold enough to be accreted. Its main purpose is to avoid gas accretion and star spawning in hot regions far from sink/star-forming regions, which can happen, e.g., if a sink leaves a galaxy. + +Let's comment on the fourth criterion, specific to our star formation scheme. This criterion restricts the swallowed mass to avoid spawning too many stars in a single time step. Swallowing too many gas particles in a time-step can lead to instabilities in the hydrodynamics, given that gas particles act as interpolation points. Also, creating many stars at once is prejudicial for two reasons. First, the stars' mass samples an IMF, but the star's metallicities do not. So, all stars end up with the same metal content. This situation does not reflect the history of metal accretion and will lead to poor galaxy properties. Second, we need to specify at runtime the maximal number of memory allocated for extra stars until the next tree rebuild. If we create more stars than this limit, the code will stop and send an error. + +Since our politics is not to arbitrarily restrict the accretion using some arbitrary mass accretion rate (in fact, the accretion must be feedback-regulated), we then lower the sink time step to swallow the remaining gas particles soon. So, instead of eating a considerable amount of mass and spawning many stars in a big time step, we swallow smaller amounts of gas/sink and create fewer stars in smaller time steps. Details about the how we reduce the timestep are given in :ref:`sink_GEAR_timesteps`. Once a gas is eligible for accretion, its properties are assigned to the sink. The sink accretes the *entire* gas particle mass and its properties are updated in the following way: @@ -167,12 +176,15 @@ Sink merging Sinks are allowed to merge if they enter one's accretion radius. We merge two sink particles if they respect a set of criteria. The criteria are similar to the gas particles, namely: +#. At least one of the sinks is not dead. A sink is considered dead if it is older than ``timestep_age_threshold_unlimited_Myr``. #. If one of the sinks falls within the other's inner accretion radius, :math:`f_{\text{acc}} r_{\text{acc}}` (:math:`0 \leq f_{\text{acc}} \leq 1`), the sinks are merged without further check. #. In the region :math:`f_{\text{acc}} r_{\text{acc}} \leq |\mathbf{x}_i| \leq r_{\text{acc}}`, then, we check: + #. The specific angular momentum is smaller than the one of a Keplerian orbit at :math:`r_{\text{acc}}`: :math:`|\mathbf{L}_{ss'}| \leq |\mathbf{L}_{\text{Kepler}}|`. #. One sink is gravitationally bound to the other: :math:`E_{\text{mec}, ss'} < 0` or :math:`E_{\text{mec}, s's} < 0`. + #. The total swallowed mass does not exceed ``n_IMF`` times the IMF mass (see the IMF sampling section), but make sure to swallow at least one particle: :math:`M_\text{swallowed} \leq n_\text{IMF} M_\text{IMF} \text{ or } M_\text{swallowed} = 0`. -We compute the total energies in the same manner as gas particles, with the difference that we do not use internal energy. Notice that we have two energies: each sink has a different potential energy since their mass can differ. +We compute the angular momenta and total energies in the same manner as gas particles, with the difference that we do not use internal energy. Notice that we have two energies: each sink has a different potential energy since their mass can differ. When sinks merge, the sink with the smallest mass merges with the sink with the largest. If the two sinks have the same mass, we check the sink ID number and add the smallest ID to the biggest one. @@ -183,9 +195,9 @@ IMF sampling :width: 400px :align: center :figclass: align-center - :alt: Initial mass function split into the continuous and discrete part. + :alt: Initial mass function split into the continuous and discrete part. - This figure shows an IMF split into two parts by :math:`m_t`: the continuous (orange) and the discrete (blue) part. + This figure shows an IMF split into two parts by :math:`m_t`: the continuous (orange) and the discrete (blue) part. The IMF mass is :math:`M_\text{IMF} = M_c + M_d`. Now remains one critical question: how are stars formed in this scheme? Simply, by sampling an IMF. In our scheme, population III stars and population II have two different IMFs. For the sake of simplicity, in the following presentation, we consider only the case of population II stars. However, this can be easily generalized to population III. @@ -197,7 +209,7 @@ Consider an IMF such as the one above. We split it into two parts at ``minimal_d As a result, we group all those low-mass stars in one stellar particle of mass ``stellar_particle_mass_Msun``. Such star particles are called "continuous", contrary to the "discrete" individual stars. With all that information, we can compute the number of stars in the continuous part of the IMF (called :math:`N_c`) and in the discrete part (called :math:`N_d`). Finally, we can compute the probabilities of each part, respectively called :math:`P_c` and :math:`P_d`. Notice that the mathematical derivation is given in the theory latex files. -Thus, the algorithm to sample the IMF and five the sink their ``target_mass`` is the following : +Thus, the algorithm to sample the IMF and determine the sink's ``target_mass`` is the following : * draw a random number :math:`\chi` from a uniform distribution in the interval :math:`(0 , \; 1 ]`; * if :math:`\chi < P_c`: ``sink.target_mass = stellar_particle_mass``; @@ -205,7 +217,9 @@ Thus, the algorithm to sample the IMF and five the sink their ``target_mass`` is We have assumed that we have a function ``sample_IMF_high()`` that correctly samples the IMF in the discrete part. -Now, what happens to the sink? After a first sink forms, we give it a target mass with the algorithm outlined above. The sink then swallows gas particles (see the task graph at the top of the page) and finally spawns stars. While the sink possesses enough mass, we can continue to choose a new target mass. When the sink does have enough mass, the algorithm stops for this timestep. The next timestep, the sink may accrete gas and spawn stars again. If the sink never reaches the target mass, then it cannot spawn stars. In practice, however, sink particles could accumulate enough pass to spawn individual (Pop III) stars with masses 240 :math:`M_\odot` and more! +Now, what happens to the sink? After a sink forms, we give it a target mass with the abovementioned algorithm. The sink then swallows gas particles (see the task graph at the top of the page) and spawns stars. *While the sink possesses enough mass*, we can continue to choose a new target mass. When the sink does have enough mass, the algorithm stops for this timestep. The next timestep, the sink may accrete gas and spawn stars again. The sink cannot spawn stars if it never reaches the target mass. In practice, sink particles could accumulate enough pass to spawn individual (Pop III) stars with masses 240 :math:`M_\odot` and more! + +For low-resolution simulations (:math:`m_\text{gas} > 100 \; M_\odot`), we also add a minimal sink mass constraint: the sink can spawn a star if ``m_sink > target_mass`` *and* ``m_sink - target_mass >= minimal_mass``. In low-resolution simulation, when a gas particle turns into a sink, the latter can have enough mass to spawn stars, depending on the sink stars and IMF parameters. As a result, the sink spawns the stars and then ends up with :math:`m_\text{sink} \ll m_\text{gas}`. Such a situation is detrimental for two reasons: 1) the sink mass is so low that gas can seldom be bound to it and thus stops spawning stars and 2) the sink can get kicked away by gravitational interactions due to the high mass difference. The parameter controlling the sink's minimal mass is ``GEARSink:sink_minimal_mass_Msun``. As explained at the beginning of this section, GEAR uses two IMFs for the population of II and III stars. The latter are called the first stars in the code. How does a sink decide which IMF to draw the target mass from? We define a threshold metallicity, ``GEARFeedback:imf_transition_metallicity`` that determines the first stars' maximal metallicity. When the sink particle's metallicity exceeds this threshold, it uses the population II IMF, defined in ``GEARFeedback:yields_table``. diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/README b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/README index 3a9be3bc2242fe6a8ec497fdf87615249dff90a0..452eb5ee93856176add2ad7100039a86166dbf5c 100644 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/README +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/README @@ -10,7 +10,7 @@ be downloaded when launching the run script (see below). To run this example with the GEAR model, SWIFT must be configured with the following options: -./configure --with-chemistry=GEAR_10 --with-cooling=grackle_0 --with-pressure-floor=GEAR --with-stars=GEAR --with-star-formation=GEAR --with-feedback=GEAR --with-grackle=${GRACKLE_ROOT} +./configure --with-subgrid=GEAR --with-grackle=${GRACKLE_ROOT} To start the simulation with the GEAR model: @@ -18,7 +18,7 @@ To start the simulation with the GEAR model: # Sink particles -To configure this example with GEAR model + GEAR sink particles, add --with-sink=GEAR. +Sink particles are already included in the option `--with-subgrid=GEAR`. To start the simulation with the sink particles, type: diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/params.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/params.yml index 02f942cc01b175ea08b60f7945933689d18e776e..69ce818141c0be0b74bd6c611b7d65915e1a38de 100644 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/params.yml +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/params.yml @@ -95,6 +95,14 @@ GEARSink: sink_formation_overlapping_sink_criterion: 1 # (Optional) Activate the overlapping sink check for sink formation. (Default: 1) disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0) + # Timesteps parameters + CFL_condition: 0.5 # Courant-Friedrich-Levy condition for time integration. + timestep_age_threshold_unlimited_Myr: 100. # (Optional) Age above which sinks have no time-step restriction any more (in Mega-years). Defaults to 0. + timestep_age_threshold_Myr: 25. # (Optional) Age at which sink switch from young to old for time-stepping purposes (in Mega-years). Defaults to FLT_MAX. + max_timestep_young_Myr: 1. # (Optional) Maximal time-step length of young sinks (in Mega-years). Defaults to FLT_MAX. + max_timestep_old_Myr: 5. # (Optional) Maximal time-step length of old sinks (in Mega-years). Defaults to FLT_MAX. + n_IMF: 2 # (Optional) Number of times the IMF mass can be swallowed in a single timestep. (Default: FLTM_MAX) + GEARPressureFloor: jeans_factor: 10 diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/makeDisk.py b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/makeDisk.py index cec23105e59edc4e1a6c551a80b722b4eb580d05..70f5c2d9c49293210f185b936b4d264ba82fd462 100755 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/makeDisk.py +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/makeDisk.py @@ -588,14 +588,14 @@ nb.rename("galaxy_multi_component.hdf5") nb.write() #%% -#Add the StellarParticleType attribute to the dataset +# Add the StellarParticleType attribute to the dataset import h5py as h5 import numpy as np nb_star = nb.select("stars") N_star = np.sum(nb_star.npart) -star_tpe = 2 # Single population stars -star_type = np.ones(N_star)*star_tpe +star_tpe = 2 # Single population stars +star_type = np.ones(N_star) * star_tpe -with h5.File("galaxy_multi_component.hdf5", "r+") as f: +with h5.File("galaxy_multi_component.hdf5", "r+") as f: f["PartType4"].create_dataset("StellarParticleType", data=star_type) diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml index 4a6647cc2a23c7e13c8237919f9ddb6c6f72f9f5..ce4671125474f8e6862056e65695676ec4418a09 100644 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml @@ -105,6 +105,14 @@ GEARSink: sink_formation_overlapping_sink_criterion: 1 # (Optional) Activate the overlapping sink check for sink formation. (Default: 1) disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0) + # Timesteps parameters + CFL_condition: 0.5 # Courant-Friedrich-Levy condition for time integration. + timestep_age_threshold_unlimited_Myr: 100. # (Optional) Age above which sinks have no time-step restriction any more (in Mega-years). Defaults to 0. + timestep_age_threshold_Myr: 25. # (Optional) Age at which sink switch from young to old for time-stepping purposes (in Mega-years). Defaults to FLT_MAX. + max_timestep_young_Myr: 0.5 # (Optional) Maximal time-step length of young sinks (in Mega-years). Defaults to FLT_MAX. + max_timestep_old_Myr: 1. # (Optional) Maximal time-step length of old sinks (in Mega-years). Defaults to FLT_MAX. + n_IMF: 2. # (Optional) Number of times the IMF mass can be swallowed in a single timestep. (Default: FLTM_MAX) + GEARFeedback: supernovae_energy_erg: 1e51 supernovae_efficiency: 0.1 diff --git a/examples/SinkParticles/HomogeneousBox/params.yml b/examples/SinkParticles/HomogeneousBox/params.yml index ccce508e35133218f4b1cea538c0693e07e11cb1..a805594d057141a4ed49c443f62b70a12f09a9b6 100755 --- a/examples/SinkParticles/HomogeneousBox/params.yml +++ b/examples/SinkParticles/HomogeneousBox/params.yml @@ -19,7 +19,7 @@ Gravity: # Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.) TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). - time_end: 0.282 #500 Myr # The end time of the simulation (in internal units). + time_end: 0.5 #0.282 #500 Myr # The end time of the simulation (in internal units). dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). dt_max: 1e-1 # The maximal time-step size of the simulation (in internal units). @@ -52,6 +52,7 @@ SPH: resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 57Ngbs with the Wendland C2 kernel). CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. h_max: 5 + h_min_ratio: 0.1 minimal_temperature: 1 @@ -71,7 +72,7 @@ GrackleCooling: maximal_density_Hpcm3: -1 # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter. GEARChemistry: - initial_metallicity: -5 + initial_metallicity: 1e-6 GEARFeedback: supernovae_energy_erg: 1e51 # supernovae energy, used only for SNIa @@ -85,15 +86,15 @@ GEARFeedback: GEARSink: cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units. - f_acc: 0.1 - temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. - density_threshold_Hpcm3: 1e1 # Minimum gas density (in g/cm3) required to form a sink particle. - maximal_density_threshold_Hpcm3: 1e5 # If the gas density exceeds this value (in g/cm3), a sink forms regardless of temperature if all other criteria are passed. + f_acc: 1e-2 + temperature_threshold_K: 1000 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. + density_threshold_Hpcm3: 1e1 # Minimum gas density (in g/cm3) required to form a sink particle. + maximal_density_threshold_Hpcm3: 1e5 # If the gas density exceeds this value (in g/cm3), a sink forms regardless of temperature if all other criteria are passed. stellar_particle_mass_Msun: 50 # Mass of the stellar particle representing the low mass stars, in solar mass minimal_discrete_mass_Msun: 8 # Minimal mass of stars represented by discrete particles, in solar mass stellar_particle_mass_first_stars_Msun: 50 # Mass of the stellar particle representing the low mass stars, in solar mass minimal_discrete_mass_first_stars_Msun: 8 # Minimal mass of stars represented by discrete particles, in solar mass - star_spawning_sigma_factor: 0.5 # Factor to rescale the velocity dispersion of the stars when they are spawned. (Default: 0.2) + star_spawning_sigma_factor: 0.5 # Factor to rescale the velocity dispersion of the stars when they are spawned. (Default: 0.2) sink_formation_contracting_gas_criterion: 1 # (Optional) Activate the contracting gas check for sink formation. (Default: 1) sink_formation_smoothing_length_criterion: 1 # (Optional) Activate the smoothing length check for sink formation. (Default: 1) sink_formation_jeans_instability_criterion: 1 # (Optional) Activate the two Jeans instability checks for sink formation. (Default: 1) @@ -101,3 +102,16 @@ GEARSink: sink_formation_overlapping_sink_criterion: 1 # (Optional) Activate the overlapping sink check for sink formation. (Default: 1) disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0) + # Timesteps parameters + CFL_condition: 0.5 # Courant-Friedrich-Levy condition for time integration. + timestep_age_threshold_unlimited_Myr: 100. # (Optional) Age above which sinks have no time-step restriction any more (in Mega-years). Defaults to 0. + timestep_age_threshold_Myr: 25. # (Optional) Age at which sink switch from young to old for time-stepping purposes (in Mega-years). Defaults to FLT_MAX. + max_timestep_young_Myr: 0.5 # (Optional) Maximal time-step length of young sinks (in Mega-years). Defaults to FLT_MAX. + max_timestep_old_Myr: 1. # (Optional) Maximal time-step length of old sinks (in Mega-years). Defaults to FLT_MAX. + n_IMF: 2. # (Optional) Number of times the IMF mass can be swallowed in a single timestep. (Default: FLTM_MAX) + +Stars: + timestep_age_threshold_unlimited_Myr: 30. # (Optional) Age above which sinks have no time-step restriction any more (in Mega-years). Defaults to 0. + timestep_age_threshold_Myr: 10. # (Optional) Age at which sink switch from young to old for time-stepping purposes (in Mega-years). Defaults to FLT_MAX. + max_timestep_young_Myr: 1 # (Optional) Maximal time-step length of young sinks (in Mega-years). Defaults to FLT_MAX. + max_timestep_old_Myr: 5 # (Optional) Maximal time-step length of old sinks (in Mega-years). Defaults to FLT_MAX. diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/README b/examples/SinkParticles/HomogeneousBoxSinkParticles/README new file mode 100644 index 0000000000000000000000000000000000000000..8669089494c529ed3e2fbe4c9ebc8c084574eeea --- /dev/null +++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/README @@ -0,0 +1,24 @@ +# Intro +This example is a modified version of the HomogeneousBox in SWIFT that allows to add sink particles to the box. By default, there is no gravity. + +The default level is 5 (N_particle = 2^(3*level)). It is quick. If you want to try the example with higher levels, we recommend using HPC facilities. + +*This example is mainly used for sink MPI debugging.* + +# Configure + +To run this example with GEAR model, + +./configure --with-chemistry=GEAR_10 --with-cooling=grackle_0 --with-stars=GEAR --with-star-formation=GEAR --with-feedback=GEAR --with-sink=GEAR --with-kernel=wendland-C2 --with-grackle=path/to/grackle --enable-debug --enable-debugging-checks + +and then + +make -j + +# ICs +The run.sh script calls `makeIC.py' script with default values. You can experiment by changing the ICs. Run `python3 makeIC.py --help` to get the list of parameters. + +# Run +Type `run.sh` (or `n_ranks=4 mpi_run.sh`), and let's go! + +You can provide parameters to the running scripts, have a look inside them. diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/getChemistryTable.sh b/examples/SinkParticles/HomogeneousBoxSinkParticles/getChemistryTable.sh new file mode 100755 index 0000000000000000000000000000000000000000..b10fd23964158ee7a38d352dbd0ddd9beb7bdd77 --- /dev/null +++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/getChemistryTable.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5 diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/getGrackleCoolingTable.sh b/examples/SinkParticles/HomogeneousBoxSinkParticles/getGrackleCoolingTable.sh new file mode 100755 index 0000000000000000000000000000000000000000..e3eb106240709c80151a48625567d2cd78e5f568 --- /dev/null +++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/getGrackleCoolingTable.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5 diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py b/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py new file mode 100755 index 0000000000000000000000000000000000000000..6ad3be507b4fd12d8bef3c0af4a1b674cb950a27 --- /dev/null +++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/makeIC.py @@ -0,0 +1,267 @@ +################################################################################ +# This file is part of SWIFT. +# Copyright (c) 2022 Yves Revaz (yves.revaz@epfl.ch) +# 2024 Darwin Roduit (darwin.roduit@epfl.ch) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. +# +################################################################################ + +import h5py +import numpy as np +import argparse +from astropy import units +from astropy import constants + + +class store_as_array(argparse._StoreAction): + """Provides numpy array as argparse arguments.""" + + def __call__(self, parser, namespace, values, option_string=None): + values = np.array(values) + return super().__call__(parser, namespace, values, option_string) + + +def parse_options(): + + usage = "usage: %prog [options] file" + parser = argparse.ArgumentParser(description=usage) + + parser.add_argument( + "--lJ", + action="store", + dest="lJ", + type=float, + default=0.250, + help="Jeans wavelength in box size unit", + ) + + parser.add_argument( + "--rho", + action="store", + dest="rho", + type=float, + default=0.1, + help="Mean gas density in atom/cm3", + ) + + parser.add_argument( + "--mass", + action="store", + dest="mass", + type=float, + default=50, + help="Gas particle mass in solar mass", + ) + + parser.add_argument( + "--level", + action="store", + dest="level", + type=int, + default=5, + help="Resolution level: N = (2**l)**3", + ) + + parser.add_argument( + "--sink_mass", + action="store", + type=float, + default=50, + help="Sink particles mass in solar mass", + ) + + parser.add_argument( + "--sinks_vel", + action=store_as_array, + nargs=3, + type=float, + default=np.array([10, 0, 0]), + help="Sink particle velocity. All sinks get the same velocity", + ) + + parser.add_argument( + "--sink_pos", + action=store_as_array, + nargs=3, + type=float, + default=np.array([0, 0, 0]), + help="Sink particle position. Only use it to place one sink particle.", + ) + + parser.add_argument( + "--n_sink", + action="store", + type=int, + default=10, + help="Number of sink particles in the box", + ) + + parser.add_argument( + "-o", + action="store", + dest="outputfilename", + type=str, + default="box.hdf5", + help="output filename", + ) + + options = parser.parse_args() + return options + + +######################################## +# main +######################################## + +opt = parse_options() + +# define standard units +UnitMass_in_cgs = 1.988409870698051e43 # 10^10 M_sun in grams +UnitLength_in_cgs = 3.0856775814913673e21 # kpc in centimeters +UnitVelocity_in_cgs = 1e5 # km/s in centimeters per second +UnitCurrent_in_cgs = 1 # Amperes +UnitTemp_in_cgs = 1 # Kelvin +UnitTime_in_cgs = UnitLength_in_cgs / UnitVelocity_in_cgs + +UnitMass = UnitMass_in_cgs * units.g +UnitLength = UnitLength_in_cgs * units.cm +UnitTime = UnitTime_in_cgs * units.s +UnitVelocity = UnitVelocity_in_cgs * units.cm / units.s + +np.random.seed(1) + +# Number of particles +N = (2 ** opt.level) ** 3 # number of particles + +# Mean density +rho = opt.rho # atom/cc +rho = rho * constants.m_p / units.cm ** 3 + +# Gas particle mass +m = opt.mass # in solar mass +m = m * units.Msun + +# Gas mass in the box +M = N * m + +# Size of the box +L = (M / rho) ** (1 / 3.0) + +# Jeans wavelength in box size unit +lJ = opt.lJ +lJ = lJ * L + +# Gravitational constant +G = constants.G + +# Jeans wave number +kJ = 2 * np.pi / lJ +# Velocity dispersion +sigma = np.sqrt(4 * np.pi * G * rho) / kJ + + +print("Boxsize : {}".format(L.to(units.kpc))) +print("Number of particles : {}".format(N)) +print("Equivalent velocity dispertion : {}".format(sigma.to(units.m / units.s))) + +# Convert to code units +m = m.to(UnitMass).value +L = L.to(UnitLength).value +rho = rho.to(UnitMass / UnitLength ** 3).value +sigma = sigma.to(UnitVelocity).value + +# Generate the particles +pos = np.random.random([N, 3]) * np.array([L, L, L]) +vel = np.zeros([N, 3]) +mass = np.ones(N) * m +u = np.ones(N) * sigma ** 2 +ids = np.arange(N) +h = np.ones(N) * 3 * L / N ** (1 / 3.0) +rho = np.ones(N) * rho + +print("Inter-particle distance (code unit) : {}".format(L / N ** (1 / 3.0))) + + +##################### +# Now, take care of the sink +##################### +N_sink = opt.n_sink +m_sink = opt.sink_mass * units.M_sun +m_sink = m_sink.to(UnitMass).value # Convert the sink mass to internal units + + +if N_sink == 1: + pos_sink = np.reshape(opt.sinks_vel, (N_sink, 3)) +else: + pos_sink = np.random.random([N_sink, 3]) * np.array([L, L, L]) + +if opt.sinks_vel is not None: + vel_sink = np.tile( + opt.sinks_vel, (N_sink, 1) + ) # Duplicate the velocity for all sinks +else: + np.zeros([N_sink, 3]) + +mass_sink = np.ones(N_sink) * m_sink +ids_sink = np.arange(N, N + N_sink) + +##################### +# Finally write the ICs in the file +##################### + +# File +fileOutput = h5py.File(opt.outputfilename, "w") +print("{} saved.".format(opt.outputfilename)) + +# Header +grp = fileOutput.create_group("/Header") +grp.attrs["BoxSize"] = [L, L, L] +grp.attrs["NumPart_Total"] = [N, 0, 0, N_sink, 0, 0] +grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] +grp.attrs["NumPart_ThisFile"] = [N, 0, 0, N_sink, 0, 0] +grp.attrs["Time"] = 0.0 +grp.attrs["NumFileOutputsPerSnapshot"] = 1 +grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0] +grp.attrs["Dimension"] = 3 + + +# Units +grp = fileOutput.create_group("/Units") +grp.attrs["Unit length in cgs (U_L)"] = UnitLength_in_cgs +grp.attrs["Unit mass in cgs (U_M)"] = UnitMass_in_cgs +grp.attrs["Unit time in cgs (U_t)"] = UnitTime_in_cgs +grp.attrs["Unit current in cgs (U_I)"] = UnitCurrent_in_cgs +grp.attrs["Unit temperature in cgs (U_T)"] = UnitTemp_in_cgs + + +# Write Gas particle group +grp = fileOutput.create_group("/PartType0") +grp.create_dataset("Coordinates", data=pos, dtype="d") +grp.create_dataset("Velocities", data=vel, dtype="f") +grp.create_dataset("Masses", data=mass, dtype="f") +grp.create_dataset("SmoothingLength", data=h, dtype="f") +grp.create_dataset("InternalEnergy", data=u, dtype="f") +grp.create_dataset("ParticleIDs", data=ids, dtype="L") +grp.create_dataset("Densities", data=rho, dtype="f") + + +# Write sink particle group +grp = fileOutput.create_group("/PartType3") +grp.create_dataset("Coordinates", data=pos_sink, dtype="d") +grp.create_dataset("Velocities", data=vel_sink, dtype="f") +grp.create_dataset("Masses", data=mass_sink, dtype="f") +grp.create_dataset("ParticleIDs", data=ids_sink, dtype="L") +fileOutput.close() diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/mpi_run.sh b/examples/SinkParticles/HomogeneousBoxSinkParticles/mpi_run.sh new file mode 100755 index 0000000000000000000000000000000000000000..8b884691d68c92902a68b62826c4a7dcfaa9dc91 --- /dev/null +++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/mpi_run.sh @@ -0,0 +1,51 @@ +#!/bin/bash + +# make run.sh fail if a subcommand fails +set -e + +n_ranks=${n_ranks:=2} #Number of MPI ranks +n_threads=${n_threads:=1} #Number of threads to use +level=${level:=5} #Number of particles = 2^(3*level) +jeans_length=${jeans_length:=0.250} #Jeans wavelenght in unit of the boxsize +n_sinks=${n_sinks:=10} #Number of sinks + +# Remove the ICs +if [ -e ICs_homogeneous_box.hdf5 ] +then + rm ICs_homogeneous_box.hdf5 +fi + +#Create the ICs if they do not exist +if [ ! -e ICs_homogeneous_box.hdf5 ] +then + echo "Generating initial conditions to run the example..." + python3 makeIC.py --level $level -o ICs_homogeneous_box.hdf5 --lJ $jeans_length --n_sink $n_sinks --sink_pos 0 0 0 --sinks_vel 10 10 0 +fi + +# Get the Grackle cooling table +if [ ! -e CloudyData_UVB=HM2012.h5 ] +then + echo "Fetching the Cloudy tables required by Grackle..." + ./getGrackleCoolingTable.sh +fi + +if [ ! -e POPIIsw.h5 ] +then + echo "Fetching the chemistry tables..." + ./getChemistryTable.sh +fi + +# Create output directory +DIR=snap #First test of units conversion +if [ -d "$DIR" ]; +then + echo "$DIR directory exists. Its content will be removed." + rm -r $DIR +else + echo "$DIR directory does not exists. It will be created." + mkdir $DIR +fi + +printf "Running simulation..." + +mpirun -n $n_ranks ../../../swift_mpi --pin --hydro --sinks --stars --external-gravity --feedback --threads=$n_threads --verbose 1 params.yml 2>&1 | tee output.log diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml b/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml new file mode 100755 index 0000000000000000000000000000000000000000..90bd18a419bad00d3191df710d337c351b058d6a --- /dev/null +++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml @@ -0,0 +1,104 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.988409870698051e+43 # 10^10 solar masses + UnitLength_in_cgs: 3.0856775814913673e+21 # 1 kpc + UnitVelocity_in_cgs: 1e5 # km/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Parameters for the self-gravity scheme +Gravity: + MAC: adaptive # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'. + epsilon_fmm: 0.001 # Tolerance parameter for the adaptive multipole acceptance criterion. + theta_cr: 0.7 # Opening angle for the purely gemoetric criterion. + eta: 0.025 # Constant dimensionless multiplier for time integration. + theta: 0.7 # Opening angle (Multipole acceptance criterion). + max_physical_baryon_softening: 0.005 # Physical softening length (in internal units) + mesh_side_length: 32 + +# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.) +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 0.282 # The end time of the simulation (in internal units). + dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). + dt_max: 2e-2 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + subdir: snap + basename: snapshot # Common part of the name of output files + time_first: 0. #230 Myr # (Optional) Time of the first output if non-cosmological time-integration (in internal units) + delta_time: 10e-3 # Time difference between consecutive outputs (in internal units) + +Scheduler: + cell_extra_gparts: 100 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation. + cell_extra_sinks: 100 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation. + cell_extra_sparts: 100 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation. + max_top_level_cells: 3 # + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 5e-3 # Time between statistics output + time_first: 0. # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units) + +# Parameters related to the initial conditions +InitialConditions: + file_name: ICs_homogeneous_box.hdf5 + periodic: 1 # Are we running with periodic ICs? + shift: [0.0,0.0,0.0] + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 57Ngbs with the Wendland C2 kernel). + CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. + h_max: 5 + minimal_temperature: 1 + + +# Cooling with Grackle 3.0 +GrackleCooling: + cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository) + with_UV_background: 0 # Enable or not the UV background + redshift: -1 # Redshift to use (-1 means time based redshift) + with_metal_cooling: 1 # Enable or not the metal cooling + provide_volumetric_heating_rates: 0 # User provide volumetric heating rates + provide_specific_heating_rates: 0 # User provide specific heating rates + self_shielding_method: -1 # Grackle (<= 3) or Gear self shielding method + self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding + max_steps: 1000 + convergence_limit: 1e-2 + thermal_time_myr: 5 + maximal_density_Hpcm3: -1 # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter. + +GEARChemistry: + initial_metallicity: 0 + +GEARFeedback: + supernovae_energy_erg: 1e51 # supernovae energy, used only for SNIa + supernovae_efficiency: 0.1 # supernovae energy efficiency, used for both SNIa and SNII + yields_table: POPIIsw.h5 + yields_table_first_stars: POPIIsw.h5 + discrete_yields: 1 + imf_transition_metallicity: -5 # Maximal metallicity ([Fe/H]) for a first star (0 to deactivate). + elements: [Fe, Mg, O, C, Al, Ca, Ba, Zn, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N). + +GEARSink: + cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units. + f_acc: 0.1 + temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_g_per_cm3 <= density <= maximal_density_threshold_g_per_cm3. + density_threshold_Hpcm3: 1e1 # Minimum gas density (in g/cm3) required to form a sink particle. + maximal_density_threshold_Hpcm3: 1e5 # If the gas density exceeds this value (in g/cm3), a sink forms regardless of temperature if all other criteria are passed. + stellar_particle_mass_Msun: 50 # Mass of the stellar particle representing the low mass stars, in solar mass + minimal_discrete_mass_Msun: 8 # Minimal mass of stars represented by discrete particles, in solar mass + stellar_particle_mass_first_stars_Msun: 50 # Mass of the stellar particle representing the low mass stars, in solar mass + minimal_discrete_mass_first_stars_Msun: 8 # Minimal mass of stars represented by discrete particles, in solar mass + star_spawning_sigma_factor: 0.5 # Factor to rescale the velocity dispersion of the stars when they are spawned. (Default: 0.2) + sink_formation_contracting_gas_criterion: 1 # (Optional) Activate the contracting gas check for sink formation. (Default: 1) + sink_formation_smoothing_length_criterion: 1 # (Optional) Activate the smoothing length check for sink formation. (Default: 1) + sink_formation_jeans_instability_criterion: 1 # (Optional) Activate the two Jeans instability checks for sink formation. (Default: 1) + sink_formation_bound_state_criterion: 1 # (Optional) Activate the bound state check for sink formation. (Default: 1) + sink_formation_overlapping_sink_criterion: 1 # (Optional) Activate the overlapping sink check for sink formation. (Default: 1) + disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0) + + # Timesteps parameters + CFL_condition: 0.5 # Courant-Friedrich-Levy condition for time integration. diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/run.sh b/examples/SinkParticles/HomogeneousBoxSinkParticles/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..9db866d75f792434d5fe2d93f9e6668299737dd2 --- /dev/null +++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/run.sh @@ -0,0 +1,55 @@ +#!/bin/bash + +# make run.sh fail if a subcommand fails +set -e + +n_threads=${n_threads:=8} #Number of threads to use +level=${level:=5} #Number of particles = 2^(3*level) +jeans_length=${jeans_length:=0.250} #Jeans wavelenght in unit of the boxsize +gas_density=${gas_density:=0.1} #Gas density in atom/cm^3 +gas_particle_mass=${gas_particle_mass:=50} #Mass of the gas particles +n_sinks=${n_sinks:=10} #Number of sinks + +# Remove the ICs +if [ -e ICs_homogeneous_box.hdf5 ] +then + rm ICs_homogeneous_box.hdf5 +fi + +#Create the ICs if they do not exist +if [ ! -e ICs_homogeneous_box.hdf5 ] +then + echo "Generating initial conditions to run the example..." + python3 makeIC.py --level $level -o ICs_homogeneous_box.hdf5 --lJ $jeans_length --n_sink $n_sinks --sink_pos 0 0 0 --sinks_vel 10 10 0 +fi + + +# Get the Grackle cooling table +if [ ! -e CloudyData_UVB=HM2012.h5 ] +then + echo "Fetching the Cloudy tables required by Grackle..." + ./getGrackleCoolingTable.sh +fi + + +if [ ! -e POPIIsw.h5 ] +then + echo "Fetching the chemistry tables..." + ./getChemistryTable.sh +fi + + +# Create output directory +DIR=snap #First test of units conversion +if [ -d "$DIR" ]; +then + echo "$DIR directory exists. Its content will be removed." + rm -r $DIR +else + echo "$DIR directory does not exists. It will be created." + mkdir $DIR +fi + +printf "Running simulation..." + +../../../swift --hydro --sinks --stars --external-gravity --feedback --threads=$n_threads params.yml 2>&1 | tee output.log diff --git a/examples/SinkParticles/PlummerSphere/params.yml b/examples/SinkParticles/PlummerSphere/params.yml index fe3e28d10b8dfe542ee53cef34cc2f07549eac05..436773e49e99d9f51caad8e8549e884e7e15b1d9 100755 --- a/examples/SinkParticles/PlummerSphere/params.yml +++ b/examples/SinkParticles/PlummerSphere/params.yml @@ -100,3 +100,4 @@ GEARSink: sink_formation_bound_state_criterion: 0 # (Optional) Activate the bound state check for sink formation. (Default: 1) sink_formation_overlapping_sink_criterion: 0 # (Optional) Activate the overlapping sink check for sink formation. (Default: 1) disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0) + CFL_condition: 0.5 # Courant-Friedrich-Levy condition for time integration. diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index f4bf5b38233d1fc5ff053a99e385addf1011624e..5feced7e25074b95b87490ea72d3ce6b08a34d15 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -878,11 +878,12 @@ DefaultSink: # GEAR sink particles GEARSink: - cut_off_radius: 1e-3 # Cut off radius of the sink particles (in internal units). This parameter should be adapted with the resolution. + cut_off_radius: 1e-3 # Cut off radius of the sink particles (in internal units). This parameter should be adapted with the resolution. f_acc: 0.1 # (Optional) Fraction of the cut_off_radius that determines if a gas particle should be swallowed wihtout additional check. It has to respect 0 <= f_acc <= 1. (Default: 0.1) - temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. - density_threshold_Hpcm3: 1e3 # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle. - maximal_density_threshold_Hpcm3: 1e5 # If the gas density exceeds this value (in Hydrogen atoms/cm3), a sink forms regardless of temperature if all other criteria are passed. (Default: FLT_MAX) + temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. + density_threshold_Hpcm3: 1e3 # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle. + maximal_density_threshold_Hpcm3: 1e5 # If the gas density exceeds this value (in Hydrogen atoms/cm3), a sink forms regardless of temperature if all other criteria are passed. (Default: FLT_MAX) + sink_minimal_mass_Msun: 0. # (Optional) Sink minimal mass in Msun. This parameter prevents m_sink << m_gas in low resolution simulations. (Default: 0.0) stellar_particle_mass_Msun: 20 # Mass of the stellar particle representing the low mass stars (continuous IMF sampling) (in solar mass) minimal_discrete_mass_Msun: 8 # Minimal mass of stars represented by discrete particles (in solar mass) stellar_particle_mass_first_stars_Msun: 20 # Mass of the stellar particle representing the low mass stars (continuous IMF sampling) (in solar mass). First stars @@ -893,7 +894,14 @@ GEARSink: sink_formation_jeans_instability_criterion: 1 # (Optional) Activate the two Jeans instability checks for sink formation. (Default: 1) sink_formation_bound_state_criterion: 1 # (Optional) Activate the bound state check for sink formation. (Default: 1) sink_formation_overlapping_sink_criterion: 1 # (Optional) Activate the overlapping sink check for sink formation. (Default: 1) - disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0) + disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0) + CFL_condition: 0.5 # Courant-Friedrich-Levy condition for time integration. + timestep_age_threshold_unlimited_Myr: 100. # (Optional) Age above which sinks no longer have time-step restrictions, except for 2-body encounters involving another young/old sink and gravity (in Mega-years). (Default: FLT_MAX) + timestep_age_threshold_Myr: 25. # (Optional) Age at which sink switch from young to old for time-stepping purposes (in Mega-years). (Default: FLT_MAX) + max_timestep_young_Myr: 1.0 # (Optional) Maximal time-step length of young sinks (in Mega-years). (Default: FLT_MAX) + max_timestep_old_Myr: 5.0 # (Optional) Maximal time-step length of old sinks (in Mega-years). (Default: FLT_MAX) + n_IMF: 2 # (Optional) Number of times the IMF mass can be swallowed in a single timestep. (Default: FLTM_MAX) + tolerance_SF_timestep: 0.5 # (Optional) Tolerance parameter for SF timestep constraint. (Default: 0.5) # Parameters related to radiative transfer --------------------------------------- diff --git a/src/Makefile.am b/src/Makefile.am index 429485fed58b4245f57c4142c6d5ae6fb74f6ba4..3d2a06d80fad524373762fd13be2392b7af6014e 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -492,7 +492,8 @@ nobase_noinst_HEADERS += pressure_floor/GEAR/pressure_floor_struct.h pressure_fl nobase_noinst_HEADERS += pressure_floor/GEAR/pressure_floor_debug.h pressure_floor/none/pressure_floor_debug.h nobase_noinst_HEADERS += sink/Default/sink.h sink/Default/sink_io.h sink/Default/sink_part.h sink/Default/sink_properties.h nobase_noinst_HEADERS += sink/Default/sink_iact.h sink/Default/sink_struct.h sink/Default/sink_debug.h -nobase_noinst_HEADERS += sink/GEAR/sink.h sink/GEAR/sink_io.h sink/GEAR/sink_part.h sink/GEAR/sink_properties.h +nobase_noinst_HEADERS += sink/GEAR/sink.h sink/GEAR/sink_io.h sink/GEAR/sink_part.h +nobase_noinst_HEADERS += sink/GEAR/sink_properties.h sink/GEAR/sink_setters.h sink/GEAR/sink_getters.h nobase_noinst_HEADERS += sink/GEAR/sink_iact.h sink/GEAR/sink_struct.h sink/GEAR/sink_debug.h nobase_noinst_HEADERS += neutrino.h neutrino_properties.h neutrino_io.h nobase_noinst_HEADERS += neutrino/Default/neutrino.h neutrino/Default/relativity.h neutrino/Default/fermi_dirac.h diff --git a/src/cell.c b/src/cell.c index 4f58054830f31a43dd96f8539262af0fc0b5ceb6..faadb774325f5bdf38ddc1c43742ed261c84a45f 100644 --- a/src/cell.c +++ b/src/cell.c @@ -293,6 +293,39 @@ int cell_link_bparts(struct cell *c, struct bpart *bparts) { return c->black_holes.count; } +/** + * @brief Link the cells recursively to the given #sink array. + * + * @param c The #cell. + * @param sinks The #sink array. + * + * @return The number of particles linked. + */ +int cell_link_sinks(struct cell *c, struct sink *sinks) { +#ifdef SWIFT_DEBUG_CHECKS + if (c->nodeID == engine_rank) + error("Linking foreign particles in a local cell!"); + + if (c->sinks.parts != NULL) + error("Linking sparts into a cell that was already linked"); +#endif + + c->sinks.parts = sinks; + c->sinks.parts_rebuild = sinks; + + /* Fill the progeny recursively, depth-first. */ + if (c->split) { + int offset = 0; + for (int k = 0; k < 8; k++) { + if (c->progeny[k] != NULL) + offset += cell_link_sinks(c->progeny[k], &sinks[offset]); + } + } + + /* Return the total number of linked particles. */ + return c->sinks.count; +} + /** * @brief Recurse down foreign cells until reaching one with hydro * tasks; then trigger the linking of the #part array from that @@ -413,6 +446,7 @@ void cell_unlink_foreign_particles(struct cell *c) { c->hydro.parts = NULL; c->stars.parts = NULL; c->black_holes.parts = NULL; + c->sinks.parts = NULL; if (c->split) { for (int k = 0; k < 8; k++) { diff --git a/src/cell.h b/src/cell.h index 2bb2d6976964461d797069a24c9bad9f9e71cdf3..9cb04b32c4a81affa50ea6584c1679ca756414dc 100644 --- a/src/cell.h +++ b/src/cell.h @@ -280,6 +280,15 @@ struct pcell_step { float dx_max_part; } black_holes; + struct { + + /*! Minimal integer end-of-timestep in this cell (sinks) */ + integertime_t ti_end_min; + + /*! Maximal distance any #part has travelled since last rebuild */ + float dx_max_part; + } sinks; + struct { /*! Minimal integer end-of-timestep in this cell (rt) */ @@ -562,6 +571,7 @@ int cell_link_parts(struct cell *c, struct part *parts); int cell_link_gparts(struct cell *c, struct gpart *gparts); int cell_link_sparts(struct cell *c, struct spart *sparts); int cell_link_bparts(struct cell *c, struct bpart *bparts); +int cell_link_sinks(struct cell *c, struct sink *sinks); int cell_link_foreign_parts(struct cell *c, struct part *parts); int cell_link_foreign_gparts(struct cell *c, struct gpart *gparts); void cell_unlink_foreign_particles(struct cell *c); diff --git a/src/cell_pack.c b/src/cell_pack.c index 21bf9ad122c788f228e811a3dd7024d39fa5cafb..83f8c4cab9d9c9de3cddfc906f789c3e0d2b2384 100644 --- a/src/cell_pack.c +++ b/src/cell_pack.c @@ -303,6 +303,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c, temp->hydro.count = 0; temp->grav.count = 0; temp->stars.count = 0; + temp->sinks.count = 0; temp->loc[0] = c->loc[0]; temp->loc[1] = c->loc[1]; temp->loc[2] = c->loc[2]; @@ -319,6 +320,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c, temp->hydro.dx_max_sort = 0.f; temp->stars.dx_max_part = 0.f; temp->stars.dx_max_sort = 0.f; + temp->sinks.dx_max_part = 0.f; temp->black_holes.dx_max_part = 0.f; temp->nodeID = c->nodeID; temp->parent = c; @@ -457,6 +459,9 @@ int cell_pack_end_step(const struct cell *c, struct pcell_step *pcells) { pcells[0].black_holes.ti_end_min = c->black_holes.ti_end_min; pcells[0].black_holes.dx_max_part = c->black_holes.dx_max_part; + pcells[0].sinks.ti_end_min = c->sinks.ti_end_min; + pcells[0].sinks.dx_max_part = c->sinks.dx_max_part; + /* Fill in the progeny, depth-first recursion. */ int count = 1; for (int k = 0; k < 8; k++) @@ -501,6 +506,9 @@ int cell_unpack_end_step(struct cell *c, const struct pcell_step *pcells) { c->black_holes.ti_end_min = pcells[0].black_holes.ti_end_min; c->black_holes.dx_max_part = pcells[0].black_holes.dx_max_part; + c->sinks.ti_end_min = pcells[0].sinks.ti_end_min; + c->sinks.dx_max_part = pcells[0].sinks.dx_max_part; + /* Fill in the progeny, depth-first recursion. */ int count = 1; for (int k = 0; k < 8; k++) diff --git a/src/cell_sinks.h b/src/cell_sinks.h index 6bd2f9a9afed9847001a9a106b6f2dac9a1b433d..9043ad191552657c7183b35ec7d443207f4af384 100644 --- a/src/cell_sinks.h +++ b/src/cell_sinks.h @@ -42,6 +42,9 @@ struct cell_sinks { /*! Pointer to the #sink data. */ struct sink *parts; + /*! Pointer to the #spart data at rebuild time. */ + struct sink *parts_rebuild; + /*! Linked list of the tasks computing this cell's sink swallow. */ struct link *swallow; diff --git a/src/cell_split.c b/src/cell_split.c index ee669ae17ac05eb8de1c8724c4656d8946d04c33..df0fc2572e0b5fa1c33164c07b7b1c7453a3804e 100644 --- a/src/cell_split.c +++ b/src/cell_split.c @@ -384,6 +384,7 @@ void cell_split(struct cell *c, const ptrdiff_t parts_offset, c->progeny[k]->sinks.count = bucket_count[k]; c->progeny[k]->sinks.count_total = c->progeny[k]->sinks.count; c->progeny[k]->sinks.parts = &c->sinks.parts[bucket_offset[k]]; + c->progeny[k]->sinks.parts_rebuild = c->progeny[k]->sinks.parts; } /* Finally, do the same song and dance for the gparts. */ diff --git a/src/cell_unskip.c b/src/cell_unskip.c index 3ce7d791a712e5dcacbc7794af08e81513fa5fdd..adfd6eb6ef1512c845a46674681cc6c1edea0bc3 100644 --- a/src/cell_unskip.c +++ b/src/cell_unskip.c @@ -3001,7 +3001,7 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) { if (cell_need_rebuild_for_sinks_pair(ci, cj)) rebuild = 1; if (cell_need_rebuild_for_sinks_pair(cj, ci)) rebuild = 1; -#ifdef WITH_MPI +#if defined(WITH_MPI) && !defined(SWIFT_DEBUG_CHECKS) error("TODO"); #endif } diff --git a/src/cooling/grackle/cooling.c b/src/cooling/grackle/cooling.c index 205a7b4f96ca663fbf0ee540febb38d4f5479488..be7d417259654ccdf370fb644de1db012f5f683a 100644 --- a/src/cooling/grackle/cooling.c +++ b/src/cooling/grackle/cooling.c @@ -68,6 +68,18 @@ double cooling_get_physical_density( const struct part* p, const struct cosmology* cosmo, const struct cooling_function_data* cooling); +/** + * @brief Record the time when cooling was switched off for a particle. + * + * @param p #part data. + * @param xp Pointer to the #xpart data. + * @param time The time when the cooling was switched off. + */ +INLINE void cooling_set_part_time_cooling_off(struct part* p, struct xpart* xp, + const double time) { + xp->cooling_data.time_last_event = time; +} + /** * @brief Common operations performed on the cooling function at a * given time-step or redshift. diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h index 92f6c921634bb62e58626effbcdc1019b2ecc75f..e5849a421eff0763217b09c8bcf15734471a1b94 100644 --- a/src/cooling/grackle/cooling.h +++ b/src/cooling/grackle/cooling.h @@ -105,6 +105,16 @@ INLINE static float cooling_get_subgrid_density(const struct part* p, return -1.f; } +/** + * @brief Record the time when cooling was switched off for a particle. + * + * @param p #part data. + * @param xp Pointer to the #xpart data. + * @param time The time when the cooling was switched off. + */ +void cooling_set_part_time_cooling_off(struct part* p, struct xpart* xp, + const double time); + float cooling_get_radiated_energy(const struct xpart* restrict xp); void cooling_print_backend(const struct cooling_function_data* cooling); diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h index dbb5e9872bf4e197f372c3f4a121c1beaaa775f9..d7955c4e4b1902796b855012e3f454f69e0eb70e 100644 --- a/src/cooling/none/cooling.h +++ b/src/cooling/none/cooling.h @@ -277,6 +277,19 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy( return 0.f; } +/** + * @brief Record the time when cooling was switched off for a particle. + * + * In none model, this function does nothing. + * + * @param p #part data. + * @param xp Pointer to the #xpart data. + * @param time The time when the cooling was switched off. + */ +INLINE static void cooling_set_part_time_cooling_off(struct part* p, + struct xpart* xp, + const double time) {} + /** * @brief Split the coolong content of a particle into n pieces * diff --git a/src/cosmology.h b/src/cosmology.h index a5fa2f32e41dfd5c4c9d1c141f83ecad489302d7..ef5adf491f9c86f4cd05b45d66f2034801e3c680 100644 --- a/src/cosmology.h +++ b/src/cosmology.h @@ -113,7 +113,7 @@ struct cosmology { /*! Scale-factor at the previous time-step */ double a_old; - /*! Redshit at the previous time-step */ + /*! Redshift at the previous time-step */ double z_old; /*------------------------------------------------------------------ */ diff --git a/src/engine.c b/src/engine.c index 3b992186dcdf2e9907accda1598cc084f4156603..1707dcf6300f1224e9fc8ccac7b769206653c64f 100644 --- a/src/engine.c +++ b/src/engine.c @@ -781,13 +781,14 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) { e->policy & (engine_policy_hydro | engine_policy_grid_hydro); const int with_stars = e->policy & engine_policy_stars; const int with_black_holes = e->policy & engine_policy_black_holes; + const int with_sinks = e->policy & engine_policy_sinks; struct space *s = e->s; ticks tic = getticks(); /* Count the number of particles we need to import and re-allocate the buffer if needed. */ size_t count_parts_in = 0, count_gparts_in = 0, count_sparts_in = 0, - count_bparts_in = 0; + count_bparts_in = 0, count_sinks_in = 0; for (int k = 0; k < nr_proxies; k++) { for (int j = 0; j < e->proxies[k].nr_cells_in; j++) { @@ -806,6 +807,11 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) { /* For black holes, we just use the numbers in the top-level cells */ count_bparts_in += e->proxies[k].cells_in[j]->black_holes.count; + + /* For sinks, we just use the numbers in the top-level cells + some + extra space */ + count_sinks_in += + e->proxies[k].cells_in[j]->sinks.count + space_extra_sinks; } } @@ -871,28 +877,49 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) { error("Failed to allocate foreign bpart data."); } + /* Allocate space for the foreign particles we will receive */ + size_t old_size_sinks_foreign = s->size_sinks_foreign; + if (!fof && count_sinks_in > s->size_sinks_foreign) { + if (s->sinks_foreign != NULL) swift_free("sinks_foreign", s->sinks_foreign); + s->size_sinks_foreign = engine_foreign_alloc_margin * count_sinks_in; + if (swift_memalign("sinks_foreign", (void **)&s->sinks_foreign, sink_align, + sizeof(struct sink) * s->size_sinks_foreign) != 0) + error("Failed to allocate foreign sink data."); + bzero(s->sinks_foreign, s->size_sinks_foreign * sizeof(struct sink)); + + /* Note: If you ever see a sink particle with id = -666, the following + lines is the ones that sets the ID to this value. */ + for (size_t i = 0; i < s->size_sinks_foreign; ++i) { + s->sinks_foreign[i].time_bin = time_bin_not_created; + s->sinks_foreign[i].id = -666; + } + } + if (e->verbose) { message( - "Allocating %zd/%zd/%zd/%zd foreign part/gpart/spart/bpart " - "(%zd/%zd/%zd/%zd MB)", + "Allocating %zd/%zd/%zd/%zd/%zd foreign part/gpart/spart/bpart/sink " + "(%zd/%zd/%zd/%zd/%zd MB)", s->size_parts_foreign, s->size_gparts_foreign, s->size_sparts_foreign, - s->size_bparts_foreign, + s->size_bparts_foreign, s->size_sinks_foreign, s->size_parts_foreign * sizeof(struct part) / (1024 * 1024), s->size_gparts_foreign * sizeof(struct gpart) / (1024 * 1024), s->size_sparts_foreign * sizeof(struct spart) / (1024 * 1024), - s->size_bparts_foreign * sizeof(struct bpart) / (1024 * 1024)); + s->size_bparts_foreign * sizeof(struct bpart) / (1024 * 1024), + s->size_sinks_foreign * sizeof(struct sink) / (1024 * 1024)); if ((s->size_parts_foreign - old_size_parts_foreign) > 0 || (s->size_gparts_foreign - old_size_gparts_foreign) > 0 || (s->size_sparts_foreign - old_size_sparts_foreign) > 0 || - (s->size_bparts_foreign - old_size_bparts_foreign) > 0) { + (s->size_bparts_foreign - old_size_bparts_foreign) > 0 || + (s->size_sinks_foreign - old_size_sinks_foreign) > 0) { message( - "Re-allocations %zd/%zd/%zd/%zd part/gpart/spart/bpart " - "(%zd/%zd/%zd/%zd MB)", + "Re-allocations %zd/%zd/%zd/%zd/%zd part/gpart/spart/bpart/sink " + "(%zd/%zd/%zd/%zd/%zd MB)", (s->size_parts_foreign - old_size_parts_foreign), (s->size_gparts_foreign - old_size_gparts_foreign), (s->size_sparts_foreign - old_size_sparts_foreign), (s->size_bparts_foreign - old_size_bparts_foreign), + (s->size_sinks_foreign - old_size_sinks_foreign), (s->size_parts_foreign - old_size_parts_foreign) * sizeof(struct part) / (1024 * 1024), (s->size_gparts_foreign - old_size_gparts_foreign) * @@ -900,7 +927,9 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) { (s->size_sparts_foreign - old_size_sparts_foreign) * sizeof(struct spart) / (1024 * 1024), (s->size_bparts_foreign - old_size_bparts_foreign) * - sizeof(struct bpart) / (1024 * 1024)); + sizeof(struct bpart) / (1024 * 1024), + (s->size_sinks_foreign - old_size_sinks_foreign) * + sizeof(struct sink) / (1024 * 1024)); } } @@ -909,6 +938,7 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) { struct gpart *gparts = s->gparts_foreign; struct spart *sparts = s->sparts_foreign; struct bpart *bparts = s->bparts_foreign; + struct sink *sinks = s->sinks_foreign; for (int k = 0; k < nr_proxies; k++) { for (int j = 0; j < e->proxies[k].nr_cells_in; j++) { @@ -940,6 +970,14 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) { cell_link_bparts(e->proxies[k].cells_in[j], bparts); bparts = &bparts[e->proxies[k].cells_in[j]->black_holes.count]; } + + if (!fof && with_sinks) { + + /* For sinks, we just use the numbers in the top-level cells */ + cell_link_sinks(e->proxies[k].cells_in[j], sinks); + sinks = + &sinks[e->proxies[k].cells_in[j]->sinks.count + space_extra_sinks]; + } } } @@ -948,6 +986,7 @@ void engine_allocate_foreign_particles(struct engine *e, const int fof) { s->nr_gparts_foreign = gparts - s->gparts_foreign; s->nr_sparts_foreign = sparts - s->sparts_foreign; s->nr_bparts_foreign = bparts - s->bparts_foreign; + s->nr_sinks_foreign = sinks - s->sinks_foreign; if (e->verbose) message("Recursively linking foreign arrays took %.3f %s.", @@ -1135,8 +1174,15 @@ int engine_estimate_nr_tasks(const struct engine *e) { #endif } if (e->policy & engine_policy_sinks) { - /* 1 drift, 2 kicks, 1 time-step, 1 sink formation */ - n1 += 5; + /* 1 drift, 2 kicks, 1 time-step, 1 sink formation | 5 + density: 1 self + 13 pairs | 14 + swallow: 1 self + 13 pairs | 14 + do_gas_swallow: 1 self + 13 pairs | 14 + do_sink_swallow: 1 self + 13 pairs | 14 + ghosts: density_ghost, sink_ghost_1, sink_ghost_2 | 3 + implicit: sink_in, sink_out | 2 */ + n1 += 66; + n2 += 3; if (e->policy & engine_policy_stars) { /* 1 star formation */ n1 += 1; @@ -1205,7 +1251,8 @@ int engine_estimate_nr_tasks(const struct engine *e) { struct cell *c = &e->s->cells_top[k]; /* Any cells with particles will have tasks (local & foreign). */ - int nparts = c->hydro.count + c->grav.count + c->stars.count; + int nparts = c->hydro.count + c->grav.count + c->stars.count + + c->black_holes.count + c->sinks.count; if (nparts > 0) { ntop++; ncells++; diff --git a/src/engine.h b/src/engine.h index 50e98011c29ff3e2cd731ebce5fcce6b8c5df085..81c6fa07109833dcabc906e435bfb6edc1fe9a0e 100644 --- a/src/engine.h +++ b/src/engine.h @@ -749,7 +749,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, size_t *Ngpart, const size_t offset_sparts, const int *ind_spart, size_t *Nspart, const size_t offset_bparts, const int *ind_bpart, - size_t *Nbpart); + size_t *Nbpart, const size_t offset_sinks, + const int *ind_sink, size_t *Nsink); void engine_rebuild(struct engine *e, int redistributed, int clean_h_values); void engine_repartition(struct engine *e); void engine_repartition_trigger(struct engine *e); diff --git a/src/engine_config.c b/src/engine_config.c index 5e6c4eb98ccc427e76f771070632c082ec9683d0..74fc0d20d6dd12fc90c17d09505c092f4da148c9 100644 --- a/src/engine_config.c +++ b/src/engine_config.c @@ -950,6 +950,8 @@ void engine_config(int restart, int fof, struct engine *e, params, "Scheduler:cell_extra_gparts", space_extra_gparts_default); space_extra_bparts = parser_get_opt_param_int( params, "Scheduler:cell_extra_bparts", space_extra_bparts_default); + space_extra_sinks = parser_get_opt_param_int( + params, "Scheduler:cell_extra_sinks", space_extra_sinks_default); /* Do we want any spare particles for on the fly creation? This condition should be the same than in space.c */ diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c index 5592e384ec9d4fc2a1609ad5744ea1b46397f8ad..e18cf26e5a1073b6f68123875d43e1f64968a908 100644 --- a/src/engine_maketasks.c +++ b/src/engine_maketasks.c @@ -365,14 +365,13 @@ void engine_addtasks_send_stars(struct engine *e, struct cell *ci, struct cell *cj, struct task *t_density, struct task *t_prep2, struct task *t_sf_counts, const int with_star_formation) { -#ifdef SWIFT_DEBUG_CHECKS +#ifdef WITH_MPI +#if !defined(SWIFT_DEBUG_CHECKS) if (e->policy & engine_policy_sinks && e->policy & engine_policy_stars) { - error("TODO"); + error("TODO: Star formation sink over MPI"); } #endif -#ifdef WITH_MPI - struct link *l = NULL; struct scheduler *s = &e->sched; const int nodeID = cj->nodeID; @@ -916,13 +915,13 @@ void engine_addtasks_recv_stars(struct engine *e, struct cell *c, struct task *t_sf_counts, struct task *const tend, const int with_star_formation) { -#ifdef SWIFT_DEBUG_CHECKS +#ifdef WITH_MPI +#if !defined(SWIFT_DEBUG_CHECKS) if (e->policy & engine_policy_sinks && e->policy & engine_policy_stars) { - error("TODO"); + error("TODO: Star formation sink over MPI"); } #endif -#ifdef WITH_MPI struct scheduler *s = &e->sched; /* Early abort (are we below the level where tasks are)? */ @@ -4378,9 +4377,9 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements, const int with_rt = (e->policy & engine_policy_rt); struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data; -#ifdef SWIFT_DEBUG_CHECKS +#if defined(WITH_MPI) && !defined(SWIFT_DEBUG_CHECKS) if (e->policy & engine_policy_sinks) { - error("TODO"); + error("TODO: Sink MPI tasks are not implemented yet!"); } #endif @@ -4451,9 +4450,9 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements, const int with_rt = (e->policy & engine_policy_rt); struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data; -#ifdef SWIFT_DEBUG_CHECKS +#if defined(WITH_MPI) && !defined(SWIFT_DEBUG_CHECKS) if (e->policy & engine_policy_sinks) { - error("TODO"); + error("TODO: Sink MPI tasks are not implemented yet!"); } #endif diff --git a/src/engine_redistribute.c b/src/engine_redistribute.c index bb8b7b6ead70796990c7ecb0311c2a65e4959bf1..cf22b146c4040d1a62c25f75b85a9975271de037 100644 --- a/src/engine_redistribute.c +++ b/src/engine_redistribute.c @@ -319,6 +319,14 @@ void ENGINE_REDISTRIBUTE_DEST_MAPPER(gpart); */ void ENGINE_REDISTRIBUTE_DEST_MAPPER(bpart); +/** + * @brief Accumulate the counts of sink particles per cell. + * Threadpool helper for accumulating the counts of particles per cell. + * + * sink version. + */ +void ENGINE_REDISTRIBUTE_DEST_MAPPER(sink); + #endif /* redist_mapper_data */ #ifdef WITH_MPI /* savelink_mapper_data */ @@ -399,12 +407,22 @@ void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 1); void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 0); #endif +/** + * @brief Save position of sink-gpart links. + * Threadpool helper for accumulating the counts of particles per cell. + */ +#ifdef SWIFT_DEBUG_CHECKS +void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(sink, 1); +#else +void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(sink, 0); +#endif + #endif /* savelink_mapper_data */ #ifdef WITH_MPI /* relink_mapper_data */ -/* Support for relinking parts, gparts, sparts and bparts after moving between - * nodes. */ +/* Support for relinking parts, gparts, sparts, bparts and sinks after moving + * between nodes. */ struct relink_mapper_data { int nodeID; int nr_nodes; @@ -412,6 +430,7 @@ struct relink_mapper_data { int *s_counts; int *g_counts; int *b_counts; + int *sink_counts; struct space *s; }; @@ -435,6 +454,7 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements, int *g_counts = mydata->g_counts; int *s_counts = mydata->s_counts; int *b_counts = mydata->b_counts; + int *sink_counts = mydata->sink_counts; struct space *s = mydata->s; for (int i = 0; i < num_elements; i++) { @@ -446,12 +466,14 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements, size_t offset_gparts = 0; size_t offset_sparts = 0; size_t offset_bparts = 0; + size_t offset_sinks = 0; for (int n = 0; n < node; n++) { int ind_recv = n * nr_nodes + nodeID; offset_parts += counts[ind_recv]; offset_gparts += g_counts[ind_recv]; offset_sparts += s_counts[ind_recv]; offset_bparts += b_counts[ind_recv]; + offset_sinks += sink_counts[ind_recv]; } /* Number of gparts sent from this node. */ @@ -493,6 +515,17 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements, s->gparts[k].id_or_neg_offset = -partner_index; s->bparts[partner_index].gpart = &s->gparts[k]; } + + /* Does this gpart have a sink partner ? */ + else if (s->gparts[k].type == swift_type_sink) { + + const ptrdiff_t partner_index = + offset_sinks - s->gparts[k].id_or_neg_offset; + + /* Re-link */ + s->gparts[k].id_or_neg_offset = -partner_index; + s->sinks[partner_index].gpart = &s->gparts[k]; + } } } } @@ -519,13 +552,6 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements, void engine_redistribute(struct engine *e) { #ifdef WITH_MPI -#ifdef SWIFT_DEBUG_CHECKS - const int nr_sinks_new = 0; -#endif - if (e->policy & engine_policy_sinks) { - error("Not implemented yet"); - } - const int nr_nodes = e->nr_nodes; const int nodeID = e->nodeID; struct space *s = e->s; @@ -536,12 +562,14 @@ void engine_redistribute(struct engine *e) { struct gpart *gparts = s->gparts; struct spart *sparts = s->sparts; struct bpart *bparts = s->bparts; + struct sink *sinks = s->sinks; ticks tic = getticks(); size_t nr_parts = s->nr_parts; size_t nr_gparts = s->nr_gparts; size_t nr_sparts = s->nr_sparts; size_t nr_bparts = s->nr_bparts; + size_t nr_sinks = s->nr_sinks; /* Start by moving inhibited particles to the end of the arrays */ for (size_t k = 0; k < nr_parts; /* void */) { @@ -609,6 +637,27 @@ void engine_redistribute(struct engine *e) { } } + /* Now move inhibited sink particles to the end of the arrays */ + for (size_t k = 0; k < nr_sinks; /* void */) { + if (sinks[k].time_bin == time_bin_inhibited || + sinks[k].time_bin == time_bin_not_created) { + nr_sinks -= 1; + + /* Swap the particle */ + memswap(&s->sinks[k], &s->sinks[nr_sinks], sizeof(struct sink)); + + /* Swap the link with the gpart */ + if (s->sinks[k].gpart != NULL) { + s->sinks[k].gpart->id_or_neg_offset = -k; + } + if (s->sinks[nr_sinks].gpart != NULL) { + s->sinks[nr_sinks].gpart->id_or_neg_offset = -nr_sinks; + } + } else { + k++; + } + } + /* Finally do the same with the gravity particles */ for (size_t k = 0; k < nr_gparts; /* void */) { if (gparts[k].time_bin == time_bin_inhibited || @@ -626,6 +675,8 @@ void engine_redistribute(struct engine *e) { s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k]; } else if (s->gparts[k].type == swift_type_black_hole) { s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k]; + } else if (s->gparts[k].type == swift_type_sink) { + s->sinks[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k]; } if (s->gparts[nr_gparts].type == swift_type_gas) { @@ -637,6 +688,9 @@ void engine_redistribute(struct engine *e) { } else if (s->gparts[nr_gparts].type == swift_type_black_hole) { s->bparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart = &s->gparts[nr_gparts]; + } else if (s->gparts[nr_gparts].type == swift_type_sink) { + s->sinks[-s->gparts[nr_sinks].id_or_neg_offset].gpart = + &s->gparts[nr_gparts]; } } else { k++; @@ -857,6 +911,73 @@ void engine_redistribute(struct engine *e) { } swift_free("b_dest", b_dest); + /* Get destination of each sink-particle */ + int *sink_counts; + if ((sink_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL) + error("Failed to allocate sink_counts temporary buffer."); + + int *sink_dest; + if ((sink_dest = (int *)swift_malloc("sink_dest", sizeof(int) * nr_sinks)) == + NULL) + error("Failed to allocate sink_dest temporary buffer."); + + redist_data.counts = sink_counts; + redist_data.dest = sink_dest; + redist_data.base = (void *)sinks; + + threadpool_map(&e->threadpool, engine_redistribute_dest_mapper_sink, sinks, + nr_sinks, sizeof(struct sink), threadpool_auto_chunk_size, + &redist_data); + + /* Sort the particles according to their cell index. */ + if (nr_sinks > 0) + space_sinks_sort(s->sinks, sink_dest, &sink_counts[nodeID * nr_nodes], + nr_nodes, 0); + +#ifdef SWIFT_DEBUG_CHECKS + /* Verify that the sink have been sorted correctly. */ + for (size_t k = 0; k < nr_sinks; k++) { + const struct sink *sink = &s->sinks[k]; + + if (sink->time_bin == time_bin_inhibited) + error("Inhibited particle found after sorting!"); + + if (sink->time_bin == time_bin_not_created) + error("Inhibited particle found after sorting!"); + + /* New cell index */ + const int new_cid = + cell_getid(s->cdim, sink->x[0] * s->iwidth[0], + sink->x[1] * s->iwidth[1], sink->x[2] * s->iwidth[2]); + + /* New cell of this sink */ + const struct cell *c = &s->cells_top[new_cid]; + const int new_node = c->nodeID; + + if (sink_dest[k] != new_node) + error("sink's new node index not matching sorted index."); + + if (sink->x[0] < c->loc[0] || sink->x[0] > c->loc[0] + c->width[0] || + sink->x[1] < c->loc[1] || sink->x[1] > c->loc[1] + c->width[1] || + sink->x[2] < c->loc[2] || sink->x[2] > c->loc[2] + c->width[2]) + error("sink not sorted into the right top-level cell!"); + } +#endif + + /* We need to re-link the gpart partners of sinks. */ + if (nr_sinks > 0) { + + struct savelink_mapper_data savelink_data; + savelink_data.nr_nodes = nr_nodes; + savelink_data.counts = sink_counts; + savelink_data.parts = (void *)sinks; + savelink_data.nodeID = nodeID; + threadpool_map(&e->threadpool, engine_redistribute_savelink_mapper_sink, + nodes, nr_nodes, sizeof(int), threadpool_auto_chunk_size, + &savelink_data); + } + swift_free("sink_dest", sink_dest); + /* Get destination of each g-particle */ int *g_counts; if ((g_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL) @@ -932,22 +1053,30 @@ void engine_redistribute(struct engine *e) { MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS) error("Failed to allreduce bparticle transfer counts."); + /* Get all the sink_counts from all the nodes. */ + if (MPI_Allreduce(MPI_IN_PLACE, sink_counts, nr_nodes * nr_nodes, MPI_INT, + MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS) + error("Failed to allreduce sink particle transfer counts."); + /* Report how many particles will be moved. */ if (e->verbose) { if (e->nodeID == 0) { - size_t total = 0, g_total = 0, s_total = 0, b_total = 0; - size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0, b_unmoved = 0; + size_t total = 0, g_total = 0, s_total = 0, b_total = 0, sink_total = 0; + size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0, b_unmoved = 0, + sink_unmoved = 0; for (int p = 0, r = 0; p < nr_nodes; p++) { for (int n = 0; n < nr_nodes; n++) { total += counts[r]; g_total += g_counts[r]; s_total += s_counts[r]; b_total += b_counts[r]; + sink_total += sink_counts[r]; if (p == n) { unmoved += counts[r]; g_unmoved += g_counts[r]; s_unmoved += s_counts[r]; b_unmoved += b_counts[r]; + sink_unmoved += sink_counts[r]; } r++; } @@ -967,14 +1096,19 @@ void engine_redistribute(struct engine *e) { message("%ld of %ld (%.2f%%) of b-particles moved", b_total - b_unmoved, b_total, 100.0 * (double)(b_total - b_unmoved) / (double)b_total); + if (sink_total > 0) + message( + "%ld of %ld (%.2f%%) of sink-particles moved", + sink_total - sink_unmoved, sink_total, + 100.0 * (double)(sink_total - sink_unmoved) / (double)sink_total); } } - /* Now each node knows how many parts, sparts, bparts, and gparts will be - * transferred to every other node. Get the new numbers of particles for this - * node. */ + /* Now each node knows how many parts, sparts, bparts, sinks and gparts will + * be transferred to every other node. Get the new numbers of particles for + * this node. */ size_t nr_parts_new = 0, nr_gparts_new = 0, nr_sparts_new = 0, - nr_bparts_new = 0; + nr_bparts_new = 0, nr_sinks_new = 0; for (int k = 0; k < nr_nodes; k++) nr_parts_new += counts[k * nr_nodes + nodeID]; for (int k = 0; k < nr_nodes; k++) @@ -983,6 +1117,8 @@ void engine_redistribute(struct engine *e) { nr_sparts_new += s_counts[k * nr_nodes + nodeID]; for (int k = 0; k < nr_nodes; k++) nr_bparts_new += b_counts[k * nr_nodes + nodeID]; + for (int k = 0; k < nr_nodes; k++) + nr_sinks_new += sink_counts[k * nr_nodes + nodeID]; #ifdef WITH_CSDS const int initial_redistribute = e->ti_current == 0; @@ -992,6 +1128,7 @@ void engine_redistribute(struct engine *e) { size_t spart_offset = 0; size_t gpart_offset = 0; size_t bpart_offset = 0; + size_t sink_offset = 0; for (int i = 0; i < nr_nodes; i++) { const size_t c_ind = engine_rank * nr_nodes + i; @@ -1002,6 +1139,7 @@ void engine_redistribute(struct engine *e) { spart_offset += s_counts[c_ind]; gpart_offset += g_counts[c_ind]; bpart_offset += b_counts[c_ind]; + sink_offset += sink_counts[c_ind]; continue; } @@ -1023,11 +1161,17 @@ void engine_redistribute(struct engine *e) { error("TODO"); } + /* Log the sinks */ + if (sink_counts[c_ind] > 0) { + error("TODO"); + } + /* Update the counters */ part_offset += counts[c_ind]; spart_offset += s_counts[c_ind]; gpart_offset += g_counts[c_ind]; bpart_offset += b_counts[c_ind]; + sink_offset += sink_counts[c_ind]; } } #endif @@ -1081,6 +1225,15 @@ void engine_redistribute(struct engine *e) { s->nr_bparts = nr_bparts_new; s->size_bparts = engine_redistribute_alloc_margin * nr_bparts_new; + /* Sink particles. */ + new_parts = engine_do_redistribute( + "sinks", sink_counts, (char *)s->sinks, nr_sinks_new, sizeof(struct sink), + sink_align, sink_mpi_type, nr_nodes, nodeID, e->syncredist); + swift_free("sinks", s->sinks); + s->sinks = (struct sink *)new_parts; + s->nr_sinks = nr_sinks_new; + s->size_sinks = engine_redistribute_alloc_margin * nr_sinks_new; + /* All particles have now arrived. Time for some final operations on the stuff we just received */ @@ -1090,6 +1243,7 @@ void engine_redistribute(struct engine *e) { size_t spart_offset = 0; size_t gpart_offset = 0; size_t bpart_offset = 0; + size_t sink_offset = 0; for (int i = 0; i < nr_nodes; i++) { const size_t c_ind = i * nr_nodes + engine_rank; @@ -1100,6 +1254,7 @@ void engine_redistribute(struct engine *e) { spart_offset += s_counts[c_ind]; gpart_offset += g_counts[c_ind]; bpart_offset += b_counts[c_ind]; + sink_offset += sink_counts[c_ind]; continue; } @@ -1121,11 +1276,17 @@ void engine_redistribute(struct engine *e) { error("TODO"); } + /* Log the sinks */ + if (sink_counts[c_ind] > 0) { + error("TODO"); + } + /* Update the counters */ part_offset += counts[c_ind]; spart_offset += s_counts[c_ind]; gpart_offset += g_counts[c_ind]; bpart_offset += b_counts[c_ind]; + sink_offset += sink_counts[c_ind]; } } #endif @@ -1139,6 +1300,7 @@ void engine_redistribute(struct engine *e) { relink_data.g_counts = g_counts; relink_data.s_counts = s_counts; relink_data.b_counts = b_counts; + relink_data.sink_counts = sink_counts; relink_data.nodeID = nodeID; relink_data.nr_nodes = nr_nodes; @@ -1151,6 +1313,7 @@ void engine_redistribute(struct engine *e) { free(g_counts); free(s_counts); free(b_counts); + free(sink_counts); #ifdef SWIFT_DEBUG_CHECKS /* Verify that all parts are in the right place. */ @@ -1186,6 +1349,15 @@ void engine_redistribute(struct engine *e) { error("Received b-particle (%zu) that does not belong here (nodeID=%i).", k, cells[cid].nodeID); } + for (size_t k = 0; k < nr_sinks_new; k++) { + const int cid = cell_getid(s->cdim, s->sinks[k].x[0] * s->iwidth[0], + s->sinks[k].x[1] * s->iwidth[1], + s->sinks[k].x[2] * s->iwidth[2]); + if (cells[cid].nodeID != nodeID) + error( + "Received sink-particle (%zu) that does not belong here (nodeID=%i).", + k, cells[cid].nodeID); + } /* Verify that the links are correct */ part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts, @@ -1200,10 +1372,11 @@ void engine_redistribute(struct engine *e) { for (int k = 0; k < nr_cells; k++) if (cells[k].nodeID == nodeID) my_cells += 1; message( - "node %i now has %zu parts, %zu sparts, %zu bparts and %zu gparts in " + "node %i now has %zu parts, %zu sparts, %zu bparts, %zu sinks and %zu " + "gparts in " "%i cells.", - nodeID, nr_parts_new, nr_sparts_new, nr_bparts_new, nr_gparts_new, - my_cells); + nodeID, nr_parts_new, nr_sparts_new, nr_bparts_new, nr_sinks_new, + nr_gparts_new, my_cells); } /* Flag that we do not have any extra particles any more */ @@ -1211,6 +1384,7 @@ void engine_redistribute(struct engine *e) { s->nr_extra_gparts = 0; s->nr_extra_sparts = 0; s->nr_extra_bparts = 0; + s->nr_extra_sinks = 0; /* Flag that a redistribute has taken place */ e->step_props |= engine_step_prop_redistribute; diff --git a/src/engine_strays.c b/src/engine_strays.c index d55909c73e08c5b594d6ab30662f3b6e64d92a25..2ecd2576fcb0091534aaf64e83550eef9d82234e 100644 --- a/src/engine_strays.c +++ b/src/engine_strays.c @@ -33,6 +33,13 @@ /* Local headers. */ #include "proxy.h" +#ifdef WITH_MPI +/* Number of particle types to wait for after launching the proxies. We have + parts, xparts, gparts, sparts, bparts and sinks to exchange, hence 6 types. + */ +#define MPI_REQUEST_NUMBER_PARTICLE_TYPES 6 +#endif + /** * @brief Exchange straying particles with other nodes. * @@ -57,18 +64,22 @@ * @param ind_bpart The foreign #cell ID of each bpart. * @param Nbpart The number of stray bparts, contains the number of bparts * received on return. + * @param offset_sinks The index in the sinks array as of which the foreign + * parts reside (i.e. the current number of local #sink). + * @param ind_sink The foreign #cell ID of each sink. + * @param Nsink The number of stray sinks, contains the number of sinks + * received on return. * * Note that this function does not mess-up the linkage between parts and * gparts, i.e. the received particles have correct linkeage. */ -void engine_exchange_strays(struct engine *e, const size_t offset_parts, - const int *restrict ind_part, size_t *Npart, - const size_t offset_gparts, - const int *restrict ind_gpart, size_t *Ngpart, - const size_t offset_sparts, - const int *restrict ind_spart, size_t *Nspart, - const size_t offset_bparts, - const int *restrict ind_bpart, size_t *Nbpart) { +void engine_exchange_strays( + struct engine *e, const size_t offset_parts, const int *restrict ind_part, + size_t *Npart, const size_t offset_gparts, const int *restrict ind_gpart, + size_t *Ngpart, const size_t offset_sparts, const int *restrict ind_spart, + size_t *Nspart, const size_t offset_bparts, const int *restrict ind_bpart, + size_t *Nbpart, const size_t offset_sinks, const int *restrict ind_sink, + size_t *Nsink) { #ifdef WITH_MPI struct space *s = e->s; @@ -80,6 +91,7 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, e->proxies[k].nr_gparts_out = 0; e->proxies[k].nr_sparts_out = 0; e->proxies[k].nr_bparts_out = 0; + e->proxies[k].nr_sinks_out = 0; } /* Put the parts into the corresponding proxies. */ @@ -204,6 +216,46 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, /* Load the bpart into the proxy */ proxy_bparts_load(&e->proxies[pid], &s->bparts[offset_bparts + k], 1); +#ifdef WITH_CSDS + if (e->policy & engine_policy_csds) { + error("Not yet implemented."); + } +#endif + } + + /* Put the sinks into the corresponding proxies. */ + for (size_t k = 0; k < *Nsink; k++) { + /* Ignore the particles we want to get rid of (inhibited, ...). */ + if (ind_sink[k] == -1) continue; + + /* Get the target node and proxy ID. */ + const int node_id = e->s->cells_top[ind_sink[k]].nodeID; + if (node_id < 0 || node_id >= e->nr_nodes) + error("Bad node ID %i.", node_id); + const int pid = e->proxy_ind[node_id]; + if (pid < 0) { + error( + "Do not have a proxy for the requested nodeID %i for part with " + "id=%lld, x=[%e,%e,%e].", + node_id, s->sinks[offset_sinks + k].id, + s->sinks[offset_sinks + k].x[0], s->sinks[offset_sinks + k].x[1], + s->sinks[offset_sinks + k].x[2]); + } + + /* Re-link the associated gpart with the buffer offset of the sink. */ + if (s->sinks[offset_sinks + k].gpart != NULL) { + s->sinks[offset_sinks + k].gpart->id_or_neg_offset = + -e->proxies[pid].nr_sinks_out; + } + +#ifdef SWIFT_DEBUG_CHECKS + if (s->sinks[offset_sinks + k].time_bin == time_bin_inhibited) + error("Attempting to exchange an inhibited particle"); +#endif + + /* Load the sink into the proxy */ + proxy_sinks_load(&e->proxies[pid], &s->sinks[offset_sinks + k], 1); + #ifdef WITH_CSDS if (e->policy & engine_policy_csds) { error("Not yet implemented."); @@ -252,8 +304,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, } /* Launch the proxies. */ - MPI_Request reqs_in[5 * engine_maxproxies]; - MPI_Request reqs_out[5 * engine_maxproxies]; + MPI_Request reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * engine_maxproxies]; + MPI_Request reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * engine_maxproxies]; for (int k = 0; k < e->nr_proxies; k++) { proxy_parts_exchange_first(&e->proxies[k]); reqs_in[k] = e->proxies[k].req_parts_count_in; @@ -281,18 +333,23 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, int count_gparts_in = 0; int count_sparts_in = 0; int count_bparts_in = 0; + int count_sinks_in = 0; for (int k = 0; k < e->nr_proxies; k++) { count_parts_in += e->proxies[k].nr_parts_in; count_gparts_in += e->proxies[k].nr_gparts_in; count_sparts_in += e->proxies[k].nr_sparts_in; count_bparts_in += e->proxies[k].nr_bparts_in; + count_sinks_in += e->proxies[k].nr_sinks_in; + message("Counting entering particles, k = %i, nr_proxies = %d", k, + e->nr_proxies); } if (e->verbose) { message( - "sent out %zu/%zu/%zu/%zu parts/gparts/sparts/bparts, got %i/%i/%i/%i " + "sent out %zu/%zu/%zu/%zu/%zu parts/gparts/sparts/bparts/sinks, got " + "%i/%i/%i/%i/%i " "back.", - *Npart, *Ngpart, *Nspart, *Nbpart, count_parts_in, count_gparts_in, - count_sparts_in, count_bparts_in); + *Npart, *Ngpart, *Nspart, *Nbpart, *Nsink, count_parts_in, + count_gparts_in, count_sparts_in, count_bparts_in, count_sinks_in); } /* Reallocate the particle arrays if necessary */ @@ -356,6 +413,24 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, } } + if (offset_sinks + count_sinks_in > s->size_sinks) { + s->size_sinks = (offset_sinks + count_sinks_in) * engine_parts_size_grow; + struct sink *sinks_new = NULL; + if (swift_memalign("sinks", (void **)&sinks_new, sink_align, + sizeof(struct sink) * s->size_sinks) != 0) + error("Failed to allocate new sink data."); + memcpy(sinks_new, s->sinks, sizeof(struct sink) * offset_sinks); + swift_free("sinks", s->sinks); + s->sinks = sinks_new; + + /* Reset the links */ + for (size_t k = 0; k < offset_sinks; k++) { + if (s->sinks[k].gpart != NULL) { + s->sinks[k].gpart->id_or_neg_offset = -k; + } + } + } + if (offset_gparts + count_gparts_in > s->size_gparts) { s->size_gparts = (offset_gparts + count_gparts_in) * engine_parts_size_grow; struct gpart *gparts_new = NULL; @@ -374,6 +449,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k]; } else if (s->gparts[k].type == swift_type_black_hole) { s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k]; + } else if (s->gparts[k].type == swift_type_sink) { + s->sinks[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k]; } } } @@ -382,82 +459,113 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, int nr_in = 0, nr_out = 0; for (int k = 0; k < e->nr_proxies; k++) { if (e->proxies[k].nr_parts_in > 0) { - reqs_in[5 * k] = e->proxies[k].req_parts_in; - reqs_in[5 * k + 1] = e->proxies[k].req_xparts_in; + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] = + e->proxies[k].req_parts_in; + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] = + e->proxies[k].req_xparts_in; nr_in += 2; } else { - reqs_in[5 * k] = reqs_in[5 * k + 1] = MPI_REQUEST_NULL; + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] = + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] = MPI_REQUEST_NULL; } if (e->proxies[k].nr_gparts_in > 0) { - reqs_in[5 * k + 2] = e->proxies[k].req_gparts_in; + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = + e->proxies[k].req_gparts_in; nr_in += 1; } else { - reqs_in[5 * k + 2] = MPI_REQUEST_NULL; + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = MPI_REQUEST_NULL; } if (e->proxies[k].nr_sparts_in > 0) { - reqs_in[5 * k + 3] = e->proxies[k].req_sparts_in; + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = + e->proxies[k].req_sparts_in; nr_in += 1; } else { - reqs_in[5 * k + 3] = MPI_REQUEST_NULL; + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = MPI_REQUEST_NULL; } if (e->proxies[k].nr_bparts_in > 0) { - reqs_in[5 * k + 4] = e->proxies[k].req_bparts_in; + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = + e->proxies[k].req_bparts_in; + nr_in += 1; + } else { + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = MPI_REQUEST_NULL; + } + if (e->proxies[k].nr_sinks_in > 0) { + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = + e->proxies[k].req_sinks_in; nr_in += 1; } else { - reqs_in[5 * k + 4] = MPI_REQUEST_NULL; + reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = MPI_REQUEST_NULL; } if (e->proxies[k].nr_parts_out > 0) { - reqs_out[5 * k] = e->proxies[k].req_parts_out; - reqs_out[5 * k + 1] = e->proxies[k].req_xparts_out; + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] = + e->proxies[k].req_parts_out; + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] = + e->proxies[k].req_xparts_out; nr_out += 2; } else { - reqs_out[5 * k] = reqs_out[5 * k + 1] = MPI_REQUEST_NULL; + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] = + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] = + MPI_REQUEST_NULL; } if (e->proxies[k].nr_gparts_out > 0) { - reqs_out[5 * k + 2] = e->proxies[k].req_gparts_out; + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = + e->proxies[k].req_gparts_out; nr_out += 1; } else { - reqs_out[5 * k + 2] = MPI_REQUEST_NULL; + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = MPI_REQUEST_NULL; } if (e->proxies[k].nr_sparts_out > 0) { - reqs_out[5 * k + 3] = e->proxies[k].req_sparts_out; + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = + e->proxies[k].req_sparts_out; nr_out += 1; } else { - reqs_out[5 * k + 3] = MPI_REQUEST_NULL; + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = MPI_REQUEST_NULL; } if (e->proxies[k].nr_bparts_out > 0) { - reqs_out[5 * k + 4] = e->proxies[k].req_bparts_out; + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = + e->proxies[k].req_bparts_out; nr_out += 1; } else { - reqs_out[5 * k + 4] = MPI_REQUEST_NULL; + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = MPI_REQUEST_NULL; + } + if (e->proxies[k].nr_sinks_out > 0) { + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = + e->proxies[k].req_sinks_out; + nr_out += 1; + } else { + reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = MPI_REQUEST_NULL; } } /* Wait for each part array to come in and collect the new parts from the proxies. */ - int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0; + int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0, + count_sinks = 0; for (int k = 0; k < nr_in; k++) { int err, pid; - if ((err = MPI_Waitany(5 * e->nr_proxies, reqs_in, &pid, - MPI_STATUS_IGNORE)) != MPI_SUCCESS) { + if ((err = MPI_Waitany(MPI_REQUEST_NUMBER_PARTICLE_TYPES * e->nr_proxies, + reqs_in, &pid, MPI_STATUS_IGNORE)) != MPI_SUCCESS) { char buff[MPI_MAX_ERROR_STRING]; int res; MPI_Error_string(err, buff, &res); error("MPI_Waitany failed (%s).", buff); } if (pid == MPI_UNDEFINED) break; - // message( "request from proxy %i has arrived." , pid / 5 ); - pid = 5 * (pid / 5); + // message( "request from proxy %i has arrived." , pid / + // MPI_REQUEST_NUMBER_PARTICLE_TYPES ); + pid = MPI_REQUEST_NUMBER_PARTICLE_TYPES * + (pid / MPI_REQUEST_NUMBER_PARTICLE_TYPES); /* If all the requests for a given proxy have arrived... */ if (reqs_in[pid + 0] == MPI_REQUEST_NULL && reqs_in[pid + 1] == MPI_REQUEST_NULL && reqs_in[pid + 2] == MPI_REQUEST_NULL && reqs_in[pid + 3] == MPI_REQUEST_NULL && - reqs_in[pid + 4] == MPI_REQUEST_NULL) { + reqs_in[pid + 4] == MPI_REQUEST_NULL && + reqs_in[pid + 5] == MPI_REQUEST_NULL) { /* Copy the particle data to the part/xpart/gpart arrays. */ - struct proxy *prox = &e->proxies[pid / 5]; + struct proxy *prox = &e->proxies[pid / MPI_REQUEST_NUMBER_PARTICLE_TYPES]; memcpy(&s->parts[offset_parts + count_parts], prox->parts_in, sizeof(struct part) * prox->nr_parts_in); memcpy(&s->xparts[offset_parts + count_parts], prox->xparts_in, @@ -468,6 +576,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, sizeof(struct spart) * prox->nr_sparts_in); memcpy(&s->bparts[offset_bparts + count_bparts], prox->bparts_in, sizeof(struct bpart) * prox->nr_bparts_in); + memcpy(&s->sinks[offset_sinks + count_sinks], prox->sinks_in, + sizeof(struct sink) * prox->nr_sinks_in); #ifdef WITH_CSDS if (e->policy & engine_policy_csds) { @@ -495,6 +605,12 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, if (prox->nr_bparts_in > 0) { error("TODO"); } + + /* Log the sinks */ + if (prox->nr_sinks_in > 0) { + /* Not implemented yet */ + error("TODO"); + } } #endif /* for (int k = offset; k < offset + count; k++) @@ -522,6 +638,11 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, &s->bparts[offset_bparts + count_bparts - gp->id_or_neg_offset]; gp->id_or_neg_offset = s->bparts - bp; bp->gpart = gp; + } else if (gp->type == swift_type_sink) { + struct sink *sink = + &s->sinks[offset_sinks + count_sinks - gp->id_or_neg_offset]; + gp->id_or_neg_offset = s->sinks - sink; + sink->gpart = gp; } } @@ -530,13 +651,14 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, count_gparts += prox->nr_gparts_in; count_sparts += prox->nr_sparts_in; count_bparts += prox->nr_bparts_in; + count_sinks += prox->nr_sinks_in; } } /* Wait for all the sends to have finished too. */ if (nr_out > 0) - if (MPI_Waitall(5 * e->nr_proxies, reqs_out, MPI_STATUSES_IGNORE) != - MPI_SUCCESS) + if (MPI_Waitall(MPI_REQUEST_NUMBER_PARTICLE_TYPES * e->nr_proxies, reqs_out, + MPI_STATUSES_IGNORE) != MPI_SUCCESS) error("MPI_Waitall on sends failed."); /* Free the proxy memory */ @@ -553,6 +675,7 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts, *Ngpart = count_gparts; *Nspart = count_sparts; *Nbpart = count_bparts; + *Nsink = count_sinks; #else error("SWIFT was not compiled with MPI support."); diff --git a/src/feedback/GEAR/feedback.c b/src/feedback/GEAR/feedback.c index 46dfc6dea1ced3d51dc7eea992f1ea3c3657ad98..204477d474612d95419117d6d745844454c4330a 100644 --- a/src/feedback/GEAR/feedback.c +++ b/src/feedback/GEAR/feedback.c @@ -21,6 +21,7 @@ #include "feedback.h" /* Local includes */ +#include "cooling.h" #include "cosmology.h" #include "engine.h" #include "error.h" @@ -49,7 +50,7 @@ void feedback_update_part(struct part* p, struct xpart* xp, const struct pressure_floor_props* pressure_floor = e->pressure_floor_props; /* Turn off the cooling */ - xp->cooling_data.time_last_event = e->time; + cooling_set_part_time_cooling_off(p, xp, e->time); /* Update mass */ const float old_mass = hydro_get_mass(p); diff --git a/src/feedback/GEAR/initial_mass_function.c b/src/feedback/GEAR/initial_mass_function.c index 5c4f1f8838943dcca832575af03bfa88c820f0b5..7509ed28598cb5c41819606f1f469298ecc12c69 100644 --- a/src/feedback/GEAR/initial_mass_function.c +++ b/src/feedback/GEAR/initial_mass_function.c @@ -642,7 +642,7 @@ INLINE double initial_mass_function_sample_power_law(double min_mass, * first stars or not. You need to verify this before this function and pass the * correct values to 'minimal_discrete_mass_Msun' and 'stellar_particle_mass'. * - * Note 2: This function implicitiely assumes M_sun since the IMF data + * Note 2: This function implicitly assumes M_sun since the IMF data * structures handles the masses in M_sun. * * @param imf The #initial_mass_function. diff --git a/src/part.c b/src/part.c index 9d9c38f577f8675d7cb9367a2227bce689b892f0..f513153d61ba9b277a0b6b88ae9e4cf61c766e10 100644 --- a/src/part.c +++ b/src/part.c @@ -544,6 +544,7 @@ MPI_Datatype xpart_mpi_type; MPI_Datatype gpart_mpi_type; MPI_Datatype spart_mpi_type; MPI_Datatype bpart_mpi_type; +MPI_Datatype sink_mpi_type; /** * @brief Registers MPI particle types. @@ -581,6 +582,11 @@ void part_create_mpi_types(void) { MPI_Type_commit(&bpart_mpi_type) != MPI_SUCCESS) { error("Failed to create MPI type for bparts."); } + if (MPI_Type_contiguous(sizeof(struct sink) / sizeof(unsigned char), MPI_BYTE, + &sink_mpi_type) != MPI_SUCCESS || + MPI_Type_commit(&sink_mpi_type) != MPI_SUCCESS) { + error("Failed to create MPI type for sink."); + } } void part_free_mpi_types(void) { @@ -590,5 +596,6 @@ void part_free_mpi_types(void) { MPI_Type_free(&gpart_mpi_type); MPI_Type_free(&spart_mpi_type); MPI_Type_free(&bpart_mpi_type); + MPI_Type_free(&sink_mpi_type); } #endif diff --git a/src/part.h b/src/part.h index b8fc625885a9f329633bc39476ecb7aac2228b79..c0e87bfde7c6525b7cdc05974d56fa4bb5be5bd9 100644 --- a/src/part.h +++ b/src/part.h @@ -171,6 +171,7 @@ extern MPI_Datatype xpart_mpi_type; extern MPI_Datatype gpart_mpi_type; extern MPI_Datatype spart_mpi_type; extern MPI_Datatype bpart_mpi_type; +extern MPI_Datatype sink_mpi_type; void part_create_mpi_types(void); void part_free_mpi_types(void); diff --git a/src/proxy.c b/src/proxy.c index 22663d42a70c75b3966d0272d93a294480057b0b..3830162ada2135d28c68b3763661a18e17bcbab0 100644 --- a/src/proxy.c +++ b/src/proxy.c @@ -705,12 +705,24 @@ void proxy_parts_exchange_first(struct proxy *p) { p->buff_out[1] = p->nr_gparts_out; p->buff_out[2] = p->nr_sparts_out; p->buff_out[3] = p->nr_bparts_out; - if (MPI_Isend(p->buff_out, 4, MPI_INT, p->nodeID, - p->mynodeID * proxy_tag_shift + proxy_tag_count, MPI_COMM_WORLD, - &p->req_parts_count_out) != MPI_SUCCESS) + p->buff_out[4] = p->nr_sinks_out; + +#ifdef SWIFT_DEBUG_CHECKS + message("Number of particles out [%i , %i, %i, %i, %i]", p->nr_parts_out, + p->nr_gparts_out, p->nr_sparts_out, p->nr_bparts_out, + p->nr_sinks_out); +#endif /* SWIFT_DEBUG_CHECKS */ + + if (MPI_Isend(p->buff_out, PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES, MPI_INT, + p->nodeID, p->mynodeID * proxy_tag_shift + proxy_tag_count, + MPI_COMM_WORLD, &p->req_parts_count_out) != MPI_SUCCESS) error("Failed to isend nr of parts."); - /* message( "isent particle counts [%i, %i] from node %i to node %i." , - p->buff_out[0], p->buff_out[1], p->mynodeID , p->nodeID ); fflush(stdout); */ +#ifdef SWIFT_DEBUG_CHECKS + message("isent particle counts [%i, %i, %i, %i, %i] from node %i to node %i.", + p->buff_out[0], p->buff_out[1], p->buff_out[2], p->buff_out[3], + p->buff_out[4], p->mynodeID, p->nodeID); + fflush(stdout); +#endif /* SWIFT_DEBUG_CHECKS */ /* Send the particle buffers. */ if (p->nr_parts_out > 0) { @@ -721,20 +733,28 @@ void proxy_parts_exchange_first(struct proxy *p) { p->mynodeID * proxy_tag_shift + proxy_tag_xparts, MPI_COMM_WORLD, &p->req_xparts_out) != MPI_SUCCESS) error("Failed to isend part data."); - // message( "isent particle data (%i) to node %i." , p->nr_parts_out , - // p->nodeID ); fflush(stdout); - /*for (int k = 0; k < p->nr_parts_out; k++) +#ifdef SWIFT_DEBUG_CHECKS + message("isent particle data (%i) to node %i.", p->nr_parts_out, p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_parts_out; k++) message("sending particle %lli, x=[%.3e %.3e %.3e], h=%.3e, to node %i.", p->parts_out[k].id, p->parts_out[k].x[0], p->parts_out[k].x[1], - p->parts_out[k].x[2], p->parts_out[k].h, p->nodeID);*/ + p->parts_out[k].x[2], p->parts_out[k].h, p->nodeID); +#endif /* SWIFT_DEBUG_CHECKS */ } if (p->nr_gparts_out > 0) { if (MPI_Isend(p->gparts_out, p->nr_gparts_out, gpart_mpi_type, p->nodeID, p->mynodeID * proxy_tag_shift + proxy_tag_gparts, MPI_COMM_WORLD, &p->req_gparts_out) != MPI_SUCCESS) error("Failed to isend gpart data."); - // message( "isent gpart data (%i) to node %i." , p->nr_gparts_out , - // p->nodeID ); fflush(stdout); +#ifdef SWIFT_DEBUG_CHECKS + message("isent gpart data (%i) to node %i.", p->nr_gparts_out, p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_parts_out; k++) + message("sending gpart %lli, x=[%.3e %.3e %.3e], to node %i.", + p->gparts_out[k].id_or_neg_offset, p->gparts_out[k].x[0], + p->gparts_out[k].x[1], p->gparts_out[k].x[2], p->nodeID); +#endif /* SWIFT_DEBUG_CHECKS */ } if (p->nr_sparts_out > 0) { @@ -742,23 +762,56 @@ void proxy_parts_exchange_first(struct proxy *p) { p->mynodeID * proxy_tag_shift + proxy_tag_sparts, MPI_COMM_WORLD, &p->req_sparts_out) != MPI_SUCCESS) error("Failed to isend spart data."); - // message( "isent spart data (%i) to node %i." , p->nr_sparts_out , - // p->nodeID ); fflush(stdout); +#ifdef SWIFT_DEBUG_CHECKS + message("isent spart data (%i) to node %i.", p->nr_sparts_out, p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_sparts_out; k++) + message("sending spart %lli, x=[%.3e %.3e %.3e], h=%.3e, to node %i.", + p->sparts_out[k].id, p->sparts_out[k].x[0], p->sparts_out[k].x[1], + p->sparts_out[k].x[2], p->sparts_out[k].h, p->nodeID); +#endif /* SWIFT_DEBUG_CHECKS */ } if (p->nr_bparts_out > 0) { if (MPI_Isend(p->bparts_out, p->nr_bparts_out, bpart_mpi_type, p->nodeID, p->mynodeID * proxy_tag_shift + proxy_tag_bparts, MPI_COMM_WORLD, &p->req_bparts_out) != MPI_SUCCESS) error("Failed to isend bpart data."); - // message( "isent bpart data (%i) to node %i." , p->nr_bparts_out , - // p->nodeID ); fflush(stdout); +#ifdef SWIFT_DEBUG_CHECKS + message("isent bpart data (%i) to node %i.", p->nr_bparts_out, p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_bparts_out; k++) + message("sending bpart %lli, x=[%.3e %.3e %.3e], h=%.3e, to node %i.", + p->bparts_out[k].id, p->bparts_out[k].x[0], p->bparts_out[k].x[1], + p->bparts_out[k].x[2], p->bparts_out[k].h, p->nodeID); +#endif /* SWIFT_DEBUG_CHECKS */ + } + if (p->nr_sinks_out > 0) { + if (MPI_Isend(p->sinks_out, p->nr_sinks_out, sink_mpi_type, p->nodeID, + p->mynodeID * proxy_tag_shift + proxy_tag_sinks, + MPI_COMM_WORLD, &p->req_sinks_out) != MPI_SUCCESS) + error("Failed to isend sink data."); +#ifdef SWIFT_DEBUG_CHECKS + message("isent sink data (%i) to node %i.", p->nr_sinks_out, p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_sinks_out; k++) + message("sending sinks %lli, x=[%.3e %.3e %.3e], h=%.3e, to node %i.", + p->sinks_out[k].id, p->sinks_out[k].x[0], p->sinks_out[k].x[1], + p->sinks_out[k].x[2], p->sinks_out[k].r_cut, p->nodeID); +#endif /* SWIFT_DEBUG_CHECKS */ } /* Receive the number of particles. */ - if (MPI_Irecv(p->buff_in, 4, MPI_INT, p->nodeID, - p->nodeID * proxy_tag_shift + proxy_tag_count, MPI_COMM_WORLD, - &p->req_parts_count_in) != MPI_SUCCESS) + if (MPI_Irecv(p->buff_in, PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES, MPI_INT, + p->nodeID, p->nodeID * proxy_tag_shift + proxy_tag_count, + MPI_COMM_WORLD, &p->req_parts_count_in) != MPI_SUCCESS) error("Failed to irecv nr of parts."); +#ifdef SWIFT_DEBUG_CHECKS + message( + "irecv particle counts [%i, %i, %i, %i, %i] from node %i, I am node %i.", + p->buff_in[0], p->buff_in[1], p->buff_in[2], p->buff_in[3], p->buff_in[4], + p->nodeID, p->mynodeID); + fflush(stdout); +#endif /* SWIFT_DEBUG_CHECKS */ #else error("SWIFT was not compiled with MPI support."); @@ -774,6 +827,7 @@ void proxy_parts_exchange_second(struct proxy *p) { p->nr_gparts_in = p->buff_in[1]; p->nr_sparts_in = p->buff_in[2]; p->nr_bparts_in = p->buff_in[3]; + p->nr_sinks_in = p->buff_in[4]; /* Is there enough space in the buffers? */ if (p->nr_parts_in > p->size_parts_in) { @@ -815,6 +869,15 @@ void proxy_parts_exchange_second(struct proxy *p) { "bparts_in", sizeof(struct bpart) * p->size_bparts_in)) == NULL) error("Failed to re-allocate bparts_in buffers."); } + if (p->nr_sinks_in > p->size_sinks_in) { + do { + p->size_sinks_in *= proxy_buffgrow; + } while (p->nr_sinks_in > p->size_sinks_in); + swift_free("sinks_in", p->sinks_in); + if ((p->sinks_in = (struct sink *)swift_malloc( + "sinks_in", sizeof(struct sink) * p->size_sinks_in)) == NULL) + error("Failed to re-allocate sinks_in buffers."); + } /* Receive the particle buffers. */ if (p->nr_parts_in > 0) { @@ -825,32 +888,73 @@ void proxy_parts_exchange_second(struct proxy *p) { p->nodeID * proxy_tag_shift + proxy_tag_xparts, MPI_COMM_WORLD, &p->req_xparts_in) != MPI_SUCCESS) error("Failed to irecv part data."); - // message( "irecv particle data (%i) from node %i." , p->nr_parts_in , - // p->nodeID ); fflush(stdout); +#ifdef SWIFT_DEBUG_CHECKS + message("irecv particle data (%i) from node %i.", p->nr_parts_in, + p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_parts_in; k++) + message("receiving parts %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.", + p->parts_in[k].id, p->parts_in[k].x[0], p->parts_in[k].x[1], + p->parts_in[k].x[2], p->parts_in[k].h, p->nodeID); +#endif /* SWIFT_DEBUG_CHECKS */ } if (p->nr_gparts_in > 0) { if (MPI_Irecv(p->gparts_in, p->nr_gparts_in, gpart_mpi_type, p->nodeID, p->nodeID * proxy_tag_shift + proxy_tag_gparts, MPI_COMM_WORLD, &p->req_gparts_in) != MPI_SUCCESS) error("Failed to irecv gpart data."); - // message( "irecv gpart data (%i) from node %i." , p->nr_gparts_in , - // p->nodeID ); fflush(stdout); +#ifdef SWIFT_DEBUG_CHECKS + message("irecv gpart data (%i) from node %i.", p->nr_gparts_in, p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_gparts_in; k++) + message("receiving gparts %lli, x=[%.3e %.3e %.3e], from node %i.", + p->gparts_in[k].id_or_neg_offset, p->gparts_in[k].x[0], + p->gparts_in[k].x[1], p->gparts_in[k].x[2], p->nodeID); +#endif /* SWIFT_DEBUG_CHECKS */ } if (p->nr_sparts_in > 0) { if (MPI_Irecv(p->sparts_in, p->nr_sparts_in, spart_mpi_type, p->nodeID, p->nodeID * proxy_tag_shift + proxy_tag_sparts, MPI_COMM_WORLD, &p->req_sparts_in) != MPI_SUCCESS) error("Failed to irecv spart data."); - // message( "irecv spart data (%i) from node %i." , p->nr_sparts_in , - // p->nodeID ); fflush(stdout); +#ifdef SWIFT_DEBUG_CHECKS + message("irecv spart data (%i) from node %i.", p->nr_sparts_in, p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_sparts_in; k++) + message( + "receiving sparts %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.", + p->sparts_in[k].id, p->sparts_in[k].x[0], p->sparts_in[k].x[1], + p->sparts_in[k].x[2], p->sparts_in[k].h, p->nodeID); +#endif } if (p->nr_bparts_in > 0) { if (MPI_Irecv(p->bparts_in, p->nr_bparts_in, bpart_mpi_type, p->nodeID, p->nodeID * proxy_tag_shift + proxy_tag_bparts, MPI_COMM_WORLD, &p->req_bparts_in) != MPI_SUCCESS) error("Failed to irecv bpart data."); - // message( "irecv bpart data (%i) from node %i." , p->nr_bparts_in , - // p->nodeID ); fflush(stdout); +#ifdef SWIFT_DEBUG_CHECKS + message("irecv bpart data (%i) from node %i.", p->nr_bparts_in, p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_bparts_in; k++) + message( + "receiving bparts %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.", + p->bparts_in[k].id, p->bparts_in[k].x[0], p->bparts_in[k].x[1], + p->bparts_in[k].x[2], p->bparts_in[k].h, p->nodeID); +#endif /* SWIFT_DEBUG_CHECKS */ + } + if (p->nr_sinks_in > 0) { + if (MPI_Irecv(p->sinks_in, p->nr_sinks_in, sink_mpi_type, p->nodeID, + p->nodeID * proxy_tag_shift + proxy_tag_sinks, MPI_COMM_WORLD, + &p->req_sinks_in) != MPI_SUCCESS) + error("Failed to irecv sink data."); +#ifdef SWIFT_DEBUG_CHECKS + message("irecv sink data (%i) from node %i.", p->nr_sinks_in, p->nodeID); + fflush(stdout); + for (int k = 0; k < p->nr_sinks_in; k++) + message("receiving sinks %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.", + p->sinks_in[k].id, p->sinks_in[k].x[0], p->sinks_in[k].x[1], + p->sinks_in[k].x[2], p->sinks_in[k].r_cut, p->nodeID); +#endif /* SWIFT_DEBUG_CHECKS */ } #else @@ -987,6 +1091,36 @@ void proxy_bparts_load(struct proxy *p, const struct bpart *bparts, int N) { p->nr_bparts_out += N; } +/** + * @brief Load sinks onto a proxy for exchange. + * + * @param p The #proxy. + * @param sinks Pointer to an array of #sink to send. + * @param N The number of sinks. + */ +void proxy_sinks_load(struct proxy *p, const struct sink *sinks, int N) { + + /* Is there enough space in the buffer? */ + if (p->nr_sinks_out + N > p->size_sinks_out) { + do { + p->size_sinks_out *= proxy_buffgrow; + } while (p->nr_sinks_out + N > p->size_sinks_out); + struct sink *tp; + if ((tp = (struct sink *)swift_malloc( + "sinks_out", sizeof(struct sink) * p->size_sinks_out)) == NULL) + error("Failed to re-allocate sinks_out buffers."); + memcpy(tp, p->sinks_out, sizeof(struct sink) * p->nr_sinks_out); + swift_free("sinks_out", p->sinks_out); + p->sinks_out = tp; + } + + /* Copy the parts and xparts data to the buffer. */ + memcpy(&p->sinks_out[p->nr_sinks_out], sinks, sizeof(struct sink) * N); + + /* Increase the counters. */ + p->nr_sinks_out += N; +} + /** * @brief Frees the memory allocated for the particle proxies and sets their * size back to the initial state. @@ -1062,6 +1196,21 @@ void proxy_free_particle_buffers(struct proxy *p) { "bparts_in", sizeof(struct bpart) * p->size_bparts_in)) == NULL) error("Failed to allocate bparts_in buffers."); } + + if (p->size_sinks_out > proxy_buffinit) { + swift_free("sinks_out", p->sinks_out); + p->size_sinks_out = proxy_buffinit; + if ((p->sinks_out = (struct sink *)swift_malloc( + "sinks_out", sizeof(struct sink) * p->size_sinks_out)) == NULL) + error("Failed to allocate sinks_out buffers."); + } + if (p->size_sinks_in > proxy_buffinit) { + swift_free("sinks_in", p->sinks_in); + p->size_sinks_in = proxy_buffinit; + if ((p->sinks_in = (struct sink *)swift_malloc( + "sinks_in", sizeof(struct sink) * p->size_sinks_in)) == NULL) + error("Failed to allocate sinks_in buffers."); + } } /** @@ -1166,6 +1315,22 @@ void proxy_init(struct proxy *p, int mynodeID, int nodeID) { error("Failed to allocate bparts_out buffers."); } p->nr_bparts_out = 0; + + /* Allocate the sinks send and receive buffers, if needed. */ + if (p->sinks_in == NULL) { + p->size_sinks_in = proxy_buffinit; + if ((p->sinks_in = (struct sink *)swift_malloc( + "sinks_in", sizeof(struct sink) * p->size_sinks_in)) == NULL) + error("Failed to allocate sinks_in buffers."); + } + p->nr_sinks_in = 0; + if (p->sinks_out == NULL) { + p->size_sinks_out = proxy_buffinit; + if ((p->sinks_out = (struct sink *)swift_malloc( + "sinks_out", sizeof(struct sink) * p->size_sinks_out)) == NULL) + error("Failed to allocate sinks_out buffers."); + } + p->nr_sinks_out = 0; } /** @@ -1184,11 +1349,13 @@ void proxy_clean(struct proxy *p) { swift_free("gparts_out", p->gparts_out); swift_free("sparts_out", p->sparts_out); swift_free("bparts_out", p->bparts_out); + swift_free("sinks_out", p->sinks_out); swift_free("parts_in", p->parts_in); swift_free("xparts_in", p->xparts_in); swift_free("gparts_in", p->gparts_in); swift_free("sparts_in", p->sparts_in); swift_free("bparts_in", p->bparts_in); + swift_free("sinks_in", p->sinks_in); } /** diff --git a/src/proxy.h b/src/proxy.h index d0d2bfb582828eefbbe56d5e7eeba850fa340ffc..46a322e3a49ff53f796353e90aa84852f855acd0 100644 --- a/src/proxy.h +++ b/src/proxy.h @@ -28,15 +28,21 @@ #define proxy_buffgrow 1.5 #define proxy_buffinit 100 +/* Number of particle types to exchange with proxies in + proxy_parts_exchange_first(). We have parts, gparts, sparts, bparts and + sinks to exchange, hence 5 types. */ +#define PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES 5 + /* Proxy tag arithmetic. */ -#define proxy_tag_shift 8 +#define proxy_tag_shift 9 #define proxy_tag_count 0 #define proxy_tag_parts 1 #define proxy_tag_xparts 2 #define proxy_tag_gparts 3 #define proxy_tag_sparts 4 #define proxy_tag_bparts 5 -#define proxy_tag_cells 6 +#define proxy_tag_sinks 6 +#define proxy_tag_cells 7 /** * @brief The different reasons a cell can be in a proxy @@ -71,6 +77,7 @@ struct proxy { struct gpart *gparts_in, *gparts_out; struct spart *sparts_in, *sparts_out; struct bpart *bparts_in, *bparts_out; + struct sink *sinks_in, *sinks_out; int size_parts_in, size_parts_out; int nr_parts_in, nr_parts_out; int size_gparts_in, size_gparts_out; @@ -79,9 +86,12 @@ struct proxy { int nr_sparts_in, nr_sparts_out; int size_bparts_in, size_bparts_out; int nr_bparts_in, nr_bparts_out; + int size_sinks_in, size_sinks_out; + int nr_sinks_in, nr_sinks_out; /* Buffer to hold the incomming/outgoing particle counts. */ - int buff_out[4], buff_in[4]; + int buff_out[PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES], + buff_in[PROXY_EXCHANGE_NUMBER_PARTICLE_TYPES]; /* MPI request handles. */ #ifdef WITH_MPI @@ -91,6 +101,7 @@ struct proxy { MPI_Request req_gparts_out, req_gparts_in; MPI_Request req_sparts_out, req_sparts_in; MPI_Request req_bparts_out, req_bparts_in; + MPI_Request req_sinks_out, req_sinks_in; MPI_Request req_cells_count_out, req_cells_count_in; MPI_Request req_cells_out, req_cells_in; #endif @@ -104,6 +115,7 @@ void proxy_parts_load(struct proxy *p, const struct part *parts, void proxy_gparts_load(struct proxy *p, const struct gpart *gparts, int N); void proxy_sparts_load(struct proxy *p, const struct spart *sparts, int N); void proxy_bparts_load(struct proxy *p, const struct bpart *bparts, int N); +void proxy_sinks_load(struct proxy *p, const struct sink *sinks, int N); void proxy_free_particle_buffers(struct proxy *p); void proxy_parts_exchange_first(struct proxy *p); void proxy_parts_exchange_second(struct proxy *p); diff --git a/src/runner_doiact_functions_sinks.h b/src/runner_doiact_functions_sinks.h index c21c17812849118d1630aeb894591da8c88fdbfc..6ff95ba1deefd310b28f4e24d30bc62d23b5663c 100644 --- a/src/runner_doiact_functions_sinks.h +++ b/src/runner_doiact_functions_sinks.h @@ -346,6 +346,10 @@ void DOPAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci, */ void DOSELF1_BRANCH_SINKS(struct runner *r, struct cell *c) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif + const struct engine *restrict e = r->e; /* Anything to do here? */ @@ -371,6 +375,10 @@ void DOSELF1_BRANCH_SINKS(struct runner *r, struct cell *c) { */ void DOPAIR1_BRANCH_SINKS(struct runner *r, struct cell *ci, struct cell *cj) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif + const struct engine *restrict e = r->e; const int ci_active = cell_is_active_sinks(ci, e); @@ -417,6 +425,9 @@ void DOPAIR1_BRANCH_SINKS(struct runner *r, struct cell *ci, struct cell *cj) { */ void DOSUB_PAIR1_SINKS(struct runner *r, struct cell *ci, struct cell *cj, int timer) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif TIMER_TIC; @@ -496,6 +507,9 @@ void DOSUB_PAIR1_SINKS(struct runner *r, struct cell *ci, struct cell *cj, * @param gettimer Do we have a timer ? */ void DOSUB_SELF1_SINKS(struct runner *r, struct cell *ci, int timer) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif TIMER_TIC; diff --git a/src/runner_others.c b/src/runner_others.c index 86b2d6fcbade1dca24eefd5629fb11b333fb8772..82d58f22c46ec5aa0f6e8240e5bd2ef477981da0 100644 --- a/src/runner_others.c +++ b/src/runner_others.c @@ -198,6 +198,9 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) { */ void runner_do_star_formation_sink(struct runner *r, struct cell *c, int timer) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif struct engine *e = r->e; const struct cosmology *cosmo = e->cosmology; @@ -567,6 +570,10 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { */ void runner_do_sink_formation(struct runner *r, struct cell *c) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif + struct engine *e = r->e; const struct cosmology *cosmo = e->cosmology; const struct sink_props *sink_props = e->sink_properties; diff --git a/src/runner_sinks.c b/src/runner_sinks.c index fa614f2cbaa23217e2b61f704dbf31fc3039d0f5..cb085dc103bfb3c37285ea92eb736edbedc98ce0 100644 --- a/src/runner_sinks.c +++ b/src/runner_sinks.c @@ -210,6 +210,9 @@ void runner_do_sinks_gas_swallow(struct runner *r, struct cell *c, int timer) { */ void runner_do_sinks_gas_swallow_self(struct runner *r, struct cell *c, int timer) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif #ifdef SWIFT_DEBUG_CHECKS if (c->nodeID != r->e->nodeID) error("Running self task on foreign node"); @@ -230,6 +233,9 @@ void runner_do_sinks_gas_swallow_self(struct runner *r, struct cell *c, */ void runner_do_sinks_gas_swallow_pair(struct runner *r, struct cell *ci, struct cell *cj, int timer) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif const struct engine *e = r->e; @@ -399,6 +405,9 @@ void runner_do_sinks_sink_swallow(struct runner *r, struct cell *c, int timer) { */ void runner_do_sinks_sink_swallow_self(struct runner *r, struct cell *c, int timer) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif #ifdef SWIFT_DEBUG_CHECKS if (c->nodeID != r->e->nodeID) error("Running self task on foreign node"); @@ -419,6 +428,9 @@ void runner_do_sinks_sink_swallow_self(struct runner *r, struct cell *c, */ void runner_do_sinks_sink_swallow_pair(struct runner *r, struct cell *ci, struct cell *cj, int timer) { +#ifdef SWIFT_DEBUG_CHECKS_MPI_DOMAIN_DECOMPOSITION + return; +#endif const struct engine *e = r->e; @@ -449,6 +461,7 @@ void runner_do_prepare_part_sink_formation(struct runner *r, struct cell *c, struct xpart *restrict xp) { struct engine *e = r->e; const struct cosmology *cosmo = e->cosmology; + const int with_cosmology = e->policy & engine_policy_cosmology; const struct sink_props *sink_props = e->sink_properties; const int count = c->hydro.count; struct part *restrict parts = c->hydro.parts; @@ -478,8 +491,8 @@ void runner_do_prepare_part_sink_formation(struct runner *r, struct cell *c, struct sink *restrict si = &sinks[i]; /* Compute the quantities required to later decide to form a sink or not. */ - sink_prepare_part_sink_formation_sink_criteria(e, p, xp, si, cosmo, - sink_props); + sink_prepare_part_sink_formation_sink_criteria(e, p, xp, si, with_cosmology, + cosmo, sink_props, e->time); } /* End of sink neighbour loop */ } diff --git a/src/sink/Default/sink.h b/src/sink/Default/sink.h index f284b9d0f916b02bd7df578274953c5eb120f706..e0d3f3b4a39f0d067d22c9bc088e481d5da9849a 100644 --- a/src/sink/Default/sink.h +++ b/src/sink/Default/sink.h @@ -30,10 +30,22 @@ /** * @brief Computes the time-step of a given sink particle. * - * @param sp Pointer to the sink-particle data. + * Note: In the Default sink, no time-step limit is applied. + * + * @param sink Pointer to the sink-particle data. + * @param sink_properties Properties of the sink model. + * @param with_cosmology Are we running with cosmological time integration. + * @param cosmo The current cosmological model (used if running with + * cosmology). + * @param grav_props The current gravity properties. + * @param time The current time (used if running without cosmology). + * @param time_base The time base. */ __attribute__((always_inline)) INLINE static float sink_compute_timestep( - const struct sink* const sp) { + const struct sink* const sink, const struct sink_props* sink_properties, + const int with_cosmology, const struct cosmology* cosmo, + const struct gravity_props* grav_props, const double time, + const double time_base) { return FLT_MAX; } @@ -372,7 +384,8 @@ INLINE static void sink_prepare_part_sink_formation_gas_criteria( */ INLINE static void sink_prepare_part_sink_formation_sink_criteria( struct engine* e, struct part* restrict p, struct xpart* restrict xp, - struct sink* restrict si, const struct cosmology* cosmo, - const struct sink_props* sink_props) {} + struct sink* restrict si, const int with_cosmology, + const struct cosmology* cosmo, const struct sink_props* sink_props, + const double time) {} #endif /* SWIFT_DEFAULT_SINK_H */ diff --git a/src/sink/GEAR/sink.h b/src/sink/GEAR/sink.h index 907d3d38aeee6ab55b5409324cf8e6f4d90e7f0e..627713379e85049f13ee709242efdd6e168f9e53 100644 --- a/src/sink/GEAR/sink.h +++ b/src/sink/GEAR/sink.h @@ -34,76 +34,122 @@ #include "feedback.h" #include "minmax.h" #include "random.h" +#include "sink_getters.h" #include "sink_part.h" #include "sink_properties.h" +#include "sink_setters.h" #include "star_formation.h" /** * @brief Computes the time-step of a given sink particle. * - * @param sp Pointer to the sink-particle data. + * @param sink Pointer to the sink-particle data. + * @param sink_properties Properties of the sink model. + * @param with_cosmology Are we running with cosmological time integration. + * @param cosmo The current cosmological model (used if running with + * cosmology). + * @param grav_props The current gravity properties. + * @param time The current time (used if running without cosmology). + * @param time_base The time base. */ __attribute__((always_inline)) INLINE static float sink_compute_timestep( - const struct sink* const sp) { - - return FLT_MAX; -} - -/** - * @brief Update the target mass of the sink particle. - * - * @param sink the sink particle. - * @param sink_props the sink properties to use. - * @param phys_const the physical constants in internal units. - * @param e The #engine - * @param star_counter The star loop counter. - */ -INLINE static void sink_update_target_mass(struct sink* sink, - const struct sink_props* sink_props, - const struct engine* e, - int star_counter) { - - float random_number = random_unit_interval_part_ID_and_index( - sink->id, star_counter, e->ti_current, random_number_sink_formation); + const struct sink* const sink, const struct sink_props* sink_properties, + const int with_cosmology, const struct cosmology* cosmo, + const struct gravity_props* grav_props, const double time, + const double time_base) { + + /* Background sink particles have no time-step limits */ + if (sink->birth_time == -1.) { + return FLT_MAX; + } - const struct feedback_props* feedback_props = e->feedback_props; + /* CFL condition for sink. Notice the conversion to physical units ------- */ + const float CFL_condition = sink_properties->CFL_condition; + const double gas_v_phys[3] = { + sink->to_collect.velocity_gas[0] * cosmo->a_inv, + sink->to_collect.velocity_gas[1] * cosmo->a_inv, + sink->to_collect.velocity_gas[2] * cosmo->a_inv}; + double gas_v_norm2 = gas_v_phys[0] * gas_v_phys[0] + + gas_v_phys[1] * gas_v_phys[1] + + gas_v_phys[2] * gas_v_phys[2]; + + const double gas_c_phys = + sink->to_collect.sound_speed_gas * cosmo->a_factor_sound_speed; + const double gas_c_phys2 = gas_c_phys * gas_c_phys; + const float denominator = sqrtf(gas_c_phys2 + gas_v_norm2); + const float h_min = + cosmo->a * + min(sink->r_cut, sink->to_collect.minimal_h_gas * kernel_gamma); + float dt_cfl = 0.0; + + /* This case can happen if the sink is just born. */ + if (gas_v_norm2 == 0.0) { + dt_cfl = FLT_MAX; + } else { + dt_cfl = 2.f * CFL_condition * h_min / denominator; + } - /* Pick the correct table. (if only one table, threshold is < 0) */ - const float metal = - chemistry_get_sink_total_iron_mass_fraction_for_feedback(sink); - const float threshold = feedback_props->metallicity_max_first_stars; + /* Free fall time condition: the sink must anticipate gas collapse ------- */ + const float rho_sink = + 3.0 * sink->mass / (4.0 * M_PI * h_min * h_min * h_min); + const float dt_ff = + sqrtf(3.0 * M_PI / (32.0 * grav_props->G_Newton * rho_sink)); - /* If metal < threshold, then the sink generates first star particles. */ - const int is_first_star = metal < threshold; - const struct stellar_model* model; - double minimal_discrete_mass_Msun; + /* Compute sink-sink orbital integration timestep ------------------------ */ + float dt_2_body = 0.0; - /* Take the correct values if your are a first star or not. */ - if (!is_first_star) /* (metal >= threshold)*/ { - model = &feedback_props->stellar_model; - minimal_discrete_mass_Msun = sink_props->minimal_discrete_mass_Msun; + /* If there are no sink neighbours, then the values are FLT_MAX. Prevent + giving a NaN to the timestep */ + if ((sink->to_collect.minimal_sink_t_c == FLT_MAX) || + (sink->to_collect.minimal_sink_t_dyn == FLT_MAX)) { + dt_2_body = FLT_MAX; } else { - model = &feedback_props->stellar_model_first_stars; - minimal_discrete_mass_Msun = - sink_props->minimal_discrete_mass_first_stars_Msun; + dt_2_body = sink->to_collect.minimal_sink_t_c * + sink->to_collect.minimal_sink_t_dyn / + (sink->to_collect.minimal_sink_t_c + + sink->to_collect.minimal_sink_t_dyn); } - const struct initial_mass_function* imf = &model->imf; + /* SF - accretion timestep ------------------------------------------------*/ + /* Now, limit timestep by computing how much we restricted the sink accretion + for SF reasons compared to an unrestricted accretion. + Note: If we divide by mass_eligible_swallow, we get the relative error + compared to unrestricted swallow. */ + const float Delta_M = + sink->to_collect.mass_swallowed - sink->to_collect.mass_eligible_swallow; + + /* We want a big timestep if the error is small. */ + float dt_SF = FLT_MAX; + + /* If Delta_M < 0, then we are limiting the accretion rate by a huge factor. + To avoid biasing the SFR too much, do a small timestep to accrete the + remaining mass sooner. */ + if (sink_properties->n_IMF > 0 && Delta_M < 0) { + /* Compute an accretion rate using this Delta_M. Use the minmal timestep + based on the local gas properties. If we use the current timestep, then we + can end up with timesteps smaller and smaller until they are smaller than + the minimal engine timestep. */ + const float dt_tmp = min3(dt_cfl, dt_ff, dt_2_body); + const float M_dot = Delta_M / dt_tmp; + dt_SF = sink_properties->tolerance_SF_timestep * + sink->to_collect.mass_eligible_swallow / fabsf(M_dot); + } + + /* Sink age (in internal units) */ + double sink_age = sink_get_sink_age(sink, with_cosmology, cosmo, time); + + /* Take the minimum dt --------------------------------------------------- */ + float dt = min3(dt_cfl, dt_ff, dt_SF); - if (random_number < imf->sink_Pc) { - /* We are dealing with the continous part of the IMF. */ - sink->target_mass_Msun = imf->stellar_particle_mass_Msun; - sink->target_type = star_population_continuous_IMF; + /* What age category are we in? */ + if (sink_age > sink_properties->age_threshold_unlimited) { + return dt_2_body; + } else if (sink_age > sink_properties->age_threshold) { + dt = min3(dt, dt_2_body, sink_properties->max_time_step_old); + return dt; } else { - /* We are dealing with the discrete part of the IMF. */ - random_number = random_unit_interval_part_ID_and_index( - sink->id, star_counter + 1, e->ti_current, - random_number_sink_formation); - double m = initial_mass_function_sample_power_law( - minimal_discrete_mass_Msun, imf->mass_max, imf->exp[imf->n_parts - 1], - random_number); - sink->target_mass_Msun = m; - sink->target_type = single_star; + dt = min3(dt, dt_2_body, sink_properties->max_time_step_young); + return dt; } } @@ -144,6 +190,18 @@ __attribute__((always_inline)) INLINE static void sink_first_init_sink( /* Initialize to the mass of the sink */ sp->mass_tot_before_star_spawning = sp->mass; + + /* Init properties based on the local gas */ + sp->to_collect.minimal_h_gas = FLT_MAX; + sp->to_collect.rho_gas = 0.0; + sp->to_collect.sound_speed_gas = 0.0; + sp->to_collect.velocity_gas[0] = 0.0; + sp->to_collect.velocity_gas[1] = 0.0; + sp->to_collect.velocity_gas[2] = 0.0; + sp->to_collect.minimal_sink_t_c = FLT_MAX; + sp->to_collect.minimal_sink_t_dyn = FLT_MAX; + sp->to_collect.mass_eligible_swallow = 0.0; + sp->to_collect.mass_swallowed = sp->mass; } /** @@ -172,7 +230,15 @@ __attribute__((always_inline)) INLINE static void sink_init_part( cpd->E_rot_neighbours[0] = 0.f; cpd->E_rot_neighbours[1] = 0.f; cpd->E_rot_neighbours[2] = 0.f; - cpd->potential = 0.f; + + /* Do not reset the potential to 0. Keep the value computed at the end of the + last step. This value is used in runner_iact_nonsym_sink() and + runner_iact_sink() to check which particle is at a potential minimum. If you + set this value to 0, then we break the check. This value is used instead of + gpart->potential because: + 1) cpd->potential does not break MPI, while gpart->potential does + 2) gpart->potential is not yet computed in runner_iact_X_sink(). */ + /* cpd->potential = 0.f; */ cpd->E_mec_bound = 0.f; /* Gravitationally bound particles will have E_mec_bound < 0. This is checked before comparing any other value with this one. So no need to put @@ -192,6 +258,17 @@ __attribute__((always_inline)) INLINE static void sink_init_sink( /* Reset to the mass of the sink */ sp->mass_tot_before_star_spawning = sp->mass; + /* Init properties based on the local gas */ + sp->to_collect.minimal_h_gas = FLT_MAX; + sp->to_collect.rho_gas = 0.0; + sp->to_collect.sound_speed_gas = 0.0; + sp->to_collect.velocity_gas[0] = 0.0; + sp->to_collect.velocity_gas[1] = 0.0; + sp->to_collect.velocity_gas[2] = 0.0; + sp->to_collect.minimal_sink_t_c = FLT_MAX; + sp->to_collect.minimal_sink_t_dyn = FLT_MAX; + sp->to_collect.mass_eligible_swallow = 0.0; + sp->to_collect.mass_swallowed = sp->mass; sp->num_ngbs = 0; #ifdef DEBUG_INTERACTIONS_SINKS @@ -243,7 +320,23 @@ __attribute__((always_inline)) INLINE static void sink_kick_extra( * @param cosmo The current cosmological model. */ __attribute__((always_inline)) INLINE static void sink_end_density( - struct sink* si, const struct cosmology* cosmo) {} + struct sink* si, const struct cosmology* cosmo) { + + const float h = si->r_cut / kernel_gamma; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + + /* --- Finish the calculation by inserting the missing h factors --- */ + si->to_collect.rho_gas *= h_inv_dim; + const float rho_inv = 1.f / si->to_collect.rho_gas; + + /* For the following, we also have to undo the mass smoothing + * (N.B.: bp->velocity_gas is in BH frame, in internal units). */ + si->to_collect.sound_speed_gas *= h_inv_dim * rho_inv; + si->to_collect.velocity_gas[0] *= h_inv_dim * rho_inv; + si->to_collect.velocity_gas[1] *= h_inv_dim * rho_inv; + si->to_collect.velocity_gas[2] *= h_inv_dim * rho_inv; +} /** * @brief Sets all particle fields to sensible values when the #sink has 0 @@ -253,7 +346,19 @@ __attribute__((always_inline)) INLINE static void sink_end_density( * @param cosmo The current cosmological model. */ __attribute__((always_inline)) INLINE static void sinks_sink_has_no_neighbours( - struct sink* restrict sp, const struct cosmology* cosmo) {} + struct sink* restrict sp, const struct cosmology* cosmo) { + + warning( + "Sink particle with ID %lld treated as having no neighbours (r_cut: %g, " + "numb_ngbs: %i).", + sp->id, sp->r_cut, sp->num_ngbs); + + /* Reset problematic values */ + sp->to_collect.velocity_gas[0] = sp->v[0]; + sp->to_collect.velocity_gas[1] = sp->v[1]; + sp->to_collect.velocity_gas[2] = sp->v[2]; + sp->to_collect.minimal_h_gas = sp->r_cut / kernel_gamma; +} /** * @brief Compute the accretion rate of the sink and any quantities @@ -279,65 +384,6 @@ __attribute__((always_inline)) INLINE static void sink_prepare_swallow( const struct entropy_floor_properties* floor_props, const double time, const int with_cosmology, const double dt, const integertime_t ti_begin) {} -/** - * @brief Compute the rotational energy of the neighbouring gas particles. - * - * Note: This function must be used after having computed the rotational energy - * per components, i.e. after sink_prepare_part_sink_formation(). - * - * @param p The gas particle. - * - */ -INLINE static double sink_compute_neighbour_rotation_energy_magnitude( - const struct part* restrict p) { - double E_rot_x = p->sink_data.E_rot_neighbours[0]; - double E_rot_y = p->sink_data.E_rot_neighbours[1]; - double E_rot_z = p->sink_data.E_rot_neighbours[2]; - double E_rot = - sqrtf(E_rot_x * E_rot_x + E_rot_y * E_rot_y + E_rot_z * E_rot_z); - return E_rot; -} - -/** - * @brief Retrieve the physical velocity divergence from the gas particle. - * - * @param p The gas particles. - * - */ -INLINE static float sink_get_physical_div_v_from_part( - const struct part* restrict p) { - - float div_v = 0.0; - - /* The implementation of div_v depends on the Hydro scheme. Furthermore, some - add a Hubble flow term, some do not. We need to take care of this */ -#ifdef SPHENIX_SPH - /* SPHENIX is already including the Hubble flow. */ - div_v = hydro_get_div_v(p); -#elif GADGET2_SPH - div_v = p->density.div_v; - - /* Add the missing term */ - div_v += hydro_dimension * cosmo->H; -#elif MINIMAL_SPH - div_v = hydro_get_div_v(p); - - /* Add the missing term */ - div_v += hydro_dimension * cosmo->H; -#elif GASOLINE_SPH - /* Copy the velocity divergence */ - div_v = (1. / 3.) * (p->viscosity.velocity_gradient[0][0] + - p->viscosity.velocity_gradient[1][1] + - p->viscosity.velocity_gradient[2][2]); -#elif HOPKINS_PU_SPH - div_v = p->density.div_v; -#else -#error \ - "This scheme is not implemented. Note that Different scheme apply the Hubble flow in different places. Be careful about it." -#endif - return div_v; -} - /** * @brief Calculate if the gas has the potential of becoming * a sink. @@ -526,6 +572,10 @@ INLINE static void sink_copy_properties( /* Note, we do not need to update sp->mass_tot_before_star_spawning because it is performed within the 'sink_init_sink()' function. */ + + /* Set the birth time of the sink */ + sink_set_sink_birth_time_or_scale_factor(sink, e->time, cosmo->a, + with_cosmology); } /** @@ -675,10 +725,13 @@ __attribute__((always_inline)) INLINE static void sink_swallow_sink( /* Add the stars spawned by the swallowed sink */ spi->n_stars += spj->n_stars; - /* Update the total mass before star spawning */ + /* Update masses */ spi->mass_tot_before_star_spawning = spi->mass; + spi->to_collect.mass_eligible_swallow += + spj->to_collect.mass_eligible_swallow; + spi->to_collect.mass_swallowed += spj->to_collect.mass_swallowed; - message("sink %lld swallow sink particle %lld. New mass: %e.", spi->id, + message("sink %lld swallows sink particle %lld. New mass: %e.", spi->id, spj->id, spi->mass); } @@ -699,8 +752,21 @@ INLINE static int sink_spawn_star(struct sink* sink, const struct engine* e, const int with_cosmology, const struct phys_const* phys_const, const struct unit_system* restrict us) { - - if (sink->mass > sink->target_mass_Msun * phys_const->const_solar_mass) + /* Convenient variables in internal units */ + const float target_mass = + sink->target_mass_Msun * phys_const->const_solar_mass; + const float minimal_mass = + sink_props->sink_minimal_mass_Msun * phys_const->const_solar_mass; + + /* To spawn a star, the sink must: + 1) m_sink > target_mass, + 2) and m_sink - target_mass >= minimal_sink_mass mass. + The second condition is relevant for low resolution simulations, where a + new born sink can already spawn stars. After spawning the stars, the sink + has a mass << gas mass and can get kicked away by gravitational + interactions. Also, if the sink's mass << gas' mass, the gas is never bound + to the sink and is thus never accreted. */ + if (sink->mass > target_mass && (sink->mass - target_mass >= minimal_mass)) return 1; else return 0; @@ -1116,13 +1182,25 @@ INLINE static void sink_prepare_part_sink_formation_gas_criteria( */ INLINE static void sink_prepare_part_sink_formation_sink_criteria( struct engine* e, struct part* restrict p, struct xpart* restrict xp, - struct sink* restrict si, const struct cosmology* cosmo, - const struct sink_props* sink_props) { + struct sink* restrict si, const int with_cosmology, + const struct cosmology* cosmo, const struct sink_props* sink_props, + const double time) { + /* Do not continue if the gas cannot form sink for any reason */ if (!p->sink_data.can_form_sink) { return; } + /* Determine if the sink is dead, i.e. if its age is bigger than the + age_threshold_unlimited */ + const int sink_age = sink_get_sink_age(si, with_cosmology, cosmo, time); + char is_dead = sink_age > sink_props->age_threshold_unlimited; + + /* If the sink is dead, do not check the criteria for the si - p pair. */ + if (is_dead) { + return; + } + /* Physical accretion radius of part p */ const float r_acc_p = sink_props->cut_off_radius * cosmo->a; diff --git a/src/sink/GEAR/sink_getters.h b/src/sink/GEAR/sink_getters.h new file mode 100644 index 0000000000000000000000000000000000000000..d2f0a1738407a39aa6b8d4ba5ab1dd9fae2f1fa9 --- /dev/null +++ b/src/sink/GEAR/sink_getters.h @@ -0,0 +1,171 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Darwin Roduit (darwin.roduit@alumni.epfl.ch) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + *******************************************************************************/ +#ifndef SWIFT_GEAR_SINK_GETTERS_H +#define SWIFT_GEAR_SINK_GETTERS_H + +#include "cosmology.h" +#include "sink_part.h" + +/** + * @file src/sink/GEAR/sink_getter.h + * @brief Getter functions for GEAR sink scheme to avoid exposing + * implementation details to the outer world. Keep the code clean and lean. + */ + +/** + * @brief Get the sink age in interal units. + * + * @param sink The #sink. + * @param with_cosmology If we run with cosmology. + * @param cosmo The co Birth scale-factor of the star. + * @param with_cosmology If we run with cosmology. + */ + +__attribute__((always_inline)) INLINE double sink_get_sink_age( + const struct sink* restrict sink, const int with_cosmology, + const struct cosmology* cosmo, const double time) { + double sink_age; + if (with_cosmology) { + + /* Deal with rounding issues */ + if (sink->birth_scale_factor >= cosmo->a) { + sink_age = 0.; + } else { + sink_age = cosmology_get_delta_time_from_scale_factors( + cosmo, sink->birth_scale_factor, cosmo->a); + } + } else { + sink_age = time - sink->birth_time; + } + return sink_age; +} + +/** + * @brief Compute the rotational energy of the neighbouring gas particles. + * + * Note: This function must be used after having computed the rotational energy + * per components, i.e. after sink_prepare_part_sink_formation(). + * + * @param p The gas particle. + * + */ +INLINE static double sink_compute_neighbour_rotation_energy_magnitude( + const struct part* restrict p) { + double E_rot_x = p->sink_data.E_rot_neighbours[0]; + double E_rot_y = p->sink_data.E_rot_neighbours[1]; + double E_rot_z = p->sink_data.E_rot_neighbours[2]; + double E_rot = + sqrtf(E_rot_x * E_rot_x + E_rot_y * E_rot_y + E_rot_z * E_rot_z); + return E_rot; +} + +/** + * @brief Retrieve the physical velocity divergence from the gas particle. + * + * @param p The gas particles. + * + */ +INLINE static float sink_get_physical_div_v_from_part( + const struct part* restrict p) { + + float div_v = 0.0; + + /* The implementation of div_v depends on the Hydro scheme. Furthermore, some + add a Hubble flow term, some do not. We need to take care of this */ +#ifdef SPHENIX_SPH + /* SPHENIX is already including the Hubble flow. */ + div_v = hydro_get_div_v(p); +#elif GADGET2_SPH + div_v = p->density.div_v; + + /* Add the missing term */ + div_v += hydro_dimension * cosmo->H; +#elif MINIMAL_SPH + div_v = hydro_get_div_v(p); + + /* Add the missing term */ + div_v += hydro_dimension * cosmo->H; +#elif GASOLINE_SPH + /* Copy the velocity divergence */ + div_v = (1. / 3.) * (p->viscosity.velocity_gradient[0][0] + + p->viscosity.velocity_gradient[1][1] + + p->viscosity.velocity_gradient[2][2]); +#elif HOPKINS_PU_SPH + div_v = p->density.div_v; +#else +#error \ + "This scheme is not implemented. Note that Different scheme apply the Hubble flow in different places. Be careful about it." +#endif + return div_v; +} + +/** + * @brief Compute the angular momentum-based criterion for sink-sink + * interaction. + * + * This function calculates the angular momentum of a sink particle relative to + * another particle (sink or gas) and evaluates the Keplerian angular momentum. + * + * @param dx Comoving vector separating the two particles (pi - pj). + * @param dv_plus_H_flow Comoving relative velocity including the Hubble flow. + * @param r Comoving distance between the two particles. + * @param r_cut_i Comoving cut-off radius of particle i. + * @param mass_i Mass of particle i. + * @param cosmo The cosmological parameters and properties + * @param grav_props The gravity scheme parameters and properties + * @param L2_kepler (return) Keplerian angular momentum squared of particle i. + * @param L2_j (return) Specific angular momentum squared relative to particle + * j. + */ +__attribute__((always_inline)) INLINE static void +sink_compute_angular_momenta_criterion( + const float dx[3], const float dv_plus_H_flow[3], const float r, + const float r_cut_i, const float mass_i, const struct cosmology* cosmo, + const struct gravity_props* grav_props, float* L2_kepler, float* L2_j) { + + /* Compute the physical relative velocity between the particles */ + const float dv_physical[3] = {dv_plus_H_flow[0] * cosmo->a_inv, + dv_plus_H_flow[1] * cosmo->a_inv, + dv_plus_H_flow[2] * cosmo->a_inv}; + + /* Compute the physical distance between the particles */ + const float dx_physical[3] = {dx[0] * cosmo->a, dx[1] * cosmo->a, + dx[2] * cosmo->a}; + const float r_physical = r * cosmo->a; + + /* Momentum check------------------------------------------------------- */ + /* Relative momentum of the gas */ + const float specific_angular_momentum[3] = { + dx_physical[1] * dv_physical[2] - dx_physical[2] * dv_physical[1], + dx_physical[2] * dv_physical[0] - dx_physical[0] * dv_physical[2], + dx_physical[0] * dv_physical[1] - dx_physical[1] * dv_physical[0]}; + + *L2_j = specific_angular_momentum[0] * specific_angular_momentum[0] + + specific_angular_momentum[1] * specific_angular_momentum[1] + + specific_angular_momentum[2] * specific_angular_momentum[2]; + + /* Keplerian angular speed squared */ + const float omega_acc_2 = + grav_props->G_Newton * mass_i / (r_physical * r_physical * r_physical); + + /*Keplerian angular momentum squared */ + *L2_kepler = (r_cut_i * r_cut_i * r_cut_i * r_cut_i) * omega_acc_2; +} + +#endif /* SWIFT_GEAR_SINK_GETTERS_H */ diff --git a/src/sink/GEAR/sink_iact.h b/src/sink/GEAR/sink_iact.h index fc75ba8842bb7ceaae6d789b277a76b21206c29f..b2c0dd86f1144505a181e1304c2f0919adc84792 100644 --- a/src/sink/GEAR/sink_iact.h +++ b/src/sink/GEAR/sink_iact.h @@ -24,6 +24,7 @@ #include "gravity.h" #include "gravity_iact.h" #include "sink.h" +#include "sink_getters.h" #include "sink_properties.h" /** @@ -33,8 +34,6 @@ * In GEAR: This function deactivates the sink formation ability of #part not * at a potential minimum. * - * Note: This functions breaks MPI. - * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). * @param hi Comoving smoothing-length of particle i. @@ -58,12 +57,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_sink( /* if the distance is less than the cut off radius */ if (r < cut_off_radius) { - - /* - * NOTE: Those lines break MPI - */ - float potential_i = pi->gpart->potential; - float potential_j = pj->gpart->potential; + float potential_i = pi->sink_data.potential; + float potential_j = pj->sink_data.potential; /* prevent the particle with the largest potential to form a sink */ if (potential_i > potential_j) { @@ -107,9 +102,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink( const float r = sqrtf(r2); if (r < cut_off_radius) { - - float potential_i = pi->gpart->potential; - float potential_j = pj->gpart->potential; + float potential_i = pi->sink_data.potential; + float potential_j = pj->sink_data.potential; /* if the potential is larger * prevent the particle to form a sink */ @@ -117,12 +111,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink( } } -/** - * @brief Density interaction between two particles (non-symmetric). - * - * @param r2 Comoving square distance between the two particles. - * @param dx Comoving vector separating both particles (pi - pj). - * @param ri Comoving cut off radius of particle i. +/* @brief Density interaction between two particles (non-symmetric). * @param hj Comoving smoothing-length of particle j. * @param si First particle (sink). * @param pj Second particle (gas, not updated). @@ -141,13 +130,93 @@ runner_iact_nonsym_sinks_gas_density( const struct sink_props *sink_props, const integertime_t ti_current, const double time) { - /* Get r. */ + /* Contribution to the number of neighbours in cutoff radius */ + si->num_ngbs++; + + float wi, wi_dx; + + /* Compute the kernel function */ const float r = sqrtf(r2); + const float hi = ri / kernel_gamma; + const float hi_inv = 1.0f / hi; + const float ui = r * hi_inv; + kernel_deval(ui, &wi, &wi_dx); + + /* Neighbour gas mass */ + const float mj = hydro_get_mass(pj); + + /* Minimum smoothing length accros the neighbours */ + /* AND the sink smoothing length */ + si->to_collect.minimal_h_gas = min(hj, si->to_collect.minimal_h_gas); + + /* Contribution to the BH gas density */ + si->to_collect.rho_gas += mj * wi; + + /* Contribution to the smoothed sound speed */ + si->to_collect.sound_speed_gas += mj * wi * hydro_get_comoving_soundspeed(pj); + + /* Neighbour's (drifted) velocity in the frame of the sink + * (we don't include a Hubble term since we are interested in the + * velocity contribution at the location of the sink) */ + const float dv[3] = {pj->v[0] - si->v[0], pj->v[1] - si->v[1], + pj->v[2] - si->v[2]}; + + /* Contribution to the smoothed velocity (gas w.r.t. black hole) */ + si->to_collect.velocity_gas[0] += mj * dv[0] * wi; + si->to_collect.velocity_gas[1] += mj * dv[1] * wi; + si->to_collect.velocity_gas[2] += mj * dv[2] * wi; +} - if (r < ri) { - /* Contribution to the number of neighbours in cutoff radius */ - si->num_ngbs++; - } +/** + * @brief Update the properties of a sink particles from its sink neighbours. + * + * Warning: No symmetric interaction for timesteps since we are in + * runner_iact_nonsym_sinks_sink_swallow() --> pay attention to not break MPI. + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param ri Comoving cut off radius of particle i. + * @param rj Comoving cut off radius of particle j. + * @param si First sink particle. + * @param sj Second sink particle. + */ +__attribute__((always_inline)) INLINE static void +sink_collect_properties_from_sink(const float r2, const float dx[3], + const float ri, const float rj, + struct sink *restrict si, + struct sink *restrict sj, + const struct gravity_props *grav_props) { + + /* Neighbour's (drifted) velocity in the frame of the sink i + * (we don't include a Hubble term since we are interested in the + * velocity contribution at the location of the sink) */ + const float dv[3] = {sj->v[0] - si->v[0], sj->v[1] - si->v[1], + sj->v[2] - si->v[2]}; + const float dv_norm = sqrtf(dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]); + + /* Get the gravitional softening */ + const float eps = gravity_get_softening(si->gpart, grav_props); + const float eps2 = eps * eps; + const float eps_inv = 1.f / eps; + const float eps_inv3 = eps_inv * eps_inv * eps_inv; + + /* Compute the kernel potential and force with mass = 1.0. We multiply by + the mass below if needed. */ + float dphi_dr, pot; + runner_iact_grav_pp_full(r2, eps2, eps_inv, eps_inv3, 1.0, &dphi_dr, &pot); + + /* From Grudic et al. (2021) eq 6, we replace the plummer functionnal form + sqrt(r^2 + eps^2) by the kernel 1.0/|phi(r,H=3*eps)| */ + const float t_c = 1.0 / (fabsf(pot) * dv_norm); + si->to_collect.minimal_sink_t_c = min(t_c, si->to_collect.minimal_sink_t_c); + + /* From Grudic et al. (2021) eq 7, we replace the plummer functionnal form + (r^2 + eps^2)^{3/2} by the kernel |(d phi(r,H=3*eps)/ dr)^{-1}| */ + const float denominator = grav_props->G_Newton * (si->mass + sj->mass); + const float numerator = 1.0 / fabsf(dphi_dr); + const float t_dyn = sqrt(numerator / denominator); + si->to_collect.minimal_sink_t_dyn = + min(t_dyn, si->to_collect.minimal_sink_t_dyn); } /** @@ -155,6 +224,9 @@ runner_iact_nonsym_sinks_gas_density( * * Note: Energies are computed with physical quantities, not the comoving ones. * + * MPI note: This functions invokes the gpart. Hence, it must be performed only + * on the local node (similarly to runner_iact_nonsym_bh_bh_repos()). + * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). * @param ri Comoving cut off radius of particle i. @@ -180,6 +252,27 @@ runner_iact_nonsym_sinks_sink_swallow( const float r = sqrtf(r2); const float f_acc_r_acc_i = sink_properties->f_acc * ri; + /* Determine if the sink is dead, i.e. if its age is bigger than the + age_threshold_unlimited */ + const int si_age = sink_get_sink_age(si, with_cosmology, cosmo, time); + char si_is_dead = si_age > sink_properties->age_threshold_unlimited; + + const int sj_age = sink_get_sink_age(sj, with_cosmology, cosmo, time); + char sj_is_dead = sj_age > sink_properties->age_threshold_unlimited; + + /* Collect the properties for 2-body interactions if one is sink is alive. If + they are both dead, we do not want to restrict the timesteps for 2-body + encounters since they won't merge. */ + if (!si_is_dead || !sj_is_dead) { + sink_collect_properties_from_sink(r2, dx, ri, rj, si, sj, grav_props); + } + + /* If si is dead, do not swallow sj. However, sj can swallow si if it alive. + */ + if (si_is_dead) { + return; + } + /* If the sink j falls within f_acc*r_acc of sink i, then the lightest is accreted on the most massive without further check. Note that this is a non-symmetric interaction. So, we do not need to check @@ -213,7 +306,6 @@ runner_iact_nonsym_sinks_sink_swallow( const float v_plus_H_flow[3] = {a2H * dx[0] + dv[0], a2H * dx[1] + dv[1], a2H * dx[2] + dv[2]}; - /* Binding energy check */ /* Compute the physical relative velocity between the particles */ const float dv_physical[3] = {v_plus_H_flow[0] * cosmo->a_inv, v_plus_H_flow[1] * cosmo->a_inv, @@ -223,7 +315,21 @@ runner_iact_nonsym_sinks_sink_swallow( dv_physical[1] * dv_physical[1] + dv_physical[2] * dv_physical[2]; - /* Kinetic energy per unit mass of the gas */ + /* Momentum check------------------------------------------------------- */ + float L2_j = 0.0; /* Relative momentum of the sink j */ + float L2_kepler = 0.0; /* Keplerian angular momentum squared */ + sink_compute_angular_momenta_criterion(dx, v_plus_H_flow, r, si->r_cut, + si->mass, cosmo, grav_props, + &L2_kepler, &L2_j); + + /* To be accreted, the sink momentum should lower than the keplerian orbit + * momentum. */ + if (L2_j > L2_kepler) { + return; + } + + /* Binding energy check------------------------------------------------- */ + /* Kinetic energy per unit mass of the sink */ const float E_kin_rel = 0.5f * dv_physical_squared; /* Compute the Newtonian or softened potential the sink exherts onto the @@ -256,6 +362,25 @@ runner_iact_nonsym_sinks_sink_swallow( return; } + /* Swallowed mass threshold--------------------------------------------- */ + si->to_collect.mass_eligible_swallow += sj->mass; + + /* Maximal mass that can be swallowed within a single timestep */ + const float mass_swallow_limit = sink_properties->n_IMF * si->mass_IMF; + + /* If the mass exceeds the threshold, do not swallow. Make sure you can at + least swallow a particle to avoid running into the problem of never being + able to spawn a star. + If n_IMF <= 0, then disable this criterion */ + if (sink_properties->n_IMF > 0 && + si->to_collect.mass_swallowed >= mass_swallow_limit && + si->to_collect.mass_eligible_swallow != 0) { + return; + } + + /* Increment the swallowd mass */ + si->to_collect.mass_swallowed += sj->mass; + /* The sink with the smaller mass will be merged onto the one with the * larger mass. * To avoid rounding issues, we additionally check for IDs if the sink @@ -286,6 +411,9 @@ runner_iact_nonsym_sinks_sink_swallow( * * Note: Energies are computed with physical quantities, not the comoving ones. * + * MPI note: This functions invokes the gpart. Hence, it must be performed only + * on the local node (similarly to runner_iact_nonsym_bh_gas_repos()). + * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). * @param ri Comoving cut off radius of particle i. @@ -311,9 +439,20 @@ runner_iact_nonsym_sinks_gas_swallow( const float r = sqrtf(r2); const float f_acc_r_acc = sink_properties->f_acc * ri; + /* Determine if the sink is dead, i.e. if its age is bigger than the + age_threshold_unlimited */ + const int sink_age = sink_get_sink_age(si, with_cosmology, cosmo, time); + char is_dead = sink_age > sink_properties->age_threshold_unlimited; + + /* If si is dead, do not swallow pj. */ + if (is_dead) { + return; + } + /* If the gas falls within f_acc*r_acc, it is accreted without further check */ if (r < f_acc_r_acc) { + warning("Gas %lld within sink %lld inner accretion radius", pj->id, si->id); /* Check if a gas particle has not been already marked to be swallowed by another sink particle. */ if (pj->sink_data.swallow_id < si->id) { @@ -344,33 +483,16 @@ runner_iact_nonsym_sinks_gas_swallow( dv_physical[1] * dv_physical[1] + dv_physical[2] * dv_physical[2]; - /* Compute the physical distance between the particles */ - const float dx_physical[3] = {dx[0] * cosmo->a, dx[1] * cosmo->a, - dx[2] * cosmo->a}; - const float r_physical = r * cosmo->a; - /* Momentum check------------------------------------------------------- */ - /* Relative momentum of the gas */ - const float specific_angular_momentum_gas[3] = { - dx_physical[1] * dv_physical[2] - dx_physical[2] * dv_physical[1], - dx_physical[2] * dv_physical[0] - dx_physical[0] * dv_physical[2], - dx_physical[0] * dv_physical[1] - dx_physical[1] * dv_physical[0]}; - const float L2_gas = - specific_angular_momentum_gas[0] * specific_angular_momentum_gas[0] + - specific_angular_momentum_gas[1] * specific_angular_momentum_gas[1] + - specific_angular_momentum_gas[2] * specific_angular_momentum_gas[2]; - - /* Keplerian angular speed squared */ - const float omega_acc_2 = grav_props->G_Newton * si->mass / - (r_physical * r_physical * r_physical); - - /*Keplerian angular momentum squared */ - const float L2_acc = - (si->r_cut * si->r_cut * si->r_cut * si->r_cut) * omega_acc_2; - - /* To be accreted, the gas momentum shoulb lower than the keplerian orbit + float L2_gas_j = 0.0; /* Relative momentum of the gas */ + float L2_kepler = 0.0; /* Keplerian angular momentum squared */ + sink_compute_angular_momenta_criterion(dx, v_plus_H_flow, r, si->r_cut, + si->mass, cosmo, grav_props, + &L2_kepler, &L2_gas_j); + + /* To be accreted, the gas momentum should lower than the keplerian orbit * momentum. */ - if (L2_gas > L2_acc) { + if (L2_gas_j > L2_kepler) { return; } @@ -405,13 +527,40 @@ runner_iact_nonsym_sinks_gas_swallow( /* To be accreted, the gas must be gravitationally bound to the sink. */ if (E_mec_sink_part >= 0) return; + /* To be accreted, the gas smoothing length must be smaller than the sink + accretion radius. This is similar to AMR codes requesting the maximum + refinement level close to the sink. */ + if (sink_properties->sink_formation_smoothing_length_criterion && + (pj->h * kernel_gamma >= si->r_cut)) + return; + /* Most bound pair check------------------------------------------------ */ /* The pair gas-sink must be the most bound among all sinks */ if (E_mec_sink_part >= pj->sink_data.E_mec_bound) { return; } - /* Since this pair gas-sink is the most bounf, keep track of the + /* Swallowed mass threshold--------------------------------------------- */ + si->to_collect.mass_eligible_swallow += hydro_get_mass(pj); + + /* Maximal mass that can be swallowed within a single timestep */ + const float mass_swallow_limit = sink_properties->n_IMF * si->mass_IMF; + + /* If the mass exceeds the threshold, do not swallow. Make sure you can at + least swallow a particle to avoid running into the problem of never being + able to spawn a star. + If n_IMF <= 0, then disable this criterion */ + if (sink_properties->n_IMF > 0 && + si->to_collect.mass_swallowed >= mass_swallow_limit && + si->to_collect.mass_eligible_swallow != 0) { + return; + } + + /* Increment the swallowd mass */ + si->to_collect.mass_swallowed += hydro_get_mass(pj); + + /* --------------------------------------------------------------------- */ + /* Since this pair gas-sink is the most bound, keep track of the E_mec_bound and set the swallow_id accordingly */ pj->sink_data.E_mec_bound = E_mec_sink_part; pj->sink_data.swallow_id = si->id; diff --git a/src/sink/GEAR/sink_io.h b/src/sink/GEAR/sink_io.h index be9da875c5f27275f431124449a784af5b82aa53..293d58d50175657c155701c07b28ddd9da61ad4d 100644 --- a/src/sink/GEAR/sink_io.h +++ b/src/sink/GEAR/sink_io.h @@ -34,7 +34,7 @@ INLINE static void sink_read_particles(struct sink* sinks, struct io_props* list, int* num_fields) { /* Say how much we want to read */ - *num_fields = 4; + *num_fields = 5; /* List what we want to read */ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY, @@ -45,6 +45,8 @@ INLINE static void sink_read_particles(struct sink* sinks, sinks, mass); list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY, UNIT_CONV_NO_UNITS, sinks, id); + list[4] = io_make_input_field("BirthTime", FLOAT, 1, OPTIONAL, UNIT_CONV_MASS, + sinks, birth_time); } INLINE static void convert_sink_pos(const struct engine* e, @@ -123,7 +125,7 @@ INLINE static void sink_write_particles(const struct sink* sinks, int with_cosmology) { /* Say how much we want to write */ - *num_fields = 9; + *num_fields = 10; /* List what we want to write */ list[0] = io_make_output_field_convert_sink( @@ -170,6 +172,17 @@ INLINE static void sink_write_particles(const struct sink* sinks, /*can convert to comoving=*/0, convert_sink_swallowed_angular_momentum, "Physical swallowed angular momentum of the particles"); + if (with_cosmology) { + list[9] = io_make_physical_output_field( + "BirthScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, + birth_scale_factor, /*can convert to comoving=*/0, + "Scale-factors at which the sinks were born"); + } else { + list[9] = + io_make_output_field("BirthTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f, sinks, + birth_time, "Times at which the sinks were born"); + } + #ifdef DEBUG_INTERACTIONS_SINKS list += *num_fields; diff --git a/src/sink/GEAR/sink_part.h b/src/sink/GEAR/sink_part.h index be7f4359972d38570d944f1be8cae187fb6d83f1..80dded0ef11eec301d93af53bf309f8b2a63e8b1 100644 --- a/src/sink/GEAR/sink_part.h +++ b/src/sink/GEAR/sink_part.h @@ -55,6 +55,9 @@ struct sink { /*! Sink target mass. In Msun. */ float target_mass_Msun; + /* Mass of the IMF this sinks is currently affected to. In internal units. */ + double mass_IMF; + /*! Integer number of neighbours */ int num_ngbs; @@ -64,6 +67,44 @@ struct sink { /*! Sink target stellar type */ enum stellar_type target_type; + /*! Union for the birth time and birth scale factor */ + union { + + /*! Birth time */ + float birth_time; + + /*! Birth scale factor */ + float birth_scale_factor; + }; + + struct { + + /*! Minimal gas smoothing length */ + float minimal_h_gas; + + /*! Density of the gas surrounding the sink. */ + float rho_gas; + + /*! Smoothed sound speed of the gas surrounding the sink. */ + float sound_speed_gas; + + /*! Smoothed velocity of the gas surrounding the sink, in the frame of the + sink (internal units) */ + float velocity_gas[3]; + + /*! Minimal t_c between all sink neighbours */ + float minimal_sink_t_c; + + /*! Minimal dynamical time between all sink neighbours */ + float minimal_sink_t_dyn; + + /*! Total mass that passes all criteria before the accretion limit */ + float mass_eligible_swallow; + + /*! Swallowed mass during this timestep */ + float mass_swallowed; + } to_collect; + /*! Particle time bin */ timebin_t time_bin; diff --git a/src/sink/GEAR/sink_properties.h b/src/sink/GEAR/sink_properties.h index 251a29ed1ba28d2796aafe3ecb9726d55774e2cb..5b60b61d369fbbc3c1b0554278b766e3a93cd4a9 100644 --- a/src/sink/GEAR/sink_properties.h +++ b/src/sink/GEAR/sink_properties.h @@ -24,6 +24,18 @@ #include "feedback_properties.h" #include "parser.h" +/* Some default values for the parameters to be read in the YAML file */ +#define sink_gear_f_acc_default 0.8 +#define sink_gear_star_spawning_sigma_factor_default 0.2 +#define sink_gear_n_imf_default FLT_MAX /* No accretion restriction */ +#define sink_gear_tolerance_sf_timestep_default 0.5 + +/* Sink formation is activated */ +#define sink_gear_disable_sink_formation_default 0 + +/* By default all current implemented criteria are active */ +#define sink_gear_sink_formation_criterion_all_default 1 + /** * @brief Properties of sink in the Default model. */ @@ -67,13 +79,40 @@ struct sink_props { char sink_formation_bound_state_criterion; char sink_formation_overlapping_sink_criterion; - /* Disable sink formation? (e.g. used in sink accretion tests). Default: 0 + /*! Disable sink formation? (e.g. used in sink accretion tests). Default: 0 (keep sink formation) */ - uint8_t disable_sink_formation; + char disable_sink_formation; - /* Factor to rescale the velocity dispersion of the stars when they are + /*! Factor to rescale the velocity dispersion of the stars when they are spawned */ double star_spawning_sigma_factor; + + /*! Minimal sink mass in Msun. This prevents m_sink << m_gas in low + resolution simulations. */ + float sink_minimal_mass_Msun; + + /***************************************************************************/ + /*! Maximal time-step length of young sinks (internal units) */ + double max_time_step_young; + + /*! Maximal time-step length of old sinks (internal units) */ + double max_time_step_old; + + /*! Age threshold for the young/old transition (internal units) */ + double age_threshold; + + /*! Age threshold for the transition to unlimited time-step size (internal + * units) */ + double age_threshold_unlimited; + + /*! Time integration CFL condition factor */ + float CFL_condition; + + /*! Number of times the IMF mass can be swallowed in a single timestep */ + float n_IMF; + + /*! Tolerance parameter for SF timestep constraint */ + float tolerance_SF_timestep; }; /** @@ -175,21 +214,12 @@ INLINE static void sink_props_init(struct sink_props *sp, "ERROR: Running with sink but without feedback. GEAR sink model needs " "to be run with --sink and --feedback"); - /* Default values */ - const float default_f_acc = 0.8; - const float default_star_spawning_sigma_factor = 0.2; - const char default_disable_sink_formation = 0; /* Sink formation is - activated */ - - /* By default all current implemented criteria are active */ - const uint8_t default_sink_formation_criterion_all = 1; - /* Read the parameters from the parameter file */ sp->cut_off_radius = parser_get_param_float(params, "GEARSink:cut_off_radius"); - sp->f_acc = - parser_get_opt_param_float(params, "GEARSink:f_acc", default_f_acc); + sp->f_acc = parser_get_opt_param_float(params, "GEARSink:f_acc", + sink_gear_f_acc_default); /* Check that sp->f_acc respects 0 <= f_acc <= 1 */ if ((sp->f_acc < 0) || (sp->f_acc > 1)) { @@ -228,33 +258,66 @@ INLINE static void sink_props_init(struct sink_props *sp, sp->star_spawning_sigma_factor = parser_get_opt_param_float(params, "GEARSink:star_spawning_sigma_factor", - default_star_spawning_sigma_factor); + sink_gear_star_spawning_sigma_factor_default); + + sp->n_IMF = parser_get_opt_param_float(params, "GEARSink:n_IMF", + sink_gear_n_imf_default); + + sp->sink_minimal_mass_Msun = + parser_get_opt_param_float(params, "GEARSink:sink_minimal_mass_Msun", 0.); /* Sink formation criterion parameters (all active by default) */ sp->sink_formation_contracting_gas_criterion = parser_get_opt_param_int( params, "GEARSink:sink_formation_contracting_gas_criterion", - default_sink_formation_criterion_all); + sink_gear_sink_formation_criterion_all_default); sp->sink_formation_smoothing_length_criterion = parser_get_opt_param_int( params, "GEARSink:sink_formation_smoothing_length_criterion", - default_sink_formation_criterion_all); + sink_gear_sink_formation_criterion_all_default); sp->sink_formation_jeans_instability_criterion = parser_get_opt_param_int( params, "GEARSink:sink_formation_jeans_instability_criterion", - default_sink_formation_criterion_all); + sink_gear_sink_formation_criterion_all_default); sp->sink_formation_bound_state_criterion = parser_get_opt_param_int( params, "GEARSink:sink_formation_bound_state_criterion", - default_sink_formation_criterion_all); + sink_gear_sink_formation_criterion_all_default); sp->sink_formation_overlapping_sink_criterion = parser_get_opt_param_int( params, "GEARSink:sink_formation_overlapping_sink_criterion", - default_sink_formation_criterion_all); + sink_gear_sink_formation_criterion_all_default); /* Should we disable sink formation ? */ sp->disable_sink_formation = parser_get_opt_param_int(params, "GEARSink:disable_sink_formation", - default_disable_sink_formation); + sink_gear_disable_sink_formation_default); + + /* Maximal time-step lengths */ + const double max_time_step_young_Myr = parser_get_opt_param_float( + params, "GEARSink:max_timestep_young_Myr", FLT_MAX); + const double max_time_step_old_Myr = parser_get_opt_param_float( + params, "GEARSink:max_timestep_old_Myr", FLT_MAX); + const double age_threshold_Myr = parser_get_opt_param_float( + params, "GEARSink:timestep_age_threshold_Myr", FLT_MAX); + const double age_threshold_unlimited_Myr = parser_get_opt_param_float( + params, "GEARSink:timestep_age_threshold_unlimited_Myr", FLT_MAX); + + /* Check for consistency */ + if (age_threshold_unlimited_Myr != 0. && age_threshold_Myr != FLT_MAX) { + if (age_threshold_unlimited_Myr < age_threshold_Myr) + error( + "The age threshold for unlimited sink time-step sizes (%e Myr) is " + "smaller than the transition threshold from young to old ages (%e " + "Myr)", + age_threshold_unlimited_Myr, age_threshold_Myr); + } + + /* Timestep tolerance paramters */ + sp->CFL_condition = parser_get_param_float(params, "GEARSink:CFL_condition"); + + sp->tolerance_SF_timestep = + parser_get_opt_param_float(params, "GEARSink:tolerance_SF_timestep", + sink_gear_tolerance_sf_timestep_default); /* Apply unit change */ sp->temperature_threshold /= @@ -267,6 +330,13 @@ INLINE static void sink_props_init(struct sink_props *sp, sp->maximal_density_threshold *= m_p_cgs / units_cgs_conversion_factor(us, UNIT_CONV_DENSITY); + const double Myr_internal_units = 1e6 * phys_const->const_year; + sp->max_time_step_young = max_time_step_young_Myr * Myr_internal_units; + sp->max_time_step_old = max_time_step_old_Myr * Myr_internal_units; + sp->age_threshold = age_threshold_Myr * Myr_internal_units; + sp->age_threshold_unlimited = + age_threshold_unlimited_Myr * Myr_internal_units; + /* here, we need to differenciate between the stellar models */ struct initial_mass_function *imf; struct stellar_model *sm; @@ -290,6 +360,8 @@ INLINE static void sink_props_init(struct sink_props *sp, sp->density_threshold); message("maximal_density_threshold = %g", sp->maximal_density_threshold); + message("sink_minimal_mass (in M_sun) = %g", + sp->sink_minimal_mass_Msun); message("stellar_particle_mass (in M_sun) = %g", sp->stellar_particle_mass_Msun); message("minimal_discrete_mass (in M_sun) = %g", @@ -313,6 +385,21 @@ INLINE static void sink_props_init(struct sink_props *sp, sp->sink_formation_bound_state_criterion); message("sink_formation_overlapping_sink_criterion = %d", sp->sink_formation_overlapping_sink_criterion); + + /* Print timestep parameters information */ + message("sink max_timestep_young = %e", + sp->max_time_step_young); + message("sink max_timestep_old = %e", + sp->max_time_step_old); + message("sink age_threshold from young to old = %e", + sp->age_threshold); + message("sink age_threshold from old to unlimited = %e", + sp->age_threshold_unlimited); + message("sink C_CFL = %e", + sp->CFL_condition); + message("tolerance_SF_timestep = %e", + sp->tolerance_SF_timestep); + message("n_IMF = %e", sp->n_IMF); } /** diff --git a/src/sink/GEAR/sink_setters.h b/src/sink/GEAR/sink_setters.h new file mode 100644 index 0000000000000000000000000000000000000000..f28cea1c1c093b8558a9ce0dabac8edb44071e55 --- /dev/null +++ b/src/sink/GEAR/sink_setters.h @@ -0,0 +1,117 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Darwin Roduit (darwin.roduit@alumni.epfl.ch) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + * + *******************************************************************************/ +#ifndef SWIFT_GEAR_SINK_SETTERS_H +#define SWIFT_GEAR_SINK_SETTERS_H + +#include "sink_part.h" + +/** + * @file src/sink/GEAR/sink_setters.h + * @brief Setters functions for GEAR sink scheme to avoid exposing + * implementation details to the outer world. Keep the code clean and lean. + */ + +/** + * @brief Set the birth time/scale-factor of a #sink particle. + * + * @param sink The #sink. + * @param birth_time Birth time of the #sink. + * @param birth_scale_factor Birth scale-factor of the star. + * @param with_cosmology If we run with cosmology. + */ + +__attribute__((always_inline)) INLINE void +sink_set_sink_birth_time_or_scale_factor(struct sink* restrict sink, + const float birth_time, + const float birth_scale_factor, + const int with_cosmology) { + if (with_cosmology) { + sink->birth_scale_factor = birth_scale_factor; + } else { + sink->birth_time = birth_time; + } +} + +/** + * @brief Update the target mass of the sink particle. + * + * @param sink the #sink particle. + * @param sink_props the sink properties to use. + * @param phys_const the physical constants in internal units. + * @param e The #engine + * @param star_counter The star loop counter. + */ +INLINE static void sink_update_target_mass(struct sink* sink, + const struct sink_props* sink_props, + const struct engine* e, + int star_counter) { + + float random_number = random_unit_interval_part_ID_and_index( + sink->id, star_counter, e->ti_current, random_number_sink_formation); + + const struct feedback_props* feedback_props = e->feedback_props; + + /* Pick the correct table. (if only one table, threshold is < 0) */ + const float metal = + chemistry_get_sink_total_iron_mass_fraction_for_feedback(sink); + const float threshold = feedback_props->metallicity_max_first_stars; + + /* If metal < threshold, then the sink generates first star particles. */ + const int is_first_star = metal < threshold; + const struct stellar_model* model; + double minimal_discrete_mass_Msun; + + /* Take the correct values if your are a first star or not. */ + if (!is_first_star) /* (metal >= threshold)*/ { + model = &feedback_props->stellar_model; + minimal_discrete_mass_Msun = sink_props->minimal_discrete_mass_Msun; + } else { + model = &feedback_props->stellar_model_first_stars; + minimal_discrete_mass_Msun = + sink_props->minimal_discrete_mass_first_stars_Msun; + } + + const struct initial_mass_function* imf = &model->imf; + + if (random_number < imf->sink_Pc) { + /* We are dealing with the continous part of the IMF. */ + sink->target_mass_Msun = imf->stellar_particle_mass_Msun; + sink->target_type = star_population_continuous_IMF; + } else { + /* We are dealing with the discrete part of the IMF. */ + random_number = random_unit_interval_part_ID_and_index( + sink->id, star_counter + 1, e->ti_current, + random_number_sink_formation); + double m = initial_mass_function_sample_power_law( + minimal_discrete_mass_Msun, imf->mass_max, imf->exp[imf->n_parts - 1], + random_number); + sink->target_mass_Msun = m; + sink->target_type = single_star; + } + + /* Also store the mass of the IMF into the sink. */ + double dummy; + ; + initial_mass_function_compute_Mc_Md_Mtot(imf, &dummy, &dummy, + &sink->mass_IMF); + + /* Convert from M_sun to internal units. */ + sink->mass_IMF *= e->physical_constants->const_solar_mass; +} +#endif /* SWIFT_GEAR_SINK_SETTERS_H */ diff --git a/src/space.c b/src/space.c index 72125200b25aa0e01edd66aeed508c89e325caf9..95d029c73da5828759b39a971fb1dee12cbef0c2 100644 --- a/src/space.c +++ b/src/space.c @@ -150,6 +150,11 @@ void space_free_foreign_parts(struct space *s, const int clear_cell_pointers) { s->size_bparts_foreign = 0; s->bparts_foreign = NULL; } + if (s->sinks_foreign != NULL) { + swift_free("sinks_foreign", s->sinks_foreign); + s->size_sinks_foreign = 0; + s->sinks_foreign = NULL; + } if (clear_cell_pointers) { for (int k = 0; k < s->e->nr_proxies; k++) { for (int j = 0; j < s->e->proxies[k].nr_cells_in; j++) { @@ -2175,6 +2180,7 @@ void space_check_drift_point(struct space *s, integertime_t ti_drift, space_map_cells_pre(s, 1, cell_check_part_drift_point, &ti_drift); space_map_cells_pre(s, 1, cell_check_gpart_drift_point, &ti_drift); space_map_cells_pre(s, 1, cell_check_spart_drift_point, &ti_drift); + space_map_cells_pre(s, 1, cell_check_sink_drift_point, &ti_drift); if (multipole) space_map_cells_pre(s, 1, cell_check_multipole_drift_point, &ti_drift); #else @@ -2321,6 +2327,55 @@ void space_check_bpart_swallow_mapper(void *map_data, int nr_bparts, #endif } +/** + * @brief #threadpool mapper function for the swallow debugging check + */ +void space_check_part_sink_swallow_mapper(void *map_data, int nr_parts, + void *extra_data) { +#ifdef SWIFT_DEBUG_CHECKS + /* Unpack the data */ + struct part *restrict parts = (struct part *)map_data; + + /* Verify that all particles have been swallowed or are untouched */ + for (int k = 0; k < nr_parts; k++) { + + if (parts[k].time_bin == time_bin_inhibited) continue; + + const long long swallow_id = sink_get_part_swallow_id(&parts[k].sink_data); + + if (swallow_id != -1) + error("Particle has not been swallowed! id=%lld", parts[k].id); + } +#else + error("Calling debugging code without debugging flag activated."); +#endif +} + +/** + * @brief #threadpool mapper function for the swallow debugging check + */ +void space_check_sink_sink_swallow_mapper(void *map_data, int nr_sinks, + void *extra_data) { +#ifdef SWIFT_DEBUG_CHECKS + /* Unpack the data */ + struct sink *restrict sinks = (struct sink *)map_data; + + /* Verify that all particles have been swallowed or are untouched */ + for (int k = 0; k < nr_sinks; k++) { + + if (sinks[k].time_bin == time_bin_inhibited) continue; + + const long long swallow_id = + sink_get_sink_swallow_id(&sinks[k].merger_data); + + if (swallow_id != -1) + error("Sink particle has not been swallowed! id=%lld", sinks[k].id); + } +#else + error("Calling debugging code without debugging flag activated."); +#endif +} + /** * @brief Checks that all particles have their swallow flag in a "no swallow" * state. @@ -2339,6 +2394,16 @@ void space_check_swallow(struct space *s) { threadpool_map(&s->e->threadpool, space_check_bpart_swallow_mapper, s->bparts, s->nr_bparts, sizeof(struct bpart), threadpool_auto_chunk_size, /*extra_data=*/NULL); + + threadpool_map(&s->e->threadpool, space_check_part_sink_swallow_mapper, + s->parts, s->nr_parts, sizeof(struct part), + threadpool_auto_chunk_size, + /*extra_data=*/NULL); + + threadpool_map(&s->e->threadpool, space_check_sink_sink_swallow_mapper, + s->sinks, s->nr_sinks, sizeof(struct sink), + threadpool_auto_chunk_size, + /*extra_data=*/NULL); #else error("Calling debugging code without debugging flag activated."); #endif @@ -2437,6 +2502,7 @@ void space_clean(struct space *s) { swift_free("sparts_foreign", s->sparts_foreign); swift_free("gparts_foreign", s->gparts_foreign); swift_free("bparts_foreign", s->bparts_foreign); + swift_free("sinks_foreign", s->sinks_foreign); #endif free(s->cells_sub); free(s->multipoles_sub); @@ -2611,6 +2677,8 @@ void space_struct_restore(struct space *s, FILE *stream) { s->size_sparts_foreign = 0; s->bparts_foreign = NULL; s->size_bparts_foreign = 0; + s->sinks_foreign = NULL; + s->size_sinks_foreign = 0; #endif /* More things to read. */ @@ -2727,9 +2795,10 @@ void space_write_cell(const struct space *s, FILE *f, const struct cell *c) { /* Write line for current cell */ fprintf(f, "%lld,%lld,%i,", c->cellID, parent, c->nodeID); - fprintf(f, "%i,%i,%i,%s,%s,%g,%g,%g,%g,%g,%g, ", c->hydro.count, - c->stars.count, c->grav.count, superID, hydro_superID, c->loc[0], - c->loc[1], c->loc[2], c->width[0], c->width[1], c->width[2]); + fprintf(f, "%i,%i,%i,%i,%s,%s,%g,%g,%g,%g,%g,%g, ", c->hydro.count, + c->stars.count, c->grav.count, c->sinks.count, superID, hydro_superID, + c->loc[0], c->loc[1], c->loc[2], c->width[0], c->width[1], + c->width[2]); fprintf(f, "%g, %g, %i, %i\n", c->hydro.h_max, c->stars.h_max, c->depth, c->maxdepth); @@ -2761,7 +2830,7 @@ void space_write_cell_hierarchy(const struct space *s, int j) { if (engine_rank == 0) { fprintf(f, "name,parent,mpi_rank,"); fprintf(f, - "hydro_count,stars_count,gpart_count,super,hydro_super," + "hydro_count,stars_count,gpart_count,sinks_count,super,hydro_super," "loc1,loc2,loc3,width1,width2,width3,"); fprintf(f, "hydro_h_max,stars_h_max,depth,maxdepth\n"); diff --git a/src/space.h b/src/space.h index df58eeda4b4bca6d9909786247877cf1caafdec2..8a85277f82d6d9169f6f984d02671f870b9a2961 100644 --- a/src/space.h +++ b/src/space.h @@ -349,6 +349,10 @@ struct space { struct bpart *bparts_foreign; size_t nr_bparts_foreign, size_bparts_foreign; + /*! Buffers for sink-parts that we will receive from foreign cells. */ + struct sink *sinks_foreign; + size_t nr_sinks_foreign, size_sinks_foreign; + #endif }; diff --git a/src/space_cell_index.c b/src/space_cell_index.c index c9d7cdd07fdcce9a2333d40e9d94fa10e064a75d..c0ad4f2df00160cc341d585cf6295c2728e4795a 100644 --- a/src/space_cell_index.c +++ b/src/space_cell_index.c @@ -692,8 +692,8 @@ void space_sinks_get_cell_index_mapper(void *map_data, int nr_sinks, if (count_extra_sink) atomic_add(&data->count_extra_sink, count_extra_sink); /* Write back the minimal part mass and velocity sum */ - atomic_min_f(&s->min_spart_mass, min_mass); - atomic_add_f(&s->sum_spart_vel_norm, sum_vel_norm); + atomic_min_f(&s->min_sink_mass, min_mass); + atomic_add_f(&s->sum_sink_vel_norm, sum_vel_norm); } /** diff --git a/src/space_extras.c b/src/space_extras.c index 9d8a88952e5a3d2a670ba63f96001dc27eda26c7..64ea69647ffc1f88ce600703a931a53626dea3d5 100644 --- a/src/space_extras.c +++ b/src/space_extras.c @@ -359,6 +359,7 @@ void space_allocate_extras(struct space *s, int verbose) { bzero(&s->sinks[i], sizeof(struct sink)); s->sinks[i].time_bin = time_bin_not_created; s->sinks[i].id = -42; + sink_mark_sink_as_not_swallowed(&s->sinks[i].merger_data); } /* Put the spare particles in their correct cell */ diff --git a/src/space_first_init.c b/src/space_first_init.c index fd36583715849c662ac56ae47d28ceff4396074a..3232e445f0fb21fb1dfeae14d91786efa2ed8688 100644 --- a/src/space_first_init.c +++ b/src/space_first_init.c @@ -484,6 +484,9 @@ void space_first_init_sinks_mapper(void *restrict map_data, int count, /* Initialize the sinks */ sink_first_init_sink(&sink[k], props, e); + /* And the sink merger markers */ + sink_mark_sink_as_not_swallowed(&sink[k].merger_data); + /* Also initialize the chemistry */ chemistry_first_init_sink(chemistry, &sink[k]); diff --git a/src/space_rebuild.c b/src/space_rebuild.c index 36d784c8b34c393d4260fbd54ae5fa02b5c4c771..7d458f6858ea64ea9970a376d177c2211cc5ea0a 100644 --- a/src/space_rebuild.c +++ b/src/space_rebuild.c @@ -491,17 +491,20 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { size_t nr_gparts_exchanged = s->nr_gparts - nr_gparts; size_t nr_sparts_exchanged = s->nr_sparts - nr_sparts; size_t nr_bparts_exchanged = s->nr_bparts - nr_bparts; + size_t nr_sinks_exchanged = s->nr_sinks - nr_sinks; engine_exchange_strays(s->e, nr_parts, &h_index[nr_parts], &nr_parts_exchanged, nr_gparts, &g_index[nr_gparts], &nr_gparts_exchanged, nr_sparts, &s_index[nr_sparts], &nr_sparts_exchanged, nr_bparts, &b_index[nr_bparts], - &nr_bparts_exchanged); + &nr_bparts_exchanged, nr_sinks, + &sink_index[nr_sinks], &nr_sinks_exchanged); /* Set the new particle counts. */ s->nr_parts = nr_parts + nr_parts_exchanged; s->nr_gparts = nr_gparts + nr_gparts_exchanged; s->nr_sparts = nr_sparts + nr_sparts_exchanged; s->nr_bparts = nr_bparts + nr_bparts_exchanged; + s->nr_sinks = nr_sinks + nr_sinks_exchanged; } else { #ifdef SWIFT_DEBUG_CHECKS @@ -511,6 +514,8 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { error("Number of sparts changing after repartition"); if (s->nr_gparts != nr_gparts) error("Number of gparts changing after repartition"); + if (s->nr_sinks != nr_sinks) + error("Number of sinks changing after repartition"); #endif } @@ -521,6 +526,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { cell_spart_counts[k] = 0; cell_gpart_counts[k] = 0; cell_bpart_counts[k] = 0; + cell_sink_counts[k] = 0; } } @@ -547,16 +553,27 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { } /* Re-allocate the index array for the bparts if needed.. */ - if (s->nr_bparts + 1 > s_index_size) { + if (s->nr_bparts + 1 > b_index_size) { int *bind_new; if ((bind_new = (int *)swift_malloc( "b_index", sizeof(int) * (s->nr_bparts + 1))) == NULL) - error("Failed to allocate temporary s-particle indices."); + error("Failed to allocate temporary b-particle indices."); memcpy(bind_new, b_index, sizeof(int) * nr_bparts); swift_free("b_index", b_index); b_index = bind_new; } + /* Re-allocate the index array for the sinks if needed.. */ + if (s->nr_sinks + 1 > sink_index_size) { + int *sink_ind_new; + if ((sink_ind_new = (int *)swift_malloc( + "sink_index", sizeof(int) * (s->nr_sinks + 1))) == NULL) + error("Failed to allocate temporary sink-particle indices."); + memcpy(sink_ind_new, sink_index, sizeof(int) * nr_sinks); + swift_free("sink_index", sink_index); + sink_index = sink_ind_new; + } + const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; const double ih[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]}; @@ -602,6 +619,20 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { } nr_bparts = s->nr_bparts; + /* Assign each received sink to its cell. */ + for (size_t k = nr_sinks; k < s->nr_sinks; k++) { + const struct sink *const sink = &s->sinks[k]; + sink_index[k] = cell_getid(cdim, sink->x[0] * ih[0], sink->x[1] * ih[1], + sink->x[2] * ih[2]); + cell_sink_counts[sink_index[k]]++; +#ifdef SWIFT_DEBUG_CHECKS + if (cells_top[sink_index[k]].nodeID != local_nodeID) + error("Received sink-part that does not belong to me (nodeID=%i).", + cells_top[sink_index[k]].nodeID); +#endif + } + nr_sinks = s->nr_sinks; + #else /* WITH_MPI */ /* Update the part, spart and bpart counters */ @@ -716,14 +747,14 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { error("Inhibited particle sorted into a cell!"); /* New cell index */ - const int new_bind = + const int new_sink_ind = cell_getid(s->cdim, sink->x[0] * s->iwidth[0], sink->x[1] * s->iwidth[1], sink->x[2] * s->iwidth[2]); /* New cell of this sink */ - const struct cell *c = &s->cells_top[new_bind]; + const struct cell *c = &s->cells_top[new_sink_ind]; - if (sink_index[k] != new_bind) + if (sink_index[k] != new_sink_ind) error("sink's new cell index not matching sorted index."); if (sink->x[0] < c->loc[0] || sink->x[0] > c->loc[0] + c->width[0] || @@ -931,6 +962,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { /* Store the state at rebuild time */ c->stars.parts_rebuild = c->stars.parts; + c->sinks.parts_rebuild = c->sinks.parts; c->grav.parts_rebuild = c->grav.parts; c->hydro.count_total = c->hydro.count + space_extra_parts; diff --git a/src/space_recycle.c b/src/space_recycle.c index 1a7a42c830e3636815feea2b4ad5e1d70bdccb45..b510eb9dced0f3ee4e9eee47a0d8c434962e0c03 100644 --- a/src/space_recycle.c +++ b/src/space_recycle.c @@ -193,6 +193,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements, c->grav.parts = NULL; c->grav.parts_rebuild = NULL; c->sinks.parts = NULL; + c->sinks.parts_rebuild = NULL; c->stars.parts = NULL; c->stars.parts_rebuild = NULL; c->black_holes.parts = NULL; diff --git a/src/space_regrid.c b/src/space_regrid.c index 487fe7c0e38495ac85622156615308cd4179ec01..dc86468e335b3cc2e50d6250327f532bc64e22eb 100644 --- a/src/space_regrid.c +++ b/src/space_regrid.c @@ -64,7 +64,8 @@ void space_regrid(struct space *s, int verbose) { if (c->black_holes.h_max > h_max) { h_max = c->black_holes.h_max; } - if (c->sinks.r_cut_max > h_max) { + /* Note: For sinks, r_cut <=> h*kernel_gamma */ + if (c->sinks.r_cut_max > h_max * kernel_gamma) { h_max = c->sinks.r_cut_max / kernel_gamma; } } @@ -82,7 +83,9 @@ void space_regrid(struct space *s, int verbose) { if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) { h_max = c->black_holes.h_max; } - if (c->nodeID == engine_rank && c->sinks.r_cut_max > h_max) { + /* Note: For sinks, r_cut <=> h*kernel_gamma */ + if (c->nodeID == engine_rank && + c->sinks.r_cut_max > h_max * kernel_gamma) { h_max = c->sinks.r_cut_max / kernel_gamma; } } diff --git a/src/space_unique_id.c b/src/space_unique_id.c index bfcb733cc6daa698e2744251d7f76b5dc93068e0..481d8a362799e310422e6020361d8ea7fc1ab913 100644 --- a/src/space_unique_id.c +++ b/src/space_unique_id.c @@ -82,9 +82,10 @@ void space_update_unique_id(struct space *s) { #endif // WITH_MPI /* Compute the size of each batch. */ - const long long batch_size = (space_extra_parts + space_extra_sparts + - space_extra_gparts + space_extra_bparts) * - s->nr_cells; + const long long batch_size = + (space_extra_parts + space_extra_sparts + space_extra_gparts + + space_extra_bparts + space_extra_sinks) * + s->nr_cells; /* Get a new batch. */ if (require_new_batch) { @@ -210,9 +211,10 @@ void space_init_unique_id(struct space *s, int nr_nodes) { s->unique_id.global_next_id++; /* Compute the size of each batch. */ - const long long batch_size = (space_extra_parts + space_extra_sparts + - space_extra_gparts + space_extra_bparts) * - s->nr_cells; + const long long batch_size = + (space_extra_parts + space_extra_sparts + space_extra_gparts + + space_extra_bparts + space_extra_sinks) * + s->nr_cells; /* Compute the initial batchs (each rank has 2 batchs). */ if (s->unique_id.global_next_id > LLONG_MAX - 2 * engine_rank * batch_size) { diff --git a/src/stars/GEAR/stars_io.h b/src/stars/GEAR/stars_io.h index 71e1aad8e196fe21e3b15b801cd217bbc971aed6..e24934b6ae8695026889490675e1c9f1d3e33274 100644 --- a/src/stars/GEAR/stars_io.h +++ b/src/stars/GEAR/stars_io.h @@ -245,10 +245,6 @@ INLINE static void stars_props_init(struct stars_props *sp, sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv)); /* Maximal time-step lengths */ - const double Myr_internal_units = - 1e6 * phys_const->const_year * - units_cgs_conversion_factor(us, UNIT_CONV_TIME); - const double max_time_step_young_Myr = parser_get_opt_param_float( params, "Stars:max_timestep_young_Myr", FLT_MAX); const double max_time_step_old_Myr = @@ -269,6 +265,7 @@ INLINE static void stars_props_init(struct stars_props *sp, } /* Convert to internal units */ + const double Myr_internal_units = 1e6 * phys_const->const_year; sp->max_time_step_young = max_time_step_young_Myr * Myr_internal_units; sp->max_time_step_old = max_time_step_old_Myr * Myr_internal_units; sp->age_threshold = age_threshold_Myr * Myr_internal_units; diff --git a/src/timestep.h b/src/timestep.h index 5e5af8074632fcfaa486e26c23f2fc1e170c3748..a6002418088e1c194f31bfd2269456f9752f8497 100644 --- a/src/timestep.h +++ b/src/timestep.h @@ -397,7 +397,9 @@ __attribute__((always_inline)) INLINE static integertime_t get_sink_timestep( const struct sink *restrict sink, const struct engine *restrict e) { /* Sink time-step */ - float new_dt_sink = sink_compute_timestep(sink); + float new_dt_sink = sink_compute_timestep( + sink, e->sink_properties, (e->policy & engine_policy_cosmology), + e->cosmology, e->gravity_properties, e->time, e->time_base); /* Gravity time-step */ float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX; diff --git a/swift.c b/swift.c index b63941cd63538e3226669bfed06bdb1837c73411..41c96e0afbc4ca5ad2a119e2089397a5fcbd540e 100644 --- a/swift.c +++ b/swift.c @@ -528,11 +528,17 @@ int main(int argc, char *argv[]) { #endif #ifdef WITH_MPI +#ifdef SWIFT_DEBUG_CHECKS + if (with_sinks) { + pretime_message("Warning: sink particles are are WIP yet with MPI."); + } +#else if (with_sinks) { pretime_message("Error: sink particles are not available yet with MPI."); return 1; } -#endif +#endif /* SWIFT_DEBUG_CHECKS */ +#endif /* WITH_MPI */ if (with_sinks && with_star_formation) { pretime_message( diff --git a/tests/difffloat.py b/tests/difffloat.py index cfa16970730bc365b284a78f8a757bd459a9b741..6f60a66535b39e9aa4d06951cc088f975435d46a 100644 --- a/tests/difffloat.py +++ b/tests/difffloat.py @@ -18,6 +18,7 @@ ############################################################################## from numpy import * + del min import sys diff --git a/theory/SinkParticles/plot_imf_image.py b/theory/SinkParticles/plot_imf_image.py index 25da38cab2b804a651c5f745063e6ba101b38a52..bc1da9cccf1c3fe60af2661e4cf65a96fb949338 100644 --- a/theory/SinkParticles/plot_imf_image.py +++ b/theory/SinkParticles/plot_imf_image.py @@ -14,18 +14,17 @@ import matplotlib.pyplot as plt import matplotlib # For swift doc, use the following and not the styleheet -plt.rcParams['text.usetex'] = True -plt.rcParams['font.family'] = 'serif' -plt.rcParams['font.serif'] = ['Computer Modern Roman'] +plt.rcParams["text.usetex"] = True +plt.rcParams["font.family"] = "serif" +plt.rcParams["font.serif"] = ["Computer Modern Roman"] mmin = 0.5 mmax = 300 -matplotlib.rcParams.update({'font.size': 16}) +matplotlib.rcParams.update({"font.size": 16}) figsize = (6.4, 4.8) -fig, ax = plt.subplots(num=1, ncols=1, nrows=1, - figsize=figsize, layout="tight") +fig, ax = plt.subplots(num=1, ncols=1, nrows=1, figsize=figsize, layout="tight") ax.set_xlim([0.3, 400]) ax.set_ylim([1, 5e4]) @@ -34,34 +33,34 @@ ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter()) ax.set_xlabel("$M_{\star}$ $[M_\odot]$") ax.set_ylabel("d$N/$d$M$ [arbitrary units]") -ax.set_xscale('log') -ax.set_yscale('log') +ax.set_xscale("log") +ax.set_yscale("log") # theoretical imf s = -1.3 -bins = 10**np.linspace(np.log10(mmin), np.log10(mmax), 100) -n = 0.9*10000*bins**s +bins = 10 ** np.linspace(np.log10(mmin), np.log10(mmax), 100) +n = 0.9 * 10000 * bins ** s ax.plot(bins, n, "k--") -bins = 10**np.linspace(np.log10(mmin), np.log10(8), 100) -n = 0.9*10000*bins**s +bins = 10 ** np.linspace(np.log10(mmin), np.log10(8), 100) +n = 0.9 * 10000 * bins ** s ax.fill_between(bins, 0.1, n, color="red", alpha=0.1) -bins = 10**np.linspace(np.log10(8), np.log10(mmax), 100) -n = 0.9*10000*bins**s +bins = 10 ** np.linspace(np.log10(8), np.log10(mmax), 100) +n = 0.9 * 10000 * bins ** s ax.fill_between(bins, 0.1, n, color="blue", alpha=0.1) -ax.text(2, 1e2, r"$M_{\rm c}$", horizontalalignment='center') -ax.text(50, 2, r"$M_{\rm d}$", horizontalalignment='center') +ax.text(2, 1e2, r"$M_{\rm c}$", horizontalalignment="center") +ax.text(50, 2, r"$M_{\rm d}$", horizontalalignment="center") # Add limit -ax.vlines(x=8, ymin=0, ymax=600, color='k', linestyle='-') +ax.vlines(x=8, ymin=0, ymax=600, color="k", linestyle="-") # Add text to the vertical line -ax.text(8, 800, r"$m_{t}$", horizontalalignment='center') +ax.text(8, 800, r"$m_{t}$", horizontalalignment="center") # fig.patch.set_facecolor('none') # Remove figure background # ax.set_facecolor('none') # Remove axes background -plt.savefig('sink_imf.png', dpi=300, bbox_inches='tight') +plt.savefig("sink_imf.png", dpi=300, bbox_inches="tight") plt.close() diff --git a/tools/data/cell_hierarchy.html b/tools/data/cell_hierarchy.html index 6c9bcdb995be7970948cf9b9544ae0a59d41caf0..051451091d70d7e6bcb000b44b71f73c2779caa6 100644 --- a/tools/data/cell_hierarchy.html +++ b/tools/data/cell_hierarchy.html @@ -162,6 +162,7 @@ "MPI rank:" + n.mpi_rank + "<br/>" + "Part: " + n.hydro_count + "<br/>" + "Spart: " + n.stars_count + "<br/>" + + "Sink: " + n.sinks_count + "<br/>" + "Super: " + n.super + "<br/>" + "Super Hydro: " + n.hydro_super + "<br/>" + "Loc: " + n.loc1 + ", " + n.loc2 + ", " + n.loc3 + "<br/>" + diff --git a/tools/plot_task_dependencies.py b/tools/plot_task_dependencies.py index 17013ca2f131ddf648d81f4fa70158df65bbd610..fa9fdb7fd1e5a04eb340f9fee0801918a6f3370c 100755 --- a/tools/plot_task_dependencies.py +++ b/tools/plot_task_dependencies.py @@ -208,7 +208,12 @@ def task_is_hydro(name): return True if "_part" in name: return True - if "density" in name and "stars" not in name and "sink" not in name and "bh" not in name: + if ( + "density" in name + and "stars" not in name + and "sink" not in name + and "bh" not in name + ): return True if "rho" in name and "bpart" not in name: return True