GEAR sink: add a minimal mass to sink

doc/RTD/source/SubgridModels/GEAR/sinks/params.rst

@@ -44,7 +44,9 @@ 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 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 +58,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

doc/RTD/source/SubgridModels/GEAR/sinks/theory.rst

@@ -197,7 +197,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 +205,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``.
 

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

src/sink/GEAR/sink.h

@@ -699,8 +699,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;

src/sink/GEAR/sink_properties.h

@@ -67,13 +67,17 @@ 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;
 };
 
 /**
@@ -230,6 +234,9 @@ INLINE static void sink_props_init(struct sink_props *sp,
       parser_get_opt_param_float(params, "GEARSink:star_spawning_sigma_factor",
                                  default_star_spawning_sigma_factor);
 
+  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",
@@ -290,6 +297,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",