Commit ed8321ea authored by Loic Hausammann's avatar Loic Hausammann Committed by Matthieu Schaller
Browse files

Gear pressure floor

parent 209dc39d
...@@ -1327,6 +1327,7 @@ with_subgrid_cooling=none ...@@ -1327,6 +1327,7 @@ with_subgrid_cooling=none
with_subgrid_chemistry=none with_subgrid_chemistry=none
with_subgrid_tracers=none with_subgrid_tracers=none
with_subgrid_entropy_floor=none with_subgrid_entropy_floor=none
with_subgrid_pressure_floor=none
with_subgrid_stars=none with_subgrid_stars=none
with_subgrid_star_formation=none with_subgrid_star_formation=none
with_subgrid_feedback=none with_subgrid_feedback=none
...@@ -1338,10 +1339,9 @@ case "$with_subgrid" in ...@@ -1338,10 +1339,9 @@ case "$with_subgrid" in
none) none)
;; ;;
GEAR) GEAR)
with_subgrid_cooling=grackle with_subgrid_cooling=grackle_0
with_subgrid_chemistry=GEAR with_subgrid_chemistry=GEAR
with_subgrid_tracers=none with_subgrid_pressure_floor=GEAR
with_subgrid_entropy_floor=none
with_subgrid_stars=GEAR with_subgrid_stars=GEAR
with_subgrid_star_formation=GEAR with_subgrid_star_formation=GEAR
with_subgrid_feedback=none with_subgrid_feedback=none
...@@ -1646,7 +1646,8 @@ esac ...@@ -1646,7 +1646,8 @@ esac
# Cooling function # Cooling function
AC_ARG_WITH([cooling], AC_ARG_WITH([cooling],
[AS_HELP_STRING([--with-cooling=<function>], [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="$withval"],
[with_cooling="none"] [with_cooling="none"]
...@@ -1673,21 +1674,10 @@ case "$with_cooling" in ...@@ -1673,21 +1674,10 @@ case "$with_cooling" in
compton) compton)
AC_DEFINE([COOLING_COMPTON], [1], [Compton cooling off the CMB]) AC_DEFINE([COOLING_COMPTON], [1], [Compton cooling off the CMB])
;; ;;
grackle) grackle_*)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library]) AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [0], [Grackle chemistry network, mode 0]) primordial_chemistry=${with_cooling:8}
;; AC_DEFINE_UNQUOTED([COOLING_GRACKLE_MODE], [$primordial_chemistry], [Grackle chemistry network])
grackle1)
AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
AC_DEFINE([COOLING_GRACKLE_MODE], [1], [Grackle chemistry network, mode 1])
;;
grackle2)
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])
;; ;;
EAGLE) EAGLE)
AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model]) AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model])
...@@ -1922,6 +1912,35 @@ case "$with_entropy_floor" in ...@@ -1922,6 +1912,35 @@ case "$with_entropy_floor" in
;; ;;
esac 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 # Star formation
AC_ARG_WITH([star-formation], AC_ARG_WITH([star-formation],
[AS_HELP_STRING([--with-star-formation=<sfm>], [AS_HELP_STRING([--with-star-formation=<sfm>],
...@@ -2056,6 +2075,7 @@ AC_MSG_RESULT([ ...@@ -2056,6 +2075,7 @@ AC_MSG_RESULT([
External potential : $with_potential External potential : $with_potential
Entropy floor : $with_entropy_floor Entropy floor : $with_entropy_floor
Pressure floor : $with_pressure_floor
Cooling function : $with_cooling Cooling function : $with_cooling
Chemistry : $with_chemistry Chemistry : $with_chemistry
Tracers : $with_tracers Tracers : $with_tracers
......
...@@ -5,6 +5,29 @@ ...@@ -5,6 +5,29 @@
GEAR model 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 Cooling: Grackle
~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~
......
...@@ -99,3 +99,6 @@ GrackleCooling: ...@@ -99,3 +99,6 @@ GrackleCooling:
GearChemistry: GearChemistry:
InitialMetallicity: 0.01295 InitialMetallicity: 0.01295
GEARPressureFloor:
Jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
...@@ -777,6 +777,13 @@ int main(int argc, char *argv[]) { ...@@ -777,6 +777,13 @@ int main(int argc, char *argv[]) {
else else
bzero(&entropy_floor, sizeof(struct entropy_floor_properties)); 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 */ /* Initialise the stars properties */
if (with_stars) if (with_stars)
stars_props_init(&stars_properties, &prog_const, &us, params, stars_props_init(&stars_properties, &prog_const, &us, params,
......
...@@ -285,6 +285,11 @@ EAGLEEntropyFloor: ...@@ -285,6 +285,11 @@ EAGLEEntropyFloor:
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. 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 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 ---------------------------------------------- # Parameters related to cooling function ----------------------------------------------
# Constant du/dt cooling function # Constant du/dt cooling function
......
...@@ -79,7 +79,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c ...@@ -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 \ collectgroup.c hydro_space.c equation_of_state.c \
chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.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 \ 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) $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES)
# Include files for distribution, not installation. # Include files for distribution, not installation.
......
...@@ -52,7 +52,7 @@ static gr_float cooling_time( ...@@ -52,7 +52,7 @@ static gr_float cooling_time(
const struct unit_system* restrict us, const struct unit_system* restrict us,
const struct cosmology* restrict cosmo, const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling, const struct cooling_function_data* restrict cooling,
const struct part* restrict p, const struct xpart* restrict xp); const struct part* restrict p, struct xpart* restrict xp);
static gr_float cooling_new_energy( static gr_float cooling_new_energy(
const struct phys_const* restrict phys_const, const struct phys_const* restrict phys_const,
const struct unit_system* restrict us, const struct unit_system* restrict us,
...@@ -573,7 +573,7 @@ __attribute__((always_inline)) INLINE static gr_float cooling_time( ...@@ -573,7 +573,7 @@ __attribute__((always_inline)) INLINE static gr_float cooling_time(
const struct unit_system* restrict us, const struct unit_system* restrict us,
const struct cosmology* restrict cosmo, const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling, const struct cooling_function_data* restrict cooling,
const struct part* restrict p, const struct xpart* restrict xp) { const struct part* restrict p, struct xpart* restrict xp) {
/* set current time */ /* set current time */
code_units units = cooling->units; code_units units = cooling->units;
......
...@@ -41,6 +41,7 @@ ...@@ -41,6 +41,7 @@
#include "hydro_space.h" #include "hydro_space.h"
#include "kernel_hydro.h" #include "kernel_hydro.h"
#include "minmax.h" #include "minmax.h"
#include "pressure_floor.h"
#include "./hydro_parameters.h" #include "./hydro_parameters.h"
...@@ -109,7 +110,8 @@ hydro_get_drifted_physical_internal_energy(const struct part *restrict p, ...@@ -109,7 +110,8 @@ hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure( __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
const struct part *restrict p) { const struct part *restrict p) {
return gas_pressure_from_entropy(p->rho, p->entropy); const float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
return pressure_floor_get_pressure(p, p->rho, comoving_pressure);
} }
/** /**
...@@ -121,7 +123,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure( ...@@ -121,7 +123,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure( __attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
const struct part *restrict p, const struct cosmology *cosmo) { const struct part *restrict p, const struct cosmology *cosmo) {
return gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy); const float phys_pressure =
gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
const float phys_rho = hydro_get_physical_density(p, cosmo);
return pressure_floor_get_pressure(p, phys_rho, phys_pressure);
} }
/** /**
...@@ -388,13 +393,15 @@ hydro_set_drifted_physical_internal_energy(struct part *p, ...@@ -388,13 +393,15 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const float rho_inv = 1.f / p->rho; const float rho_inv = 1.f / p->rho;
/* Compute the pressure */ /* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the sound speed */ /* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); const float soundspeed =
gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */ /* Divide the pressure by the density squared to get the SPH term */
const float P_over_rho2 = pressure * rho_inv * rho_inv; const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
/* Update variables. */ /* Update variables. */
p->force.P_over_rho2 = P_over_rho2; p->force.P_over_rho2 = P_over_rho2;
...@@ -591,13 +598,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -591,13 +598,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float abs_div_physical_v = fabsf(div_physical_v); const float abs_div_physical_v = fabsf(div_physical_v);
/* Compute the pressure */ /* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the sound speed */ /* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); const float soundspeed =
gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */ /* Divide the pressure by the density squared to get the SPH term */
const float P_over_rho2 = pressure * rho_inv * rho_inv; const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
/* Compute the Balsara switch */ /* Compute the Balsara switch */
/* Pre-multiply in the AV factor; hydro_props are not passed to the iact /* Pre-multiply in the AV factor; hydro_props are not passed to the iact
...@@ -665,14 +674,16 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( ...@@ -665,14 +674,16 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->entropy = xp->entropy_full; p->entropy = xp->entropy_full;
/* Re-compute the pressure */ /* Re-compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the new sound speed */ /* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); const float soundspeed =
gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */ /* Divide the pressure by the density squared to get the SPH term */
const float rho_inv = 1.f / p->rho; const float rho_inv = 1.f / p->rho;
const float P_over_rho2 = pressure * rho_inv * rho_inv; const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
/* Update variables */ /* Update variables */
p->force.soundspeed = soundspeed; p->force.soundspeed = soundspeed;
...@@ -730,14 +741,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -730,14 +741,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf(w2); p->rho *= expf(w2);
/* Re-compute the pressure */ /* Re-compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the new sound speed */ /* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); const float soundspeed =
gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */ /* Divide the pressure by the density squared to get the SPH term */
const float rho_inv = 1.f / p->rho; const float rho_inv = 1.f / p->rho;
const float P_over_rho2 = pressure * rho_inv * rho_inv; const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
/* Update variables */ /* Update variables */
p->force.soundspeed = soundspeed; p->force.soundspeed = soundspeed;
...@@ -840,14 +853,16 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( ...@@ -840,14 +853,16 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
} }
/* Compute the pressure */ /* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
/* Compute the sound speed */ /* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); const float soundspeed =
gas_soundspeed_from_pressure(p->rho, comoving_pressure);
/* Divide the pressure by the density squared to get the SPH term */ /* Divide the pressure by the density squared to get the SPH term */
const float rho_inv = 1.f / p->rho; const float rho_inv = 1.f / p->rho;
const float P_over_rho2 = pressure * rho_inv * rho_inv; const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
p->force.soundspeed = soundspeed; p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = P_over_rho2; p->force.P_over_rho2 = P_over_rho2;
......
...@@ -46,6 +46,7 @@ ...@@ -46,6 +46,7 @@
#include "hydro_space.h" #include "hydro_space.h"
#include "kernel_hydro.h" #include "kernel_hydro.h"
#include "minmax.h" #include "minmax.h"
#include "pressure_floor.h"
#include "./hydro_parameters.h" #include "./hydro_parameters.h"
...@@ -221,7 +222,9 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) { ...@@ -221,7 +222,9 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) {
/* Compute the sound speed -- see theory section for justification */ /* Compute the sound speed -- see theory section for justification */
/* IDEAL GAS ONLY -- P-U does not work with generic EoS. */ /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho); const float comoving_pressure =
pressure_floor_get_pressure(p, p->rho, p->pressure_bar);
const float square_rooted = sqrtf(hydro_gamma * comoving_pressure / p->rho);
return square_rooted; return square_rooted;
} }
...@@ -236,7 +239,10 @@ __attribute__((always_inline)) INLINE static float ...@@ -236,7 +239,10 @@ __attribute__((always_inline)) INLINE static float
hydro_get_physical_soundspeed(const struct part *restrict p, hydro_get_physical_soundspeed(const struct part *restrict p,
const struct cosmology *cosmo) { const struct cosmology *cosmo) {
return cosmo->a_factor_sound_speed * p->force.soundspeed; const float phys_rho = hydro_get_physical_density(p, cosmo);
return pressure_floor_get_pressure(
p, phys_rho, cosmo->a_factor_sound_speed * p->force.soundspeed);
} }
/** /**
...@@ -640,10 +646,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -640,10 +646,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
hydro_one_over_gamma_minus_one) / hydro_one_over_gamma_minus_one) /
(1.f + common_factor * p->density.wcount_dh); (1.f + common_factor * p->density.wcount_dh);
/* Get the pressures */
const float comoving_pressure_with_floor =
pressure_floor_get_pressure(p, p->rho, p->pressure_bar);
/* Update variables. */ /* Update variables. */
p->force.f = grad_h_term; p->force.f = grad_h_term;
p->force.soundspeed = soundspeed; p->force.soundspeed = soundspeed;
p->force.balsara = balsara; p->force.balsara = balsara;
p->force.pressure_bar_with_floor = comoving_pressure_with_floor;
} }
/** /**
...@@ -755,6 +766,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -755,6 +766,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const float soundspeed = hydro_get_comoving_soundspeed(p); const float soundspeed = hydro_get_comoving_soundspeed(p);
p->force.soundspeed = soundspeed; p->force.soundspeed = soundspeed;
/* update the required variables */
const float comoving_pressure_with_floor =
pressure_floor_get_pressure(p, p->rho, p->pressure_bar);
p->force.pressure_bar_with_floor = comoving_pressure_with_floor;
} }
/** /**
......
...@@ -259,11 +259,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -259,11 +259,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Convolve with the kernel */ /* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* Compute the ratio of pressures */
const float pressure_inverse_i =
pi->force.pressure_bar_with_floor / (pi->pressure_bar * pi->pressure_bar);
const float pressure_inverse_j =
pj->force.pressure_bar_with_floor / (pj->pressure_bar * pj->pressure_bar);
/* SPH acceleration term */ /* SPH acceleration term */
const float sph_acc_term = const float sph_acc_term = pj->u * pi->u * hydro_gamma_minus_one *
pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one * hydro_gamma_minus_one *
((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) * ((f_ij * pressure_inverse_i) * wi_dr +
r_inv; (f_ji * pressure_inverse_j) * wj_dr) *
r_inv;
/* Assemble the acceleration */ /* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term; const float acc = sph_acc_term + visc_acc_term;
...@@ -278,11 +285,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -278,11 +285,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->a_hydro[2] += mi * acc * dx[2]; pj->a_hydro[2] += mi * acc * dx[2];
/* Get the time derivative for u. */ /* Get the time derivative for u. */
const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one * const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
pj->u * pi->u * (f_ij / pi->pressure_bar) * pj->u * pi->u * (f_ij * pressure_inverse_i) *
wi_dr * dvdr * r_inv; wi_dr * dvdr * r_inv;
const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one * const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one *
pi->u * pj->u * (f_ji / pj->pressure_bar) * pi->u * pj->u * (f_ji * pressure_inverse_j) *
wj_dr * dvdr * r_inv; wj_dr * dvdr * r_inv;
/* Viscosity term */ /* Viscosity term */
...@@ -386,11 +395,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -386,11 +395,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Convolve with the kernel */ /* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* Compute the ratio of pressures */
const float pressure_inverse_i =
pi->force.pressure_bar_with_floor / (pi->pressure_bar * pi->pressure_bar);
const float pressure_inverse_j =
pj->force.pressure_bar_with_floor / (pj->pressure_bar * pj->pressure_bar);
/* SPH acceleration term */ /* SPH acceleration term */
const float sph_acc_term = const float sph_acc_term = pj->u * pi->u * hydro_gamma_minus_one *
pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one * hydro_gamma_minus_one *
((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) * ((f_ij * pressure_inverse_i) * wi_dr +
r_inv; (f_ji * pressure_inverse_j) * wj_dr) *
r_inv;
/* Assemble the acceleration */ /* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term; const float acc = sph_acc_term + visc_acc_term;
...@@ -402,7 +418,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -402,7 +418,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Get the time derivative for u. */ /* Get the time derivative for u. */
const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one * const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
pj->u * pi->u * (f_ij / pi->pressure_bar) * pj->u * pi->u * (f_ij * pressure_inverse_i) *
wi_dr * dvdr * r_inv; wi_dr * dvdr * r_inv;
/* Viscosity term */ /* Viscosity term */
......
...@@ -168,6 +168,9 @@ struct part { ...@@ -168,6 +168,9 @@ struct part {
/*! Balsara switch */ /*! Balsara switch */
float balsara; float balsara;
/*! Pressure term including the pressure floor */
float pressure_bar_with_floor;