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

Added more conversion functions from comoving to physical for the thermodynamic quantities.

parent 0b4ee25b
Branches
Tags
1 merge request!509Cosmological time integration
...@@ -48,6 +48,7 @@ tests/swift_dopair_standard.dat ...@@ -48,6 +48,7 @@ tests/swift_dopair_standard.dat
tests/brute_force_perturbed.dat tests/brute_force_perturbed.dat
tests/swift_dopair_perturbed.dat tests/swift_dopair_perturbed.dat
tests/test27cells tests/test27cells
tests/test27cells_subset
tests/testPeriodicBC tests/testPeriodicBC
tests/test125cells tests/test125cells
tests/brute_force_27_standard.dat tests/brute_force_27_standard.dat
......
...@@ -137,7 +137,7 @@ void cosmology_update(struct cosmology *c, integertime_t ti_current) { ...@@ -137,7 +137,7 @@ void cosmology_update(struct cosmology *c, integertime_t ti_current) {
c->a = a; c->a = a;
c->a_inv = a_inv; c->a_inv = a_inv;
c->a3_inv = a_inv * a_inv * a_inv; c->a3_inv = a_inv * a_inv * a_inv;
c->a_factor_sig_vel = c->a_factor_sound_speed =
pow(a, -1.5 * hydro_gamma_minus_one); /* a^{3*(1-gamma)/2} */ pow(a, -1.5 * hydro_gamma_minus_one); /* a^{3*(1-gamma)/2} */
c->a_factor_grav_accel = a_inv * a_inv; /* 1 / a^2 */ c->a_factor_grav_accel = a_inv * a_inv; /* 1 / a^2 */
c->a_factor_hydro_accel = c->a_factor_hydro_accel =
...@@ -458,7 +458,7 @@ void cosmology_init_no_cosmo(const struct swift_params *params, ...@@ -458,7 +458,7 @@ void cosmology_init_no_cosmo(const struct swift_params *params,
c->a = 1.; c->a = 1.;
c->a_inv = 1.; c->a_inv = 1.;
c->a3_inv = 1.; c->a3_inv = 1.;
c->a_factor_sig_vel = 1.; c->a_factor_sound_speed = 1.;
c->a_factor_hydro_accel = 1.; c->a_factor_hydro_accel = 1.;
c->a_factor_grav_accel = 1.; c->a_factor_grav_accel = 1.;
......
...@@ -41,8 +41,8 @@ struct cosmology { ...@@ -41,8 +41,8 @@ struct cosmology {
/*! Inverse cube of the current expansion factor of the Universe */ /*! Inverse cube of the current expansion factor of the Universe */
double a3_inv; double a3_inv;
/*! Power of the scale-factor used for signal-velocity conversion */ /*! Power of the scale-factor used for sound-speed conversion */
double a_factor_sig_vel; double a_factor_sound_speed;
/*! Power of the scale-factor used for gravity accelerations */ /*! Power of the scale-factor used for gravity accelerations */
double a_factor_grav_accel; double a_factor_grav_accel;
......
...@@ -42,7 +42,7 @@ ...@@ -42,7 +42,7 @@
#include "minmax.h" #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 * @param p The particle of interest
*/ */
...@@ -53,7 +53,7 @@ hydro_get_comoving_internal_energy(const struct part *restrict p) { ...@@ -53,7 +53,7 @@ hydro_get_comoving_internal_energy(const struct part *restrict p) {
} }
/** /**
* @brief Returns the internal energy of a particle * @brief Returns the physical internal energy of a particle
* *
* @param p The particle of interest. * @param p The particle of interest.
* @param cosmo The cosmological model. * @param cosmo The cosmological model.
...@@ -66,16 +66,27 @@ hydro_get_physical_internal_energy(const struct part *restrict p, ...@@ -66,16 +66,27 @@ hydro_get_physical_internal_energy(const struct part *restrict p,
} }
/** /**
* @brief Returns the pressure of a particle * @brief Returns the comoving pressure of a particle
* *
* @param p The particle of interest * @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) { const struct part *restrict p) {
return gas_pressure_from_entropy(p->rho, p->entropy); return gas_pressure_from_entropy(p->rho, p->entropy);
} }
/**
* @brief Returns the physical pressure of a particle
*
* @param p The particle of interest
*/
__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 * cosmo->a3_inv, p->entropy);
}
/** /**
* @brief Returns the comoving entropy of a particle * @brief Returns the comoving entropy of a particle
* *
...@@ -100,16 +111,29 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_entropy( ...@@ -100,16 +111,29 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
} }
/** /**
* @brief Returns the sound speed of a particle * @brief Returns the comoving sound speed of a particle
* *
* @param p The particle of interest * @param p The particle of interest
*/ */
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed( __attribute__((always_inline)) INLINE static float
const struct part *restrict p) { hydro_get_comoving_soundspeed(const struct part *restrict p) {
return p->force.soundspeed; return p->force.soundspeed;
} }
/**
* @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 * @brief Returns the comoving density of a particle
* *
...@@ -209,7 +233,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep( ...@@ -209,7 +233,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
/* CFL condition */ /* CFL condition */
const float dt_cfl = 2.f * kernel_gamma * CFL * cosmo->a * p->h / const float dt_cfl = 2.f * kernel_gamma * CFL * cosmo->a * p->h /
(cosmo->a_factor_sig_vel * p->force.v_sig); (cosmo->a_factor_sound_speed * p->force.v_sig);
return dt_cfl; return dt_cfl;
} }
......
...@@ -31,7 +31,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle( ...@@ -31,7 +31,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, hydro_get_pressure(p), p->force.P_over_rho2, p->entropy, p->rho, hydro_get_comoving_pressure(p), p->force.P_over_rho2, p->entropy,
p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0], p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
p->density.rot_v[1], p->density.rot_v[2], p->force.balsara, p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
p->force.v_sig, p->force.h_dt, p->time_bin); p->force.v_sig, p->force.h_dt, p->time_bin);
......
...@@ -62,7 +62,7 @@ void convert_u(const struct engine* e, const struct part* p, float* ret) { ...@@ -62,7 +62,7 @@ void convert_u(const struct engine* e, const struct part* p, float* ret) {
void convert_P(const struct engine* e, const struct part* p, float* ret) { 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, void convert_part_pos(const struct engine* e, const struct part* p,
......
/******************************************************************************* /*******************************************************************************
* This file is part of SWIFT. * This file is part of SWIFT.
* Copyright (C) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk). * Copyright (C) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
...@@ -411,8 +410,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell, ...@@ -411,8 +410,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
#endif #endif
hydro_get_comoving_entropy(&main_cell->parts[pid]), hydro_get_comoving_entropy(&main_cell->parts[pid]),
hydro_get_comoving_internal_energy(&main_cell->parts[pid]), hydro_get_comoving_internal_energy(&main_cell->parts[pid]),
hydro_get_pressure(&main_cell->parts[pid]), hydro_get_comoving_pressure(&main_cell->parts[pid]),
hydro_get_soundspeed(&main_cell->parts[pid]), hydro_get_comoving_soundspeed(&main_cell->parts[pid]),
main_cell->parts[pid].a_hydro[0], main_cell->parts[pid].a_hydro[1], main_cell->parts[pid].a_hydro[0], main_cell->parts[pid].a_hydro[1],
main_cell->parts[pid].a_hydro[2], main_cell->parts[pid].force.h_dt, main_cell->parts[pid].a_hydro[2], main_cell->parts[pid].force.h_dt,
#if defined(GADGET2_SPH) #if defined(GADGET2_SPH)
......
...@@ -137,8 +137,8 @@ void dump_indv_particle_fields(char *fileName, struct part *p) { ...@@ -137,8 +137,8 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
p->density.div_v, p->density.div_v,
#endif #endif
hydro_get_comoving_entropy(p), hydro_get_comoving_internal_energy(p), hydro_get_comoving_entropy(p), hydro_get_comoving_internal_energy(p),
hydro_get_pressure(p), hydro_get_soundspeed(p), p->a_hydro[0], hydro_get_comoving_pressure(p), hydro_get_comoving_soundspeed(p),
p->a_hydro[1], p->a_hydro[2], p->force.h_dt, p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], p->force.h_dt,
#if defined(GADGET2_SPH) #if defined(GADGET2_SPH)
p->force.v_sig, p->entropy_dt, 0.f p->force.v_sig, p->entropy_dt, 0.f
#elif defined(DEFAULT_SPH) #elif defined(DEFAULT_SPH)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment