diff --git a/src/Makefile.am b/src/Makefile.am
index 99dce90211806eee062a9acf81b77ddca9fb3842..7ef359658ba071e2a29c05f3876e2f109db4f4f4 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -141,7 +141,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/GizmoMFM/hydro_slope_limiters_face.h \
hydro/GizmoMFM/hydro_slope_limiters.h \
hydro/GizmoMFM/hydro_unphysical.h \
- hydro/GizmoMFM/hydro_parameters.h \
+ hydro/GizmoMFM/hydro_parameters.h \
+ hydro/GizmoMFM/hydro_flag_variable.h \
hydro/Shadowswift/hydro_debug.h \
hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
hydro/Shadowswift/hydro_iact.h \
diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index 9fb2fdd2c1be45896f94c5d20e8e98eaf3b20fb0..5b3864a43e92b0aaac9e2a6bc99fe405e927bd83 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -25,6 +25,7 @@
#include "cosmology.h"
#include "entropy_floor.h"
#include "equation_of_state.h"
+#include "hydro_flag_variable.h"
#include "hydro_gradients.h"
#include "hydro_properties.h"
#include "hydro_space.h"
@@ -188,6 +189,8 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->geometry.centroid[0] = 0.0f;
p->geometry.centroid[1] = 0.0f;
p->geometry.centroid[2] = 0.0f;
+
+ hydro_flag_variable_init(p);
}
/**
@@ -349,6 +352,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
the geometry matrix is close to singular) */
p->density.wcount *= p->geometry.wcorr;
p->density.wcount_dh *= p->geometry.wcorr;
+
+ if (p->geometry.wcorr < 1.f) {
+ hydro_add_flag(p, GIZMO_FLAG_MORE_NEIGHBOURS);
+ }
+ if (p->geometry.wcorr <= const_gizmo_min_wcorr) {
+ hydro_add_flag(p, GIZMO_FLAG_REVERT_TO_SPH);
+ }
}
/**
@@ -529,9 +539,9 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
p->conserved.energy /= cosmo->a_factor_internal_energy;
/* initialize the entropy */
- p->conserved.entropy =
- p->conserved.energy * pow_minus_gamma_minus_one(p->rho) *
- hydro_gamma_minus_one;
+ p->conserved.entropy = p->conserved.energy *
+ pow_minus_gamma_minus_one(p->rho) *
+ hydro_gamma_minus_one;
}
/**
@@ -566,6 +576,10 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->h *= h_corr;
}
+ p->v[0] = xp->v_full[0];
+ p->v[1] = xp->v_full[1];
+ p->v[2] = xp->v_full[2];
+
/* drift the primitive variables based on the old fluxes */
if (p->conserved.mass > 0.0f) {
const float m_inv = 1.0f / p->conserved.mass;
@@ -666,31 +680,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.momentum[2] -= Pcorr * p->gradients.P[2];
#endif
- /* Apply the minimal energy limit */
- const float min_energy =
- hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
- if (p->conserved.energy < min_energy * p->conserved.mass) {
- p->conserved.energy = min_energy * p->conserved.mass;
- p->flux.energy = 0.0f;
- }
-
- // MATTHIEU: Apply the entropy floor here.
-
- gizmo_check_physical_quantities(
- "mass", "energy", p->conserved.mass, p->conserved.momentum[0],
- p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.energy);
-
-#ifdef SWIFT_DEBUG_CHECKS
- /* Note that this check will only have effect if no GIZMO_UNPHYSICAL option
- was selected. */
- if (p->conserved.energy < 0.0f) {
- error(
- "Negative energy after conserved variables update (energy: %g, "
- "denergy: %g)!",
- p->conserved.energy, p->flux.energy);
- }
-#endif
-
/* Add gravity. We only do this if we have gravity activated. */
if (p->gpart) {
/* Retrieve the current value of the gravitational acceleration from the
@@ -711,18 +700,51 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
const float dEgrav = p->conserved.mass *
sqrtf(a_grav[0] * a_grav[0] + a_grav[1] * a_grav[1] +
a_grav[2] * a_grav[2]) *
- p->h;
- if (p->conserved.energy <
- 0.001 * (p->force.Ekinmax + p->conserved.energy) ||
- p->conserved.energy < 0.001 * dEgrav) {
+ cosmo->a_inv * p->h;
+ if (p->conserved.energy < const_entropy_switch_ekin_fac *
+ (p->force.Ekinmax + p->conserved.energy) ||
+ p->conserved.energy < const_entropy_switch_grav_fac * dEgrav) {
p->conserved.energy = hydro_one_over_gamma_minus_one *
p->conserved.entropy * pow_gamma_minus_one(p->rho);
+ if (p->conserved.energy < const_entropy_switch_ekin_fac *
+ (p->force.Ekinmax + p->conserved.energy)) {
+ hydro_add_flag(p, GIZMO_FLAG_EKIN_SWITCH);
+ }
+ if (p->conserved.energy < const_entropy_switch_grav_fac * dEgrav) {
+ hydro_add_flag(p, GIZMO_FLAG_GRAVITY_SWITCH);
+ }
+ hydro_add_flag(p, GIZMO_FLAG_ENTROPY);
} else {
p->conserved.entropy = hydro_gamma_minus_one * p->conserved.energy *
pow_minus_gamma_minus_one(p->rho);
}
}
+ gizmo_check_physical_quantities(
+ "mass", "energy", p->conserved.mass, p->conserved.momentum[0],
+ p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.energy);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Note that this check will only have effect if no GIZMO_UNPHYSICAL option
+ was selected. */
+ if (p->conserved.energy < 0.0f) {
+ error(
+ "Negative energy after conserved variables update (energy: %g, "
+ "denergy: %g)!",
+ p->conserved.energy, p->flux.energy);
+ }
+#endif
+
+ /* Apply the minimal energy limit */
+ const float min_energy =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+ if (p->conserved.energy < min_energy * p->conserved.mass) {
+ p->conserved.energy = min_energy * p->conserved.mass;
+ p->conserved.entropy = hydro_gamma_minus_one * p->conserved.energy *
+ pow_minus_gamma_minus_one(p->rho);
+ p->flux.energy = 0.0f;
+ }
+
/* Set the velocities: */
/* We first set the particle velocity */
if (p->conserved.mass > 0.0f && p->rho > 0.0f) {
@@ -1124,10 +1146,9 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
p->conserved.energy = u * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic energy */
- p->conserved.energy +=
- 0.5f *
- (p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
- p->conserved.momentum[2] * p->v[2]);
+ p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->v[0] +
+ p->conserved.momentum[1] * p->v[1] +
+ p->conserved.momentum[2] * p->v[2]);
#endif
p->P = hydro_gamma_minus_one * p->rho * u;
}
@@ -1148,10 +1169,9 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
hydro_one_over_gamma_minus_one * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic energy */
- p->conserved.energy +=
- 0.5f *
- (p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
- p->conserved.momentum[2] * p->v[2]);
+ p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->v[0] +
+ p->conserved.momentum[1] * p->v[1] +
+ p->conserved.momentum[2] * p->v[2]);
#endif
p->P = S * pow_gamma(p->rho);
}
diff --git a/src/hydro/GizmoMFM/hydro_flag_variable.h b/src/hydro/GizmoMFM/hydro_flag_variable.h
new file mode 100644
index 0000000000000000000000000000000000000000..a0cd5116f2e29b341a717dcaee188396937b554c
--- /dev/null
+++ b/src/hydro/GizmoMFM/hydro_flag_variable.h
@@ -0,0 +1,102 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_HYDRO_GIZMO_MFM_FLAG_VARIABLE_H
+#define SWIFT_HYDRO_GIZMO_MFM_FLAG_VARIABLE_H
+
+#include "io_properties.h"
+
+enum gizmo_flag_variables {
+ GIZMO_FLAG_ENTROPY = 0,
+ GIZMO_FLAG_MORE_NEIGHBOURS = 1,
+ GIZMO_FLAG_REVERT_TO_SPH = 2,
+ GIZMO_FLAG_EKIN_SWITCH = 3,
+ GIZMO_FLAG_GRAVITY_SWITCH = 4
+};
+
+#if defined(GIZMO_FLAG_VARIABLE)
+
+/**
+ * @brief Initialize the flag variable.
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_flag_variable_init(
+ struct part *p) {
+ p->flag_variable = 0;
+}
+
+/**
+ * @brief Add the flag variable to the snapshot output.
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+__attribute__((always_inline)) INLINE static void hydro_write_flag_variable(
+ const struct part* parts, struct io_props* list, int* num_fields) {
+
+ list[*num_fields] = io_make_output_field("FlagVariable", INT, 1, UNIT_CONV_NO_UNITS, parts, flag_variable);
+ ++(*num_fields);
+}
+
+/**
+ * @brief Add the given flag value to the flag variable.
+ *
+ * @param p Particle.
+ * @param value Flag value to add.
+ */
+__attribute__((always_inline)) INLINE static void hydro_add_flag(
+ struct part* p, int value) {
+
+ p->flag_variable |= (1 << value);
+}
+
+#else
+
+/**
+ * @brief Initialize the flag variable.
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_flag_variable_init(
+ struct part *p) {}
+
+/**
+ * @brief Add the flag variable to the snapshot output.
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+__attribute__((always_inline)) INLINE static void hydro_write_flag_variable(
+ const struct part* parts, struct io_props* list, int* num_fields) {}
+
+/**
+ * @brief Add the given flag value to the flag variable.
+ *
+ * @param p Particle.
+ * @param value Flag value to add.
+ */
+__attribute__((always_inline)) INLINE static void hydro_add_flag(
+ struct part* p, int value) {}
+
+#endif
+
+#endif /* SWIFT_HYDRO_GIZMO_MFM_FLAG_VARIABLE_H */
diff --git a/src/hydro/GizmoMFM/hydro_io.h b/src/hydro/GizmoMFM/hydro_io.h
index 1f956edf3fdc31990c6aba254603ea69a98238eb..7f7e7473722896bfc4d3fd854977a062c1eb00dc 100644
--- a/src/hydro/GizmoMFM/hydro_io.h
+++ b/src/hydro/GizmoMFM/hydro_io.h
@@ -215,6 +215,8 @@ INLINE static void hydro_write_particles(const struct part* parts,
list[10] = io_make_output_field_convert_part("Potential", FLOAT, 1,
UNIT_CONV_POTENTIAL, parts,
xparts, convert_part_potential);
+
+ hydro_write_flag_variable(parts, list, num_fields);
}
/**
diff --git a/src/hydro/GizmoMFM/hydro_parameters.h b/src/hydro/GizmoMFM/hydro_parameters.h
index 5c345d6055be330ff7c8145c88edecc609bb9009..4a81a6bddefeedd255d636420362db76803b1d01 100644
--- a/src/hydro/GizmoMFM/hydro_parameters.h
+++ b/src/hydro/GizmoMFM/hydro_parameters.h
@@ -52,6 +52,15 @@
* Beta is defined as in e.g. Price (2010) Eqn (103) */
#define const_viscosity_beta 3.0f
+/* Prefactor for the kinetic energy condition for the entropy switch. */
+#define const_entropy_switch_ekin_fac 0.0001f
+
+/* Prefactor for the gravitational energy condition for the entropy switch. */
+#define const_entropy_switch_grav_fac 0.0001f
+
+/*! Activate this to write a diagnostic flag variable to the snapshots. */
+#define GIZMO_FLAG_VARIABLE
+
/* Structs that store the relevant variables */
/*! Artificial viscosity parameters */
@@ -110,6 +119,9 @@ static INLINE void viscosity_print(
**/
static INLINE void viscosity_print_snapshot(
hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ /* dummy to make the plot scripts happy */
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", 0.0f);
io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
}
#endif
@@ -159,4 +171,4 @@ static INLINE void diffusion_print_snapshot(
hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
#endif
-#endif /* SWIFT_GIZMOMFM_HYDRO_PARAMETERS_H */
\ No newline at end of file
+#endif /* SWIFT_GIZMOMFM_HYDRO_PARAMETERS_H */
diff --git a/src/hydro/GizmoMFM/hydro_part.h b/src/hydro/GizmoMFM/hydro_part.h
index 17bc03836c99530c2734efac86875a281d54c2c8..b1de8a5385902bdeb7e099e2744dce4472d7742a 100644
--- a/src/hydro/GizmoMFM/hydro_part.h
+++ b/src/hydro/GizmoMFM/hydro_part.h
@@ -21,6 +21,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
+#include "hydro_parameters.h"
#include "star_formation_struct.h"
#include "tracers_struct.h"
@@ -206,6 +207,11 @@ struct part {
/* Need waking-up ? */
timebin_t wakeup;
+#ifdef GIZMO_FLAG_VARIABLE
+ /* Flag variable containing diagnostic information about the particle. */
+ int flag_variable;
+#endif
+
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
diff --git a/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h b/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
index b77f76f47dda66a40526ac76eb906b23976d8f8d..1ce1635adfa834ac9531f05e80528ad1e18f0801 100644
--- a/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
+++ b/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
@@ -23,7 +23,7 @@
/*! @brief Maximum value for the extrapolated quantity as a function of the
* maximum value across the particle neighbours. */
-#define HYDRO_SLOPE_LIMITER_ALPHAMAX 2.0f
+#define HYDRO_SLOPE_LIMITER_ALPHAMAX 1.0f
/**
* @brief Initialize variables for the cell wide slope limiter
diff --git a/src/hydro/GizmoMFM/hydro_unphysical.h b/src/hydro/GizmoMFM/hydro_unphysical.h
index 81c644e6cce00e0d871aafc39140e420e042a926..9d2d98f4487dc66bfea30248cdce30ac4a576ccd 100644
--- a/src/hydro/GizmoMFM/hydro_unphysical.h
+++ b/src/hydro/GizmoMFM/hydro_unphysical.h
@@ -52,7 +52,6 @@
gizmo_check_physical_quantity(energy_name, energy); \
/* now check for vacuum and make sure we have a real vacuum */ \
if (mass == 0.f || energy == 0.f) { \
- mass = 0.f; \
momentum_x = 0.f; \
momentum_y = 0.f; \
momentum_z = 0.f; \
diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h
index caf2819b086ba51fa448108a6b9b2c9305305d77..312cba76c2029f81b7656dd9c16b24c31fe3bf49 100644
--- a/src/hydro/GizmoMFV/hydro.h
+++ b/src/hydro/GizmoMFV/hydro.h
@@ -1164,10 +1164,9 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
p->conserved.energy = u * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic 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 += 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]);
#endif
p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u;
}
@@ -1188,10 +1187,9 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
hydro_one_over_gamma_minus_one * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic 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 += 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]);
#endif
p->primitives.P = S * pow_gamma(p->primitives.rho);
}
diff --git a/src/timestep.h b/src/timestep.h
index c2b1a10fcb3b0426e7c34625d65c1fd5353d25e9..f87c8ce1c21819e157379b859ccb3b9ed1839bd9 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -155,6 +155,22 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
? fabsf(e->hydro_properties->log_max_h_change * p->h / p->force.h_dt)
: FLT_MAX;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (new_dt_grav == 0.0f) {
+ error("Zero gravity time step!");
+ }
+ if (new_dt_hydro == 0.0f) {
+ printParticle_single(p, xp);
+ error("Zero hydro time step!");
+ }
+ if (new_dt_cooling == 0.0f) {
+ error("Zero cooling time step!");
+ }
+ if (dt_h_change == 0.0f) {
+ error("Zero h_change time step!");
+ }
+#endif
+
new_dt = min(new_dt, dt_h_change);
/* Apply the maximal displacement constraint (FLT_MAX if non-cosmological)*/
@@ -163,6 +179,13 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
/* Apply cosmology correction (This is 1 if non-cosmological) */
new_dt *= e->cosmology->time_step_factor;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (new_dt == 0.0f) {
+ error("Zero timestep! (%e, %e)", e->dt_max_RMS_displacement,
+ e->cosmology->time_step_factor);
+ }
+#endif
+
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);