Skip to content
Snippets Groups Projects
Commit 66f9ca53 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into multi_softening

parents 0d1c2ebb 8c86c97b
No related branches found
No related tags found
1 merge request!884Support for multiple softening lengths in the gravity solver
Showing
with 308 additions and 85 deletions
......@@ -1327,6 +1327,7 @@ with_subgrid_cooling=none
with_subgrid_chemistry=none
with_subgrid_tracers=none
with_subgrid_entropy_floor=none
with_subgrid_pressure_floor=none
with_subgrid_stars=none
with_subgrid_star_formation=none
with_subgrid_feedback=none
......@@ -1338,10 +1339,9 @@ case "$with_subgrid" in
none)
;;
GEAR)
with_subgrid_cooling=grackle
with_subgrid_cooling=grackle_0
with_subgrid_chemistry=GEAR
with_subgrid_tracers=none
with_subgrid_entropy_floor=none
with_subgrid_pressure_floor=GEAR
with_subgrid_stars=GEAR
with_subgrid_star_formation=GEAR
with_subgrid_feedback=none
......@@ -1400,6 +1400,12 @@ AC_ARG_WITH([gravity],
[with_gravity="default"]
)
if test "$with_subgrid" = "EAGLE"; then
if test "$with_gravity" = "default"; then
with_gravity="with-potential"
fi
fi
case "$with_gravity" in
with-potential)
AC_DEFINE([POTENTIAL_GRAVITY], [1], [Gravity scheme with potential calculation])
......@@ -1643,7 +1649,8 @@ esac
# Cooling function
AC_ARG_WITH([cooling],
[AS_HELP_STRING([--with-cooling=<function>],
[cooling function @<:@none, const-du, const-lambda, EAGLE, grackle, grackle1, grackle2, grackle3 default: none@:>@]
[cooling function @<:@none, const-du, const-lambda, EAGLE, grackle_* default: none@:>@.
For Grackle, you need to provide the primordial chemistry parameter (e.g. grackle_0)]
)],
[with_cooling="$withval"],
[with_cooling="none"]
......@@ -1670,21 +1677,10 @@ case "$with_cooling" in
compton)
AC_DEFINE([COOLING_COMPTON], [1], [Compton cooling off the CMB])
;;
grackle)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [0], [Grackle chemistry network, mode 0])
;;
grackle1)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [1], [Grackle chemistry network, mode 1])
;;
grackle2)
grackle_*)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [2], [Grackle chemistry network, mode 2])
;;
grackle3)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [3], [Grackle chemistry network, mode 3])
primordial_chemistry=${with_cooling:8}
AC_DEFINE_UNQUOTED([COOLING_GRACKLE_MODE], [$primordial_chemistry], [Grackle chemistry network])
;;
EAGLE)
AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model])
......@@ -1919,6 +1915,35 @@ case "$with_entropy_floor" in
;;
esac
# Pressure floor
AC_ARG_WITH([pressure-floor],
[AS_HELP_STRING([--with-pressure-floor=<floor>],
[pressure floor @<:@none, GEAR, default: none@:>@
The hydro model needs to be compatible.]
)],
[with_pressure_floor="$withval"],
[with_pressure_floor="none"]
)
if test "$with_subgrid" != "none"; then
if test "$with_pressure_floor" != "none"; then
AC_MSG_ERROR([Cannot provide with-subgrid and with-pressure-floor together])
else
with_pressure_floor="$with_subgrid_pressure_floor"
fi
fi
case "$with_pressure_floor" in
none)
AC_DEFINE([PRESSURE_FLOOR_NONE], [1], [No pressure floor])
;;
GEAR)
AC_DEFINE([PRESSURE_FLOOR_GEAR], [1], [GEAR pressure floor])
;;
*)
AC_MSG_ERROR([Unknown pressure floor model])
;;
esac
# Star formation
AC_ARG_WITH([star-formation],
[AS_HELP_STRING([--with-star-formation=<sfm>],
......@@ -1964,6 +1989,13 @@ AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Mul
AC_PATH_PROG([GIT_CMD], [git])
AC_SUBST([GIT_CMD])
# Check that the gravity model is compatible with the subgrid
if test $with_black_holes = "EAGLE"; then
if test $with_gravity != "with-potential"; then
AC_MSG_ERROR([The EAGLE BH model needs the gravity scheme to provide potentials. The code must be compile with --with-gravity=with-potential.])
fi
fi
# Make the documentation. Add conditional to handle disable option.
DX_INIT_DOXYGEN(libswift,doc/Doxyfile,doc/)
AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
......@@ -2046,6 +2078,7 @@ AC_MSG_RESULT([
External potential : $with_potential
Entropy floor : $with_entropy_floor
Pressure floor : $with_pressure_floor
Cooling function : $with_cooling
Chemistry : $with_chemistry
Tracers : $with_tracers
......
......@@ -577,6 +577,11 @@ Black-hole creation
Black-hole accretion
~~~~~~~~~~~~~~~~~~~~
.. _EAGLE_black_hole_reposition:
Black-hole repositioning
~~~~~~~~~~~~~~~~~~~~~~~~
.. _EAGLE_black_hole_feedback:
AGN feedback
......
......@@ -5,6 +5,29 @@
GEAR model
===========
Pressure Floor
~~~~~~~~~~~~~~
In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles.
This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by:
.. math::
P_\textrm{Jeans} = \frac{\rho}{\gamma} \left( \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3} - \sigma^2 \right)
where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant,
:math:`h` the smoothing length, :math:`N_\textrm{Jeans}` is the number of particle required in order to resolve a clump and
:math:`\sigma` the velocity dispersion.
This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2 and Pressure-Energy) are currently implemented.
In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
.. math::
m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right]
and simply replace the :math:`P_i, P_j` by the pressure with the floor.
Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy).
Cooling: Grackle
~~~~~~~~~~~~~~~~
......
......@@ -170,3 +170,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
......@@ -171,3 +171,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
......@@ -171,5 +171,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
......@@ -162,3 +162,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
......@@ -169,3 +169,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
......@@ -161,3 +161,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
......@@ -171,3 +171,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
......@@ -99,3 +99,6 @@ GrackleCooling:
GearChemistry:
InitialMetallicity: 0.01295
GEARPressureFloor:
Jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
......@@ -771,6 +771,13 @@ int main(int argc, char *argv[]) {
else
bzero(&entropy_floor, sizeof(struct entropy_floor_properties));
/* Initialise the pressure floor */
if (with_hydro)
pressure_floor_init(&pressure_floor_props, &prog_const, &us,
&hydro_properties, params);
else
bzero(&pressure_floor_props, sizeof(struct pressure_floor_properties));
/* Initialise the stars properties */
if (with_stars)
stars_props_init(&stars_properties, &prog_const, &us, params,
......
......@@ -288,6 +288,11 @@ EAGLEEntropyFloor:
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
# Parameters related to pressure floors ----------------------------------------------
GEARPressureFloor:
Jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
# Parameters related to cooling function ----------------------------------------------
# Constant du/dt cooling function
......@@ -341,6 +346,11 @@ EAGLEChemistry:
# Parameters related to star formation models -----------------------------------------------
# GEAR star formation model (Revaz and Jablonka 2018)
GEARStarFormation:
star_formation_efficiency: 0.01 # star formation efficiency (c_*)
maximal_temperature: 3e4 # Upper limit to the temperature of a star forming particle
# EAGLE star formation model (Schaye and Dalla Vecchia 2008)
EAGLEStarFormation:
EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
......@@ -406,4 +416,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
......@@ -79,7 +79,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c
collectgroup.c hydro_space.c equation_of_state.c \
chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
outputlist.c velociraptor_dummy.c logger_io.c memuse.c fof.c \
hashmap.c \
hashmap.c pressure_floor.c \
$(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES)
# Include files for distribution, not installation.
......@@ -98,28 +98,28 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro.h hydro_io.h hydro_parameters.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
hydro/Minimal/hydro_parameters.h \
hydro/Minimal/hydro_parameters.h \
hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \
hydro/Default/hydro_debug.h hydro/Default/hydro_part.h \
hydro/Default/hydro_parameters.h \
hydro/Default/hydro_parameters.h \
hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h \
hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
hydro/Gadget2/hydro_parameters.h \
hydro/Gadget2/hydro_parameters.h \
hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \
hydro/PressureEntropy/hydro_parameters.h \
hydro/PressureEntropy/hydro_parameters.h \
hydro/PressureEnergy/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
hydro/PressureEnergy/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
hydro/PressureEnergy/hydro_parameters.h \
hydro/PressureEnergy/hydro_parameters.h \
hydro/PressureEnergyMorrisMonaghanAV/hydro.h hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h \
hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h \
hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h \
hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h \
hydro/AnarchyPU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
hydro/AnarchyPU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
hydro/AnarchyPU/hydro_parameters.h \
hydro/AnarchyPU/hydro_parameters.h \
hydro/AnarchyDU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
hydro/AnarchyDU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
hydro/AnarchyDU/hydro_parameters.h \
hydro/AnarchyDU/hydro_parameters.h \
hydro/GizmoMFV/hydro.h hydro/GizmoMFV/hydro_iact.h \
hydro/GizmoMFV/hydro_io.h hydro/GizmoMFV/hydro_debug.h \
hydro/GizmoMFV/hydro_part.h \
......@@ -131,7 +131,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/GizmoMFV/hydro_slope_limiters.h \
hydro/GizmoMFV/hydro_unphysical.h \
hydro/GizmoMFV/hydro_velocities.h \
hydro/GizmoMFV/hydro_parameters.h \
hydro/GizmoMFV/hydro_parameters.h \
hydro/GizmoMFM/hydro.h hydro/GizmoMFM/hydro_iact.h \
hydro/GizmoMFM/hydro_io.h hydro/GizmoMFM/hydro_debug.h \
hydro/GizmoMFM/hydro_part.h \
......@@ -142,7 +142,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/GizmoMFM/hydro_slope_limiters_face.h \
hydro/GizmoMFM/hydro_slope_limiters.h \
hydro/GizmoMFM/hydro_unphysical.h \
hydro/GizmoMFM/hydro_parameters.h \
hydro/GizmoMFM/hydro_parameters.h \
hydro/Shadowswift/hydro_debug.h \
hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
hydro/Shadowswift/hydro_iact.h \
......@@ -190,7 +190,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
cooling/const_lambda/cooling_io.h \
cooling/grackle/cooling.h cooling/grackle/cooling_struct.h \
cooling/grackle/cooling_io.h \
cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h \
cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h cooling/EAGLE/cooling_tables.h \
cooling/EAGLE/cooling_io.h cooling/EAGLE/interpolate.h cooling/EAGLE/cooling_rates.h \
chemistry/none/chemistry.h \
chemistry/none/chemistry_io.h \
......@@ -221,7 +221,10 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
black_holes/Default/black_holes_struct.h \
black_holes/EAGLE/black_holes.h black_holes/EAGLE/black_holes_io.h \
black_holes/EAGLE/black_holes_part.h black_holes/EAGLE/black_holes_iact.h \
black_holes/EAGLE/black_holes_properties.h
black_holes/EAGLE/black_holes_properties.h \
black_holes/EAGLE/black_holes_struct.h \
pressure_floor.h \
pressure_floor/GEAR/pressure_floor.h pressure_floor/none/pressure_floor.h
# Sources and flags for regular library
......
......@@ -187,6 +187,20 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
const struct phys_const* constants, const struct cosmology* cosmo,
const double dt) {}
/**
* @brief Finish the calculation of the new BH position.
*
* Nothing to do here.
*
* @param bp The black hole particle.
* @param props The properties of the black hole scheme.
* @param constants The physical constants (in internal units).
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
struct bpart* restrict bp, const struct black_holes_props* props,
const struct phys_const* constants, const struct cosmology* cosmo) {}
/**
* @brief Reset acceleration fields of a particle
*
......
......@@ -30,16 +30,15 @@
* @param pj Second particle (gas, not updated).
* @param xpj The extended data of the second particle (not updated).
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
const float hi, const float hj,
struct bpart *restrict bi,
const struct part *restrict pj,
const struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {
runner_iact_nonsym_bh_gas_density(
const float r2, const float *dx, const float hi, const float hj,
struct bpart *restrict bi, const struct part *restrict pj,
const struct xpart *restrict xpj, const struct cosmology *cosmo,
const struct gravity_props *grav_props, const integertime_t ti_current) {
float wi, wi_dx;
......@@ -80,16 +79,15 @@ runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
const float hi, const float hj,
const struct bpart *restrict bi,
struct part *restrict pj,
struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {}
runner_iact_nonsym_bh_gas_swallow(
const float r2, const float *dx, const float hi, const float hj,
const struct bpart *restrict bi, struct part *restrict pj,
struct xpart *restrict xpj, const struct cosmology *cosmo,
const struct gravity_props *grav_props, const integertime_t ti_current) {}
/**
* @brief Swallowing interaction between two BH particles (non-symmetric).
......@@ -104,7 +102,7 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
* @param bi First particle (black hole).
* @param bj Second particle (black hole)
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...)
* @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
......@@ -127,16 +125,15 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
const float hi, const float hj,
struct bpart *restrict bi,
struct part *restrict pj,
struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {
runner_iact_nonsym_bh_gas_feedback(
const float r2, const float *dx, const float hi, const float hj,
struct bpart *restrict bi, struct part *restrict pj,
struct xpart *restrict xpj, const struct cosmology *cosmo,
const struct gravity_props *grav_props, const integertime_t ti_current) {
#ifdef DEBUG_INTERACTIONS_BH
/* Update ngb counters */
if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_BH)
......
......@@ -102,6 +102,7 @@ INLINE static void convert_bpart_vel(const struct engine* e,
* @param bparts The b-particle array.
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
* @param with_cosmology Are we running a cosmological simulation?
*/
INLINE static void black_holes_write_particles(const struct bpart* bparts,
struct io_props* list,
......
......@@ -24,6 +24,7 @@
#include "black_holes_struct.h"
#include "cosmology.h"
#include "dimension.h"
#include "gravity.h"
#include "kernel_hydro.h"
#include "minmax.h"
#include "physical_constants.h"
......@@ -89,16 +90,43 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart(
bp->circular_velocity_gas[2] = 0.f;
bp->ngb_mass = 0.f;
bp->num_ngbs = 0;
bp->reposition.x[0] = -FLT_MAX;
bp->reposition.x[1] = -FLT_MAX;
bp->reposition.x[2] = -FLT_MAX;
bp->reposition.min_potential = FLT_MAX;
}
/**
* @brief Predict additional particle fields forward in time when drifting
*
* The fields do not get predicted but we move the BH to its new position
* if a new one was calculated in the repositioning loop.
*
* @param bp The particle
* @param dt_drift The drift time-step for positions.
*/
__attribute__((always_inline)) INLINE static void black_holes_predict_extra(
struct bpart* restrict bp, float dt_drift) {}
struct bpart* restrict bp, float dt_drift) {
/* Are we doing some repositioning? */
if (bp->reposition.min_potential != FLT_MAX) {
#ifdef SWIFT_DEBUG_CHECKS
if (bp->reposition.x[0] == -FLT_MAX || bp->reposition.x[1] == -FLT_MAX ||
bp->reposition.x[2] == -FLT_MAX) {
error("Something went wrong with the new repositioning position");
}
#endif
bp->x[0] = bp->reposition.x[0];
bp->x[1] = bp->reposition.x[1];
bp->x[2] = bp->reposition.x[2];
bp->gpart->x[0] = bp->reposition.x[0];
bp->gpart->x[1] = bp->reposition.x[1];
bp->gpart->x[2] = bp->reposition.x[2];
}
}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
......@@ -430,6 +458,35 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
}
}
/**
* @brief Finish the calculation of the new BH position.
*
* Here, we check that the BH should indeed be moved in the next drift.
*
* @param bp The black hole particle.
* @param props The properties of the black hole scheme.
* @param constants The physical constants (in internal units).
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
struct bpart* restrict bp, const struct black_holes_props* props,
const struct phys_const* constants, const struct cosmology* cosmo) {
const float potential = gravity_get_comoving_potential(bp->gpart);
/* Is the potential lower (i.e. the BH is at the bottom already)
* OR is the BH massive enough that we don't reposition? */
if (potential < bp->reposition.min_potential ||
bp->subgrid_mass > props->max_reposition_mass) {
/* No need to reposition */
bp->reposition.min_potential = FLT_MAX;
bp->reposition.x[0] = -FLT_MAX;
bp->reposition.x[1] = -FLT_MAX;
bp->reposition.x[2] = -FLT_MAX;
}
}
/**
* @brief Reset acceleration fields of a particle
*
......
......@@ -20,6 +20,7 @@
#define SWIFT_EAGLE_BH_IACT_H
/* Local includes */
#include "gravity.h"
#include "hydro.h"
#include "random.h"
#include "space.h"
......@@ -35,16 +36,15 @@
* @param pj Second particle (gas, not updated).
* @param xpj The extended data of the second particle (not updated).
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
const float hi, const float hj,
struct bpart *restrict bi,
const struct part *restrict pj,
const struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {
runner_iact_nonsym_bh_gas_density(
const float r2, const float *dx, const float hi, const float hj,
struct bpart *restrict bi, const struct part *restrict pj,
const struct xpart *restrict xpj, const struct cosmology *cosmo,
const struct gravity_props *grav_props, const integertime_t ti_current) {
float wi, wi_dx;
......@@ -117,16 +117,15 @@ runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
const float hi, const float hj,
const struct bpart *restrict bi,
struct part *restrict pj,
struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {
runner_iact_nonsym_bh_gas_swallow(
const float r2, const float *dx, const float hi, const float hj,
struct bpart *restrict bi, struct part *restrict pj,
struct xpart *restrict xpj, const struct cosmology *cosmo,
const struct gravity_props *grav_props, const integertime_t ti_current) {
float wi;
......@@ -140,6 +139,44 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
const float ui = r * hi_inv;
kernel_eval(ui, &wi);
/* Start by checking the repositioning criteria */
/* Note the factor 9 is taken from EAGLE. Will be turned into a parameter */
const float max_dist_repos2 =
kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv * 9.f *
grav_props->epsilon_cur2;
/* This gas neighbour is close enough that we can consider it's potential
for repositioning */
if (r2 < max_dist_repos2) {
/* Compute relative peculiar velocity between the two BHs
* Recall that in SWIFT v is (v_pec * a) */
const float delta_v[3] = {bi->v[0] - pj->v[0], bi->v[1] - pj->v[1],
bi->v[2] - pj->v[2]};
const float v2 = delta_v[0] * delta_v[0] + delta_v[1] * delta_v[1] +
delta_v[2] * delta_v[2];
const float v2_pec = v2 * cosmo->a2_inv;
/* Check the velocity criterion */
if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) {
const float potential = gravity_get_comoving_potential(pj->gpart);
/* Is the potential lower? */
if (potential < bi->reposition.min_potential) {
/* Store this as our new best */
bi->reposition.min_potential = potential;
bi->reposition.x[0] = pj->x[0];
bi->reposition.x[1] = pj->x[1];
bi->reposition.x[2] = pj->x[2];
}
}
}
/* Is the BH hungry? */
if (bi->subgrid_mass > bi->mass) {
......@@ -188,13 +225,13 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
* @param bi First particle (black hole).
* @param bj Second particle (black hole)
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...)
* @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const float hi, const float hj,
const struct bpart *restrict bi,
struct bpart *restrict bi,
struct bpart *restrict bj,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
......@@ -209,6 +246,33 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const float v2_pec = v2 * cosmo->a2_inv;
/* Note the factor 9 is taken from EAGLE. Will be turned into a parameter */
const float max_dist_repos2 =
kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv * 9.f *
grav_props->epsilon_cur2;
/* This gas neighbour is close enough that we can consider it's potential
for repositioning */
if (r2 < max_dist_repos2) {
/* Check the velocity criterion */
if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) {
const float potential = gravity_get_comoving_potential(bj->gpart);
/* Is the potential lower? */
if (potential < bi->reposition.min_potential) {
/* Store this as our new best */
bi->reposition.min_potential = potential;
bi->reposition.x[0] = bj->x[0];
bi->reposition.x[1] = bj->x[1];
bi->reposition.x[2] = bj->x[2];
}
}
}
/* Find the most massive of the two BHs */
float M = bi->subgrid_mass;
float h = hi;
......@@ -218,10 +282,11 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
}
/* Note the factor 9 is taken from EAGLE. Will be turned into a parameter */
const float max_distance2 = kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv *
grav_props->epsilon_baryon_cur *
grav_props->epsilon_baryon_cur * 9.f;
const float max_dist_merge2 =
kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv * 9.f *
grav_props->epsilon_baryon_cur *
grav_props->epsilon_baryon_cur;
const float G_Newton = grav_props->G_Newton;
......@@ -235,7 +300,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
/* Merge if gravitationally bound AND if within max distance
* Note that we use the kernel support here as the size and not just the
* smoothing length */
if (v2_pec < G_Newton * M / (kernel_gamma * h) && (r2 < max_distance2)) {
if (v2_pec < G_Newton * M / (kernel_gamma * h) && (r2 < max_dist_merge2)) {
/* This particle is swallowed by the BH with the largest ID of all the
* candidates wanting to swallow it */
......@@ -270,16 +335,15 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
* @param cosmo The cosmological model.
* @param grav_props The properties of the gravity scheme (softening, G, ...).
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
const float hi, const float hj,
const struct bpart *restrict bi,
struct part *restrict pj,
struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {
runner_iact_nonsym_bh_gas_feedback(
const float r2, const float *dx, const float hi, const float hj,
const struct bpart *restrict bi, struct part *restrict pj,
struct xpart *restrict xpj, const struct cosmology *cosmo,
const struct gravity_props *grav_props, const integertime_t ti_current) {
/* Get the heating probability */
const float prob = bi->to_distribute.AGN_heating_probability;
......
......@@ -104,6 +104,7 @@ INLINE static void convert_bpart_vel(const struct engine* e,
* @param bparts The b-particle array.
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
* @param with_cosmology Are we running a cosmological simulation?
*/
INLINE static void black_holes_write_particles(const struct bpart* bparts,
struct io_props* list,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment