Commit 6cb8470b authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added some cosmological terms to the gizmo scheme.

parent cfbd74f8
......@@ -22,6 +22,7 @@
#include "adiabatic_index.h"
#include "approx_math.h"
#include "cosmology.h"
#include "equation_of_state.h"
#include "hydro_gradients.h"
#include "hydro_properties.h"
......@@ -41,10 +42,12 @@
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
* @param hydro_properties Pointer to the hydro 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;
......@@ -63,8 +66,11 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
vmax = max(vmax, p->timestepvars.vmax);
const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
hydro_dimension_inv);
// MATTHIEU: Bert is this correct? Do we need more cosmology terms here?
const float psize =
cosmo->a * powf(p->geometry.volume / hydro_dimension_unit_sphere,
hydro_dimension_inv);
float dt = FLT_MAX;
if (vmax > 0.) {
dt = psize / vmax;
......@@ -223,9 +229,10 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* passive particles.
*
* @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;
......@@ -318,6 +325,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
}
#endif
// MATTHIEU: Bert is this correct? Do we need cosmology terms here?
float momentum[3];
momentum[0] = p->conserved.momentum[0];
momentum[1] = p->conserved.momentum[1];
......@@ -377,9 +385,11 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
*
* @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;
......@@ -420,13 +430,17 @@ __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 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) {
/* Initialize time step criterion variables */
p->timestepvars.vmax = 0.;
// MATTHIEU: Bert is this correct? Do we need cosmology terms here?
/* Set the actual velocity of the particle */
hydro_velocities_prepare_force(p, xp);
}
......@@ -510,13 +524,11 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
*
* @param p Particle to act upon.
* @param xp The extended particle data to act upon.
* @param dt The drift time-step.
* @param t0 Integer start time of the drift interval.
* @param t1 Integer end time of the drift interval.
* @param timeBase Conversion factor between integer and physical time.
* @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* p, struct xpart* xp, float dt) {
struct part* p, struct xpart* xp, float dt_drift, float dt_therm) {
#ifdef GIZMO_LLOYD_ITERATION
return;
......@@ -525,7 +537,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const float h_inv = 1.0f / 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;
float h_corr;
if (fabsf(w1) < 0.2f)
h_corr = approx_expf(w1); /* 4th order expansion of exp(w) */
......@@ -540,19 +552,20 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* drift the primitive variables based on the old fluxes */
if (p->geometry.volume > 0.) {
p->primitives.rho += p->conserved.flux.mass * dt / p->geometry.volume;
p->primitives.rho += p->conserved.flux.mass * dt_drift / p->geometry.volume;
}
if (p->conserved.mass > 0.) {
p->primitives.v[0] +=
p->conserved.flux.momentum[0] * dt / p->conserved.mass;
p->conserved.flux.momentum[0] * dt_drift / p->conserved.mass;
p->primitives.v[1] +=
p->conserved.flux.momentum[1] * dt / p->conserved.mass;
p->conserved.flux.momentum[1] * dt_drift / p->conserved.mass;
p->primitives.v[2] +=
p->conserved.flux.momentum[2] * dt / p->conserved.mass;
p->conserved.flux.momentum[2] * dt_drift / p->conserved.mass;
// MATTHIEU: Bert is this correct?
#if !defined(EOS_ISOTHERMAL_GAS)
const float u = p->conserved.energy + p->conserved.flux.energy * dt;
const float u = p->conserved.energy + p->conserved.flux.energy * dt_therm;
p->primitives.P =
hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
#endif
......@@ -576,12 +589,14 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* should only happen at the start of the simulation.
*
* @param p Particle to act upon.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {
struct part* p, const struct cosmology* cosmo) {
/* set the variables that are used to drift the primitive variables */
// MATTHIEU: Bert is this correct? Do we need cosmology terms here?
hydro_velocities_end_force(p);
}
......@@ -679,12 +694,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
}
/**
* @brief Returns the internal energy of a particle
* @brief Returns the comoving internal energy of a particle
*
* @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) {
if (p->primitives.rho > 0.)
return gas_internal_energy_from_pressure(p->primitives.rho,
......@@ -694,11 +709,25 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
}
/**
* @brief Returns the entropy of a particle
* @brief Returns the physical internal energy of a particle
*
* @param p The particle of interest.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
__attribute__((always_inline)) INLINE static float
hydro_get_physical_internal_energy(const struct part* restrict p,
const struct cosmology* cosmo) {
return cosmo->a_factor_internal_energy *
hydro_get_comoving_internal_energy(p);
}
/**
* @brief Returns the comoving entropy of a particle
*
* @param p The particle of interest.
*/
__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
const struct part* restrict p) {
if (p->primitives.rho > 0.) {
......@@ -708,13 +737,27 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
}
}
/**
* @brief Returns the physical internal energy of a particle
*
* @param p The particle of interest.
* @param cosmo The cosmological model.
*/
__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 hydro_get_comoving_entropy(p);
}
/**
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest.
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part* restrict p) {
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_soundspeed(const struct part* restrict p) {
if (p->primitives.rho > 0.)
return gas_soundspeed_from_pressure(p->primitives.rho, p->primitives.P);
......@@ -723,16 +766,41 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
}
/**
* @brief Returns the pressure 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_physical_soundspeed(const struct part* restrict p,
const struct cosmology* cosmo) {
return cosmo->a_factor_sound_speed * hydro_get_comoving_soundspeed(p);
}
/**
* @brief Returns the comoving pressure of a particle
*
* @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 p->primitives.P;
}
/**
* @brief Returns the comoving pressure of a particle
*
* @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 * p->primitives.P;
}
/**
* @brief Returns the mass of a particle
*
......@@ -749,38 +817,56 @@ __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]) {
const struct part* restrict p, const struct xpart* xp, float dt_kick_hydro,
float dt_kick_grav, float v[3]) {
if (p->conserved.mass > 0.) {
v[0] = p->primitives.v[0] +
p->conserved.flux.momentum[0] * dt / p->conserved.mass;
p->conserved.flux.momentum[0] * dt_kick_hydro / p->conserved.mass;
v[1] = p->primitives.v[1] +
p->conserved.flux.momentum[1] * dt / p->conserved.mass;
p->conserved.flux.momentum[1] * dt_kick_hydro / p->conserved.mass;
v[2] = p->primitives.v[2] +
p->conserved.flux.momentum[2] * dt / p->conserved.mass;
p->conserved.flux.momentum[2] * dt_kick_hydro / p->conserved.mass;
} else {
v[0] = p->primitives.v[0];
v[1] = p->primitives.v[1];
v[2] = p->primitives.v[2];
}
// MATTHIEU: Bert is this correct?
v[0] += xp->a_grav[0] * dt_kick_grav;
v[1] += xp->a_grav[1] * dt_kick_grav;
v[2] += xp->a_grav[2] * dt_kick_grav;
}
/**
* @brief Returns the density of a particle
* @brief Returns the comoving density of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_density(
__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
const struct part* restrict p) {
return p->primitives.rho;
}
/**
* @brief Returns the physical 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->primitives.rho;
}
/**
* @brief Modifies the thermal state of a particle to the imposed internal
* energy
......
......@@ -45,7 +45,7 @@
* @param p Particle.
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_init(
struct part* p) {}
struct part *p) {}
/**
* @brief Gradient calculations done during the neighbour loop
......@@ -58,8 +58,8 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
* @param pj Particle j.
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
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) {}
/**
* @brief Gradient calculations done during the neighbour loop: non-symmetric
......@@ -73,8 +73,9 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
* @param pj Particle j.
*/
__attribute__((always_inline)) INLINE static void
hydro_gradients_nonsym_collect(float r2, float* dx, float hi, float hj,
struct part* pi, struct part* pj) {}
hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
struct part *restrict pi,
struct part *restrict pj) {}
/**
* @brief Finalize the gradient variables after all data have been collected
......@@ -82,7 +83,7 @@ hydro_gradients_nonsym_collect(float r2, float* dx, float hi, float hj,
* @param p Particle.
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
struct part* p) {}
struct part *p) {}
#endif
......@@ -91,8 +92,8 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
* gradients_none does nothing, since all gradients are zero -- are they?).
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_predict(
struct part* pi, struct part* pj, float hi, float hj, float* dx, float r,
float* xij_i, float* Wi, float* Wj) {
struct part* restrict pi, struct part* restrict pj, float hi, float hj,
const float* dx, float r, float* xij_i, float* Wi, float* Wj) {
float dWi[5], dWj[5];
float xij_j[3];
......
......@@ -62,7 +62,8 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
* @param pj Particle j.
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
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 r = sqrtf(r2);
float xi, xj;
......@@ -299,8 +300,9 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
* @param pj Particle j.
*/
__attribute__((always_inline)) INLINE static void
hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
struct part *pi, struct part *pj) {
hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
struct part *restrict pi,
struct part *restrict pj) {
float r = sqrtf(r2);
float xi;
......
......@@ -61,7 +61,8 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
* @param pj Particle j.
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
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 wi, wi_dx, xi, hi_inv;
float wj, wj_dx, xj, hj_inv;
......@@ -163,8 +164,9 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
* @param pj Particle j.
*/
__attribute__((always_inline)) INLINE static void
hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
struct part *pi, struct part *pj) {
hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
struct part *restrict pi,
struct part *restrict pj) {
float wi, wi_dx, xi, hi_inv;
float r = sqrtf(r2);
......
......@@ -37,15 +37,19 @@
* order accurate gradient calculations and for the calculation of the interface
* surface areas.
*
* @param r2 Squared distance between particle i and particle j.
* @param dx Distance vector between the particles (dx = pi->x - pj->x).
* @param hi Smoothing length of particle i.
* @param hj Smoothing length of particle j.
* @param r2 Comoving squared distance between particle i and particle j.
* @param dx Comoving distance vector between the particles (dx = pi->x -
* pj->x).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi Particle i.
* @param pj Particle j.
* @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 r = sqrtf(r2);
float xi, xj;
......@@ -99,15 +103,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
* order accurate gradient calculations and for the calculation of the interface
* surface areas.
*
* @param r2 Squared distance between particle i and particle j.
* @param dx Distance vector between the particles (dx = pi->x - pj->x).
* @param hi Smoothing length of particle i.
* @param hj Smoothing length of particle j.
* @param r2 Comoving squared distance between particle i and particle j.
* @param dx Comoving distance vector between the particles (dx = pi->x -
* pj->x).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi Particle i.
* @param pj Particle j.
* @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 r;
float xi;
......@@ -141,15 +149,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
* This method wraps around hydro_gradients_collect, which can be an empty
* method, in which case no gradients are used.
*
* @param r2 Squared distance between particle i and particle j.
* @param dx Distance vector between the particles (dx = pi->x - pj->x).
* @param hi Smoothing length of particle i.
* @param hj Smoothing length of particle j.
* @param r2 Comoving squared distance between particle i and particle j.
* @param dx Comoving distance vector between the particles (dx = pi->x -
* pj->x).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi Particle i.
* @param pj Particle j.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_gradient(
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) {
hydro_gradients_collect(r2, dx, hi, hj, pi, pj);
}
......@@ -161,15 +173,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
* This method wraps around hydro_gradients_nonsym_collect, which can be an
* empty method, in which case no gradients are used.
*
* @param r2 Squared distance between particle i and particle j.
* @param dx Distance vector between the particles (dx = pi->x - pj->x).
* @param hi Smoothing length of particle i.
* @param hj Smoothing length of particle j.
* @param r2 Comoving squared distance between particle i and particle j.
* @param dx Comoving distance vector between the particles (dx = pi->x -
* pj->x).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi Particle i.
* @param pj Particle j.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
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) {
hydro_gradients_nonsym_collect(r2, dx, hi, hj, pi, pj);
}
......@@ -192,16 +208,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
* This method also calculates the maximal velocity used to calculate the time
* step.
*
* @param r2 Squared distance between particle i and particle j.
* @param dx Distance vector between the particles (dx = pi->x - pj->x).
* @param hi Smoothing length of particle i.
* @param hj Smoothing length of particle j.
* @param r2 Comoving squared distance between particle i and particle j.
* @param dx Comoving distance vector between the particles (dx = pi->x -
* pj->x).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi Particle i.
* @param pj Particle j.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
int mode) {
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj, int mode, float a, float H) {
float r = sqrtf(r2);
float xi, xj;
......@@ -457,17 +476,21 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
*
* This method calls runner_iact_fluxes_common with mode 1.
*
* @param r2 Squared distance between particle i and particle j.
* @param dx Distance vector between the particles (dx = pi->x - pj->x).
* @param hi Smoothing length of particle i.
* @param hj Smoothing length of particle j.
* @param r2 Comoving squared distance between particle i and particle j.
* @param dx Comoving distance vector between the particles (dx = pi->x -
* pj->x).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi Particle i.
* @param pj Particle j.
* @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) {
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1);
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1, a, H);
}
/**
......@@ -476,17 +499,21 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
*
* This method calls runner_iact_fluxes_common with mode 0.
*
* @param r2 Squared distance between particle i and particle j.
* @param dx Distance vector between the particles (dx = pi->x - pj->x).
* @param hi Smoothing length of particle i.
* @param hj Smoothing length of particle j.
* @param r2 Comoving squared distance between particle i and particle j.
* @param dx Comoving distance vector between the particles (dx = pi->x -
* pj->x).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi Particle i.
* @param pj Particle j.
* @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,
struct part *restrict pj, float a, float H) {
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0);
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
}
#endif /* SWIFT_GIZMO_HYDRO_IACT_H */
......@@ -74,7 +74,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
*/
void convert_u(const struct engine* e, const struct part* p, float* ret) {
ret[0] = hydro_get_internal_energy(p);
ret[0] = hydro_get_comoving_internal_energy(p);
}