Commit 8fad2212 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Cleaned up 1D Voronoi algorithm some more. Added SHADOWFAX_TOTAL_ENERGY flag...

Cleaned up 1D Voronoi algorithm some more. Added SHADOWFAX_TOTAL_ENERGY flag and support for an isothermal EoS.
parent 249885fd
......@@ -60,6 +60,8 @@
/* Options to control SHADOWFAX_SPH */
/* This option disables cell movement */
//#define SHADOWFAX_FIX_CELLS
/* This option evolves the total energy instead of the thermal energy */
//#define SHADOWFAX_TOTAL_ENERGY
/* Source terms */
#define SOURCETERMS_NONE
......
......@@ -69,6 +69,21 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
return entropy * pow_gamma(density);
}
/**
* @brief Returns the entropy given density and pressure.
*
* Computes \f$A = \frac{P}{\rho^\gamma}\f$.
*
* @param density The density \f$\rho\f$.
* @param pressure The pressure \f$P\f$.
* @return The entropy \f$A\f$.
*/
__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
float density, float pressure) {
return pressure * pow_minus_gamma(density);
}
/**
* @brief Returns the sound speed given density and entropy
*
......@@ -111,6 +126,20 @@ gas_pressure_from_internal_energy(float density, float u) {
return hydro_gamma_minus_one * u * density;
}
/**
* @brief Returns the internal energy given density and pressure.
*
* Computes \f$u = \frac{1}{\gamma - 1}\frac{P}{\rho}\f$.
*
* @param density The density \f$\rho\f$.
* @param pressure The pressure \f$P\f$.
* @return The internal energy \f$u\f$.
*/
__attribute__((always_inline)) INLINE static float
gas_internal_energy_from_pressure(float density, float pressure) {
return hydro_one_over_gamma_minus_one * pressure / density;
}
/**
* @brief Returns the sound speed given density and internal energy
*
......@@ -172,6 +201,22 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
return hydro_gamma_minus_one * const_isothermal_internal_energy * density;
}
/**
* @brief Returns the entropy given density and pressure.
*
* Computes \f$A = \frac{P}{\rho^\gamma}\f$.
*
* @param density The density \f$\rho\f$.
* @param pressure The pressure \f$P\f$ (ignored).
* @return The entropy \f$A\f$.
*/
__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
float density, float pressure) {
return hydro_gamma_minus_one * const_isothermal_internal_energy *
pow_minus_gamma_minus_one(density);
}
/**
* @brief Returns the sound speed given density and entropy
*
......@@ -219,6 +264,20 @@ gas_pressure_from_internal_energy(float density, float u) {
return hydro_gamma_minus_one * const_isothermal_internal_energy * density;
}
/**
* @brief Returns the internal energy given density and pressure.
*
* Just returns the constant internal energy.
*
* @param density The density \f$\rho\f$ (ignored).
* @param pressure The pressure \f$P\f$ (ignored).
* @return The internal energy \f$u\f$ (which is constant).
*/
__attribute__((always_inline)) INLINE static float
gas_internal_energy_from_pressure(float density, float pressure) {
return const_isothermal_energy;
}
/**
* @brief Returns the sound speed given density and internal energy
*
......
......@@ -39,7 +39,13 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
float R = get_radius_dimension_sphere(p->cell.volume);
return CFL_condition * R / fabsf(p->timestepvars.vmax);
if (p->timestepvars.vmax == 0.) {
/* vmax can be zero in vacuum. We force the time step to become the maximal
time step in this case */
return FLT_MAX;
} else {
return CFL_condition * R / fabsf(p->timestepvars.vmax);
}
}
/**
......@@ -91,7 +97,17 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->conserved.momentum[1] = mass * p->primitives.v[1];
p->conserved.momentum[2] = mass * p->primitives.v[2];
#ifdef EOS_ISOTHERMAL_GAS
p->conserved.energy = mass * const_isothermal_internal_energy;
#else
p->conserved.energy *= mass;
#endif
#ifdef SHADOWFAX_TOTAL_ENERGY
p->conserved.energy += 0.5f * (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
#if defined(SHADOWFAX_FIX_CELLS)
p->v[0] = 0.;
......@@ -99,6 +115,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->v[2] = 0.;
#endif
/* set the initial velocity of the cells */
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
......@@ -117,7 +134,6 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount = 0.0f;
p->density.wcount_dh = 0.0f;
p->rho = 0.0f;
voronoi_cell_init(&p->cell, p->x);
......@@ -165,13 +181,22 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Iteration not succesful: we force h to become 1.1*hnew */
p->density.wcount = 0.0f;
p->density.wcount_dh = 1.0f / (1.1f * hnew - p->h);
return;
}
volume = p->cell.volume;
#ifdef SWIFT_DEBUG_CHECKS
/* the last condition checks for NaN: a NaN value always evaluates to false,
even when checked against itself */
if (volume == 0. || volume == INFINITY || volume != volume) {
error("Invalid value for cell volume (%g)!", volume);
}
#endif
/* compute primitive variables */
/* eqns (3)-(5) */
m = p->conserved.mass;
if (m) {
if (m > 0.) {
momentum[0] = p->conserved.momentum[0];
momentum[1] = p->conserved.momentum[1];
momentum[2] = p->conserved.momentum[2];
......@@ -179,9 +204,36 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->primitives.v[0] = momentum[0] / m;
p->primitives.v[1] = momentum[1] / m;
p->primitives.v[2] = momentum[2] / m;
energy = p->conserved.energy;
p->primitives.P = hydro_gamma_minus_one * energy / volume;
#ifdef SHADOWFAX_TOTAL_ENERGY
energy -= 0.5f * (momentum[0] * p->primitives.v[0] +
momentum[1] * p->primitives.v[1] +
momentum[2] * p->primitives.v[2]);
#endif
energy /= m;
p->primitives.P =
gas_pressure_from_internal_energy(p->primitives.rho, energy);
} else {
p->primitives.rho = 0.;
p->primitives.v[0] = 0.;
p->primitives.v[1] = 0.;
p->primitives.v[2] = 0.;
p->primitives.P = 0.;
}
#ifdef SWIFT_DEBUG_CHECKS
if (p->primitives.rho < 0.) {
error("Negative density!");
}
if (p->primitives.P < 0.) {
error("Negative pressure!");
}
#endif
}
/**
......@@ -197,15 +249,10 @@ __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 ti_current Current integer time.
* @param timeBase Conversion factor between integer time and physical time.
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp) {
/* Set the physical time step */
// p->force.dt = (p->ti_end - p->ti_begin) * timeBase;
/* Initialize time step criterion variables */
p->timestepvars.vmax = 0.0f;
......@@ -292,15 +339,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, float dt) {}
/**
* @brief Set the particle acceleration after the flux loop
*
* We use the new conserved variables to calculate the new velocity of the
* particle, and use that to derive the change of the velocity over the particle
* time step. We also add a correction to the velocity which steers the
* generator towards the centroid of the cell.
*
* If the particle time step is zero, we set the accelerations to zero. This
* should only happen at the start of the simulation.
* @brief Set the particle acceleration after the flux loop.
*
* @param p Particle to act upon.
*/
......@@ -318,6 +357,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt) {
float vcell[3];
/* Update the conserved variables. We do this here and not in the kick,
......@@ -328,6 +368,18 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.momentum[2] += p->conserved.flux.momentum[2];
p->conserved.energy += p->conserved.flux.energy;
#ifdef EOS_ISOTHERMAL_GAS
/* reset the thermal energy */
p->conserved.energy = p->conserved.mass * const_isothermal_internal_energy;
#ifdef SHADOWFAX_TOTAL_ENERGY
p->conserved.energy += 0.5f * (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
#endif
/* reset fluxes */
/* we can only do this here, since we need to keep the fluxes for inactive
particles */
......@@ -337,11 +389,17 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.flux.momentum[2] = 0.0f;
p->conserved.flux.energy = 0.0f;
/* We want the cell velocity to be as close as possible to the fluid
velocity */
vcell[0] = p->conserved.momentum[0] / p->conserved.mass;
vcell[1] = p->conserved.momentum[1] / p->conserved.mass;
vcell[2] = p->conserved.momentum[2] / p->conserved.mass;
if (p->conserved.mass > 0.) {
/* We want the cell velocity to be as close as possible to the fluid
velocity */
vcell[0] = p->conserved.momentum[0] / p->conserved.mass;
vcell[1] = p->conserved.momentum[1] / p->conserved.mass;
vcell[2] = p->conserved.momentum[2] / p->conserved.mass;
} else {
vcell[0] = 0.;
vcell[1] = 0.;
vcell[2] = 0.;
}
/* To prevent stupid things like cell crossovers or generators that move
outside their cell, we steer the motion of the cell somewhat */
......@@ -399,7 +457,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part* restrict p) {
return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
if (p->primitives.rho > 0.) {
return gas_internal_energy_from_pressure(p->primitives.rho,
p->primitives.P);
} else {
return 0.;
}
}
/**
......@@ -411,7 +474,11 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part* restrict p) {
return p->primitives.P / pow_gamma(p->primitives.rho);
if (p->primitives.rho > 0.) {
return gas_entropy_from_pressure(p->primitives.rho, p->primitives.P);
} else {
return 0.;
}
}
/**
......@@ -423,7 +490,11 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part* restrict p) {
return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
if (p->primitives.rho > 0.) {
return gas_soundspeed_from_pressure(p->primitives.rho, p->primitives.P);
} else {
return 0.;
}
}
/**
......@@ -473,7 +544,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_density(
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
struct part* restrict p, float u) {
p->conserved.energy = u * p->conserved.mass;
if (p->primitives.rho > 0.) {
p->conserved.energy = u * p->conserved.mass;
#ifdef SHADOWFAX_TOTAL_ENERGY
p->conserved.energy +=
0.5f * (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 = gas_pressure_from_internal_energy(p->primitives.rho, u);
}
}
/**
......@@ -488,6 +570,18 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
__attribute__((always_inline)) INLINE static void hydro_set_entropy(
struct part* restrict p, float S) {
p->conserved.energy = gas_internal_energy_from_entropy(p->primitives.rho, S) *
p->conserved.mass;
if (p->primitives.rho > 0.) {
p->conserved.energy =
gas_internal_energy_from_entropy(p->primitives.rho, S) *
p->conserved.mass;
#ifdef SHADOWFAX_TOTAL_ENERGY
p->conserved.energy +=
0.5f * (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 = gas_pressure_from_entropy(p->primitives.rho, S);
}
}
......@@ -45,8 +45,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"vmax=%.3e}, "
"density={"
"wcount_dh=%.3e, "
"wcount=%.3e}, "
"mass=%.3e\n",
"wcount=%.3e}",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
p->a_hydro[1], p->a_hydro[2], p->h, p->primitives.v[0],
p->primitives.v[1], p->primitives.v[2], p->primitives.rho,
......@@ -66,5 +65,5 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
p->primitives.limiter.maxr, p->conserved.momentum[0],
p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
p->conserved.energy, p->timestepvars.vmax, p->density.wcount_dh,
p->density.wcount, p->mass);
p->density.wcount);
}
......@@ -168,17 +168,21 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
dtj = pj->force.dt;
/* calculate the maximal signal velocity */
if (Wi[0] && Wj[0]) {
vmax =
sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
} else {
vmax = 0.0f;
vmax = 0.0f;
if (Wi[0] > 0.) {
vmax += gas_soundspeed_from_pressure(Wi[0], Wi[4]);
}
if (Wj[0] > 0.) {
vmax += gas_soundspeed_from_pressure(Wj[0], Wj[4]);
}
dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
(Wi[3] - Wj[3]) * dx[2];
if (dvdotdx > 0.) {
vmax -= dvdotdx / r;
}
pi->timestepvars.vmax = fmaxf(pi->timestepvars.vmax, vmax);
if (mode == 1) {
pj->timestepvars.vmax = fmaxf(pj->timestepvars.vmax, vmax);
......@@ -187,8 +191,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* The flux will be exchanged using the smallest time step of the two
* particles */
mindt = fminf(dti, dtj);
dti = mindt;
dtj = mindt;
/* compute the normal vector of the interface */
for (k = 0; k < 3; ++k) {
......@@ -260,19 +262,21 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* Update conserved variables */
/* eqn. (16) */
pi->conserved.flux.mass -= dti * A * totflux[0];
pi->conserved.flux.momentum[0] -= dti * A * totflux[1];
pi->conserved.flux.momentum[1] -= dti * A * totflux[2];
pi->conserved.flux.momentum[2] -= dti * A * totflux[3];
pi->conserved.flux.energy -= dti * A * totflux[4];
pi->conserved.flux.mass -= mindt * A * totflux[0];
pi->conserved.flux.momentum[0] -= mindt * A * totflux[1];
pi->conserved.flux.momentum[1] -= mindt * A * totflux[2];
pi->conserved.flux.momentum[2] -= mindt * A * totflux[3];
pi->conserved.flux.energy -= mindt * A * totflux[4];
#ifndef SHADOWFAX_TOTAL_ENERGY
float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
pi->primitives.v[1] * pi->primitives.v[1] +
pi->primitives.v[2] * pi->primitives.v[2]);
pi->conserved.flux.energy += dti * A * totflux[1] * pi->primitives.v[0];
pi->conserved.flux.energy += dti * A * totflux[2] * pi->primitives.v[1];
pi->conserved.flux.energy += dti * A * totflux[3] * pi->primitives.v[2];
pi->conserved.flux.energy -= dti * A * totflux[0] * ekin;
pi->conserved.flux.energy += mindt * A * totflux[1] * pi->primitives.v[0];
pi->conserved.flux.energy += mindt * A * totflux[2] * pi->primitives.v[1];
pi->conserved.flux.energy += mindt * A * totflux[3] * pi->primitives.v[2];
pi->conserved.flux.energy -= mindt * A * totflux[0] * ekin;
#endif
/* here is how it works:
Mode will only be 1 if both particles are ACTIVE and they are in the same
......@@ -287,19 +291,21 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
*/
if (mode == 1 || pj->force.active == 0) {
pj->conserved.flux.mass += dtj * A * totflux[0];
pj->conserved.flux.momentum[0] += dtj * A * totflux[1];
pj->conserved.flux.momentum[1] += dtj * A * totflux[2];
pj->conserved.flux.momentum[2] += dtj * A * totflux[3];
pj->conserved.flux.energy += dtj * A * totflux[4];
pj->conserved.flux.mass += mindt * A * totflux[0];
pj->conserved.flux.momentum[0] += mindt * A * totflux[1];
pj->conserved.flux.momentum[1] += mindt * A * totflux[2];
pj->conserved.flux.momentum[2] += mindt * A * totflux[3];
pj->conserved.flux.energy += mindt * A * totflux[4];
#ifndef SHADOWFAX_TOTAL_ENERGY
ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
pj->primitives.v[1] * pj->primitives.v[1] +
pj->primitives.v[2] * pj->primitives.v[2]);
pj->conserved.flux.energy -= dtj * A * totflux[1] * pj->primitives.v[0];
pj->conserved.flux.energy -= dtj * A * totflux[2] * pj->primitives.v[1];
pj->conserved.flux.energy -= dtj * A * totflux[3] * pj->primitives.v[2];
pj->conserved.flux.energy += dtj * A * totflux[0] * ekin;
pj->conserved.flux.energy -= mindt * A * totflux[1] * pj->primitives.v[0];
pj->conserved.flux.energy -= mindt * A * totflux[2] * pj->primitives.v[1];
pj->conserved.flux.energy -= mindt * A * totflux[3] * pj->primitives.v[2];
pj->conserved.flux.energy += mindt * A * totflux[0] * ekin;
#endif
}
}
......
......@@ -18,6 +18,7 @@
******************************************************************************/
#include "adiabatic_index.h"
#include "equation_of_state.h"
#include "hydro_gradients.h"
#include "hydro_slope_limiters.h"
#include "io_properties.h"
......@@ -63,7 +64,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
* @return Internal energy of the particle
*/
float convert_u(struct engine* e, struct part* p) {
return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
if (p->primitives.rho > 0.) {
return gas_internal_energy_from_pressure(p->primitives.rho,
p->primitives.P);
} else {
return 0.;
}
}
/**
......@@ -74,7 +80,11 @@ float convert_u(struct engine* e, struct part* p) {
* @return Entropic function of the particle
*/
float convert_A(struct engine* e, struct part* p) {
return p->primitives.P / pow_gamma(p->primitives.rho);
if (p->primitives.rho > 0.) {
return gas_entropy_from_pressure(p->primitives.rho, p->primitives.P);
} else {
return 0.;
}
}
/**
......@@ -85,13 +95,21 @@ float convert_A(struct engine* e, struct part* p) {
* @return Total energy of the particle
*/
float convert_Etot(struct engine* e, struct part* p) {
float momentum2;
momentum2 = p->conserved.momentum[0] * p->conserved.momentum[0] +
p->conserved.momentum[1] * p->conserved.momentum[1] +
p->conserved.momentum[2] * p->conserved.momentum[2];
return p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
#ifdef SHADOWFAX_TOTAL_ENERGY
return p->conserved.energy;
#else
if (p->conserved.mass > 0.) {
float momentum2;
momentum2 = p->conserved.momentum[0] * p->conserved.momentum[0] +
p->conserved.momentum[1] * p->conserved.momentum[1] +
p->conserved.momentum[2] * p->conserved.momentum[2];
return p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
} else {
return 0.;
}
#endif
}
/**
......
......@@ -31,9 +31,6 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
/* Old density. */
float omega;
/* Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
......@@ -42,6 +39,12 @@ struct xpart {
/* Data of a single particle. */
struct part {
/* Particle ID. */
long long id;
/* Associated gravitas. */
struct gpart *gpart;
/* Particle position. */
double x[3];
......@@ -163,18 +166,6 @@ struct part {
} force;
/* Particle mass (this field is also part of the conserved quantities...). */
float mass;
/* Particle ID. */
long long id;
/* Associated gravitas. */
struct gpart *gpart;
/* Variables needed for the code to compile (should be removed/replaced). */
float rho;
/* Time-step length */
timebin_t time_bin;
......
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