diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 025779f17496e7bf30fdf12353c4381c7d6292ce..d70d58c6ba508ba4282ac9dd32565478afb40692 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_SHADOWSWIFT_HYDRO_H
+#define SWIFT_SHADOWSWIFT_HYDRO_H
#include <float.h>
#include "adiabatic_index.h"
@@ -41,15 +43,23 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float CFL_condition = hydro_properties->CFL_condition;
- float R = get_radius_dimension_sphere(p->cell.volume);
-
- 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);
+ float vrel[3];
+ vrel[0] = p->primitives.v[0] - xp->v_full[0];
+ vrel[1] = p->primitives.v[1] - xp->v_full[1];
+ vrel[2] = p->primitives.v[2] - xp->v_full[2];
+ float vmax =
+ sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
+ sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+ vmax = max(vmax, p->timestepvars.vmax);
+
+ const float psize =
+ cosmo->a *
+ powf(p->cell.volume / hydro_dimension_unit_sphere, hydro_dimension_inv);
+ float dt = FLT_MAX;
+ if (vmax > 0.) {
+ dt = psize / vmax;
}
+ return CFL_condition * dt;
}
/**
@@ -102,7 +112,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->conserved.momentum[2] = mass * p->primitives.v[2];
#ifdef EOS_ISOTHERMAL_GAS
- p->conserved.energy = mass * const_isothermal_internal_energy;
+ p->conserved.energy = mass * gas_internal_energy_from_entropy(0.f, 0.f);
#else
p->conserved.energy *= mass;
#endif
@@ -117,12 +127,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->v[0] = 0.;
p->v[1] = 0.;
p->v[2] = 0.;
+#else
+ p->v[0] = p->primitives.v[0];
+ p->v[1] = p->primitives.v[1];
+ p->v[2] = p->primitives.v[2];
#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];
+
+ /* ignore accelerations present in the initial condition */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
}
/**
@@ -137,7 +156,8 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part* p, const struct hydro_space* hs) {
- p->density.wcount = 0.0f;
+ /* make sure we don't enter the no neighbour case in runner.c */
+ p->density.wcount = 1.0f;
p->density.wcount_dh = 0.0f;
voronoi_cell_init(&p->cell, p->x, hs->anchor, hs->side);
@@ -180,12 +200,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
if (hnew < p->h) {
/* Iteration succesful: we accept, but manually set h to a smaller value
for the next time step */
- p->density.wcount = 1.0f;
+ const float hinvdim = pow_dimension(1.0f / p->h);
+ p->density.wcount = hinvdim;
p->h = 1.1f * hnew;
} else {
/* 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);
+ p->density.wcount_dh = p->h / (1.1f * hnew - p->h);
return;
}
volume = p->cell.volume;
@@ -286,8 +307,48 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.v_full[0] = xp->v_full[0];
p->force.v_full[1] = xp->v_full[1];
p->force.v_full[2] = xp->v_full[2];
+
+ p->conserved.flux.mass = 0.0f;
+ p->conserved.flux.momentum[0] = 0.0f;
+ p->conserved.flux.momentum[1] = 0.0f;
+ p->conserved.flux.momentum[2] = 0.0f;
+ p->conserved.flux.energy = 0.0f;
+}
+
+/**
+ * @brief Prepare a particle for the gradient calculation.
+ *
+ * This function is called after the density loop and before the gradient loop.
+ *
+ * We use it to set the physical timestep for the particle and to copy the
+ * actual velocities, which we need to boost our interfaces during the flux
+ * calculation. We also initialize the variables used for the time step
+ * calculation.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
+ struct part* restrict p, struct xpart* restrict xp,
+ const struct cosmology* cosmo) {
+
+ /* Initialize time step criterion variables */
+ p->timestepvars.vmax = 0.;
}
+/**
+ * @brief Resets the variables that are required for a gradient calculation.
+ *
+ * This function is called after hydro_prepare_gradient.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_gradient(
+ struct part* restrict p) {}
+
/**
* @brief Finishes the gradient calculation.
*
@@ -385,53 +446,36 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {
- float vcell[3];
-
/* Update the conserved variables. We do this here and not in the kick,
since we need the updated variables below. */
- p->conserved.mass += p->conserved.flux.mass;
- p->conserved.momentum[0] += p->conserved.flux.momentum[0];
- p->conserved.momentum[1] += p->conserved.flux.momentum[1];
- p->conserved.momentum[2] += p->conserved.flux.momentum[2];
- p->conserved.energy += p->conserved.flux.energy;
+ p->conserved.mass += p->conserved.flux.mass * dt;
+ p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt;
+ p->conserved.momentum[1] += p->conserved.flux.momentum[1] * dt;
+ p->conserved.momentum[2] += p->conserved.flux.momentum[2] * dt;
#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]);
+ p->conserved.energy =
+ p->conserved.mass * gas_internal_energy_from_entropy(0.f, 0.f);
+#else
+ p->conserved.energy += p->conserved.flux.energy * dt;
#endif
-#endif
+#if defined(SHADOWFAX_FIX_CELLS)
+ p->v[0] = 0.0f;
+ p->v[1] = 0.0f;
+ p->v[2] = 0.0f;
+#else
+ if (p->conserved.mass > 0.0f && p->primitives.rho > 0.0f) {
- /* reset fluxes */
- /* we can only do this here, since we need to keep the fluxes for inactive
- particles */
- p->conserved.flux.mass = 0.0f;
- p->conserved.flux.momentum[0] = 0.0f;
- p->conserved.flux.momentum[1] = 0.0f;
- p->conserved.flux.momentum[2] = 0.0f;
- p->conserved.flux.energy = 0.0f;
+ const float inverse_mass = 1.f / 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.;
- }
+ /* Normal case: set particle velocity to fluid velocity. */
+ p->v[0] = p->conserved.momentum[0] * inverse_mass;
+ p->v[1] = p->conserved.momentum[1] * inverse_mass;
+ p->v[2] = p->conserved.momentum[2] * inverse_mass;
#ifdef SHADOWFAX_STEER_CELL_MOTION
- /* To prevent stupid things like cell crossovers or generators that move
- outside their cell, we steer the motion of the cell somewhat */
- if (p->primitives.rho) {
float centroid[3], d[3];
float volume, csnd, R, vfac, fac, dnrm;
voronoi_get_centroid(&p->cell, centroid);
@@ -449,32 +493,29 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
} else {
vfac = csnd / dnrm;
}
- } else {
- vfac = 0.0f;
+ p->v[0] += vfac * d[0];
+ p->v[1] += vfac * d[1];
+ p->v[2] += vfac * d[2];
}
- vcell[0] += vfac * d[0];
- vcell[1] += vfac * d[1];
- vcell[2] += vfac * d[2];
- }
#endif
-#if defined(SHADOWFAX_FIX_CELLS)
- xp->v_full[0] = 0.;
- xp->v_full[1] = 0.;
- xp->v_full[2] = 0.;
+ } else {
+ p->v[0] = 0.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+ }
+#endif
- p->v[0] = 0.;
- p->v[1] = 0.;
- p->v[2] = 0.;
-#else
- xp->v_full[0] = vcell[0];
- xp->v_full[1] = vcell[1];
- xp->v_full[2] = vcell[2];
+ /* Now make sure all velocity variables are up to date. */
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
- p->v[0] = xp->v_full[0];
- p->v[1] = xp->v_full[1];
- p->v[2] = xp->v_full[2];
-#endif
+ if (p->gpart) {
+ p->gpart->v_full[0] = p->v[0];
+ p->gpart->v_full[1] = p->v[1];
+ p->gpart->v_full[2] = p->v[2];
+ }
}
/**
@@ -799,3 +840,5 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_density(
return cosmo->a3_inv * p->primitives.rho;
}
+
+#endif /* SWIFT_SHADOWSWIFT_HYDRO_H */
diff --git a/src/hydro/Shadowswift/hydro_gradients.h b/src/hydro/Shadowswift/hydro_gradients.h
index 285d889a1a6e10662a06979f69290aabd4206059..4e7a9911d8d4fc586fe7a56687dd4c4ae9ec8de2 100644
--- a/src/hydro/Shadowswift/hydro_gradients.h
+++ b/src/hydro/Shadowswift/hydro_gradients.h
@@ -86,7 +86,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_predict(
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 r, float* xij_i, float* Wi, float* Wj) {
float dWi[5], dWj[5];
float xij_j[3];
@@ -132,59 +132,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
- /* time */
- dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
- Wi[2] * pi->primitives.gradients.rho[1] +
- Wi[3] * pi->primitives.gradients.rho[2] +
- Wi[0] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
- dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
- Wi[2] * pi->primitives.gradients.v[0][1] +
- Wi[3] * pi->primitives.gradients.v[0][2] +
- pi->primitives.gradients.P[0] / Wi[0]);
- dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
- Wi[2] * pi->primitives.gradients.v[1][1] +
- Wi[3] * pi->primitives.gradients.v[1][2] +
- pi->primitives.gradients.P[1] / Wi[0]);
- dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
- Wi[2] * pi->primitives.gradients.v[2][1] +
- Wi[3] * pi->primitives.gradients.v[2][2] +
- pi->primitives.gradients.P[2] / Wi[0]);
- dWi[4] -=
- 0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
- Wi[2] * pi->primitives.gradients.P[1] +
- Wi[3] * pi->primitives.gradients.P[2] +
- hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
-
- dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
- Wj[2] * pj->primitives.gradients.rho[1] +
- Wj[3] * pj->primitives.gradients.rho[2] +
- Wj[0] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
- dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
- Wj[2] * pj->primitives.gradients.v[0][1] +
- Wj[3] * pj->primitives.gradients.v[0][2] +
- pj->primitives.gradients.P[0] / Wj[0]);
- dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
- Wj[2] * pj->primitives.gradients.v[1][1] +
- Wj[3] * pj->primitives.gradients.v[1][2] +
- pj->primitives.gradients.P[1] / Wj[0]);
- dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
- Wj[2] * pj->primitives.gradients.v[2][1] +
- Wj[3] * pj->primitives.gradients.v[2][2] +
- pj->primitives.gradients.P[2] / Wj[0]);
- dWj[4] -=
- 0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
- Wj[2] * pj->primitives.gradients.P[1] +
- Wj[3] * pj->primitives.gradients.P[2] +
- hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
-
Wi[0] += dWi[0];
Wi[1] += dWi[1];
Wi[2] += dWi[2];
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index 9ac1debf3184c25603412867c41c62a1131345f3..eda8e3759d9e08dac8073ebed9fb36dd0c5b99f6 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -143,7 +143,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float vmax, dvdotdx;
float vi[3], vj[3], vij[3];
float Wi[5], Wj[5];
- float dti, dtj, mindt;
float n_unit[3];
A = voronoi_get_face(&pi->cell, pj->id, xij_i);
@@ -168,9 +167,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[3] = pj->primitives.v[2];
Wj[4] = pj->primitives.P;
- dti = pi->force.dt;
- dtj = pj->force.dt;
-
/* calculate the maximal signal velocity */
vmax = 0.0f;
if (Wi[0] > 0.) {
@@ -192,10 +188,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pj->timestepvars.vmax = fmaxf(pj->timestepvars.vmax, vmax);
}
- /* The flux will be exchanged using the smallest time step of the two
- * particles */
- mindt = fminf(dti, dtj);
-
/* compute the normal vector of the interface */
for (k = 0; k < 3; ++k) {
n_unit[k] = -dx[k] / r;
@@ -219,13 +211,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[2] -= vij[1];
Wj[3] -= vij[2];
- hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj, mindt);
+ hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
/* we don't need to rotate, we can use the unit vector in the Riemann problem
* itself (see GIZMO) */
if (Wi[0] < 0.0f || Wj[0] < 0.0f || Wi[4] < 0.0f || Wj[4] < 0.0f) {
- printf("mindt: %g\n", mindt);
printf("WL: %g %g %g %g %g\n", pi->primitives.rho, pi->primitives.v[0],
pi->primitives.v[1], pi->primitives.v[2], pi->primitives.P);
#ifdef USE_GRADIENTS
@@ -266,20 +257,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* Update conserved variables */
/* eqn. (16) */
- 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];
+ pi->conserved.flux.mass -= A * totflux[0];
+ pi->conserved.flux.momentum[0] -= A * totflux[1];
+ pi->conserved.flux.momentum[1] -= A * totflux[2];
+ pi->conserved.flux.momentum[2] -= A * totflux[3];
+ pi->conserved.flux.energy -= 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 += 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;
+ pi->conserved.flux.energy += A * totflux[1] * pi->primitives.v[0];
+ pi->conserved.flux.energy += A * totflux[2] * pi->primitives.v[1];
+ pi->conserved.flux.energy += A * totflux[3] * pi->primitives.v[2];
+ pi->conserved.flux.energy -= A * totflux[0] * ekin;
#endif
/* here is how it works:
@@ -295,20 +286,20 @@ __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 += 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];
+ pj->conserved.flux.mass += A * totflux[0];
+ pj->conserved.flux.momentum[0] += A * totflux[1];
+ pj->conserved.flux.momentum[1] += A * totflux[2];
+ pj->conserved.flux.momentum[2] += A * totflux[3];
+ pj->conserved.flux.energy += 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 -= 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;
+ pj->conserved.flux.energy -= A * totflux[1] * pj->primitives.v[0];
+ pj->conserved.flux.energy -= A * totflux[2] * pj->primitives.v[1];
+ pj->conserved.flux.energy -= A * totflux[3] * pj->primitives.v[2];
+ pj->conserved.flux.energy += A * totflux[0] * ekin;
#endif
}
}
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index c5448f77353e1859c1f8853394bbefbe26d0a3a9..5540d8c33c6263213b35697cfed2b6b5b07de2b7 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -72,6 +72,7 @@ void hydro_props_init(struct hydro_props *p,
/* change the meaning of target_neighbours and delta_neighbours */
p->target_neighbours = 1.0f;
p->delta_neighbours = 0.0f;
+ p->eta_neighbours = 1.0f;
#endif
/* Maximal smoothing length */
diff --git a/tests/test27cells.c b/tests/test27cells.c
index ada1b782cfff3866bf26937391007947e9c9a175..f62c169486250ba940c22f21ba1556cca4060c1a 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -217,8 +217,20 @@ void clean_up(struct cell *ci) {
* @brief Initializes all particles field to be ready for a density calculation
*/
void zero_particle_fields(struct cell *c) {
+#ifdef SHADOWFAX_SPH
+ struct hydro_space hs;
+ hs.anchor[0] = 0.;
+ hs.anchor[1] = 0.;
+ hs.anchor[2] = 0.;
+ hs.side[0] = 1.;
+ hs.side[1] = 1.;
+ hs.side[2] = 1.;
+ struct hydro_space *hspointer = &hs;
+#else
+ struct hydro_space *hspointer = NULL;
+#endif
for (int pid = 0; pid < c->count; pid++) {
- hydro_init_part(&c->parts[pid], NULL);
+ hydro_init_part(&c->parts[pid], hspointer);
}
}