Commit 32ff7194 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'fix_minimal_sph' into 'master'

Fix the minimal version of SPH

Ok, I have now fixed the minimalistic implementation of SPH. That's the one that should be used as a template for other people to implement their stuff. It's overly documented and now works on all the test-cases we have in the examples directory.  

Don't know if you want to look at this or not @pdraper given that it does not affect our default setup at all.

Solves #148.

See merge request !218
parents 954b572d 38398a02
...@@ -107,9 +107,10 @@ void writeSPHflavour(hid_t h_grpsph) { ...@@ -107,9 +107,10 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */ /* Viscosity and thermal conduction */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model", writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
"(No treatment) Legacy Gadget-2 as in Springel (2005)"); "(No treatment) as in Springel (2005)");
writeAttribute_s(h_grpsph, "Viscosity Model", writeAttribute_s(
"Legacy Gadget-2 as in Springel (2005)"); 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 alpha", const_viscosity_alpha);
writeAttribute_f(h_grpsph, "Viscosity beta", 3.f); writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
} }
......
...@@ -16,10 +16,29 @@ ...@@ -16,10 +16,29 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * 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 "adiabatic_index.h"
#include "approx_math.h" #include "dimension.h"
#include "equation_of_state.h" #include "equation_of_state.h"
#include "hydro_properties.h"
#include "kernel_hydro.h"
/** /**
* @brief Returns the internal energy of a particle * @brief Returns the internal energy of a particle
...@@ -34,7 +53,7 @@ ...@@ -34,7 +53,7 @@
__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, 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( ...@@ -46,7 +65,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
__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, 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( ...@@ -62,7 +83,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
__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, 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( ...@@ -74,7 +97,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
__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, 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( ...@@ -150,7 +175,6 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->v_full[0] = p->v[0]; xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1]; xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2]; 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( ...@@ -164,6 +188,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
*/ */
__attribute__((always_inline)) INLINE static void hydro_init_part( __attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p) { struct part *restrict p) {
p->density.wcount = 0.f; p->density.wcount = 0.f;
p->density.wcount_dh = 0.f; p->density.wcount_dh = 0.f;
p->rho = 0.f; p->rho = 0.f;
...@@ -227,7 +252,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -227,7 +252,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current, struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) { 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( ...@@ -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, struct part *restrict p, const struct xpart *restrict xp, int t0, int t1,
double timeBase) { double timeBase) {
p->u = xp->u_full; /* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
/* Need to recompute the pressure as well */ p->force.pressure = hydro_get_pressure(p, dt_entr);
p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
} }
/** /**
...@@ -304,11 +331,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -304,11 +331,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt, struct part *restrict p, struct xpart *restrict xp, float dt,
float half_dt) { float half_dt) {
/* Kick in momentum space */ /* Do not decrease the energy by more than a factor of 2*/
xp->u_full += p->u_dt * dt; 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 */ /* Do not 'overcool' when timestep increases */
p->u = xp->u_full - half_dt * p->u_dt; 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( ...@@ -323,3 +354,5 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/ */
__attribute__((always_inline)) INLINE static void hydro_convert_quantities( __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {} struct part *restrict p) {}
#endif /* SWIFT_MINIMAL_HYDRO_H */
...@@ -16,18 +16,36 @@ ...@@ -16,18 +16,36 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * 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( __attribute__((always_inline)) INLINE static void hydro_debug_particle(
const struct part* p, const struct xpart* xp) { const struct part* p, const struct xpart* xp) {
printf( printf(
"x=[%.3e,%.3e,%.3e], " "x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.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" "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
"h=%.3e, dh/dt=%.3e " "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, "
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%d, t_end=%d\n", "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], 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],
xp->u_full, p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt,
p->force.h_dt, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
p->ti_begin, p->ti_end); p->ti_end);
} }
#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
...@@ -16,18 +16,25 @@ ...@@ -16,18 +16,25 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * along with this program. If not, see <http://www.gnu.org/licenses/>.
* *
******************************************************************************/ ******************************************************************************/
#ifndef SWIFT_RUNNER_IACT_MINIMAL_H #ifndef SWIFT_MINIMAL_HYDRO_IACT_H
#define SWIFT_RUNNER_IACT_MINIMAL_H #define SWIFT_MINIMAL_HYDRO_IACT_H
#include "adiabatic_index.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 * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* thermal conduction terms are implemented. * \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 * @brief Density loop
*/ */
...@@ -115,6 +122,8 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj, ...@@ -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( __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, 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 = sqrtf(r2);
const float r_inv = 1.0f / r; const float r_inv = 1.0f / r;
...@@ -143,34 +152,60 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -143,34 +152,60 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float wj_dr = hjd_inv * wj_dx; const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */ /* Compute gradient terms */
const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh; const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh; const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */ /* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] + (pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2]; (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 ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(hydro_gamma * pressurej / rhoj); 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 */ /* 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 ! */ /* Use the force Luke ! */
pi->a_hydro[0] -= mj * sph_term * dx[0]; pi->a_hydro[0] -= mj * acc * dx[0];
pi->a_hydro[1] -= mj * sph_term * dx[1]; pi->a_hydro[1] -= mj * acc * dx[1];
pi->a_hydro[2] -= mj * sph_term * dx[2]; pi->a_hydro[2] -= mj * acc * dx[2];
pj->a_hydro[0] += mi * sph_term * dx[0]; pj->a_hydro[0] += mi * acc * dx[0];
pj->a_hydro[1] += mi * sph_term * dx[1]; pj->a_hydro[1] += mi * acc * dx[1];
pj->a_hydro[2] += mi * sph_term * dx[2]; pj->a_hydro[2] += mi * acc * dx[2];
/* Get the time derivative for u. */ /* 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;
pj->u_dt += P_over_rho_j * mi * dvdr * r_inv * wj_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. */ /* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr; 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( ...@@ -198,6 +233,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_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) { 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 = sqrtf(r2);
const float r_inv = 1.0f / r; const float r_inv = 1.0f / r;
...@@ -226,29 +263,53 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -226,29 +263,53 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float wj_dr = hjd_inv * wj_dx; const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */ /* Compute gradient terms */
const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh; const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh; const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */ /* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] + (pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2]; (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 ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(hydro_gamma * pressurej / rhoj); 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 */ /* 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 ! */ /* Use the force Luke ! */
pi->a_hydro[0] -= mj * sph_term * dx[0]; pi->a_hydro[0] -= mj * acc * dx[0];
pi->a_hydro[1] -= mj * sph_term * dx[1]; pi->a_hydro[1] -= mj * acc * dx[1];
pi->a_hydro[2] -= mj * sph_term * dx[2]; pi->a_hydro[2] -= mj * acc * dx[2];
/* Get the time derivative for u. */ /* 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. */ /* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr; 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( ...@@ -268,4 +329,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
"function does not exist yet!"); "function does not exist yet!");
} }
#endif /* SWIFT_RUNNER_IACT_MINIMAL_H */ #endif /* SWIFT_MINIMAL_HYDRO_IACT_H */
...@@ -16,7 +16,25 @@ ...@@ -16,7 +16,25 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * 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 "io_properties.h"
#include "kernel_hydro.h" #include "kernel_hydro.h"
...@@ -51,6 +69,16 @@ void hydro_read_particles(struct part* parts, struct io_props* list, ...@@ -51,6 +69,16 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho); 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 * @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, ...@@ -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, void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) { int* num_fields) {
*num_fields = 8; *num_fields = 10;
/* List what we want to write */ /* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 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, ...@@ -80,6 +108,11 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_ACCELERATION, parts, a_hydro); UNIT_CONV_ACCELERATION, parts, a_hydro);
list[7] = list[7] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho); 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) { ...@@ -90,8 +123,9 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */ /* Viscosity and thermal conduction */
/* Nothing in this minimal model... */ /* Nothing in this minimal model... */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No model"); writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
writeAttribute_s(h_grpsph, "Viscosity Model", "No model"); writeAttribute_s(h_grpsph, "Viscosity Model",
"Minimal treatment as in Monaghan (1992)");
/* Time integration properties */ /* Time integration properties */
writeAttribute_f(h_grpsph, "Maximal Delta u change over dt", writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
...@@ -104,3 +138,5 @@ void writeSPHflavour(hid_t h_grpsph) { ...@@ -104,3 +138,5 @@ void writeSPHflavour(hid_t h_grpsph) {
* @return 1 if entropy is in 'internal energy', 0 otherwise. * @return 1 if entropy is in 'internal energy', 0 otherwise.
*/ */
int writeEntropyFlag() { return 0; } int writeEntropyFlag() { return 0; }
#endif /* SWIFT_MINIMAL_HYDRO_IO_H */