Commit 63ada53a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'shadowswift' into 'master'

Shadowswift repairs

See merge request !551
parents 26483c3f 93e2c2a4
......@@ -20,6 +20,7 @@
#include <float.h>
#include "adiabatic_index.h"
#include "approx_math.h"
#include "cosmology.h"
#include "equation_of_state.h"
#include "hydro_gradients.h"
#include "hydro_properties.h"
......@@ -35,7 +36,8 @@
*/
__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;
......@@ -158,7 +160,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
struct part* restrict p) {
struct part* restrict p, const struct cosmology* cosmo) {
float volume;
float m, momentum[3], energy;
......@@ -246,7 +248,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
* @param xp The extended particle data to act upon
*/
__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;
......@@ -273,7 +276,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
* @param xp The extended particle data to act upon.
*/
__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.0f;
......@@ -346,7 +350,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @param xp The extended particle data to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p, struct xpart* xp) {}
struct part* p, struct xpart* xp, const struct cosmology* cosmo) {}
/**
* @brief Extra operations to be done during the drift
......@@ -358,7 +362,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
* @param dt The drift time-step.
*/
__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) {}
/**
* @brief Set the particle acceleration after the flux loop.
......@@ -366,7 +370,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @param p Particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {}
struct part* p, const struct cosmology* cosmo) {}
/**
* @brief Extra operations done during the kick
......@@ -378,7 +382,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param dt Physical time step.
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt) {
struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {
float vcell[3];
......@@ -553,8 +558,13 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
* @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]) {
v[0] = p->v[0];
v[1] = p->v[1];
v[2] = p->v[2];
}
/**
* @brief Returns the density of a particle
......@@ -621,3 +631,171 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
p->primitives.P = gas_pressure_from_entropy(p->primitives.rho, S);
}
}
/**
* @brief Sets the mass of a particle
*
* @param p The particle of interest
* @param m The mass to set.
*/
__attribute__((always_inline)) INLINE static void hydro_set_mass(
struct part* restrict p, float m) {
p->conserved.mass = m;
}
/**
* @brief Overwrite the initial internal energy of a particle.
*
* Note that in the cases where the thermodynamic variable is not
* internal energy but gets converted later, we must overwrite that
* field. The conversion to the actual variable happens later after
* the initial fake time-step.
*
* @param p The #part to write to.
* @param u_init The new initial internal energy.
*/
__attribute__((always_inline)) INLINE static void
hydro_set_init_internal_energy(struct part* p, float u_init) {
p->conserved.energy = u_init * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic energy */
p->conserved.energy += 0.5f * p->conserved.mass *
(p->conserved.momentum[0] * p->primitives.v[0] +
p->conserved.momentum[1] * p->primitives.v[1] +
p->conserved.momentum[2] * p->primitives.v[2]);
#endif
p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
}
/**
* @brief Returns the comoving internal energy of a particle
*
* @param p The particle of interest.
*/
__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,
p->primitives.P);
else
return 0.;
}
/**
* @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.) {
return gas_entropy_from_pressure(p->primitives.rho, p->primitives.P);
} else {
return 0.;
}
}
/**
* @brief Returns the 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) {
if (p->primitives.rho > 0.)
return gas_soundspeed_from_pressure(p->primitives.rho, p->primitives.P);
else
return 0.;
}
/**
* @brief Returns the comoving pressure of a particle
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
const struct part* restrict p) {
return p->primitives.P;
}
/**
* @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->primitives.rho;
}
/**
* @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 *
hydro_get_comoving_internal_energy(p);
}
/**
* @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 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.
* @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 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;
}
......@@ -51,8 +51,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
......@@ -66,8 +66,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,
const struct part* restrict pj) {}
/**
* @brief Finalize the gradient variables after all data have been collected
......@@ -84,8 +85,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, float mindt) {
struct part* pi, struct part* pj, float hi, float hj, const float* dx,
float r, float* xij_i, float* Wi, float* Wj, float mindt) {
float dWi[5], dWj[5];
float xij_j[3];
......
......@@ -65,7 +65,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
* @param grad Current value of the gradient for the quantity (is updated).
*/
__attribute__((always_inline)) INLINE void hydro_gradients_single_quantity(
float qL, float qR, float *cLR, float *xLR, float rLR, float A,
float qL, float qR, float *cLR, const float *xLR, float rLR, float A,
float *grad) {
grad[0] += A * ((qR - qL) * cLR[0] / rLR - 0.5f * (qL + qR) * xLR[0] / rLR);
......@@ -84,7 +84,8 @@ __attribute__((always_inline)) INLINE void hydro_gradients_single_quantity(
* @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 *pi,
struct part *pj) {
float A, midpoint[3];
......@@ -148,8 +149,8 @@ __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 *pi, const struct part *pj) {
float A, midpoint[3];
......
......@@ -35,7 +35,8 @@
* @param pj Particle j.
*/
__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 mindx[3];
......@@ -59,7 +60,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
* @param pj Particle j.
*/
__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) {
voronoi_cell_interact(&pi->cell, dx, pj->id);
}
......@@ -78,7 +80,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
* @param pj Particle j.
*/
__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);
}
......@@ -98,7 +101,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
* @param pj Particle j.
*/
__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);
}
......@@ -129,8 +133,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
* @param pj Particle j.
*/
__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);
int k;
......@@ -322,9 +326,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
* @param pj Particle j.
*/
__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);
}
/**
......@@ -341,7 +346,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
* @param pj Particle j.
*/
__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);
}
......@@ -64,7 +64,8 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
* @param p Particle.
* @return Internal energy of the particle
*/
void convert_u(const struct engine* e, const struct part* p, float* ret) {
void convert_u(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
ret[0] = hydro_get_internal_energy(p);
}
......@@ -75,7 +76,8 @@ void convert_u(const struct engine* e, const struct part* p, float* ret) {
* @param p Particle.
* @return Entropic function of the particle
*/
void convert_A(const struct engine* e, const struct part* p, float* ret) {
void convert_A(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
ret[0] = hydro_get_entropy(p);
}
......@@ -86,7 +88,8 @@ void convert_A(const struct engine* e, const struct part* p, float* ret) {
* @param p Particle.
* @return Total energy of the particle
*/
void convert_Etot(const struct engine* e, const struct part* p, float* ret) {
void convert_Etot(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
#ifdef SHADOWFAX_TOTAL_ENERGY
return p->conserved.energy;
#else
......@@ -105,7 +108,7 @@ void convert_Etot(const struct engine* e, const struct part* p, float* ret) {
}
void convert_part_pos(const struct engine* e, const struct part* p,
double* ret) {
const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
......@@ -125,14 +128,15 @@ void convert_part_pos(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
struct io_props* list, int* num_fields) {
*num_fields = 13;
/* List what we want to write */
list[0] = io_make_output_field_convert_part(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
UNIT_CONV_LENGTH, parts, xparts,
convert_part_pos);
list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts,
primitives.v);
list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
......@@ -141,7 +145,7 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
parts, h);
list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
parts, convert_u);
parts, xparts, convert_u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
list[6] = io_make_output_field("Acceleration", FLOAT, 3,
......@@ -153,11 +157,11 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
list[9] = io_make_output_field("GradDensity", FLOAT, 3, UNIT_CONV_DENSITY,
parts, primitives.gradients.rho);
list[10] = io_make_output_field_convert_part(
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, convert_A);
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, xparts, convert_A);
list[11] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
parts, primitives.P);
list[12] = io_make_output_field_convert_part(
"TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, convert_Etot);
"TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, xparts, convert_Etot);
}
/**
......
......@@ -35,6 +35,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;
......
......@@ -50,7 +50,8 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
* @param r Distance between particle i and particle j.
*/
__attribute__((always_inline)) INLINE static void
hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
hydro_slope_limit_cell_collect(struct part* pi, const struct part* pj,
float r) {
/* basic slope limiter: collect the maximal and the minimal value for the
* primitive variables among the ngbs */
......
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