Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
SWIFT
SWIFTsim
Commits
fa3dab63
Commit
fa3dab63
authored
Aug 09, 2019
by
Matthieu Schaller
Browse files
Merge branch 'gear_pressure_floor' into 'master'
Gear pressure floor See merge request
!811
parents
209dc39d
ed8321ea
Changes
17
Hide whitespace changes
Inline
Side-by-side
configure.ac
View file @
fa3dab63
...
...
@@ -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
...
...
@@ -1646,7 +1646,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"]
...
...
@@ -1673,21 +1674,10 @@ case "$with_cooling" in
compton)
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_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)
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])
...
...
@@ -1922,6 +1912,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>],
...
...
@@ -2056,6 +2075,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
...
...
doc/RTD/source/SubgridModels/GEAR/index.rst
View file @
fa3dab63
...
...
@@ -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
~~~~~~~~~~~~~~~~
...
...
examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
View file @
fa3dab63
...
...
@@ -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.
examples/main.c
View file @
fa3dab63
...
...
@@ -777,6 +777,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
,
...
...
examples/parameter_example.yml
View file @
fa3dab63
...
...
@@ -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_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
...
...
src/Makefile.am
View file @
fa3dab63
...
...
@@ -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.
...
...
src/cooling/grackle/cooling.h
View file @
fa3dab63
...
...
@@ -52,7 +52,7 @@ static gr_float cooling_time(
const
struct
unit_system
*
restrict
us
,
const
struct
cosmology
*
restrict
cosmo
,
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
(
const
struct
phys_const
*
restrict
phys_const
,
const
struct
unit_system
*
restrict
us
,
...
...
@@ -573,7 +573,7 @@ __attribute__((always_inline)) INLINE static gr_float cooling_time(
const
struct
unit_system
*
restrict
us
,
const
struct
cosmology
*
restrict
cosmo
,
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 */
code_units
units
=
cooling
->
units
;
...
...
src/hydro/Gadget2/hydro.h
View file @
fa3dab63
...
...
@@ -41,6 +41,7 @@
#include
"hydro_space.h"
#include
"kernel_hydro.h"
#include
"minmax.h"
#include
"pressure_floor.h"
#include
"./hydro_parameters.h"
...
...
@@ -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
(
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(
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_physical_pressure
(
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,
const
float
rho_inv
=
1
.
f
/
p
->
rho
;
/* 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 */
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 */
const
float
P_over_rho2
=
pressure
*
rho_inv
*
rho_inv
;
const
float
P_over_rho2
=
comoving_
pressure
*
rho_inv
*
rho_inv
;
/* Update variables. */
p
->
force
.
P_over_rho2
=
P_over_rho2
;
...
...
@@ -591,13 +598,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const
float
abs_div_physical_v
=
fabsf
(
div_physical_v
);
/* 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 */
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 */
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 */
/* 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(
p
->
entropy
=
xp
->
entropy_full
;
/* 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 */
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 */
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 */
p
->
force
.
soundspeed
=
soundspeed
;
...
...
@@ -730,14 +741,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p
->
rho
*=
expf
(
w2
);
/* 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 */
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 */
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 */
p
->
force
.
soundspeed
=
soundspeed
;
...
...
@@ -840,14 +853,16 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
}
/* 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 */
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 */
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
.
P_over_rho2
=
P_over_rho2
;
...
...
src/hydro/PressureEnergy/hydro.h
View file @
fa3dab63
...
...
@@ -46,6 +46,7 @@
#include
"hydro_space.h"
#include
"kernel_hydro.h"
#include
"minmax.h"
#include
"pressure_floor.h"
#include
"./hydro_parameters.h"
...
...
@@ -221,7 +222,9 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) {
/* Compute the sound speed -- see theory section for justification */
/* 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
;
}
...
...
@@ -236,7 +239,10 @@ __attribute__((always_inline)) INLINE static float
hydro_get_physical_soundspeed
(
const
struct
part
*
restrict
p
,
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(
hydro_one_over_gamma_minus_one
)
/
(
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. */
p
->
force
.
f
=
grad_h_term
;
p
->
force
.
soundspeed
=
soundspeed
;
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(
const
float
soundspeed
=
hydro_get_comoving_soundspeed
(
p
);
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
;
}
/**
...
...
src/hydro/PressureEnergy/hydro_iact.h
View file @
fa3dab63
...
...
@@ -259,11 +259,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Convolve with the kernel */
const
float
visc_acc_term
=
0
.
5
f
*
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 */
const
float
sph_acc_term
=
pj
->
u
*
pi
->
u
*
hydro_gamma_minus_one
*
hydro_gamma_minus_one
*
((
f_ij
/
pi
->
pressure_bar
)
*
wi_dr
+
(
f_ji
/
pj
->
pressure_bar
)
*
wj_dr
)
*
r_inv
;
const
float
sph_acc_term
=
pj
->
u
*
pi
->
u
*
hydro_gamma_minus_one
*
hydro_gamma_minus_one
*
((
f_ij
*
pressure_inverse_i
)
*
wi_dr
+
(
f_ji
*
pressure_inverse_j
)
*
wj_dr
)
*
r_inv
;
/* Assemble the acceleration */
const
float
acc
=
sph_acc_term
+
visc_acc_term
;
...
...
@@ -278,11 +285,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj
->
a_hydro
[
2
]
+=
mi
*
acc
*
dx
[
2
];
/* Get the time derivative for u. */
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
;
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
;
/* Viscosity term */
...
...
@@ -386,11 +395,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Convolve with the kernel */
const
float
visc_acc_term
=
0
.
5
f
*
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 */
const
float
sph_acc_term
=
pj
->
u
*
pi
->
u
*
hydro_gamma_minus_one
*
hydro_gamma_minus_one
*
((
f_ij
/
pi
->
pressure_bar
)
*
wi_dr
+
(
f_ji
/
pj
->
pressure_bar
)
*
wj_dr
)
*
r_inv
;
const
float
sph_acc_term
=
pj
->
u
*
pi
->
u
*
hydro_gamma_minus_one
*
hydro_gamma_minus_one
*
((
f_ij
*
pressure_inverse_i
)
*
wi_dr
+
(
f_ji
*
pressure_inverse_j
)
*
wj_dr
)
*
r_inv
;
/* Assemble the acceleration */
const
float
acc
=
sph_acc_term
+
visc_acc_term
;
...
...
@@ -402,7 +418,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Get the time derivative for u. */
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
;
/* Viscosity term */
...
...
src/hydro/PressureEnergy/hydro_part.h
View file @
fa3dab63
...
...
@@ -168,6 +168,9 @@ struct part {
/*! Balsara switch */
float
balsara
;
/*! Pressure term including the pressure floor */
float
pressure_bar_with_floor
;
}
force
;
};
...
...
src/hydro_properties.c
View file @
fa3dab63
...
...
@@ -34,6 +34,7 @@
#include
"hydro.h"
#include
"kernel_hydro.h"
#include
"parser.h"
#include
"pressure_floor.h"
#include
"units.h"
#define hydro_props_default_max_iterations 30
...
...
@@ -186,6 +187,9 @@ void hydro_props_print(const struct hydro_props *p) {
/* Print equation of state first */
eos_print
(
&
eos
);
/* Then the pressure floor */
pressure_floor_print
(
&
pressure_floor_props
);
/* Now SPH */
message
(
"Hydrodynamic scheme: %s in %dD."
,
SPH_IMPLEMENTATION
,
(
int
)
hydro_dimension
);
...
...
@@ -238,6 +242,7 @@ void hydro_props_print(const struct hydro_props *p) {
void
hydro_props_print_snapshot
(
hid_t
h_grpsph
,
const
struct
hydro_props
*
p
)
{
eos_print_snapshot
(
h_grpsph
,
&
eos
);
pressure_floor_print_snapshot
(
h_grpsph
);
io_write_attribute_i
(
h_grpsph
,
"Dimension"
,
(
int
)
hydro_dimension
);
io_write_attribute_s
(
h_grpsph
,
"Scheme"
,
SPH_IMPLEMENTATION
);
...
...
src/pressure_floor.c
0 → 100644
View file @
fa3dab63
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Loic Hausammann (loic.hausammann@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/>.
*
******************************************************************************/
/* This object's header. */
#include
"pressure_floor.h"
/* Pressure floor for the physics model. */
struct
pressure_floor_properties
pressure_floor_props
;
src/pressure_floor.h
0 → 100644
View file @
fa3dab63
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
*
* 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_PRESSURE_FLOOR_H
#define SWIFT_PRESSURE_FLOOR_H
/**
* @file src/pressure_floor.h
* @brief Branches between the different pressure floor models
*/
/* Config parameters. */
#include
"../config.h"
/* Local includes */
#include
"common_io.h"
#include
"cosmology.h"
#include
"error.h"
#include
"inline.h"
extern
struct
pressure_floor_properties
pressure_floor_props
;
/* Check if pressure floor is implemented in hydro */
#ifndef PRESSURE_FLOOR_NONE
#if defined(GADGET2_SPH) || defined(HOPKINS_PU_SPH)
/* Implemented */
#else
#error Pressure floor not implemented with this hydro scheme
#endif
#endif
/* Import the right pressure floor definition */
#if defined(PRESSURE_FLOOR_NONE)
#include
"./pressure_floor/none/pressure_floor.h"
#elif defined(PRESSURE_FLOOR_GEAR)
#include
"./pressure_floor/GEAR/pressure_floor.h"
#endif
#endif
/* SWIFT_PRESSURE_FLOOR_H */
src/pressure_floor/GEAR/pressure_floor.h
0 → 100644
View file @
fa3dab63
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
*
* 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_PRESSURE_FLOOR_GEAR_H
#define SWIFT_PRESSURE_FLOOR_GEAR_H
#include
"adiabatic_index.h"
#include
"cosmology.h"
#include
"equation_of_state.h"
#include
"hydro_properties.h"
#include
"parser.h"
#include
"part.h"
#include
"units.h"
/**
* @file src/pressure_floor/GEAR/pressure_floor.h
* @brief Pressure floor used in the GEAR model
*/
/**
* @brief Properties of the pressure floor in the GEAR model.
*/
struct
pressure_floor_properties
{
/*! Jeans factor */
float
n_jeans
;
/*! The constants in internal units (4 G N_jeans^(2/3) / PI) */
float
constants
;
};
/**
* @brief Compute the physical pressure floor of a given #part.
*
* Note that the particle is not updated!!
*
* The inputs can be either in physical or comoving coordinates (the output is
* in the same coordinates).
*
* @param p The #part.
* @param rho The physical or comoving density.
* @param pressure The physical or comoving pressure without any pressure floor.
*
* @return The physical or comoving pressure with the floor.
*/