diff --git a/src/cell_grid.c b/src/cell_grid.c
index 0908b478ba195f87593517d003d0262cbe469d20..302d60f51a99d8dcbddc48a6a75e176ba628243f 100644
--- a/src/cell_grid.c
+++ b/src/cell_grid.c
@@ -141,14 +141,14 @@ void cell_grid_update_self_completeness(struct cell *c, int force) {
c->grid.self_completeness = grid_complete;
c->grid.complete = 1;
} else {
-#ifdef SHADOWSWIFT_RELAXED_COMPLETENESS
- if (kernel_gamma * c->hydro.h_max < 0.5 * c->dmin &&
- c->hydro.dx_max_part < 0.333 * c->dmin) {
- c->grid.self_completeness = grid_complete;
- c->grid.complete = 1;
- return;
- }
-#endif
+//#ifdef SHADOWSWIFT_RELAXED_COMPLETENESS
+// if (kernel_gamma * c->hydro.h_max < 0.5 * c->dmin &&
+// c->hydro.dx_max_part < 0.333 * c->dmin) {
+// c->grid.self_completeness = grid_complete;
+// c->grid.complete = 1;
+// return;
+// }
+//#endif
c->grid.self_completeness = grid_incomplete;
c->grid.complete = 0;
}
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index b6771c95459d46f641070cda00169a9584fb2392..e9e8304fbec0773195146cc6244638c1da8aa19e 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -396,12 +396,6 @@ hydro_convert_conserved_to_primitive(
/* Check for vacuum (or near limits of floating point precision), which give
* rise to precision errors resulting in NaNs... */
-#ifdef SWIFT_DEBUG_CHECKS
- if (Q[0] < 0. || Q[4] < 0.)
- warning(
- "Negative mass or energy after applying fluxes! Q[0] = %E, Q[4] = %E",
- Q[0], Q[4]);
-#endif
const double epsilon = 16. * FLT_MIN;
if (Q[0] < epsilon || Q[4] < epsilon) {
for (int k = 0; k < 6; k++) {
@@ -620,6 +614,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
if (p->timestepvars.last_kick == KICK1) {
/* I.e. we are in kick2 (end of timestep), since the dt_therm > 0. */
+ float Q[6];
+ hydro_part_get_conserved_variables(p, Q);
+ float dE_springel;
/* Add gravity. We only do this if we have gravity activated. */
if (p->gpart) {
/* Retrieve the current value of the gravitational acceleration from the
@@ -640,9 +637,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* apply both half kicks to the momentum */
/* Note that this also affects the particle movement, as the velocity for
the particles is set after this. */
- p->conserved.momentum[0] += grav_kick[0];
- p->conserved.momentum[1] += grav_kick[1];
- p->conserved.momentum[2] += grav_kick[2];
+ Q[1] += grav_kick[0];
+ Q[2] += grav_kick[1];
+ Q[3] += grav_kick[2];
/* Extra hydrodynamic energy due to gravity kick, see eq. 94 in Springel
* (2010) or eq. 62 in theory/Cosmology/cosmology.pdf. */
@@ -655,19 +652,17 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
dt_grav_corr1 = p->gravity.dt;
dt_grav_corr2 = dt_grav;
}
- float dE_springel = hydro_gravity_energy_update_term(
+ dE_springel = hydro_gravity_energy_update_term(
dt_grav_corr1, dt_grav_corr2, xp->a_grav, a_grav, p->gravity.mflux,
p->v_part_full, grav_kick);
- p->conserved.energy += dE_springel;
+ Q[4] += dE_springel;
}
#ifdef SWIFT_DEBUG_CHECKS
assert(p->flux.dt >= 0.0f);
#endif
- float Q[6];
- hydro_part_get_conserved_variables(p, Q);
if (p->flux.dt > 0.0f) {
/* Apply the fluxes */
/* We are in kick2 of a normal timestep (not the very beginning of the
@@ -706,6 +701,20 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
#endif
}
+ #ifdef SHADOWSWIFT_WARNINGS
+ if (Q[0] < 0. || Q[4] < 0.)
+ warning(
+ "Negative mass or energy after applying fluxes! "
+ "\n\tupdated mass = %E, updated energy = %E, "
+ "\n\told mass = %E, old energy = %E, "
+ "\n\tmflux = %E, energyflux = %E, dE_grav = %E"
+ "\n\tdensity = %E, pressure = %E"
+ "\n\tp->x = (%E, %E, %E), p->v = (%E, %E, %E)",
+ Q[0], Q[4], p->conserved.mass, p->conserved.energy, p->flux.mass,
+ p->flux.energy, dE_springel, p->rho, p->P, p->x[0], p->x[1], p->x[2],
+ p->v[0], p->v[1], p->v[2]);
+ #endif
+
/* Compute primitive quantities. Note that this may also modify the vector
* of conserved quantities (e.g. for cold flows).
* This also applies entropy floor/minimal internal energy limit. */
diff --git a/src/hydro/Shadowswift/hydro_gravity.h b/src/hydro/Shadowswift/hydro_gravity.h
index 140cf6cc5ee7b8fff11bf2f9974943ed92225354..a00aebcf2099e16e56f2c1de7c950e0b5ed55a9a 100644
--- a/src/hydro/Shadowswift/hydro_gravity.h
+++ b/src/hydro/Shadowswift/hydro_gravity.h
@@ -68,6 +68,19 @@ hydro_gravity_mass_update_term(const float mass_flux, const float dt) {
/**
* @brief Applies the gravitational work term at the face between pi and pj to
* both particles.
+ *
+ * @param pi, pj The particles
+ * @param shift Shift to appie to pj's coordinates, if wrapping
+ * around simulation box
+ * @param Whalf The state of primitive vectors at the interface (output from
+ * riemann solver, in reference frame of interface)
+ * @param vij Interface velocity
+ * @param cij Interface centroid
+ * @param n_unit Interface normal. Note that this is also the direction from
+ * pi->x to pj->x (voronoi faces are perpendicular to the corresponding Delaunay
+ * edges).
+ * @param area Interface area
+ * @param dt the timestep over which fluxes are exchanged. (currently unused)
*/
__attribute__((always_inline)) INLINE static void
hydro_grav_work_from_half_state(struct part* pi, struct part* pj,
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index 84e25f47fb92a9308a967815396e4724b7dc1a24..0570e550bd47f450d58b3efbddd8ba336a9ea074 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -253,6 +253,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_flux_exchange(
totflux[3] *= surface_area;
totflux[4] *= surface_area;
totflux[5] = surface_area * entropy_flux;
+#ifdef SHADOWSWIFT_FLUX_LIMITER
+ hydro_part_positivity_limiter_fluxes(pi, pj, n_unit, vij, surface_area,
+ hydro->epsilon_rho, hydro->epsilon_P,
+ totflux);
+#endif
hydro_grav_work_from_half_state(pi, pj, shift, Whalf, vij, centroid, n_unit,
surface_area, min_dt);
diff --git a/src/hydro/Shadowswift/hydro_iact_boundary.h b/src/hydro/Shadowswift/hydro_iact_boundary.h
index b302beec2febf4daf8ce6137a55a204f28d6ea7e..8173117079c9a94abc175a7051d49bdc8ac471cc 100644
--- a/src/hydro/Shadowswift/hydro_iact_boundary.h
+++ b/src/hydro/Shadowswift/hydro_iact_boundary.h
@@ -394,7 +394,7 @@ runner_iact_boundary_flux_exchange(struct part *p, struct part *p_boundary,
#ifdef SHADOWSWIFT_EXACT_GRAV_WORK
double shift[3] = {0., 0., 0.};
- hydro_grav_work_from_half_state(p, p_boundary, shift, Whalf, vij, centroid,
+ hydro_grav_work_from_half_state(p, p_boundary, shift, WL, vij, centroid,
n_unit, surface_area, p->flux.dt);
#else
hydro_grav_work_from_mass_flux(p, p_boundary, dx, totflux[0], p->flux.dt);