diff --git a/src/hydro/Gizmo/MFM/hydro.h b/src/hydro/Gizmo/MFM/hydro.h
index bfbee24307c4c719ed596d4a460a1801757a4a59..1f81722abe644083f7b8ada58d5220ff09de0a07 100644
--- a/src/hydro/Gizmo/MFM/hydro.h
+++ b/src/hydro/Gizmo/MFM/hydro.h
@@ -36,343 +36,6 @@
#include <float.h>
-/**
- * @brief Finishes the volume calculation.
- *
- * Multiplies the density and number of neighbours by the appropiate constants
- * and adds the self-contribution term. Calculates the volume and uses it to
- * update the primitive variables (based on the conserved variables). The latter
- * should only be done for active particles. This is okay, since this method is
- * only called for active particles.
- *
- * Multiplies the components of the matrix E with the appropriate constants and
- * inverts it. Initializes the variables used during the gradient loop. This
- * cannot be done in hydro_prepare_force, since that method is called for all
- * particles, and not just the active ones. If we would initialize the
- * variables there, gradients for passive particles would be zero, while we
- * actually use the old gradients in the flux calculation between active and
- * passive particles.
- *
- * @param p The particle to act upon.
- * @param cosmo The cosmological model.
- */
-__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part* restrict p, const struct cosmology* cosmo) {
-
- /* Some smoothing length multiples. */
- const float h = p->h;
- const float ih = 1.0f / h;
- const float ihdim = pow_dimension(ih);
- const float ihdim_plus_one = ihdim * ih;
-
- /* Final operation on the density. */
- p->density.wcount += kernel_root;
- p->density.wcount *= ihdim;
-
- p->density.wcount_dh -= hydro_dimension * kernel_root;
- p->density.wcount_dh *= ihdim_plus_one;
-
- /* Final operation on the geometry. */
- /* we multiply with the smoothing kernel normalization ih3 and calculate the
- * volume */
- const float volume_inv = ihdim * (p->geometry.volume + kernel_root);
- const float volume = 1.0f / volume_inv;
- p->geometry.volume = volume;
-
- /* we multiply with the smoothing kernel normalization */
- p->geometry.matrix_E[0][0] = ihdim * p->geometry.matrix_E[0][0];
- p->geometry.matrix_E[0][1] = ihdim * p->geometry.matrix_E[0][1];
- p->geometry.matrix_E[0][2] = ihdim * p->geometry.matrix_E[0][2];
- p->geometry.matrix_E[1][0] = ihdim * p->geometry.matrix_E[1][0];
- p->geometry.matrix_E[1][1] = ihdim * p->geometry.matrix_E[1][1];
- p->geometry.matrix_E[1][2] = ihdim * p->geometry.matrix_E[1][2];
- p->geometry.matrix_E[2][0] = ihdim * p->geometry.matrix_E[2][0];
- p->geometry.matrix_E[2][1] = ihdim * p->geometry.matrix_E[2][1];
- p->geometry.matrix_E[2][2] = ihdim * p->geometry.matrix_E[2][2];
-
- p->geometry.centroid[0] *= kernel_norm;
- p->geometry.centroid[1] *= kernel_norm;
- p->geometry.centroid[2] *= kernel_norm;
-
- const float wcount_inv = 1.0f / p->density.wcount;
- p->geometry.centroid[0] *= wcount_inv;
- p->geometry.centroid[1] *= wcount_inv;
- p->geometry.centroid[2] *= wcount_inv;
-
- /* Check the condition number to see if we have a stable geometry. */
- float condition_number_E = 0.0f;
- int i, j;
- for (i = 0; i < 3; ++i) {
- for (j = 0; j < 3; ++j) {
- condition_number_E +=
- p->geometry.matrix_E[i][j] * p->geometry.matrix_E[i][j];
- }
- }
-
- invert_dimension_by_dimension_matrix(p->geometry.matrix_E);
-
- float condition_number_Einv = 0.0f;
- for (i = 0; i < 3; ++i) {
- for (j = 0; j < 3; ++j) {
- condition_number_Einv +=
- p->geometry.matrix_E[i][j] * p->geometry.matrix_E[i][j];
- }
- }
-
- float condition_number =
- hydro_dimension_inv * sqrtf(condition_number_E * condition_number_Einv);
-
- if (condition_number > const_gizmo_max_condition_number &&
- p->geometry.wcorr > const_gizmo_min_wcorr) {
-#ifdef GIZMO_PATHOLOGICAL_ERROR
- error("Condition number larger than %g (%g)!",
- const_gizmo_max_condition_number, condition_number);
-#endif
-#ifdef GIZMO_PATHOLOGICAL_WARNING
- message("Condition number too large: %g (> %g, p->id: %llu)!",
- condition_number, const_gizmo_max_condition_number, p->id);
-#endif
- /* add a correction to the number of neighbours for this particle */
- p->geometry.wcorr *= const_gizmo_w_correction_factor;
- }
-
- /* compute primitive variables */
- /* eqns (3)-(5) */
- const float m = p->conserved.mass;
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (m < 0.0f) {
- error("Mass is negative!");
- }
-
- if (volume == 0.0f) {
- error("Volume is 0!");
- }
-#endif
-
- // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
- float momentum[3];
- momentum[0] = p->conserved.momentum[0];
- momentum[1] = p->conserved.momentum[1];
- momentum[2] = p->conserved.momentum[2];
- p->rho = m * volume_inv;
- if (m == 0.0f) {
- p->v[0] = 0.0f;
- p->v[1] = 0.0f;
- p->v[2] = 0.0f;
- } else {
- const float m_inv = 1.0f / m;
- p->v[0] = momentum[0] * m_inv;
- p->v[1] = momentum[1] * m_inv;
- p->v[2] = momentum[2] * m_inv;
- }
-
-#ifdef EOS_ISOTHERMAL_GAS
- /* although the pressure is not formally used anywhere if an isothermal eos
- has been selected, we still make sure it is set to the correct value */
- p->P = gas_pressure_from_internal_energy(p->rho, 0.0f);
-#else
-
- float energy = p->conserved.energy;
-
-#ifdef GIZMO_TOTAL_ENERGY
- /* subtract the kinetic energy; we want the thermal energy */
- energy -= 0.5f * (momentum[0] * p->v[0] + momentum[1] * p->v[1] +
- momentum[2] * p->v[2]);
-#endif
-
- /* energy contains the total thermal energy, we want the specific energy.
- this is why we divide by the volume, and not by the density */
- p->P = hydro_gamma_minus_one * energy * volume_inv;
-#endif
-
- /* sanity checks */
- gizmo_check_physical_quantities("density", "pressure", p->rho, p->v[0],
- p->v[1], p->v[2], p->P);
-
- /* Add a correction factor to wcount (to force a neighbour number increase if
- the geometry matrix is close to singular) */
- p->density.wcount *= p->geometry.wcorr;
- p->density.wcount_dh *= p->geometry.wcorr;
-}
-
-/**
- * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
- *
- * @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_part_has_no_neighbours(
- struct part* restrict p, struct xpart* restrict xp,
- const struct cosmology* cosmo) {
-
- /* Some smoothing length multiples. */
- const float h = p->h;
- const float h_inv = 1.0f / h; /* 1/h */
- const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
-
- /* Re-set problematic values */
- p->density.wcount = kernel_root * h_inv_dim;
- p->density.wcount_dh = 0.0f;
- p->geometry.volume = 1.0f;
- p->geometry.matrix_E[0][0] = 1.0f;
- p->geometry.matrix_E[0][1] = 0.0f;
- p->geometry.matrix_E[0][2] = 0.0f;
- p->geometry.matrix_E[1][0] = 0.0f;
- p->geometry.matrix_E[1][1] = 1.0f;
- p->geometry.matrix_E[1][2] = 0.0f;
- p->geometry.matrix_E[2][0] = 0.0f;
- p->geometry.matrix_E[2][1] = 0.0f;
- p->geometry.matrix_E[2][2] = 1.0f;
- /* centroid is relative w.r.t. particle position */
- /* by setting the centroid to 0.0f, we make sure no velocity correction is
- applied */
- p->geometry.centroid[0] = 0.0f;
- p->geometry.centroid[1] = 0.0f;
- p->geometry.centroid[2] = 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.
- * @param hydro_props Hydrodynamic properties.
- */
-__attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
- struct part* restrict p, struct xpart* restrict xp,
- const struct cosmology* cosmo, const struct hydro_props* hydro_props) {
-
- /* Initialize time step criterion variables */
- p->timestepvars.vmax = 0.0f;
-
- hydro_gradients_init(p);
-
- // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
-}
-
-/**
- * @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.
- *
- * Just a wrapper around hydro_gradients_finalize, which can be an empty method,
- * in which case no gradients are used.
- *
- * This method also initializes the force loop variables.
- *
- * @param p The particle to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_end_gradient(
- struct part* p) {
-
- hydro_gradients_finalize(p);
-
-#ifdef GIZMO_LLOYD_ITERATION
- /* reset the gradients to zero, as we don't want them */
- hydro_gradients_init(p);
-#endif
-}
-
-/**
- * @brief Prepare a particle for the force calculation.
- *
- * This function is called in the ghost task to convert some quantities coming
- * from the density loop over neighbours into quantities ready to be used in the
- * force loop over neighbours. Quantities are typically read from the density
- * sub-structure and written to the force sub-structure.
- * Examples of calculations done here include the calculation of viscosity term
- * constants, thermal conduction terms, hydro conversions, etc.
- *
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- * @param cosmo The current cosmological model.
- * @param hydro_props Hydrodynamic properties.
- * @param dt_alpha The time-step used to evolve non-cosmological quantities such
- * as the artificial viscosity.
- */
-__attribute__((always_inline)) INLINE static void hydro_prepare_force(
- struct part* restrict p, struct xpart* restrict xp,
- const struct cosmology* cosmo, const struct hydro_props* hydro_props,
- const float dt_alpha) {
-
- /* Initialise values that are used in the force loop */
- p->flux.momentum[0] = 0.0f;
- p->flux.momentum[1] = 0.0f;
- p->flux.momentum[2] = 0.0f;
- p->flux.energy = 0.0f;
-}
-
-/**
- * @brief Reset acceleration fields of a particle
- *
- * This is actually not necessary for GIZMO, since we just set the accelerations
- * after the flux calculation.
- *
- * @param p The particle to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
- struct part* p) {
-
- /* Reset the acceleration. */
- p->a_hydro[0] = 0.0f;
- p->a_hydro[1] = 0.0f;
- p->a_hydro[2] = 0.0f;
-
- /* Reset the time derivatives. */
- p->force.h_dt = 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.
- * @param cosmo The cosmological model.
- */
-__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
- struct part* restrict p, const struct xpart* restrict xp,
- const struct cosmology* cosmo) {
- // MATTHIEU: Do we need something here?
-}
-
-/**
- * @brief Converts the hydrodynamic variables from the initial condition file to
- * conserved variables that can be used during the integration
- *
- * We no longer do this, as the mass needs to be provided in the initial
- * condition file, and the mass alone is enough to initialize all conserved
- * variables. This is now done in hydro_first_init_part.
- *
- * @param p The particle to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p, struct xpart* xp, const struct cosmology* cosmo,
- const struct hydro_props* hydro_props) {
-
- p->conserved.energy /= cosmo->a_factor_internal_energy;
-}
-
/**
* @brief Extra operations to be done during the drift
*
diff --git a/src/hydro/Gizmo/MFM/hydro_getters.h b/src/hydro/Gizmo/MFM/hydro_getters.h
index 9ae88b199607c5b4bfa7d3aeed24ae688cebdbee..c7353797002552566666f11179afb418af9185e7 100644
--- a/src/hydro/Gizmo/MFM/hydro_getters.h
+++ b/src/hydro/Gizmo/MFM/hydro_getters.h
@@ -40,6 +40,17 @@ hydro_part_get_primitive_variables(const struct part* restrict p, float* W) {
W[4] = p->P;
}
+/**
+ * @brief Get the degenerate case correction factor for the given particle.
+ *
+ * @param p Particle.
+ * @return Degenerate case correction factor.
+ */
+__attribute__((always_inline)) INLINE static float hydro_part_get_wcorr(
+ const struct part* restrict p) {
+ return p->geometry.wcorr;
+}
+
/**
* @brief Get the gradients of the primitive variables for the given particle.
*
diff --git a/src/hydro/Gizmo/MFM/hydro_setters.h b/src/hydro/Gizmo/MFM/hydro_setters.h
index dc19cb1b9f2b8f4950d907c2902ebf1bf3f6d0d3..a503833daccdce24124dfa25c64575230e3e5327 100644
--- a/src/hydro/Gizmo/MFM/hydro_setters.h
+++ b/src/hydro/Gizmo/MFM/hydro_setters.h
@@ -67,6 +67,20 @@ __attribute__((always_inline)) INLINE static void hydro_part_set_wcorr(
p->geometry.wcorr = wcorr;
}
+/**
+ * @brief Reset the fluxes for the given particle.
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_reset_fluxes(
+ struct part* restrict p) {
+
+ p->flux.momentum[0] = 0.0f;
+ p->flux.momentum[1] = 0.0f;
+ p->flux.momentum[2] = 0.0f;
+ p->flux.energy = 0.0f;
+}
+
/**
* @brief Update the fluxes for the particle with the given contributions,
* assuming the particle is to the left of the interparticle interface.
diff --git a/src/hydro/Gizmo/MFV/hydro.h b/src/hydro/Gizmo/MFV/hydro.h
index f6673fbba3dca6ca6f6f6be4cb845a0950f43540..30e06c9f771c38431e4bdf25fcc303d5469815da 100644
--- a/src/hydro/Gizmo/MFV/hydro.h
+++ b/src/hydro/Gizmo/MFV/hydro.h
@@ -37,357 +37,6 @@
#include <float.h>
-/**
- * @brief Finishes the volume calculation.
- *
- * Multiplies the density and number of neighbours by the appropiate constants
- * and adds the self-contribution term. Calculates the volume and uses it to
- * update the primitive variables (based on the conserved variables). The latter
- * should only be done for active particles. This is okay, since this method is
- * only called for active particles.
- *
- * Multiplies the components of the matrix E with the appropriate constants and
- * inverts it. Initializes the variables used during the gradient loop. This
- * cannot be done in hydro_prepare_force, since that method is called for all
- * particles, and not just the active ones. If we would initialize the
- * variables there, gradients for passive particles would be zero, while we
- * actually use the old gradients in the flux calculation between active and
- * passive particles.
- *
- * @param p The particle to act upon.
- * @param cosmo The cosmological model.
- */
-__attribute__((always_inline)) INLINE static void hydro_end_density(
- struct part* restrict p, const struct cosmology* cosmo) {
-
- /* Some smoothing length multiples. */
- const float h = p->h;
- const float ih = 1.0f / h;
- const float ihdim = pow_dimension(ih);
- const float ihdim_plus_one = ihdim * ih;
-
- /* Final operation on the density. */
- p->density.wcount += kernel_root;
- p->density.wcount *= ihdim;
-
- p->density.wcount_dh -= hydro_dimension * kernel_root;
- p->density.wcount_dh *= ihdim_plus_one;
-
- /* Final operation on the geometry. */
- /* we multiply with the smoothing kernel normalization ih3 and calculate the
- * volume */
- const float volume = 1.f / (ihdim * (p->geometry.volume + kernel_root));
- p->geometry.volume = volume;
-
- /* we multiply with the smoothing kernel normalization */
- p->geometry.matrix_E[0][0] = ihdim * p->geometry.matrix_E[0][0];
- p->geometry.matrix_E[0][1] = ihdim * p->geometry.matrix_E[0][1];
- p->geometry.matrix_E[0][2] = ihdim * p->geometry.matrix_E[0][2];
- p->geometry.matrix_E[1][0] = ihdim * p->geometry.matrix_E[1][0];
- p->geometry.matrix_E[1][1] = ihdim * p->geometry.matrix_E[1][1];
- p->geometry.matrix_E[1][2] = ihdim * p->geometry.matrix_E[1][2];
- p->geometry.matrix_E[2][0] = ihdim * p->geometry.matrix_E[2][0];
- p->geometry.matrix_E[2][1] = ihdim * p->geometry.matrix_E[2][1];
- p->geometry.matrix_E[2][2] = ihdim * p->geometry.matrix_E[2][2];
-
- p->geometry.centroid[0] *= kernel_norm;
- p->geometry.centroid[1] *= kernel_norm;
- p->geometry.centroid[2] *= kernel_norm;
-
- p->geometry.centroid[0] /= p->density.wcount;
- p->geometry.centroid[1] /= p->density.wcount;
- p->geometry.centroid[2] /= p->density.wcount;
-
- /* Check the condition number to see if we have a stable geometry. */
- float condition_number_E = 0.0f;
- int i, j;
- for (i = 0; i < 3; ++i) {
- for (j = 0; j < 3; ++j) {
- condition_number_E +=
- p->geometry.matrix_E[i][j] * p->geometry.matrix_E[i][j];
- }
- }
-
- invert_dimension_by_dimension_matrix(p->geometry.matrix_E);
-
- float condition_number_Einv = 0.0f;
- for (i = 0; i < 3; ++i) {
- for (j = 0; j < 3; ++j) {
- condition_number_Einv +=
- p->geometry.matrix_E[i][j] * p->geometry.matrix_E[i][j];
- }
- }
-
- float condition_number =
- hydro_dimension_inv * sqrtf(condition_number_E * condition_number_Einv);
-
- if (condition_number > const_gizmo_max_condition_number &&
- p->density.wcorr > const_gizmo_min_wcorr) {
-#ifdef GIZMO_PATHOLOGICAL_ERROR
- error("Condition number larger than %g (%g)!",
- const_gizmo_max_condition_number, condition_number);
-#endif
-#ifdef GIZMO_PATHOLOGICAL_WARNING
- message("Condition number too large: %g (> %g, p->id: %llu)!",
- condition_number, const_gizmo_max_condition_number, p->id);
-#endif
- /* add a correction to the number of neighbours for this particle */
- p->density.wcorr *= const_gizmo_w_correction_factor;
- }
-
- hydro_gradients_init(p);
-
- /* compute primitive variables */
- /* eqns (3)-(5) */
- const float m = p->conserved.mass;
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (m < 0.) {
- error("Mass is negative!");
- }
-
- if (volume == 0.) {
- error("Volume is 0!");
- }
-#endif
-
- // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
- float momentum[3];
- momentum[0] = p->conserved.momentum[0];
- momentum[1] = p->conserved.momentum[1];
- momentum[2] = p->conserved.momentum[2];
- p->primitives.rho = m / volume;
- if (m == 0.) {
- p->primitives.v[0] = 0.;
- p->primitives.v[1] = 0.;
- p->primitives.v[2] = 0.;
- } else {
- p->primitives.v[0] = momentum[0] / m;
- p->primitives.v[1] = momentum[1] / m;
- p->primitives.v[2] = momentum[2] / m;
- }
-
-#ifdef EOS_ISOTHERMAL_GAS
- /* although the pressure is not formally used anywhere if an isothermal eos
- has been selected, we still make sure it is set to the correct value */
- p->primitives.P = gas_pressure_from_internal_energy(p->primitives.rho, 0.);
-#else
-
- float energy = p->conserved.energy;
-
-#ifdef GIZMO_TOTAL_ENERGY
- /* subtract the kinetic energy; we want the thermal energy */
- energy -= 0.5f * (momentum[0] * p->primitives.v[0] +
- momentum[1] * p->primitives.v[1] +
- momentum[2] * p->primitives.v[2]);
-#endif
-
- /* energy contains the total thermal energy, we want the specific energy.
- this is why we divide by the volume, and not by the density */
- p->primitives.P = hydro_gamma_minus_one * energy / volume;
-#endif
-
- /* sanity checks */
- gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
- p->primitives.v[0], p->primitives.v[1],
- p->primitives.v[2], p->primitives.P);
-
-#ifdef GIZMO_LLOYD_ITERATION
- /* overwrite primitive variables to make sure they still have safe values */
- p->primitives.rho = 1.;
- p->primitives.v[0] = 0.;
- p->primitives.v[1] = 0.;
- p->primitives.v[2] = 0.;
- p->primitives.P = 1.;
-#endif
-
- /* Add a correction factor to wcount (to force a neighbour number increase if
- the geometry matrix is close to singular) */
- p->density.wcount *= p->density.wcorr;
- p->density.wcount_dh *= p->density.wcorr;
-}
-
-/**
- * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
- *
- * @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_part_has_no_neighbours(
- struct part* restrict p, struct xpart* restrict xp,
- const struct cosmology* cosmo) {
-
- /* Some smoothing length multiples. */
- const float h = p->h;
- const float h_inv = 1.0f / h; /* 1/h */
- const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
-
- /* Re-set problematic values */
- p->density.wcount = kernel_root * h_inv_dim;
- p->density.wcount_dh = 0.f;
- p->geometry.volume = 1.0f;
- p->geometry.matrix_E[0][0] = 1.0f;
- p->geometry.matrix_E[0][1] = 0.0f;
- p->geometry.matrix_E[0][2] = 0.0f;
- p->geometry.matrix_E[1][0] = 0.0f;
- p->geometry.matrix_E[1][1] = 1.0f;
- p->geometry.matrix_E[1][2] = 0.0f;
- p->geometry.matrix_E[2][0] = 0.0f;
- p->geometry.matrix_E[2][1] = 0.0f;
- p->geometry.matrix_E[2][2] = 1.0f;
- /* centroid is relative w.r.t. particle position */
- /* by setting the centroid to 0.0f, we make sure no velocity correction is
- applied */
- p->geometry.centroid[0] = 0.0f;
- p->geometry.centroid[1] = 0.0f;
- p->geometry.centroid[2] = 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.
- * @param hydro_props Hydrodynamic properties.
- */
-__attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
- struct part* restrict p, struct xpart* restrict xp,
- const struct cosmology* cosmo, const struct hydro_props* hydro_props) {
-
- /* Initialize time step criterion variables */
- p->timestepvars.vmax = 0.;
-
- // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
-
- /* Set the actual velocity of the particle */
- hydro_velocities_prepare_force(p, xp);
-}
-
-/**
- * @brief Finishes the gradient calculation.
- *
- * Just a wrapper around hydro_gradients_finalize, which can be an empty method,
- * in which case no gradients are used.
- *
- * @param p The particle to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_end_gradient(
- struct part* p) {
-
- hydro_gradients_finalize(p);
-
-#ifdef GIZMO_LLOYD_ITERATION
- /* reset the gradients to zero, as we don't want them */
- hydro_gradients_init(p);
-#endif
-}
-
-/**
- * @brief Prepare a particle for the force calculation.
- *
- * This function is called in the ghost task to convert some quantities coming
- * from the density loop over neighbours into quantities ready to be used in the
- * force loop over neighbours. Quantities are typically read from the density
- * sub-structure and written to the force sub-structure.
- * Examples of calculations done here include the calculation of viscosity term
- * constants, thermal conduction terms, hydro conversions, etc.
- *
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- * @param cosmo The current cosmological model.
- * @param hydro_props Hydrodynamic properties.
- * @param dt_alpha The time-step used to evolve non-cosmological quantities such
- * as the artificial viscosity.
- */
-__attribute__((always_inline)) INLINE static void hydro_prepare_force(
- struct part* restrict p, struct xpart* restrict xp,
- const struct cosmology* cosmo, const struct hydro_props* hydro_props,
- const float dt_alpha) {
-
- /* Initialise values that are used in the force loop */
- p->gravity.mflux[0] = 0.0f;
- p->gravity.mflux[1] = 0.0f;
- p->gravity.mflux[2] = 0.0f;
-
- 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 Reset acceleration fields of a particle
- *
- * This is actually not necessary for GIZMO, since we just set the accelerations
- * after the flux calculation.
- *
- * @param p The particle to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
- struct part* p) {
-
- /* Reset the acceleration. */
- p->a_hydro[0] = 0.0f;
- p->a_hydro[1] = 0.0f;
- p->a_hydro[2] = 0.0f;
-
- /* Reset the time derivatives. */
- p->force.h_dt = 0.0f;
-}
-
-/**
- * @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 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.
- * @param cosmo The cosmological model.
- */
-__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
- struct part* restrict p, const struct xpart* restrict xp,
- const struct cosmology* cosmo) {
- // MATTHIEU: Apply the entropy floor here.
-}
-
-/**
- * @brief Converts the hydrodynamic variables from the initial condition file to
- * conserved variables that can be used during the integration
- *
- * We no longer do this, as the mass needs to be provided in the initial
- * condition file, and the mass alone is enough to initialize all conserved
- * variables. This is now done in hydro_first_init_part.
- *
- * @param p The particle to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p, struct xpart* xp, const struct cosmology* cosmo,
- const struct hydro_props* hydro_props) {
-
- p->conserved.energy /= cosmo->a_factor_internal_energy;
-}
-
/**
* @brief Extra operations to be done during the drift
*
diff --git a/src/hydro/Gizmo/MFV/hydro_getters.h b/src/hydro/Gizmo/MFV/hydro_getters.h
index c66b675d673cd81809427bddc4343851dff85867..94aca11276ef2f67bea930f780663c4546ee578a 100644
--- a/src/hydro/Gizmo/MFV/hydro_getters.h
+++ b/src/hydro/Gizmo/MFV/hydro_getters.h
@@ -40,6 +40,17 @@ hydro_part_get_primitive_variables(const struct part* restrict p, float* W) {
W[4] = p->primitives.P;
}
+/**
+ * @brief Get the degenerate case correction factor for the given particle.
+ *
+ * @param p Particle.
+ * @return Degenerate case correction factor.
+ */
+__attribute__((always_inline)) INLINE static float hydro_part_get_wcorr(
+ const struct part* restrict p) {
+ return p->density.wcorr;
+}
+
/**
* @brief Get the gradients of the primitive variables for the given particle.
*
diff --git a/src/hydro/Gizmo/MFV/hydro_setters.h b/src/hydro/Gizmo/MFV/hydro_setters.h
index 9e5cc7d6202d3800526828e25b2cf8232380f991..770336ea63ee325b71fdbd5c2648bc9e2ff9189c 100644
--- a/src/hydro/Gizmo/MFV/hydro_setters.h
+++ b/src/hydro/Gizmo/MFV/hydro_setters.h
@@ -67,6 +67,25 @@ __attribute__((always_inline)) INLINE static void hydro_part_set_wcorr(
p->density.wcorr = wcorr;
}
+/**
+ * @brief Reset the fluxes for the given particle.
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_reset_fluxes(
+ struct part* restrict p) {
+
+ 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;
+
+ p->gravity.mflux[0] = 0.0f;
+ p->gravity.mflux[1] = 0.0f;
+ p->gravity.mflux[2] = 0.0f;
+}
+
/**
* @brief Update the fluxes for the particle with the given contributions,
* assuming the particle is to the left of the interparticle interface.
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 3f729e5eb6dd0aafbf741ad38c4a3502427a5931..4b59f498acba5dbc010273d3693e3eb0456ca3d7 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -22,6 +22,7 @@
//#define GIZMO_LLOYD_ITERATION
#include "hydro_getters.h"
+#include "hydro_gradients.h"
#include "hydro_setters.h"
#include "hydro_space.h"
@@ -213,6 +214,354 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->geometry.centroid[2] = 0.0f;
}
+/**
+ * @brief Finishes the volume calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and adds the self-contribution term. Calculates the volume and uses it to
+ * update the primitive variables (based on the conserved variables). The latter
+ * should only be done for active particles. This is okay, since this method is
+ * only called for active particles.
+ *
+ * Multiplies the components of the matrix E with the appropriate constants and
+ * inverts it. Initializes the variables used during the gradient loop. This
+ * cannot be done in hydro_prepare_force, since that method is called for all
+ * particles, and not just the active ones. If we would initialize the
+ * variables there, gradients for passive particles would be zero, while we
+ * actually use the old gradients in the flux calculation between active and
+ * passive particles.
+ *
+ * @param p The particle to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part* restrict p, const struct cosmology* cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float ih = 1.0f / h;
+ const float ihdim = pow_dimension(ih);
+ const float ihdim_plus_one = ihdim * ih;
+
+ /* Final operation on the density. */
+ p->density.wcount += kernel_root;
+ p->density.wcount *= ihdim;
+
+ p->density.wcount_dh -= hydro_dimension * kernel_root;
+ p->density.wcount_dh *= ihdim_plus_one;
+
+ /* Final operation on the geometry. */
+ /* we multiply with the smoothing kernel normalization ih3 and calculate the
+ * volume */
+ const float volume_inv = ihdim * (p->geometry.volume + kernel_root);
+ const float volume = 1.0f / volume_inv;
+ p->geometry.volume = volume;
+
+ /* we multiply with the smoothing kernel normalization */
+ p->geometry.matrix_E[0][0] = ihdim * p->geometry.matrix_E[0][0];
+ p->geometry.matrix_E[0][1] = ihdim * p->geometry.matrix_E[0][1];
+ p->geometry.matrix_E[0][2] = ihdim * p->geometry.matrix_E[0][2];
+ p->geometry.matrix_E[1][0] = ihdim * p->geometry.matrix_E[1][0];
+ p->geometry.matrix_E[1][1] = ihdim * p->geometry.matrix_E[1][1];
+ p->geometry.matrix_E[1][2] = ihdim * p->geometry.matrix_E[1][2];
+ p->geometry.matrix_E[2][0] = ihdim * p->geometry.matrix_E[2][0];
+ p->geometry.matrix_E[2][1] = ihdim * p->geometry.matrix_E[2][1];
+ p->geometry.matrix_E[2][2] = ihdim * p->geometry.matrix_E[2][2];
+
+ p->geometry.centroid[0] *= kernel_norm;
+ p->geometry.centroid[1] *= kernel_norm;
+ p->geometry.centroid[2] *= kernel_norm;
+
+ const float wcount_inv = 1.0f / p->density.wcount;
+ p->geometry.centroid[0] *= wcount_inv;
+ p->geometry.centroid[1] *= wcount_inv;
+ p->geometry.centroid[2] *= wcount_inv;
+
+ /* Check the condition number to see if we have a stable geometry. */
+ float condition_number_E = 0.0f;
+ int i, j;
+ for (i = 0; i < 3; ++i) {
+ for (j = 0; j < 3; ++j) {
+ condition_number_E +=
+ p->geometry.matrix_E[i][j] * p->geometry.matrix_E[i][j];
+ }
+ }
+
+ invert_dimension_by_dimension_matrix(p->geometry.matrix_E);
+
+ float condition_number_Einv = 0.0f;
+ for (i = 0; i < 3; ++i) {
+ for (j = 0; j < 3; ++j) {
+ condition_number_Einv +=
+ p->geometry.matrix_E[i][j] * p->geometry.matrix_E[i][j];
+ }
+ }
+
+ const float condition_number =
+ hydro_dimension_inv * sqrtf(condition_number_E * condition_number_Einv);
+
+ if (condition_number > const_gizmo_max_condition_number &&
+ hydro_part_get_wcorr(p) > const_gizmo_min_wcorr) {
+#ifdef GIZMO_PATHOLOGICAL_ERROR
+ error("Condition number larger than %g (%g)!",
+ const_gizmo_max_condition_number, condition_number);
+#endif
+#ifdef GIZMO_PATHOLOGICAL_WARNING
+ message("Condition number too large: %g (> %g, p->id: %llu)!",
+ condition_number, const_gizmo_max_condition_number, p->id);
+#endif
+ /* add a correction to the number of neighbours for this particle */
+ hydro_part_set_wcorr(
+ p, const_gizmo_w_correction_factor * hydro_part_get_wcorr(p));
+ }
+
+ /* need to figure out whether or not to move this to hydro_prepare_gradient */
+ hydro_gradients_init(p);
+
+ /* compute primitive variables */
+ /* eqns (3)-(5) */
+ const float Q[5] = {p->conserved.mass, p->conserved.momentum[0],
+ p->conserved.momentum[1], p->conserved.momentum[2],
+ p->conserved.energy};
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (Q[0] < 0.) {
+ error("Mass is negative!");
+ }
+
+ if (volume == 0.) {
+ error("Volume is 0!");
+ }
+#endif
+
+ float W[5];
+
+ W[0] = Q[0] * volume_inv;
+ if (Q[0] == 0.0f) {
+ W[1] = 0.;
+ W[2] = 0.;
+ W[3] = 0.;
+ } else {
+ const float m_inv = 1.0f / Q[0];
+ W[1] = Q[1] * m_inv;
+ W[2] = Q[2] * m_inv;
+ W[3] = Q[3] * m_inv;
+ }
+
+#ifdef EOS_ISOTHERMAL_GAS
+ /* although the pressure is not formally used anywhere if an isothermal eos
+ has been selected, we still make sure it is set to the correct value */
+ W[4] = gas_pressure_from_internal_energy(W[0], 0.0f);
+#else
+
+#ifdef GIZMO_TOTAL_ENERGY
+ /* subtract the kinetic energy; we want the thermal energy */
+ Q[4] -= 0.5f * (Q[1] * W[1] + Q[2] * W[2] + Q[3] * W[3]);
+#endif
+
+ /* energy contains the total thermal energy, we want the specific energy.
+ this is why we divide by the volume, and not by the density */
+ W[4] = hydro_gamma_minus_one * Q[4] * volume_inv;
+#endif
+
+ /* sanity checks */
+ gizmo_check_physical_quantities("density", "pressure", W[0], W[1], W[2], W[3],
+ W[4]);
+
+#ifdef GIZMO_LLOYD_ITERATION
+ /* overwrite primitive variables to make sure they still have safe values */
+ W[0] = 1.0f;
+ W[1] = 0.0f;
+ W[2] = 0.0f;
+ W[3] = 0.0f;
+ W[4] = 1.0f;
+#endif
+
+ hydro_part_set_primitive_variables(p, W);
+
+ /* Add a correction factor to wcount (to force a neighbour number increase if
+ the geometry matrix is close to singular) */
+ p->density.wcount *= hydro_part_get_wcorr(p);
+ p->density.wcount_dh *= hydro_part_get_wcorr(p);
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * @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_part_has_no_neighbours(
+ struct part* restrict p, struct xpart* restrict xp,
+ const struct cosmology* cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+ /* Re-set problematic values */
+ p->density.wcount = kernel_root * h_inv_dim;
+ p->density.wcount_dh = 0.f;
+ p->geometry.volume = 1.0f;
+ p->geometry.matrix_E[0][0] = 1.0f;
+ p->geometry.matrix_E[0][1] = 0.0f;
+ p->geometry.matrix_E[0][2] = 0.0f;
+ p->geometry.matrix_E[1][0] = 0.0f;
+ p->geometry.matrix_E[1][1] = 1.0f;
+ p->geometry.matrix_E[1][2] = 0.0f;
+ p->geometry.matrix_E[2][0] = 0.0f;
+ p->geometry.matrix_E[2][1] = 0.0f;
+ p->geometry.matrix_E[2][2] = 1.0f;
+ /* centroid is relative w.r.t. particle position */
+ /* by setting the centroid to 0.0f, we make sure no velocity correction is
+ applied */
+ p->geometry.centroid[0] = 0.0f;
+ p->geometry.centroid[1] = 0.0f;
+ p->geometry.centroid[2] = 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.
+ * @param hydro_props Hydrodynamic properties.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
+ struct part* restrict p, struct xpart* restrict xp,
+ const struct cosmology* cosmo, const struct hydro_props* hydro_props) {
+
+ /* Initialize time step criterion variables */
+ p->timestepvars.vmax = 0.;
+
+ /* not sure if we need to do this here or in end_density() */
+ hydro_gradients_init(p);
+
+#if defined(GIZMO_MFV_SPH)
+ /* Set the actual velocity of the particle */
+ hydro_velocities_prepare_force(p, xp);
+#endif
+}
+
+/**
+ * @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.
+ *
+ * Just a wrapper around hydro_gradients_finalize, which can be an empty method,
+ * in which case no gradients are used.
+ *
+ * This method also initializes the force loop variables.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_gradient(
+ struct part* p) {
+
+ hydro_gradients_finalize(p);
+
+#ifdef GIZMO_LLOYD_ITERATION
+ /* reset the gradients to zero, as we don't want them */
+ hydro_gradients_init(p);
+#endif
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * This function is called in the ghost task to convert some quantities coming
+ * from the density loop over neighbours into quantities ready to be used in the
+ * force loop over neighbours. Quantities are typically read from the density
+ * sub-structure and written to the force sub-structure.
+ * Examples of calculations done here include the calculation of viscosity term
+ * constants, thermal conduction terms, hydro conversions, etc.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties.
+ * @param dt_alpha The time-step used to evolve non-cosmological quantities such
+ * as the artificial viscosity.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part* restrict p, struct xpart* restrict xp,
+ const struct cosmology* cosmo, const struct hydro_props* hydro_props,
+ const float dt_alpha) {
+
+ hydro_part_reset_fluxes(p);
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * This is actually not necessary for GIZMO, since we just set the accelerations
+ * after the flux calculation.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part* p) {
+
+ /* Reset the acceleration. */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ /* Reset the time derivatives. */
+ p->force.h_dt = 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.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+ struct part* restrict p, const struct xpart* restrict xp,
+ const struct cosmology* cosmo) {
+ // MATTHIEU: Apply the entropy floor here.
+}
+
+/**
+ * @brief Converts the hydrodynamic variables from the initial condition file to
+ * conserved variables that can be used during the integration
+ *
+ * We no longer do this, as the mass needs to be provided in the initial
+ * condition file, and the mass alone is enough to initialize all conserved
+ * variables. This is now done in hydro_first_init_part.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part* p, struct xpart* xp, const struct cosmology* cosmo,
+ const struct hydro_props* hydro_props) {
+
+ p->conserved.energy /= cosmo->a_factor_internal_energy;
+}
+
#if defined(GIZMO_MFV_SPH)
#include "MFV/hydro.h"
#define SPH_IMPLEMENTATION "GIZMO MFV (Hopkins 2015)"