Skip to content
Snippets Groups Projects
Commit cfbd74f8 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added cosmological terms to the default hydro scheme.

parent 31701feb
No related branches found
No related tags found
1 merge request!509Cosmological time integration
Hide whitespace changes
Inline Side-by-side
src/hydro/Default/hydro.h
+ 130
− 48
View file @ cfbd74f8
......@@ -21,71 +21,134 @@
#include "adiabatic_index.h"
#include "approx_math.h"
#include "cosmology.h"
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
#include <float.h>
/**
* @brief Returns the internal energy of a particle
* @brief Returns the comoving internal energy of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__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
*
* @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 cosmo->a_factor_internal_energy * p->u;
}
/**
* @brief Returns the comoving pressure of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__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
*
* @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
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__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
*
* @param p The particle of interest
* @param dt Time since the last kick
* @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 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
*
......@@ -102,16 +165,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;
}
/**
......@@ -146,17 +213,19 @@ __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 constants used in the scheme
* @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;
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);
/* Limit change in u */
const float dt_u_change =
......@@ -187,6 +256,7 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p, const struct hydro_space *hs) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
......@@ -204,9 +274,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;
......@@ -226,14 +297,15 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
const float irho = 1.f / p->rho;
const float a_inv2 = cosmo->a2_inv;
/* Finish calculation of the velocity curl components */
p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * irho;
p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * irho;
p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * irho;
/* Finish calculation of the velocity divergence */
p->density.div_v *= h_inv_dim_plus_one * irho;
p->density.div_v *= h_inv_dim_plus_one * a_inv2 * irho;
}
/**
......@@ -241,9 +313,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;
......@@ -268,10 +342,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 time The current time
* @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 = cosmo->a_factor_mu;
/* Some smoothing length multiples. */
const float h = p->h;
......@@ -296,7 +373,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.P_over_rho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
/* Balsara switch */
p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * h_inv);
p->force.balsara =
normDiv_v / (normDiv_v + normRot_v + 0.0001f * fac_mu * fc * h_inv);
/* Viscosity parameter decay time */
/* const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
......@@ -349,6 +427,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->v[0] = xp->v_full[0];
p->v[1] = xp->v_full[1];
p->v[2] = xp->v_full[2];
p->u = xp->u_full;
}
/**
......@@ -356,19 +436,18 @@ __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 t0 The time at the start of the drift
* @param t1 The time at the end of the drift
* @param timeBase The minimal time-step size
* @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, struct xpart *restrict xp, float dt) {
struct part *restrict p, struct xpart *restrict xp, float dt_drift,
float dt_therm) {
float u, w;
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
......@@ -382,7 +461,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf(w2);
/* Predict internal energy */
w = p->force.u_dt / p->u * dt;
w = p->force.u_dt / p->u * dt_therm;
if (fabsf(w) < 0.2f)
u = p->u *= approx_expf(w);
else
......@@ -398,9 +477,10 @@ __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;
}
......@@ -409,11 +489,10 @@ __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) {}
/**
* @brief Converts hydro quantity of a particle at the start of a run
......@@ -441,6 +520,9 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
xp->a_grav[0] = 0.f;
xp->a_grav[1] = 0.f;
xp->a_grav[2] = 0.f;
xp->u_full = p->u;
hydro_reset_acceleration(p);
......
src/hydro/Default/hydro_iact.h
+ 64
− 16
View file @ cfbd74f8
......@@ -26,22 +26,32 @@
* @brief SPH interaction functions following the Gadget-2 version of SPH.
*
* The interactions computed here are the ones presented in the Gadget-2 paper
*and use the same
* and use the same
* numerical coefficients as the Gadget-2 code. When used with the Spline-3
*kernel, the results
* kernel, the results
* should be equivalent to the ones obtained with Gadget-2 up to the rounding
*errors and interactions
* errors and interactions
* missed by the Gadget-2 tree-code neighbours search.
*
* The code uses internal energy instead of entropy as a thermodynamical
*variable.
* variable.
*/
/**
* @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 r = sqrtf(r2), ri = 1.0f / r;
float xi, xj;
......@@ -97,10 +107,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 r, ri;
float xi;
......@@ -145,10 +165,20 @@ __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 r = sqrtf(r2), ri = 1.0f / r;
float xi, xj;
......@@ -160,6 +190,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
float f;
int k;
/* Cosmological factors entering the EoMs */
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
/* Get some values in local variables. */
mi = pi->mass;
mj = pj->mass;
......@@ -184,12 +218,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute dv dot r. */
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;
dvdr *= ri;
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
omega_ij = min(dvdr, 0.f);
omega_ij = min(fac_mu * dvdr, 0.f);
/* Compute signal velocity */
v_sig = pi->force.soundspeed + pj->force.soundspeed - 2.0f * omega_ij;
......@@ -240,10 +274,20 @@ __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 r = sqrtf(r2), ri = 1.0f / r;
float xi, xj;
......@@ -255,6 +299,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float f;
int k;
/* Cosmological factors entering the EoMs */
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
/* Get some values in local variables. */
// mi = pi->mass;
mj = pj->mass;
......@@ -279,12 +327,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute dv dot r. */
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;
dvdr *= ri;
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
omega_ij = min(dvdr, 0.f);
omega_ij = min(fac_mu * dvdr, 0.f);
/* Compute signal velocity */
v_sig = pi->force.soundspeed + pj->force.soundspeed - 2.0f * omega_ij;
......
src/hydro/Default/hydro_part.h
+ 3
− 0
View file @ cfbd74f8
......@@ -34,6 +34,9 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
/* Gravitational acceleration at the last full step. */
float a_grav[3];
/* Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment