Commit d2e3ba8a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added cosmological terms to the Pressure-Entropy SPH scheme.

parent dd882b93
......@@ -99,7 +99,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
}
/**
* @brief Returns the entropy of a particle
* @brief Returns the physical entropy of a particle
*
* @param p The particle of interest.
* @param cosmo The cosmological model.
......
......@@ -33,6 +33,7 @@
#include "adiabatic_index.h"
#include "approx_math.h"
#include "cosmology.h"
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
......@@ -41,60 +42,123 @@
#include "minmax.h"
/**
* @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) {
return gas_internal_energy_from_entropy(p->rho_bar, p->entropy);
}
/**
* @brief Returns the pressure 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_physical_internal_energy(const struct part *restrict p,
const struct cosmology *cosmo) {
return gas_internal_energy_from_entropy(p->rho_bar * cosmo->a3_inv,
p->entropy);
}
/**
* @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 gas_pressure_from_entropy(p->rho_bar, p->entropy);
}
/**
* @brief Returns the entropy of a particle
* @brief Returns the physical pressure of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
__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_bar * cosmo->a3_inv, p->entropy);
}
/**
* @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) {
return p->entropy;
}
/**
* @brief Returns the sound speed of a particle
* @brief Returns the physical entropy 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 p->entropy;
}
/**
* @brief Returns the comoving 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) {
return p->force.soundspeed;
}
/**
* @brief Returns the physical 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_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_density(
__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
const struct part *restrict p) {
return p->rho;
}
/**
* @brief Returns the physical density of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
const struct part *restrict p, const struct cosmology *cosmo) {
return p->rho * cosmo->a3_inv;
}
/**
* @brief Returns the mass of a particle
*
......@@ -111,16 +175,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;
}
/**
......@@ -153,19 +221,21 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
/**
* @brief Computes the hydro time-step of a given particle
*
* @param p Pointer to the particle data
* @param xp Pointer to the extended particle data
*
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
* @param hydro_properties The constants used in the scheme.
* @param cosmo The current 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;
const float CFL = 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 * cosmo->a * p->h /
(cosmo->a_factor_sound_speed * p->force.v_sig);
return dt_cfl;
}
......@@ -212,9 +282,10 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* and add the self-contribution term.
*
* @param p The particle to act upon
* @param cosmo The current 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;
......@@ -240,15 +311,16 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.wcount_dh *= h_inv_dim_plus_one;
const float rho_inv = 1.f / p->rho;
const float a_inv2 = cosmo->a2_inv;
const float entropy_minus_one_over_gamma = 1.f / p->entropy_one_over_gamma;
/* Final operation on the weighted density */
p->rho_bar *= entropy_minus_one_over_gamma;
/* Finish calculation of the velocity curl components */
p->density.rot_v[0] *= h_inv_dim_plus_one * rho_inv;
p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
p->density.rot_v[2] *= h_inv_dim_plus_one * rho_inv;
p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
/* Finish calculation of the velocity divergence */
p->density.div_v *= h_inv_dim_plus_one * rho_inv;
......@@ -259,9 +331,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 current 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;
......@@ -288,11 +362,13 @@ __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) {
const float fac_mu = 1.f; /* Will change with cosmological integration */
const float fac_mu = cosmo->a_factor_mu;
/* Compute the norm of the curl */
const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
......@@ -310,7 +386,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the Balsara switch */
const float balsara =
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
/* Divide the pressure by the density squared to get the SPH term */
const float rho_bar_inv = 1.f / p->rho_bar;
......@@ -380,15 +456,17 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
*
* @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.2f)
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
else
......@@ -405,7 +483,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
}
/* Predict the entropy */
p->entropy += p->entropy_dt * dt;
p->entropy += p->entropy_dt * dt_therm;
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
......@@ -429,14 +507,15 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* Multiplies the forces and accelerationsby the appropiate constants
*
* @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;
p->entropy_dt =
0.5f * gas_entropy_from_internal_energy(p->rho_bar, p->entropy_dt);
p->entropy_dt = 0.5f * cosmo->a2_inv *
gas_entropy_from_internal_energy(p->rho_bar, p->entropy_dt);
}
/**
......@@ -444,17 +523,16 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
*
* @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 half_dt The half time-step for this kick
* @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 entropy (temperature) by more than a factor of 2*/
if (dt > 0. && p->entropy_dt * dt < -0.5f * xp->entropy_full) {
p->entropy_dt = -0.5f * xp->entropy_full / dt;
if (dt_therm > 0. && p->entropy_dt * dt_therm < -0.5f * xp->entropy_full) {
p->entropy_dt = -0.5f * xp->entropy_full / dt_therm;
}
xp->entropy_full += p->entropy_dt * dt;
xp->entropy_full += p->entropy_dt * dt_therm;
/* Compute the pressure */
const float pressure =
......
......@@ -40,10 +40,10 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
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->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
p->rho, p->rho_bar, hydro_get_pressure(p), p->density.pressure_dh,
p->force.P_over_rho2, p->entropy, p->entropy_one_over_gamma,
p->entropy_dt, p->force.soundspeed, p->force.v_sig, p->force.h_dt,
p->time_bin);
p->rho, p->rho_bar, hydro_get_comoving_pressure(p),
p->density.pressure_dh, p->force.P_over_rho2, p->entropy,
p->entropy_one_over_gamma, p->entropy_dt, p->force.soundspeed,
p->force.v_sig, p->force.h_dt, p->time_bin);
}
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_DEBUG_H */
......@@ -21,7 +21,7 @@
/**
* @file PressureEntropy/hydro_iact.h
* @brief Pressure-Entropy implementation of SPH (Neighbour loop equations)
* @brief Pressure-Entropy implementation of SPH (Particle interactions)
*
* The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension.
......@@ -31,10 +31,20 @@
*/
/**
* @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, wi_dx;
float wj, wj_dx;
......@@ -111,10 +121,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;
float dv[3], curlvr[3];
......@@ -164,14 +184,26 @@ __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) {
float wi, wj, wi_dx, wj_dx;
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.0f / r;
......@@ -215,7 +247,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;
/* Balsara term */
const float balsara_i = pi->force.balsara;
......@@ -265,14 +297,26 @@ __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) {
float wi, wj, wi_dx, wj_dx;
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.0f / r;
......@@ -316,7 +360,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;
/* Balsara term */
const float balsara_i = pi->force.balsara;
......
......@@ -69,12 +69,12 @@ 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);
}
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,
......
......@@ -45,6 +45,9 @@ struct xpart {
/*! Velocity at the last full step. */
float v_full[3];
/*! Gravitational acceleration at the last full step. */
float a_grav[3];
/*! Entropy at the last full step. */
float entropy_full;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment