Skip to content
GitLab
Commit 2f6f0fd4 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' of gitlab.cosma.dur.ac.uk:swift/swiftsim

parents e87ca4f1 32ff7194
Hide whitespace changes
Inline Side-by-side
src/hydro/Gadget2/hydro_io.h
View file @ 2f6f0fd4
......@@ -107,9 +107,10 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
"(No treatment) Legacy Gadget-2 as in Springel (2005)");
writeAttribute_s(h_grpsph, "Viscosity Model",
"Legacy Gadget-2 as in Springel (2005)");
"(No treatment) as in Springel (2005)");
writeAttribute_s(
h_grpsph, "Viscosity Model",
"as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
}
......
src/hydro/Minimal/hydro.h
View file @ 2f6f0fd4
......@@ -16,10 +16,29 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_MINIMAL_HYDRO_H
#define SWIFT_MINIMAL_HYDRO_H
/**
* @file Minimal/hydro.h
* @brief Minimal conservative implementation of SPH (Non-neighbour loop
* equations)
*
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
* Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
* pp. 759-794.
*/
#include "adiabatic_index.h"
#include "approx_math.h"
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "kernel_hydro.h"
/**
* @brief Returns the internal energy of a particle
......@@ -34,7 +53,7 @@
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part *restrict p, float dt) {
return p->u;
return p->u + p->u_dt * dt;
}
/**
......@@ -46,7 +65,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
return gas_pressure_from_internal_energy(p->rho, p->u);
const float u = p->u + p->u_dt * dt;
return gas_pressure_from_internal_energy(p->rho, u);
}
/**
......@@ -62,7 +83,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part *restrict p, float dt) {
return gas_entropy_from_internal_energy(p->rho, p->u);
const float u = p->u + p->u_dt * dt;
return gas_entropy_from_internal_energy(p->rho, u);
}
/**
......@@ -74,7 +97,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) {
return gas_soundspeed_from_internal_energy(p->rho, p->u);
const float u = p->u + p->u_dt * dt;
return gas_soundspeed_from_internal_energy(p->rho, u);
}
/**
......@@ -150,7 +175,6 @@ __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->u_full = p->u;
}
/**
......@@ -164,6 +188,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
......@@ -227,7 +252,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) {
p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, half_dt);
p->force.pressure = pressure;
}
/**
......@@ -268,10 +296,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part *restrict p, const struct xpart *restrict xp, int t0, int t1,
double timeBase) {
p->u = xp->u_full;
/* Need to recompute the pressure as well */
p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
/* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
p->force.pressure = hydro_get_pressure(p, dt_entr);
}
/**
......@@ -304,11 +331,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt,
float half_dt) {
/* Kick in momentum space */
xp->u_full += p->u_dt * dt;
/* Do not decrease the energy by more than a factor of 2*/
const float u_change = p->u_dt * dt;
if (u_change > -0.5f * p->u)
p->u += u_change;
else
p->u *= 0.5f;
/* Get the predicted internal energy */
p->u = xp->u_full - half_dt * p->u_dt;
/* Do not 'overcool' when timestep increases */
if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt;
}
/**
......@@ -323,3 +354,5 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {}
#endif /* SWIFT_MINIMAL_HYDRO_H */
src/hydro/Minimal/hydro_debug.h
View file @ 2f6f0fd4
......@@ -16,18 +16,36 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_MINIMAL_HYDRO_DEBUG_H
#define SWIFT_MINIMAL_HYDRO_DEBUG_H
/**
* @file Minimal/hydro_debug.h
* @brief Minimal conservative implementation of SPH (Debugging routines)
*
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
* Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
* pp. 759-794.
*/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
const struct part* p, const struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
"u_full=%.3e, u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
"h=%.3e, dh/dt=%.3e "
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%d, t_end=%d\n",
"u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
"h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, "
"t_begin=%d, t_end=%d\n",
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->u_full, p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h,
p->force.h_dt, (int)p->density.wcount, p->mass, p->rho_dh, p->rho,
p->ti_begin, p->ti_end);
p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt,
(int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
p->ti_end);
}
#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
src/hydro/Minimal/hydro_iact.h
View file @ 2f6f0fd4
......@@ -16,18 +16,25 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_RUNNER_IACT_MINIMAL_H
#define SWIFT_RUNNER_IACT_MINIMAL_H
#include "adiabatic_index.h"
#ifndef SWIFT_MINIMAL_HYDRO_IACT_H
#define SWIFT_MINIMAL_HYDRO_IACT_H
/**
* @brief Minimal conservative implementation of SPH
* @file Minimal/hydro_iact.h
* @brief Minimal conservative implementation of SPH (Neighbour loop equations)
*
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* The thermal variable is the internal energy (u). No viscosity nor
* thermal conduction terms are implemented.
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
* Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
* pp. 759-794.
*/
#include "adiabatic_index.h"
/**
* @brief Density loop
*/
......@@ -115,6 +122,8 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
__attribute__((always_inline)) INLINE static void runner_iact_force(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
const float fac_mu = 1.f; /* Will change with cosmological integration */
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
......@@ -143,34 +152,60 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */
const float 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];
/* Compute sound speeds */
/* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
const float v_sig = ci + cj;
const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* SPH acceleration term */
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
const float sph_acc_term =
(P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
/* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term;
/* Use the force Luke ! */
pi->a_hydro[0] -= mj * sph_term * dx[0];
pi->a_hydro[1] -= mj * sph_term * dx[1];
pi->a_hydro[2] -= mj * sph_term * dx[2];
pi->a_hydro[0] -= mj * acc * dx[0];
pi->a_hydro[1] -= mj * acc * dx[1];
pi->a_hydro[2] -= mj * acc * dx[2];
pj->a_hydro[0] += mi * sph_term * dx[0];
pj->a_hydro[1] += mi * sph_term * dx[1];
pj->a_hydro[2] += mi * sph_term * dx[2];
pj->a_hydro[0] += mi * acc * dx[0];
pj->a_hydro[1] += mi * acc * dx[1];
pj->a_hydro[2] += mi * acc * dx[2];
/* Get the time derivative for u. */
pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
pj->u_dt += P_over_rho_j * mi * dvdr * r_inv * wj_dr;
const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
const float sph_du_term_j = P_over_rho2_j * dvdr * r_inv * wj_dr;
/* Viscosity term */
const float visc_du_term = 0.5f * visc_acc_term * dvdr;
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term;
const float du_dt_j = sph_du_term_j + visc_du_term;
/* Internal energy time derivatibe */
pi->u_dt += du_dt_i * mj;
pj->u_dt += du_dt_j * mi;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
......@@ -198,6 +233,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
const float fac_mu = 1.f; /* Will change with cosmological integration */
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
......@@ -226,29 +263,53 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */
const float 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];
/* Compute sound speeds */
/* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
const float v_sig = ci + cj;
const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* SPH acceleration term */
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
const float sph_acc_term =
(P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
/* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term;
/* Use the force Luke ! */
pi->a_hydro[0] -= mj * sph_term * dx[0];
pi->a_hydro[1] -= mj * sph_term * dx[1];
pi->a_hydro[2] -= mj * sph_term * dx[2];
pi->a_hydro[0] -= mj * acc * dx[0];
pi->a_hydro[1] -= mj * acc * dx[1];
pi->a_hydro[2] -= mj * acc * dx[2];
/* Get the time derivative for u. */
pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
/* Viscosity term */
const float visc_du_term = 0.5f * visc_acc_term * dvdr;
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term;
/* Internal energy time derivatibe */
pi->u_dt += du_dt_i * mj;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
......@@ -268,4 +329,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
"function does not exist yet!");
}
#endif /* SWIFT_RUNNER_IACT_MINIMAL_H */
#endif /* SWIFT_MINIMAL_HYDRO_IACT_H */
src/hydro/Minimal/hydro_io.h
View file @ 2f6f0fd4
......@@ -16,7 +16,25 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_MINIMAL_HYDRO_IO_H
#define SWIFT_MINIMAL_HYDRO_IO_H
/**
* @file Minimal/hydro_io.h
* @brief Minimal conservative implementation of SPH (i/o routines)
*
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
* Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
* pp. 759-794.
*/
#include "adiabatic_index.h"
#include "hydro.h"
#include "io_properties.h"
#include "kernel_hydro.h"
......@@ -51,6 +69,16 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
float convert_S(struct engine* e, struct part* p) {
return hydro_get_entropy(p, 0);
}
float convert_P(struct engine* e, struct part* p) {
return hydro_get_pressure(p, 0);
}
/**
* @brief Specifies which particle fields to write to a dataset
*
......@@ -61,7 +89,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
*num_fields = 8;
*num_fields = 10;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
......@@ -80,6 +108,11 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_ACCELERATION, parts, a_hydro);
list[7] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
UNIT_CONV_ENTROPY_PER_UNIT_MASS,
parts, rho, convert_S);
list[9] = io_make_output_field_convert_part(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
}
/**
......@@ -90,8 +123,9 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */
/* Nothing in this minimal model... */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No model");
writeAttribute_s(h_grpsph, "Viscosity Model", "No model");
writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
writeAttribute_s(h_grpsph, "Viscosity Model",
"Minimal treatment as in Monaghan (1992)");
/* Time integration properties */
writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
......@@ -104,3 +138,5 @@ void writeSPHflavour(hid_t h_grpsph) {
* @return 1 if entropy is in 'internal energy', 0 otherwise.
*/
int writeEntropyFlag() { return 0; }
#endif /* SWIFT_MINIMAL_HYDRO_IO_H */
src/hydro/Minimal/hydro_part.h
View file @ 2f6f0fd4
......@@ -16,6 +16,22 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_MINIMAL_HYDRO_PART_H
#define SWIFT_MINIMAL_HYDRO_PART_H
/**
* @file Minimal/hydro_part.h
* @brief Minimal conservative implementation of SPH (Particle definition)
*
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
* Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
* pp. 759-794.
*/
/**
* @brief Particle fields not needed during the SPH loops over neighbours.
......@@ -31,8 +47,6 @@ struct xpart {
float v_full[3]; /*!< Velocity at the last full step. */
float u_full; /*!< Thermal energy at the last full step. */
} __attribute__((aligned(xpart_align)));
/**
......@@ -107,3 +121,5 @@ struct part {
struct gpart* gpart; /*!< Pointer to corresponding gravity part. */
} __attribute__((aligned(part_align)));
#endif /* SWIFT_MINIMAL_HYDRO_PART_H */
src/hydro_properties.c
View file @ 2f6f0fd4
......@@ -70,7 +70,7 @@ void hydro_props_print(const struct hydro_props *p) {
message(
"Hydrodynamic integration: Max change of volume: %.2f "
"(max|dlog(h)/dt|=%f).",
powf(expf(p->log_max_h_change), hydro_dimension), p->log_max_h_change);
pow_dimension(expf(p->log_max_h_change)), p->log_max_h_change);
if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
message("Maximal iterations in ghost task set to %d (default is %d)",
......@@ -90,8 +90,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
writeAttribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
writeAttribute_f(h_grpsph, "Volume log(max(delta h))", p->log_max_h_change);
writeAttribute_f(h_grpsph, "Volume max change time-step",
powf(expf(p->log_max_h_change), 3.f));
writeAttribute_f(h_grpsph, "Max ghost iterations",
pow_dimension(expf(p->log_max_h_change)));
writeAttribute_i(h_grpsph, "Max ghost iterations",
p->max_smoothing_iterations);
}
#endif
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