Commit 8ac9c614 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Updated Minimal-SPH to the new scheme

parent 6264dfea
......@@ -210,26 +210,6 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
/**
* @brief Initialises the particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions or assignments between the particle
* and extended particle fields.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part *restrict p, struct xpart *restrict xp) {
p->ti_begin = 0;
p->ti_end = 0;
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
}
/**
* @brief Prepares a particle for the density calculation.
*
......@@ -292,16 +272,12 @@ __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 The current time (on the timeline)
* @param timeBase The minimal time-step size
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) {
struct part *restrict p, struct xpart *restrict xp) {
/* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, half_dt);
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
......@@ -339,6 +315,25 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
p->force.v_sig = 0.0f;
}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
*
* @param p The particle.
* @param xp The extended data of this particle.
*/
__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
struct part *restrict p, const struct xpart *restrict xp) {
/* Re-set the predicted velocities */
p->v[0] = xp->v_full[0];
p->v[1] = xp->v_full[1];
p->v[2] = xp->v_full[2];
/* Re-set the entropy */
p->u = xp->u_full;
}
/**
* @brief Predict additional particle fields forward in time when drifting
*
......@@ -348,13 +343,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
* @param p The particle.
* @param xp The extended data of the particle.
* @param dt The drift time-step.
* @param t0 The time at the start of the drift (on the timeline).
* @param t1 The time at the end of the drift (on the timeline).
* @param timeBase The minimal time-step size.
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part *restrict p, const struct xpart *restrict xp, float dt, int t0,
int t1, double timeBase) {
struct part *restrict p, const struct xpart *restrict xp, float dt) {
const float h_inv = 1.f / p->h;
......@@ -372,9 +363,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
else
p->rho *= expf(w2);
/* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, dt_entr);
/* Predict the internal energy */
p->u += p->u_dt * dt;
/* Compute the new pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
/* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
......@@ -407,24 +400,19 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt The time-step for this kick
* @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt,
float half_dt) {
struct part *restrict p, struct xpart *restrict xp, float 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;
if (u_change > -0.5f * xp->u_full)
xp->u_full += u_change;
else
p->u *= 0.5f;
/* 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;
xp->u_full *= 0.5f;
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full);
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
......@@ -442,9 +430,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* This can be used to convert internal energy into entropy for instance.
*
* @param p The particle to act upon
* @param xp The extended particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {
struct part *restrict p, struct xpart *restrict xp) {
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
......@@ -456,4 +445,27 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
p->force.soundspeed = soundspeed;
}
/**
* @brief Initialises the particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions or assignments between the particle
* and extended particle fields.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part *restrict p, struct xpart *restrict xp) {
p->time_bin = 0;
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;
hydro_reset_acceleration(p);
hydro_init_part(p);
}
#endif /* SWIFT_MINIMAL_HYDRO_H */
......@@ -40,12 +40,11 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.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",
"time_bin=%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],
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->density.rho_dh, p->rho, p->ti_begin,
p->ti_end);
(int)p->density.wcount, p->mass, p->density.rho_dh, p->rho, p->time_bin);
}
#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
......@@ -49,6 +49,9 @@ struct xpart {
/*! Velocity at the last full step. */
float v_full[3];
/*! Internal energy at the last full step. */
float u_full;
/*! Additional data used to record cooling information */
struct cooling_xpart_data cooling_data;
......@@ -63,6 +66,12 @@ struct xpart {
*/
struct part {
/*! Particle unique ID. */
long long id;
/*! Pointer to corresponding gravity part. */
struct gpart* gpart;
/*! Particle position. */
double x[3];
......@@ -78,12 +87,6 @@ struct part {
/*! Particle smoothing length. */
float h;
/*! Time at the beginning of time-step. */
int ti_begin;
/*! Time at the end of time-step. */
int ti_end;
/*! Particle internal energy. */
float u;
......@@ -143,11 +146,18 @@ struct part {
} force;
};
/*! Particle unique ID. */
long long id;
/*! Time-step length */
timebin_t time_bin;
/*! Pointer to corresponding gravity part. */
struct gpart* gpart;
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
integertime_t ti_drift;
/* Time of the last kick */
integertime_t ti_kick;
#endif
} SWIFT_STRUCT_ALIGN;
......
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