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
83c92087
Commit
83c92087
authored
Feb 27, 2018
by
Matthieu Schaller
Browse files
Added cosmological factors to the Minimal SPH scheme.
parent
a5a448db
Changes
5
Hide whitespace changes
Inline
Side-by-side
src/hydro/Minimal/hydro.h
View file @
83c92087
...
...
@@ -35,6 +35,7 @@
#include
"adiabatic_index.h"
#include
"approx_math.h"
#include
"cosmology.h"
#include
"dimension.h"
#include
"equation_of_state.h"
#include
"hydro_properties.h"
...
...
@@ -43,7 +44,7 @@
#include
"minmax.h"
/**
* @brief Returns the internal energy of a particle
* @brief Returns the
comoving
internal energy of a particle
*
* For implementations where the main thermodynamic variable
* is not internal energy, this function computes the internal
...
...
@@ -51,25 +52,61 @@
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_internal_energy
(
const
struct
part
*
restrict
p
)
{
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_comoving_internal_energy
(
const
struct
part
*
restrict
p
)
{
return
p
->
u
;
}
/**
* @brief Returns the pressure of a particle
* @brief Returns the physical internal energy of a particle
*
* For implementations where the main thermodynamic variable
* is not internal energy, this function computes the internal
* energy from the thermodynamic variable and converts it to
* physical coordinates.
*
* @param p The particle of interest.
* @param cosmo The cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_physical_internal_energy
(
const
struct
part
*
restrict
p
,
const
struct
cosmology
*
cosmo
)
{
return
p
->
u
*
cosmo
->
a_factor_internal_energy
;
}
/**
* @brief Returns the comoving pressure of a particle
*
* Computes the pressure based on the particle's properties.
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_pressure
(
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_
comoving_
pressure
(
const
struct
part
*
restrict
p
)
{
return
gas_pressure_from_internal_energy
(
p
->
rho
,
p
->
u
);
}
/**
* @brief Returns the entropy of a particle
* @brief Returns the physical pressure of a particle
*
* Computes the pressure based on the particle's properties and
* convert it to physical coordinates.
*
* @param p The particle of interest
* @param cosmo The cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_physical_pressure
(
const
struct
part
*
restrict
p
,
const
struct
cosmology
*
cosmo
)
{
return
cosmo
->
a_factor_pressure
*
gas_pressure_from_internal_energy
(
p
->
rho
,
p
->
u
);
}
/**
* @brief Returns the comoving entropy of a particle
*
* For implementations where the main thermodynamic variable
* is not entropy, this function computes the entropy from
...
...
@@ -77,34 +114,78 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_entropy
(
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_
comoving_
entropy
(
const
struct
part
*
restrict
p
)
{
return
gas_entropy_from_internal_energy
(
p
->
rho
,
p
->
u
);
}
/**
* @brief Returns the sound speed of a particle
* @brief Returns the physical entropy of a particle
*
* For implementations where the main thermodynamic variable
* is not entropy, this function computes the entropy from
* the thermodynamic variable and converts it to
* physical coordinates.
*
* @param p The particle of interest
* @param cosmo The cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_soundspeed
(
const
struct
part
*
restrict
p
)
{
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_physical_entropy
(
const
struct
part
*
restrict
p
,
const
struct
cosmology
*
cosmo
)
{
/* Note: no cosmological conversion required here with our choice of
* coordinates. */
return
gas_entropy_from_internal_energy
(
p
->
rho
,
p
->
u
);
}
/**
* @brief Returns the comoving sound speed of a particle
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_comoving_soundspeed
(
const
struct
part
*
restrict
p
)
{
return
p
->
force
.
soundspeed
;
}
/**
* @brief Returns the
density
of a particle
* @brief Returns the
physical sound speed
of a particle
*
* @param p The particle of interest
* @param cosmo The cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_density
(
__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
;
}
/**
* @brief Returns the comoving density of a particle
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_comoving_density
(
const
struct
part
*
restrict
p
)
{
return
p
->
rho
;
}
/**
* @brief Returns the comoving density of a particle.
*
* @param p The particle of interest
* @param cosmo The cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_physical_density
(
const
struct
part
*
restrict
p
,
const
struct
cosmology
*
cosmo
)
{
return
cosmo
->
a3_inv
*
p
->
rho
;
}
/**
* @brief Returns the mass of a particle
*
...
...
@@ -121,16 +202,20 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
*
* @param p The particle of interest
* @param xp The extended data of the particle.
* @param dt The time since the last kick.
* @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
* @param dt_kick_grav The time (for gravity accelerations) since the last kick.
* @param v (return) The velocities at the current time.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_get_drifted_velocities
(
const
struct
part
*
restrict
p
,
const
struct
xpart
*
xp
,
float
dt
,
float
v
[
3
])
{
v
[
0
]
=
xp
->
v_full
[
0
]
+
p
->
a_hydro
[
0
]
*
dt
;
v
[
1
]
=
xp
->
v_full
[
1
]
+
p
->
a_hydro
[
1
]
*
dt
;
v
[
2
]
=
xp
->
v_full
[
2
]
+
p
->
a_hydro
[
2
]
*
dt
;
const
struct
part
*
restrict
p
,
const
struct
xpart
*
xp
,
float
dt_kick_hydro
,
float
dt_kick_grav
,
float
v
[
3
])
{
v
[
0
]
=
xp
->
v_full
[
0
]
+
p
->
a_hydro
[
0
]
*
dt_kick_hydro
+
xp
->
a_grav
[
0
]
*
dt_kick_grav
;
v
[
1
]
=
xp
->
v_full
[
1
]
+
p
->
a_hydro
[
1
]
*
dt_kick_hydro
+
xp
->
a_grav
[
1
]
*
dt_kick_grav
;
v
[
2
]
=
xp
->
v_full
[
2
]
+
p
->
a_hydro
[
2
]
*
dt_kick_hydro
+
xp
->
a_grav
[
2
]
*
dt_kick_grav
;
}
/**
...
...
@@ -168,17 +253,18 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
* @param p Pointer to the particle data
* @param xp Pointer to the extended particle data
* @param hydro_properties The SPH parameters
*
*
@param cosmo The cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_compute_timestep
(
const
struct
part
*
restrict
p
,
const
struct
xpart
*
restrict
xp
,
const
struct
hydro_props
*
restrict
hydro_properties
)
{
const
struct
hydro_props
*
restrict
hydro_properties
,
const
struct
cosmology
*
restrict
cosmo
)
{
const
float
CFL_condition
=
hydro_properties
->
CFL_condition
;
/* CFL condition */
const
float
dt_cfl
=
2
.
f
*
kernel_gamma
*
CFL_condition
*
p
->
h
/
p
->
force
.
v_sig
;
const
float
dt_cfl
=
2
.
f
*
kernel_gamma
*
CFL_condition
*
cosmo
->
a
*
p
->
h
/
(
cosmo
->
a_factor_sound_speed
*
p
->
force
.
v_sig
)
;
return
dt_cfl
;
}
...
...
@@ -218,13 +304,15 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* Multiplies the density and number of neighbours by the appropiate constants
* and add the self-contribution term.
* Additional quantities such as velocity gradients will also get the final
*terms
* added to them here.
* terms added to them here.
*
* Also adds/multiplies the cosmological terms if need be.
*
* @param p The particle to act upon
* @param cosmo The cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_end_density
(
struct
part
*
restrict
p
)
{
struct
part
*
restrict
p
,
const
struct
cosmology
*
cosmo
)
{
/* Some smoothing length multiples. */
const
float
h
=
p
->
h
;
...
...
@@ -248,11 +336,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/**
* @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
*
* In the desperate case where a particle has no neighbours (likely because
* of the h_max ceiling), set the particle fields to something sensible to avoid
* NaNs in the next calculations.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param cosmo The cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_part_has_no_neighbours
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
)
{
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
,
const
struct
cosmology
*
cosmo
)
{
/* Some smoothing length multiples. */
const
float
h
=
p
->
h
;
...
...
@@ -278,9 +372,11 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param cosmo The current cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_prepare_force
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
)
{
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
,
const
struct
cosmology
*
cosmo
)
{
/* Compute the pressure */
const
float
pressure
=
gas_pressure_from_internal_energy
(
p
->
rho
,
p
->
u
);
...
...
@@ -346,17 +442,22 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* Additional hydrodynamic quantites are drifted forward in time here. These
* include thermal quantities (thermal energy or total energy or entropy, ...).
*
* Note the different time-step sizes used for the different quantities as they
* include cosmological factors.
*
* @param p The particle.
* @param xp The extended data of the particle.
* @param dt The drift time-step.
* @param dt_drift The drift time-step for positions.
* @param dt_therm The drift time-step for thermal quantities.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_predict_extra
(
struct
part
*
restrict
p
,
const
struct
xpart
*
restrict
xp
,
float
dt
)
{
struct
part
*
restrict
p
,
const
struct
xpart
*
restrict
xp
,
float
dt_drift
,
float
dt_therm
)
{
const
float
h_inv
=
1
.
f
/
p
->
h
;
/* Predict smoothing length */
const
float
w1
=
p
->
force
.
h_dt
*
h_inv
*
dt
;
const
float
w1
=
p
->
force
.
h_dt
*
h_inv
*
dt
_drift
;
if
(
fabsf
(
w1
)
<
0
.
2
f
)
p
->
h
*=
approx_expf
(
w1
);
/* 4th order expansion of exp(w) */
else
...
...
@@ -370,7 +471,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p
->
rho
*=
expf
(
w2
);
/* Predict the internal energy */
p
->
u
+=
p
->
u_dt
*
dt
;
p
->
u
+=
p
->
u_dt
*
dt
_therm
;
/* Compute the new pressure */
const
float
pressure
=
gas_pressure_from_internal_energy
(
p
->
rho
,
p
->
u
);
...
...
@@ -386,13 +487,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @brief Finishes the force calculation.
*
* Multiplies the force and accelerations by the appropiate constants
* and add the self-contribution term. In most cases, there is
nothing
* and add the self-contribution term. In most cases, there is
little
* to do here.
*
* Cosmological terms are also added/multiplied here.
*
* @param p The particle to act upon
* @param cosmo The current cosmological model.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_end_force
(
struct
part
*
restrict
p
)
{
struct
part
*
restrict
p
,
const
struct
cosmology
*
cosmo
)
{
p
->
force
.
h_dt
*=
p
->
h
*
hydro_dimension_inv
;
}
...
...
@@ -403,18 +507,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* Additional hydrodynamic quantites are kicked forward in time here. These
* include thermal quantities (thermal energy or total energy or entropy, ...).
*
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt The time-step for this kick
* @param p The particle to act upon
.
* @param xp The particle extended data to act upon
.
* @param dt
_therm
The time-step for this kick
(for thermodynamic quantities).
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_kick_extra
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
,
float
dt
)
{
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
,
float
dt
_therm
)
{
/* Do not decrease the energy by more than a factor of 2*/
if
(
dt
>
0
.
&&
p
->
u_dt
*
dt
<
-
0
.
5
f
*
xp
->
u_full
)
{
p
->
u_dt
=
-
0
.
5
f
*
xp
->
u_full
/
dt
;
if
(
dt
_therm
>
0
.
&&
p
->
u_dt
*
dt
_therm
<
-
0
.
5
f
*
xp
->
u_full
)
{
p
->
u_dt
=
-
0
.
5
f
*
xp
->
u_full
/
dt
_therm
;
}
xp
->
u_full
+=
p
->
u_dt
*
dt
;
xp
->
u_full
+=
p
->
u_dt
*
dt
_therm
;
/* Compute the pressure */
const
float
pressure
=
gas_pressure_from_internal_energy
(
p
->
rho
,
xp
->
u_full
);
...
...
src/hydro/Minimal/hydro_debug.h
View file @
83c92087
...
...
@@ -43,8 +43,9 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"time_bin=%d
\n
"
,
p
->
x
[
0
],
p
->
x
[
1
],
p
->
x
[
2
],
p
->
v
[
0
],
p
->
v
[
1
],
p
->
v
[
2
],
xp
->
v_full
[
0
],
xp
->
v_full
[
1
],
xp
->
v_full
[
2
],
p
->
a_hydro
[
0
],
p
->
a_hydro
[
1
],
p
->
a_hydro
[
2
],
p
->
u
,
p
->
u_dt
,
p
->
force
.
v_sig
,
p
->
force
.
pressure
,
p
->
h
,
p
->
force
.
h_dt
,
(
int
)
p
->
density
.
wcount
,
p
->
mass
,
p
->
density
.
rho_dh
,
p
->
rho
,
p
->
time_bin
);
p
->
u
,
p
->
u_dt
,
p
->
force
.
v_sig
,
hydro_get_comoving_pressure
(
p
),
p
->
h
,
p
->
force
.
h_dt
,
(
int
)
p
->
density
.
wcount
,
p
->
mass
,
p
->
density
.
rho_dh
,
p
->
rho
,
p
->
time_bin
);
}
#endif
/* SWIFT_MINIMAL_HYDRO_DEBUG_H */
src/hydro/Minimal/hydro_iact.h
View file @
83c92087
...
...
@@ -36,10 +36,20 @@
#include
"minmax.h"
/**
* @brief Density loop
* @brief Density interaction between two particles.
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi First particle.
* @param pj Second particle.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_density
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
restrict
pi
,
struct
part
*
restrict
pj
,
float
a
,
float
H
)
{
float
wi
,
wj
,
wi_dx
,
wj_dx
;
...
...
@@ -71,10 +81,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
}
/**
* @brief Density loop (non-symmetric version)
* @brief Density interaction between two particles (non-symmetric).
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi First particle.
* @param pj Second particle (not updated).
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_nonsym_density
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
restrict
pi
,
const
struct
part
*
restrict
pj
,
float
a
,
float
H
)
{
float
wi
,
wi_dx
;
...
...
@@ -95,12 +115,24 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
}
/**
* @brief Force loop
* @brief Force interaction between two particles.
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi First particle.
* @param pj Second particle.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_force
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
restrict
pi
,
struct
part
*
restrict
pj
,
float
a
,
float
H
)
{
const
float
fac_mu
=
1
.
f
;
/* Will change with cosmological integration */
/* Cosmological factors entering the EoMs */
const
float
fac_mu
=
pow_three_gamma_minus_five_over_two
(
a
);
const
float
a2_Hubble
=
a
*
a
*
H
;
const
float
r
=
sqrtf
(
r2
);
const
float
r_inv
=
1
.
0
f
/
r
;
...
...
@@ -136,7 +168,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute dv dot r. */
const
float
dvdr
=
(
pi
->
v
[
0
]
-
pj
->
v
[
0
])
*
dx
[
0
]
+
(
pi
->
v
[
1
]
-
pj
->
v
[
1
])
*
dx
[
1
]
+
(
pi
->
v
[
2
]
-
pj
->
v
[
2
])
*
dx
[
2
];
(
pi
->
v
[
2
]
-
pj
->
v
[
2
])
*
dx
[
2
]
+
a2_Hubble
*
r2
;
/* Are the particles moving towards each others ? */
const
float
omega_ij
=
min
(
dvdr
,
0
.
f
);
...
...
@@ -195,12 +227,24 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
}
/**
* @brief Force loop (non-symmetric version)
* @brief Force interaction between two particles (non-symmetric).
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi First particle.
* @param pj Second particle (not updated).
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_nonsym_force
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
restrict
pi
,
const
struct
part
*
restrict
pj
,
float
a
,
float
H
)
{
const
float
fac_mu
=
1
.
f
;
/* Will change with cosmological integration */
/* Cosmological factors entering the EoMs */
const
float
fac_mu
=
pow_three_gamma_minus_five_over_two
(
a
);
const
float
a2_Hubble
=
a
*
a
*
H
;
const
float
r
=
sqrtf
(
r2
);
const
float
r_inv
=
1
.
0
f
/
r
;
...
...
@@ -236,7 +280,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute dv dot r. */
const
float
dvdr
=
(
pi
->
v
[
0
]
-
pj
->
v
[
0
])
*
dx
[
0
]
+
(
pi
->
v
[
1
]
-
pj
->
v
[
1
])
*
dx
[
1
]
+
(
pi
->
v
[
2
]
-
pj
->
v
[
2
])
*
dx
[
2
];
(
pi
->
v
[
2
]
-
pj
->
v
[
2
])
*
dx
[
2
]
+
a2_Hubble
*
r2
;
/* Are the particles moving towards each others ? */
const
float
omega_ij
=
min
(
dvdr
,
0
.
f
);
...
...
src/hydro/Minimal/hydro_io.h
View file @
83c92087
...
...
@@ -71,12 +71,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
void
convert_S
(
const
struct
engine
*
e
,
const
struct
part
*
p
,
float
*
ret
)
{
ret
[
0
]
=
hydro_get_entropy
(
p
);
ret
[
0
]
=
hydro_get_
comoving_
entropy
(
p
);
}
void
convert_P
(
const
struct
engine
*
e
,
const
struct
part
*
p
,
float
*
ret
)
{
ret
[
0
]
=
hydro_get_pressure
(
p
);
ret
[
0
]
=
hydro_get_
comoving_
pressure
(
p
);
}
void
convert_part_pos
(
const
struct
engine
*
e
,
const
struct
part
*
p
,
...
...
src/hydro/Minimal/hydro_part.h
View file @
83c92087
...
...
@@ -53,6 +53,9 @@ struct xpart {
/*! Velocity at the last full step. */
float
v_full
[
3
];
/*! Gravitational acceleration at the last full step. */
float
a_grav
[
3
];
/*! Internal energy at the last full step. */
float
u_full
;
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment