From 332e9d44074ce458a0b97f656d4cf2483dc7d2fd Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@gmail.com>
Date: Thu, 16 May 2019 11:54:27 +0100
Subject: [PATCH] Added optional flag variable to detect usage of entropy
switch and degeneracy handling code. Added parameters for entropy switch.
Made sure mass of MFM particles is never set to zero to prevent RMS
displacement timestep from dropping to zero.
---
src/Makefile.am | 3 +-
src/hydro/GizmoMFM/hydro.h | 100 ++++++++++-------
src/hydro/GizmoMFM/hydro_flag_variable.h | 102 ++++++++++++++++++
src/hydro/GizmoMFM/hydro_io.h | 2 +
src/hydro/GizmoMFM/hydro_parameters.h | 14 ++-
src/hydro/GizmoMFM/hydro_part.h | 6 ++
.../GizmoMFM/hydro_slope_limiters_cell.h | 2 +-
src/hydro/GizmoMFM/hydro_unphysical.h | 1 -
src/hydro/GizmoMFV/hydro.h | 14 ++-
src/timestep.h | 23 ++++
10 files changed, 215 insertions(+), 52 deletions(-)
create mode 100644 src/hydro/GizmoMFM/hydro_flag_variable.h
diff --git a/src/Makefile.am b/src/Makefile.am
index 99dce90211..7ef359658b 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 9fb2fdd2c1..5b3864a43e 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 0000000000..a0cd5116f2
--- /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 1f956edf3f..7f7e747372 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 5c345d6055..4a81a6bdde 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 17bc03836c..b1de8a5385 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 b77f76f47d..1ce1635adf 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 81c644e6cc..9d2d98f448 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 caf2819b08..312cba76c2 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 c2b1a10fcb..f87c8ce1c2 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);
--
GitLab