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

Cooling also fine with Pressure-Entropy

parent be7b197a
...@@ -76,7 +76,6 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy( ...@@ -76,7 +76,6 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
* @brief Returns the sound speed of a particle * @brief Returns the sound speed of a particle
* *
* @param p The particle of interest * @param p The particle of interest
* @param dt Time since the last kick
*/ */
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed( __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p) { const struct part *restrict p) {
......
...@@ -43,50 +43,42 @@ ...@@ -43,50 +43,42 @@
* @brief Returns the internal energy of a particle * @brief Returns the internal energy of a particle
* *
* @param p The particle of interest * @param p The particle of interest
* @param dt Time since the last kick
*/ */
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy( __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part *restrict p, float dt) { const struct part *restrict p) {
const float entropy = p->entropy + p->entropy_dt * dt;
return gas_internal_energy_from_entropy(p->rho_bar, entropy); return gas_internal_energy_from_entropy(p->rho_bar, p->entropy);
} }
/** /**
* @brief Returns the pressure of a particle * @brief Returns the pressure of a particle
* *
* @param p The particle of interest * @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_pressure(
const struct part *restrict p, float dt) { const struct part *restrict p) {
const float entropy = p->entropy + p->entropy_dt * dt;
return gas_pressure_from_entropy(p->rho_bar, entropy); return gas_pressure_from_entropy(p->rho_bar, p->entropy);
} }
/** /**
* @brief Returns the entropy of a particle * @brief Returns the entropy of a particle
* *
* @param p The particle of interest * @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_entropy(
const struct part *restrict p, float dt) { const struct part *restrict p) {
return p->entropy + p->entropy_dt * dt; return p->entropy;
} }
/** /**
* @brief Returns the sound speed of a particle * @brief Returns the sound speed of a particle
* *
* @param p The particle of interest * @param p The particle of interest
* @param dt Time since the last kick
*/ */
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed( __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) { const struct part *restrict p) {
return p->force.soundspeed; return p->force.soundspeed;
} }
...@@ -114,72 +106,30 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass( ...@@ -114,72 +106,30 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
} }
/** /**
* @brief Modifies the thermal state of a particle to the imposed internal * @brief Returns the time derivative of internal energy of a particle
* energy
* *
* This overwrites the current state of the particle but does *not* change its * We assume a constant density.
* time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
* updated.
* *
* @param p The particle * @param p The particle of interest
* @param u The new internal energy
*/ */
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy( __attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
struct part *restrict p, float u) { const struct part *restrict p) {
p->entropy = gas_entropy_from_internal_energy(p->rho_bar, u);
p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
/* Compute the sound speed from the pressure*/
const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
/* Update the signal velocity */
const float v_sig_old = p->force.v_sig;
const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
const float v_sig = max(v_sig_old, v_sig_new);
const float rho_bar_inv = 1.f / p->rho_bar;
p->force.soundspeed = soundspeed; return gas_internal_energy_from_entropy(p->rho_bar, p->entropy_dt);
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
p->force.v_sig = v_sig;
} }
/** /**
* @brief Modifies the thermal state of a particle to the imposed entropy * @brief Returns the time derivative of internal energy of a particle
* *
* This overwrites the current state of the particle but does *not* change its * We assume a constant density.
* time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
* updated.
* *
* @param p The particle * @param p The particle of interest.
* @param S The new entropy * @param du_dt The new time derivative of the internal energy.
*/ */
__attribute__((always_inline)) INLINE static void hydro_set_entropy( __attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
struct part *restrict p, float S) { struct part *restrict p, float du_dt) {
p->entropy = S; p->entropy_dt = gas_entropy_from_internal_energy(p->rho_bar, du_dt);
p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
/* Compute the sound speed from the pressure*/
const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
/* Update the signal velocity */
const float v_sig_old = p->force.v_sig;
const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
const float v_sig = max(v_sig_old, v_sig_new);
const float rho_bar_inv = 1.f / p->rho_bar;
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
p->force.v_sig = v_sig;
} }
/** /**
......
...@@ -40,7 +40,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle( ...@@ -40,7 +40,7 @@ __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], 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], 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->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
p->rho, p->rho_bar, hydro_get_pressure(p, 0.), p->density.pressure_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->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->entropy_dt, p->force.soundspeed, p->force.v_sig, p->force.h_dt,
p->time_bin); p->time_bin);
......
...@@ -69,12 +69,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list, ...@@ -69,12 +69,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
float convert_u(struct engine* e, struct part* p) { float convert_u(struct engine* e, struct part* p) {
return hydro_get_internal_energy(p, 0); return hydro_get_internal_energy(p);
} }
float convert_P(struct engine* e, struct part* p) { float convert_P(struct engine* e, struct part* p) {
return hydro_get_pressure(p, 0); return hydro_get_pressure(p);
} }
/** /**
......
...@@ -32,7 +32,7 @@ typedef long long integertime_t; ...@@ -32,7 +32,7 @@ typedef long long integertime_t;
typedef char timebin_t; typedef char timebin_t;
/*! The number of time bins */ /*! The number of time bins */
#define num_time_bins 50 #define num_time_bins 56
/*! The maximal number of timesteps in a simulation */ /*! The maximal number of timesteps in a simulation */
#define max_nr_timesteps (1LL << (num_time_bins + 1)) #define max_nr_timesteps (1LL << (num_time_bins + 1))
......
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