diff --git a/src/engine.c b/src/engine.c
index 584cf20146418ecdba67f786a1bd1e7165e3caa0..a8ec2d1ddaf9220c7b2553fb223c968c53768eec 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1827,7 +1827,11 @@ void engine_maketasks(struct engine *e) {
is the number of cells (s->tot_cells) times the number of neighbours (27)
times the number of interaction types (2, density and force). */
if (e->links != NULL) free(e->links);
+#ifdef EXTRA_HYDRO_LOOP
+ e->size_links = s->tot_cells * 27 * 3;
+#else
e->size_links = s->tot_cells * 27 * 2;
+#endif
if ((e->links = malloc(sizeof(struct link) * e->size_links)) == NULL)
error("Failed to allocate cell-task links.");
e->nr_links = 0;
@@ -2716,7 +2720,7 @@ void engine_step(struct engine *e) {
mask |= 1 << task_type_sub_pair;
mask |= 1 << task_type_ghost;
mask |= 1 << task_type_extra_ghost; /* Adding unnecessary things to the mask
- does not harm */
+ does not harm */
submask |= 1 << task_subtype_density;
submask |= 1 << task_subtype_gradient;
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 7bc75059e680e96585263e951f1a0e12fcd57fad..3dc8b62918c9df34cbfbbd0c7f4f34a7c524b185 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -21,14 +21,12 @@
#include "adiabatic_index.h"
#include "hydro_gradients.h"
-#define THERMAL_ENERGY
-
/**
* @brief Computes the hydro time-step of a given particle
*
- * @param p Pointer to the particle data
- * @param xp Pointer to the extended particle data
- *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ * @param hydro_properties Pointer to the hydro parameters.
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const struct part* restrict p, const struct xpart* restrict xp,
@@ -45,6 +43,12 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* This function is called only once just after the ICs have been
* read in to do some conversions.
*
+ * In this case, we copy the particle velocities into the corresponding
+ * primitive variable field. We do this because the particle velocities in GIZMO
+ * can be independent of the actual fluid velocity. The latter is stored as a
+ * primitive variable and integrated using the linear momentum, a conserved
+ * variable.
+ *
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
@@ -59,14 +63,13 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
/**
* @brief Prepares a particle for the volume calculation.
*
+ * Simply makes sure all necessary variables are initialized to zero.
+ *
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part* p) {
- /* We need to do this before we reset the volume */
- hydro_gradients_init_density_loop(p);
-
p->density.wcount = 0.0f;
p->density.wcount_dh = 0.0f;
p->geometry.volume = 0.0f;
@@ -82,12 +85,24 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
}
/**
- * @brief Finishes the density calculation.
+ * @brief Finishes the volume calculation.
*
* Multiplies the density and number of neighbours by the appropiate constants
- * and add the self-contribution term.
+ * 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.
*
- * @param p The particle to act upon
+ * 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 The current physical time.
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
struct part* restrict p, float time) {
@@ -107,10 +122,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
float volume;
float m, momentum[3], energy;
-#ifndef THERMAL_ENERGY
- float momentum2;
-#endif
-
/* Final operation on the geometry. */
/* we multiply with the smoothing kernel normalization ih3 and calculate the
* volume */
@@ -129,7 +140,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
invert_dimension_by_dimension_matrix(p->geometry.matrix_E);
- hydro_gradients_prepare_force_loop(p, ih, volume);
+ hydro_gradients_init(p);
/* compute primitive variables */
/* eqns (3)-(5) */
@@ -138,31 +149,30 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
momentum[0] = p->conserved.momentum[0];
momentum[1] = p->conserved.momentum[1];
momentum[2] = p->conserved.momentum[2];
-#ifndef THERMAL_ENERGY
- momentum2 = (momentum[0] * momentum[0] + momentum[1] * momentum[1] +
- momentum[2] * momentum[2]);
-#endif
p->primitives.rho = m / volume;
p->primitives.v[0] = momentum[0] / m;
p->primitives.v[1] = momentum[1] / m;
p->primitives.v[2] = momentum[2] / m;
energy = p->conserved.energy;
-#ifndef THERMAL_ENERGY
- p->primitives.P =
- hydro_gamma_minus_one * (energy - 0.5 * momentum2 / m) / volume;
-#else
p->primitives.P = hydro_gamma_minus_one * energy / volume;
-#endif
}
}
/**
- * @brief Prepare a particle for the force calculation.
+ * @brief Prepare a particle for the gradient calculation.
*
- * Computes viscosity term, conduction term and smoothing length gradient terms.
+ * The name of this method is confusing, as this method is really called after
+ * the density loop and before the gradient loop.
*
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
+ * 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 ti_current Current integer time.
+ * @param timeBase Conversion factor between integer time and physical time.
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp, int ti_current,
@@ -170,6 +180,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Set the physical time step */
p->force.dt = (p->ti_end - p->ti_begin) * timeBase;
+
+ /* Initialize time step criterion variables */
+ p->timestepvars.vmax = 0.0f;
+
/* Set the actual velocity of the particle */
p->force.v_full[0] = xp->v_full[0];
p->force.v_full[1] = xp->v_full[1];
@@ -179,166 +193,24 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/**
* @brief Finishes the gradient calculation.
*
- * @param p The particle to act upon
+ * 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) {
-#ifndef SPH_GRADIENTS
- float h, ih, ih2, ih3;
-#ifdef SLOPE_LIMITER
- float gradrho[3], gradv[3][3], gradP[3];
- float gradtrue, gradmax, gradmin, alpha;
-#endif
-
- /* add kernel normalization to gradients */
- h = p->h;
- ih = 1.0f / h;
- ih2 = ih * ih;
- ih3 = ih * ih2;
-
- p->primitives.gradients.rho[0] *= ih3;
- p->primitives.gradients.rho[1] *= ih3;
- p->primitives.gradients.rho[2] *= ih3;
-
- p->primitives.gradients.v[0][0] *= ih3;
- p->primitives.gradients.v[0][1] *= ih3;
- p->primitives.gradients.v[0][2] *= ih3;
- p->primitives.gradients.v[1][0] *= ih3;
- p->primitives.gradients.v[1][1] *= ih3;
- p->primitives.gradients.v[1][2] *= ih3;
- p->primitives.gradients.v[2][0] *= ih3;
- p->primitives.gradients.v[2][1] *= ih3;
- p->primitives.gradients.v[2][2] *= ih3;
-
- p->primitives.gradients.P[0] *= ih3;
- p->primitives.gradients.P[1] *= ih3;
- p->primitives.gradients.P[2] *= ih3;
-
-/* slope limiter */
-#ifdef SLOPE_LIMITER
- gradrho[0] = p->primitives.gradients.rho[0];
- gradrho[1] = p->primitives.gradients.rho[1];
- gradrho[2] = p->primitives.gradients.rho[2];
-
- gradv[0][0] = p->primitives.gradients.v[0][0];
- gradv[0][1] = p->primitives.gradients.v[0][1];
- gradv[0][2] = p->primitives.gradients.v[0][2];
-
- gradv[1][0] = p->primitives.gradients.v[1][0];
- gradv[1][1] = p->primitives.gradients.v[1][1];
- gradv[1][2] = p->primitives.gradients.v[1][2];
-
- gradv[2][0] = p->primitives.gradients.v[2][0];
- gradv[2][1] = p->primitives.gradients.v[2][1];
- gradv[2][2] = p->primitives.gradients.v[2][2];
-
- gradP[0] = p->primitives.gradients.P[0];
- gradP[1] = p->primitives.gradients.P[1];
- gradP[2] = p->primitives.gradients.P[2];
-
- gradtrue = gradrho[0] * gradrho[0] + gradrho[1] * gradrho[1] +
- gradrho[2] * gradrho[2];
- /* gradtrue might be zero. In this case, there is no gradient and we don't
- need to slope limit anything... */
- if (gradtrue) {
- gradtrue = sqrtf(gradtrue);
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
- gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
- /* gradmin and gradmax might be negative if the value of the current
- particle is larger/smaller than all neighbouring values */
- gradmax = fabs(gradmax);
- gradmin = fabs(gradmin);
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.rho[0] *= alpha;
- p->primitives.gradients.rho[1] *= alpha;
- p->primitives.gradients.rho[2] *= alpha;
- }
-
- gradtrue = gradv[0][0] * gradv[0][0] + gradv[0][1] * gradv[0][1] +
- gradv[0][2] * gradv[0][2];
- if (gradtrue) {
- gradtrue = sqrtf(gradtrue);
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
- gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
- gradmax = fabs(gradmax);
- gradmin = fabs(gradmin);
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.v[0][0] *= alpha;
- p->primitives.gradients.v[0][1] *= alpha;
- p->primitives.gradients.v[0][2] *= alpha;
- }
-
- gradtrue = gradv[1][0] * gradv[1][0] + gradv[1][1] * gradv[1][1] +
- gradv[1][2] * gradv[1][2];
- if (gradtrue) {
- gradtrue = sqrtf(gradtrue);
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
- gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
- gradmax = fabs(gradmax);
- gradmin = fabs(gradmin);
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.v[1][0] *= alpha;
- p->primitives.gradients.v[1][1] *= alpha;
- p->primitives.gradients.v[1][2] *= alpha;
- }
-
- gradtrue = gradv[2][0] * gradv[2][0] + gradv[2][1] * gradv[2][1] +
- gradv[2][2] * gradv[2][2];
- if (gradtrue) {
- gradtrue = sqrtf(gradtrue);
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
- gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
- gradmax = fabs(gradmax);
- gradmin = fabs(gradmin);
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.v[2][0] *= alpha;
- p->primitives.gradients.v[2][1] *= alpha;
- p->primitives.gradients.v[2][2] *= alpha;
- }
-
- gradtrue = gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2];
- if (gradtrue) {
- gradtrue = sqrtf(gradtrue);
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.P[1] - p->primitives.P;
- gradmin = p->primitives.P - p->primitives.limiter.P[0];
- gradmax = fabs(gradmax);
- gradmin = fabs(gradmin);
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.P[0] *= alpha;
- p->primitives.gradients.P[1] *= alpha;
- p->primitives.gradients.P[2] *= alpha;
- }
-#endif // SLOPE_LIMITER
-
-#endif // SPH_GRADIENTS
-}
-
-/**
- * @brief Prepare a particle for the fluxes calculation.
- *
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_prepare_fluxes(
- struct part* p, struct xpart* xp) {
-
- /* initialize variables used for timestep calculation */
- p->timestepvars.vmax = 0.0f;
+ hydro_gradients_finalize(p);
}
/**
* @brief Reset acceleration fields of a particle
*
- * Resets all hydro acceleration and time derivative fields in preparation
- * for the sums taking place in the variaous force tasks
+ * This is actually not necessary for GIZMO, since we just set the accelerations
+ * after the flux calculation.
*
- * @param p The particle to act upon
+ * @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
struct part* p) {
@@ -350,11 +222,20 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
}
/**
- * @brief Converts hydro quantity of a particle
+ * @brief Converts the hydrodynamic variables from the initial condition file to
+ * conserved variables that can be used during the integration
*
- * Requires the volume to be known
+ * Requires the volume to be known.
*
- * @param p The particle to act upon
+ * The initial condition file contains a mixture of primitive and conserved
+ * variables. Mass is a conserved variable, and we just copy the particle
+ * mass into the corresponding conserved quantity. We need the volume to
+ * also derive a density, which is then used to convert the internal energy
+ * to a pressure. However, we do not actually use these variables anymore.
+ * We do need to initialize the linear momentum, based on the mass and the
+ * velocity of the particle.
+ *
+ * @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p) {
@@ -362,9 +243,6 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
float volume;
float m;
float momentum[3];
-#ifndef THERMAL_ENERGY
- float momentum2;
-#endif
volume = p->geometry.volume;
p->conserved.mass = m = p->mass;
@@ -376,19 +254,35 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
p->conserved.momentum[0] = momentum[0] = m * p->primitives.v[0];
p->conserved.momentum[1] = momentum[1] = m * p->primitives.v[1];
p->conserved.momentum[2] = momentum[2] = m * p->primitives.v[2];
-#ifndef THERMAL_ENERGY
- momentum2 = momentum[0] * momentum[0] + momentum[1] * momentum[1] +
- momentum[2] * momentum[2];
- p->conserved.energy =
- p->primitives.P / hydro_gamma_minus_one * volume + 0.5 * momentum2 / m;
-#else
p->conserved.energy = p->primitives.P / hydro_gamma_minus_one * volume;
-#endif
}
+/**
+ * @brief Extra operations to be done during the drift
+ *
+ * Not used for GIZMO.
+ *
+ * @param p Particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param t0 Integer start time of the drift interval.
+ * @param t1 Integer end time of the drift interval.
+ * @param timeBase Conversion factor between integer and physical time.
+ */
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {}
+/**
+ * @brief Set the particle acceleration after the flux loop
+ *
+ * We use the new conserved variables to calculate the new velocity of the
+ * particle, and use that to derive the change of the velocity over the particle
+ * time step.
+ *
+ * If the particle time step is zero, we set the accelerations to zero. This
+ * should only happen at the start of the simulation.
+ *
+ * @param p Particle to act upon.
+ */
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {
@@ -414,17 +308,25 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
}
}
+/**
+ * @brief Extra operations done during the kick
+ *
+ * Not used for GIZMO.
+ *
+ * @param p Particle to act upon.
+ * @param xp Extended particle data to act upon.
+ * @param dt Physical time step.
+ * @param half_dt Half the physical time step.
+ */
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part* p, struct xpart* xp, float dt, float half_dt) {
-
- /* Nothing needs to be done in this case */
-}
+ struct part* p, struct xpart* xp, float dt, float half_dt) {}
/**
* @brief Returns the internal energy of a particle
*
- * @param p The particle of interest
- * @param dt Time since the last kick
+ * @param p The particle of interest.
+ * @param dt Time since the last kick.
+ * @return Internal energy of the particle.
*/
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part* restrict p, float dt) {
@@ -435,8 +337,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
/**
* @brief Returns the entropy of a particle
*
- * @param p The particle of interest
- * @param dt Time since the last kick
+ * @param p The particle of interest.
+ * @param dt Time since the last kick.
+ * @return Entropy of the particle.
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part* restrict p, float dt) {
@@ -447,8 +350,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
/**
* @brief Returns the sound speed of a particle
*
- * @param p The particle of interest
- * @param dt Time since the last kick
+ * @param p The particle of interest.
+ * @param dt Time since the last kick.
+ * @param Sound speed of the particle.
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part* restrict p, float dt) {
@@ -461,6 +365,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
*
* @param p The particle of interest
* @param dt Time since the last kick
+ * @param Pressure of the particle.
*/
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part* restrict p, float dt) {
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index ff1e39d68bc2f9b02b80743e497c270fd97337fa..5a97b014288e7a1837dd78f47a7228ff7c9e7038 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -20,43 +20,70 @@
#ifndef SWIFT_HYDRO_GRADIENTS_H
#define SWIFT_HYDRO_GRADIENTS_H
-#define SPH_GRADIENTS
+//#define SPH_GRADIENTS
+#define GIZMO_GRADIENTS
#include "hydro_slope_limiters.h"
#if defined(SPH_GRADIENTS)
+
+#define HYDRO_GRADIENT_IMPLEMENTATION "SPH gradients (Price 2012)"
#include "hydro_gradients_sph.h"
+
#elif defined(GIZMO_GRADIENTS)
+
+#define HYDRO_GRADIENT_IMPLEMENTATION "GIZMO gradients (Hopkins 2015)"
#include "hydro_gradients_gizmo.h"
+
#else
/* No gradients. Perfectly acceptable, but we have to provide empty functions */
+#define HYDRO_GRADIENT_IMPLEMENTATION "No gradients (first order scheme)"
/**
- * @brief Initialize variables before the density loop
+ * @brief Initialize gradient variables
+ *
+ * @param p Particle.
*/
-__attribute__((always_inline)) INLINE static void
-hydro_gradients_init_density_loop(struct part *p) {}
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+ struct part* p) {}
/**
- * @brief Gradient calculations done during the density loop
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
*/
-__attribute__((always_inline)) INLINE static void hydro_gradients_density_loop(
- struct part *pi, struct part *pj, float wi_dx, float wj_dx, float *dx,
- float r, int mode) {}
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, float* dx, float hi, float hj, struct part* pi, struct part* pj) {
+}
/**
- * @brief Calculations done before the force loop
+ * @brief Gradient calculations done during the neighbour loop: non-symmetric
+ * version
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
*/
-__attribute__((always_inline)) INLINE static void
-hydro_gradients_prepare_force_loop(struct part *p, float ih2, float volume) {}
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, float* dx, float hi, float hj, struct part* pi, struct part* pj) {
+}
/**
- * @brief Gradient calculations done during the gradient loop
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
*/
-__attribute__((always_inline)) INLINE static void hydro_gradients_gradient_loop(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
- int mode) {}
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+ struct part* p) {}
#endif
diff --git a/src/hydro/Gizmo/hydro_gradients_gizmo.h b/src/hydro/Gizmo/hydro_gradients_gizmo.h
index 97bbcf98d76eca730ac49a5dcfbfae7b2f336c2c..aa6e4406b94e7a5cafcd0ca556162476003477de 100644
--- a/src/hydro/Gizmo/hydro_gradients_gizmo.h
+++ b/src/hydro/Gizmo/hydro_gradients_gizmo.h
@@ -18,20 +18,12 @@
******************************************************************************/
/**
- * @brief Initialize variables before the density loop
- */
-__attribute__((always_inline)) INLINE static void
-hydro_gradients_init_density_loop(struct part *p) {}
-
-/**
- * @brief Gradient calculations done during the density loop
+ * @brief Initialize gradient variables
+ *
+ * @param p Particle.
*/
-__attribute__((always_inline)) INLINE static void hydro_gradients_density_loop(
- struct part *pi, struct part *pj, float wi_dx, float wj_dx, float *dx,
- float r, int mode) {}
-
-__attribute__((always_inline)) INLINE static void
-hydro_gradients_prepare_force_loop(struct part *p, float ih2, volume) {
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+ struct part *p) {
p->primitives.gradients.rho[0] = 0.0f;
p->primitives.gradients.rho[1] = 0.0f;
@@ -53,26 +45,21 @@ hydro_gradients_prepare_force_loop(struct part *p, float ih2, volume) {
p->primitives.gradients.P[1] = 0.0f;
p->primitives.gradients.P[2] = 0.0f;
- p->primitives.limiter.rho[0] = FLT_MAX;
- p->primitives.limiter.rho[1] = -FLT_MAX;
- p->primitives.limiter.v[0][0] = FLT_MAX;
- p->primitives.limiter.v[0][1] = -FLT_MAX;
- p->primitives.limiter.v[1][0] = FLT_MAX;
- p->primitives.limiter.v[1][1] = -FLT_MAX;
- p->primitives.limiter.v[2][0] = FLT_MAX;
- p->primitives.limiter.v[2][1] = -FLT_MAX;
- p->primitives.limiter.P[0] = FLT_MAX;
- p->primitives.limiter.P[1] = -FLT_MAX;
-
- p->primitives.limiter.maxr = -FLT_MAX;
+ hydro_slope_limit_cell_init(p);
}
/**
- * @brief Gradient calculations done during the gradient loop
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
*/
-__attribute__((always_inline)) INLINE static void hydro_gradients_gradient_loop(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
- int mode) {
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
float r = sqrtf(r2);
float xi, xj;
@@ -157,116 +144,198 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_gradient_loop(
(Wi[4] - Wj[4]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- /* basic slope limiter: collect the maximal and the minimal value for the
- * primitive variables among the ngbs */
- pi->primitives.limiter.rho[0] =
- fmin(pj->primitives.rho, pi->primitives.limiter.rho[0]);
- pi->primitives.limiter.rho[1] =
- fmax(pj->primitives.rho, pi->primitives.limiter.rho[1]);
-
- pi->primitives.limiter.v[0][0] =
- fmin(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
- pi->primitives.limiter.v[0][1] =
- fmax(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
- pi->primitives.limiter.v[1][0] =
- fmin(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
- pi->primitives.limiter.v[1][1] =
- fmax(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
- pi->primitives.limiter.v[2][0] =
- fmin(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
- pi->primitives.limiter.v[2][1] =
- fmax(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
-
- pi->primitives.limiter.P[0] =
- fmin(pj->primitives.P, pi->primitives.limiter.P[0]);
- pi->primitives.limiter.P[1] =
- fmax(pj->primitives.P, pi->primitives.limiter.P[1]);
-
- pi->primitives.limiter.maxr = fmax(r, pi->primitives.limiter.maxr);
-
- if (mode == 1) {
- /* Compute kernel of pj. */
- hj_inv = 1.0 / hj;
- xj = r * hj_inv;
- kernel_deval(xj, &wj, &wj_dx);
-
- /* Compute gradients for pj */
- /* there is no sign difference w.r.t. eqn. (6) because dx is now what we
- * want
- * it to be */
- pj->primitives.gradients.rho[0] +=
- (Wi[0] - Wj[0]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.rho[1] +=
- (Wi[0] - Wj[0]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.rho[2] +=
- (Wi[0] - Wj[0]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
- pj->primitives.gradients.v[0][0] +=
- (Wi[1] - Wj[1]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.v[0][1] +=
- (Wi[1] - Wj[1]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.v[0][2] +=
- (Wi[1] - Wj[1]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->primitives.gradients.v[1][0] +=
- (Wi[2] - Wj[2]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.v[1][1] +=
- (Wi[2] - Wj[2]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.v[1][2] +=
- (Wi[2] - Wj[2]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->primitives.gradients.v[2][0] +=
- (Wi[3] - Wj[3]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.v[2][1] +=
- (Wi[3] - Wj[3]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.v[2][2] +=
- (Wi[3] - Wj[3]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
- pj->primitives.gradients.P[0] +=
- (Wi[4] - Wj[4]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.P[1] +=
- (Wi[4] - Wj[4]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.P[2] +=
- (Wi[4] - Wj[4]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
- /* basic slope limiter: collect the maximal and the minimal value for the
- * primitive variables among the ngbs */
- pj->primitives.limiter.rho[0] =
- fmin(pi->primitives.rho, pj->primitives.limiter.rho[0]);
- pj->primitives.limiter.rho[1] =
- fmax(pi->primitives.rho, pj->primitives.limiter.rho[1]);
-
- pj->primitives.limiter.v[0][0] =
- fmin(pi->primitives.v[0], pj->primitives.limiter.v[0][0]);
- pj->primitives.limiter.v[0][1] =
- fmax(pi->primitives.v[0], pj->primitives.limiter.v[0][1]);
- pj->primitives.limiter.v[1][0] =
- fmin(pi->primitives.v[1], pj->primitives.limiter.v[1][0]);
- pj->primitives.limiter.v[1][1] =
- fmax(pi->primitives.v[1], pj->primitives.limiter.v[1][1]);
- pj->primitives.limiter.v[2][0] =
- fmin(pi->primitives.v[2], pj->primitives.limiter.v[2][0]);
- pj->primitives.limiter.v[2][1] =
- fmax(pi->primitives.v[2], pj->primitives.limiter.v[2][1]);
-
- pj->primitives.limiter.P[0] =
- fmin(pi->primitives.P, pj->primitives.limiter.P[0]);
- pj->primitives.limiter.P[1] =
- fmax(pi->primitives.P, pj->primitives.limiter.P[1]);
-
- pj->primitives.limiter.maxr = fmax(r, pj->primitives.limiter.maxr);
+ hydro_slope_limit_cell_collect(pi, pj, r);
+
+ /* Compute kernel of pj. */
+ hj_inv = 1.0 / hj;
+ xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* Compute gradients for pj */
+ /* there is no sign difference w.r.t. eqn. (6) because dx is now what we
+ * want
+ * it to be */
+ pj->primitives.gradients.rho[0] +=
+ (Wi[0] - Wj[0]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.rho[1] +=
+ (Wi[0] - Wj[0]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.rho[2] +=
+ (Wi[0] - Wj[0]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ pj->primitives.gradients.v[0][0] +=
+ (Wi[1] - Wj[1]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.v[0][1] +=
+ (Wi[1] - Wj[1]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.v[0][2] +=
+ (Wi[1] - Wj[1]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+ pj->primitives.gradients.v[1][0] +=
+ (Wi[2] - Wj[2]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.v[1][1] +=
+ (Wi[2] - Wj[2]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.v[1][2] +=
+ (Wi[2] - Wj[2]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+ pj->primitives.gradients.v[2][0] +=
+ (Wi[3] - Wj[3]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.v[2][1] +=
+ (Wi[3] - Wj[3]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.v[2][2] +=
+ (Wi[3] - Wj[3]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ pj->primitives.gradients.P[0] +=
+ (Wi[4] - Wj[4]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.P[1] +=
+ (Wi[4] - Wj[4]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.P[2] +=
+ (Wi[4] - Wj[4]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ hydro_slope_limit_cell_collect(pj, pi, r);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
+ struct part *pi, struct part *pj) {
+
+ float r = sqrtf(r2);
+ float xi;
+ float hi_inv;
+ float wi, wi_dx;
+ int k, l;
+ float Bi[3][3];
+ float Wi[5], Wj[5];
+
+ /* Initialize local variables */
+ for (k = 0; k < 3; k++) {
+ for (l = 0; l < 3; l++) {
+ Bi[k][l] = pi->geometry.matrix_E[k][l];
+ }
}
+ Wi[0] = pi->primitives.rho;
+ Wi[1] = pi->primitives.v[0];
+ Wi[2] = pi->primitives.v[1];
+ Wi[3] = pi->primitives.v[2];
+ Wi[4] = pi->primitives.P;
+ Wj[0] = pj->primitives.rho;
+ Wj[1] = pj->primitives.v[0];
+ Wj[2] = pj->primitives.v[1];
+ Wj[3] = pj->primitives.v[2];
+ Wj[4] = pj->primitives.P;
+
+ /* Compute kernel of pi. */
+ hi_inv = 1.0 / hi;
+ xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* Compute gradients for pi */
+ /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+ * definition of dx */
+ pi->primitives.gradients.rho[0] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.rho[1] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.rho[2] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.v[0][0] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[0][1] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[0][2] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[1][0] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[1][1] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[1][2] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[2][0] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[2][1] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[2][2] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.P[0] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.P[1] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.P[2] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+}
+
+/**
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+ struct part *p) {
+
+ float h, ih;
+
+ /* add kernel normalization to gradients */
+ h = p->h;
+ ih = 1.0f / h;
+ const float ihdim = pow_dimension(ih);
+
+ p->primitives.gradients.rho[0] *= ihdim;
+ p->primitives.gradients.rho[1] *= ihdim;
+ p->primitives.gradients.rho[2] *= ihdim;
+
+ p->primitives.gradients.v[0][0] *= ihdim;
+ p->primitives.gradients.v[0][1] *= ihdim;
+ p->primitives.gradients.v[0][2] *= ihdim;
+ p->primitives.gradients.v[1][0] *= ihdim;
+ p->primitives.gradients.v[1][1] *= ihdim;
+ p->primitives.gradients.v[1][2] *= ihdim;
+ p->primitives.gradients.v[2][0] *= ihdim;
+ p->primitives.gradients.v[2][1] *= ihdim;
+ p->primitives.gradients.v[2][2] *= ihdim;
+
+ p->primitives.gradients.P[0] *= ihdim;
+ p->primitives.gradients.P[1] *= ihdim;
+ p->primitives.gradients.P[2] *= ihdim;
+
+ hydro_slope_limit_cell(p);
}
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/Gizmo/hydro_gradients_sph.h
index f8dc7ea00458dd11c45aa0a340e6139376a286ae..f635faecea549f7da280ade9b944021a5e4aeb4c 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/src/hydro/Gizmo/hydro_gradients_sph.h
@@ -18,10 +18,12 @@
******************************************************************************/
/**
- * @brief Initialize variables before the density loop
+ * @brief Initialize gradient variables
+ *
+ * @param p Particle.
*/
-__attribute__((always_inline)) INLINE static void
-hydro_gradients_init_density_loop(struct part *p) {
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+ struct part *p) {
p->primitives.gradients.rho[0] = 0.0f;
p->primitives.gradients.rho[1] = 0.0f;
@@ -46,11 +48,25 @@ hydro_gradients_init_density_loop(struct part *p) {
}
/**
- * @brief Gradient calculations done during the density loop
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
*/
-__attribute__((always_inline)) INLINE static void hydro_gradients_density_loop(
- struct part *pi, struct part *pj, float wi_dx, float wj_dx, float *dx,
- float r, int mode) {
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float wi, wi_dx, xi, hi_inv;
+ float wj, wj_dx, xj, hj_inv;
+ float r = sqrtf(r2);
+
+ hi_inv = 1.0f / hi;
+ xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
/* very basic gradient estimate */
pi->primitives.gradients.rho[0] -=
@@ -90,54 +106,123 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_density_loop(
hydro_slope_limit_cell_collect(pi, pj, r);
- if (mode == 1) {
- /* signs are the same as before, since we swap i and j twice */
- pj->primitives.gradients.rho[0] -=
- wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
- pj->primitives.gradients.rho[1] -=
- wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
- pj->primitives.gradients.rho[2] -=
- wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
-
- pj->primitives.gradients.v[0][0] -=
- wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
- pj->primitives.gradients.v[0][1] -=
- wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
- pj->primitives.gradients.v[0][2] -=
- wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
-
- pj->primitives.gradients.v[1][0] -=
- wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
- pj->primitives.gradients.v[1][1] -=
- wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
- pj->primitives.gradients.v[1][2] -=
- wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
- pj->primitives.gradients.v[2][0] -=
- wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
- pj->primitives.gradients.v[2][1] -=
- wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
- pj->primitives.gradients.v[2][2] -=
- wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
-
- pj->primitives.gradients.P[0] -=
- wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
- pj->primitives.gradients.P[1] -=
- wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
- pj->primitives.gradients.P[2] -=
- wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
- hydro_slope_limit_cell_collect(pj, pi, r);
- }
+ hj_inv = 1.0f / hj;
+ xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* signs are the same as before, since we swap i and j twice */
+ pj->primitives.gradients.rho[0] -=
+ wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pj->primitives.gradients.rho[1] -=
+ wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pj->primitives.gradients.rho[2] -=
+ wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+ pj->primitives.gradients.v[0][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pj->primitives.gradients.v[0][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pj->primitives.gradients.v[0][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+
+ pj->primitives.gradients.v[1][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pj->primitives.gradients.v[1][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pj->primitives.gradients.v[1][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pj->primitives.gradients.v[2][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pj->primitives.gradients.v[2][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pj->primitives.gradients.v[2][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+ pj->primitives.gradients.P[0] -=
+ wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+ pj->primitives.gradients.P[1] -=
+ wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+ pj->primitives.gradients.P[2] -=
+ wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+ hydro_slope_limit_cell_collect(pj, pi, r);
}
/**
- * @brief Calculations done before the force loop
+ * @brief Gradient calculations done during the neighbour loop: non-symmetric
+ * version
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
*/
__attribute__((always_inline)) INLINE static void
-hydro_gradients_prepare_force_loop(struct part *p, float ih, float volume) {
+hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
+ struct part *pi, struct part *pj) {
+
+ float wi, wi_dx, xi, hi_inv;
+ float r = sqrtf(r2);
+
+ hi_inv = 1.0f / hi;
+ xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* very basic gradient estimate */
+ pi->primitives.gradients.rho[0] -=
+ wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[1] -=
+ wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[2] -=
+ wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+ pi->primitives.gradients.v[0][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[1][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+
+ pi->primitives.gradients.v[2][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+ pi->primitives.gradients.P[0] -=
+ wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[1] -=
+ wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[2] -=
+ wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+}
+
+/**
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+ struct part *p) {
+
+ const float h = p->h;
+ const float ih = 1.0f / h;
const float ihdimp1 = pow_dimension_plus_one(ih);
+ float volume = p->geometry.volume;
+
/* finalize gradients by multiplying with volume */
p->primitives.gradients.rho[0] *= ihdimp1 * volume;
p->primitives.gradients.rho[1] *= ihdimp1 * volume;
@@ -161,10 +246,3 @@ hydro_gradients_prepare_force_loop(struct part *p, float ih, float volume) {
hydro_slope_limit_cell(p);
}
-
-/**
- * @brief Gradient calculations done during the gradient loop
- */
-__attribute__((always_inline)) INLINE static void hydro_gradients_gradient_loop(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
- int mode) {}
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index f7f48997194feb27ff2d19a20b5ca01b4d04fc88..1d5fdf15322767432c87db8d35f63b99695961bb 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -23,11 +23,23 @@
#include "hydro_gradients.h"
#include "riemann.h"
-//#define USE_GRADIENTS
-//#define PER_FACE_LIMITER
-/* #define PRINT_ID 0 */
-
-/* this corresponds to task_subtype_hydro_loop1 */
+/**
+ * @brief Calculate the volume interaction between particle i and particle j
+ *
+ * The volume is in essence the same as the weighted number of neighbours in a
+ * classical SPH density calculation.
+ *
+ * We also calculate the components of the matrix E, which is used for second
+ * order accurate gradient calculations and for the calculation of the interface
+ * surface areas.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
__attribute__((always_inline)) INLINE static void runner_iact_density(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
@@ -62,11 +74,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->geometry.volume += wj;
for (k = 0; k < 3; k++)
for (l = 0; l < 3; l++) pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
-
- hydro_gradients_density_loop(pi, pj, wi_dx, wj_dx, dx, r, 1);
}
-/* this corresponds to task_subtype_hydro_loop1 */
+/**
+ * @brief Calculate the volume interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * The volume is in essence the same as the weighted number of neighbours in a
+ * classical SPH density calculation.
+ *
+ * We also calculate the components of the matrix E, which is used for second
+ * order accurate gradient calculations and for the calculation of the interface
+ * surface areas.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
@@ -90,22 +117,72 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->geometry.volume += wi;
for (k = 0; k < 3; k++)
for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
-
- hydro_gradients_density_loop(pi, pj, wi_dx, 0.0f, dx, r, 0);
}
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j
+ *
+ * This method wraps around hydro_gradients_collect, which can be an empty
+ * method, in which case no gradients are used.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
__attribute__((always_inline)) INLINE static void runner_iact_gradient(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
- hydro_gradients_gradient_loop(r2, dx, hi, hj, pi, pj, 1);
+ hydro_gradients_collect(r2, dx, hi, hj, pi, pj);
}
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * This method wraps around hydro_gradients_nonsym_collect, which can be an
+ * empty method, in which case no gradients are used.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
- hydro_gradients_gradient_loop(r2, dx, hi, hj, pi, pj, 0);
+ hydro_gradients_nonsym_collect(r2, dx, hi, hj, pi, pj);
}
+/**
+ * @brief Common part of the flux calculation between particle i and j
+ *
+ * Since the only difference between the symmetric and non-symmetric version
+ * of the flux calculation is in the update of the conserved variables at the
+ * very end (which is not done for particle j if mode is 0 and particle j is
+ * active), both runner_iact_force and runner_iact_nonsym_force call this
+ * method, with an appropriate mode.
+ *
+ * This method calculates the surface area of the interface between particle i
+ * and particle j, as well as the interface position and velocity. These are
+ * then used to reconstruct and predict the primitive variables, which are then
+ * fed to a Riemann solver that calculates a flux. This flux is used to update
+ * the conserved variables of particle i or both particles.
+ *
+ * This method also calculates the maximal velocity used to calculate the time
+ * step.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
int mode) {
@@ -327,21 +404,45 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
}
}
-/* this corresponds to task_subtype_fluxes */
+/**
+ * @brief Flux calculation between particle i and particle j
+ *
+ * This method calls runner_iact_fluxes_common with mode 1.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
__attribute__((always_inline)) INLINE static void runner_iact_force(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1);
}
-/* this corresponds to task_subtype_fluxes */
+/**
+ * @brief Flux calculation between particle i and particle j: non-symmetric
+ * version
+ *
+ * This method calls runner_iact_fluxes_common with mode 0.
+ *
+ * @param r2 Squared distance between particle i and particle j.
+ * @param dx Distance vector between the particles (dx = pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0);
}
-/* PLACEHOLDERS */
+//// EMPTY VECTORIZED VERSIONS (gradients methods are missing...)
+
__attribute__((always_inline)) INLINE static void runner_iact_vec_density(
float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
struct part **pj) {}
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index c94451aad943e7603f3393a400fe5864e624d16b..6c4f64b2576c76e81fc7d432a4d4a6f6ecb5189c 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -18,6 +18,8 @@
******************************************************************************/
#include "adiabatic_index.h"
+#include "hydro_gradients.h"
+#include "hydro_slope_limiters.h"
#include "io_properties.h"
/**
@@ -106,7 +108,17 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
* @brief Writes the current model of SPH to the file
* @param h_grpsph The HDF5 group in which to write
*/
-void writeSPHflavour(hid_t h_grpsph) {}
+void writeSPHflavour(hid_t h_grpsph) {
+ /* Gradient information */
+ writeAttribute_s(h_grpsph, "Gradient reconstruction model",
+ HYDRO_GRADIENT_IMPLEMENTATION);
+
+ /* Slope limiter information */
+ writeAttribute_s(h_grpsph, "Cell wide slope limiter model",
+ HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION);
+ writeAttribute_s(h_grpsph, "Piecewise slope limiter model",
+ HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
+}
/**
* @brief Are we writing entropy in the internal energy field ?
diff --git a/src/hydro/Gizmo/hydro_slope_limiters.h b/src/hydro/Gizmo/hydro_slope_limiters.h
index 357ba2a84eb0a0c340906ddfa9d20e6a4153a405..e1f841411a65fc5fce78ecab07e272c9c6df475b 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters.h
@@ -23,10 +23,32 @@
#define PER_FACE_LIMITER
#define CELL_WIDE_LIMITER
+#include "dimension.h"
+#include "kernel_hydro.h"
+
#ifdef PER_FACE_LIMITER
+
+#define HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION \
+ "GIZMO piecewise slope limiter (Hopkins 2015)"
#include "hydro_slope_limiters_face.h"
+
#else
+#define HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION "No piecewise slope limiter"
+
+/**
+ * @brief Slope limit the slopes at the interface between two particles
+ *
+ * @param Wi Hydrodynamic variables of particle i.
+ * @param Wj Hydrodynamic variables of particle j.
+ * @param dWi Difference between the hydrodynamic variables of particle i at the
+ * position of particle i and at the interface position.
+ * @param dWj Difference between the hydrodynamic variables of particle j at the
+ * position of particle j and at the interface position.
+ * @param xij_i Relative position vector of the interface w.r.t. particle i.
+ * @param xij_j Relative position vector of the interface w.r.t. partilce j.
+ * @param r Distance between particle i and particle j.
+ */
__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
float r) {}
@@ -34,15 +56,39 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
#endif
#ifdef CELL_WIDE_LIMITER
+
+#define HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION \
+ "Cell wide slope limiter (Springel 2010)"
#include "hydro_slope_limiters_cell.h"
+
#else
+#define HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION "No cell wide slope limiter"
+
+/**
+ * @brief Initialize variables for the cell wide slope limiter
+ *
+ * @param p Particle.
+ */
__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
struct part *p) {}
+/**
+ * @brief Collect information for the cell wide slope limiter during the
+ * neighbour loop
+ *
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param r Distance between particle i and particle j.
+ */
__attribute__((always_inline)) INLINE static void
hydro_slope_limit_cell_collect(struct part *pi, struct part *pj, float r) {}
+/**
+ * @brief Slope limit cell gradients
+ *
+ * @param p Particle.
+ */
__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
struct part *p) {}
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_cell.h b/src/hydro/Gizmo/hydro_slope_limiters_cell.h
index 32df7c53625eea90f1c7f97fcf36c6c2c3f9be4e..aa99b43721f669f47a7888a5da0b1933ca1ebd62 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_cell.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters_cell.h
@@ -17,6 +17,13 @@
*
******************************************************************************/
+#include <float.h>
+
+/**
+ * @brief Initialize variables for the cell wide slope limiter
+ *
+ * @param p Particle.
+ */
__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
struct part* p) {
@@ -34,6 +41,14 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
p->primitives.limiter.maxr = -FLT_MAX;
}
+/**
+ * @brief Collect information for the cell wide slope limiter during the
+ * neighbour loop
+ *
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param r Distance between particle i and particle j.
+ */
__attribute__((always_inline)) INLINE static void
hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
@@ -65,6 +80,11 @@ hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
pi->primitives.limiter.maxr = fmax(r, pi->primitives.limiter.maxr);
}
+/**
+ * @brief Slope limit cell gradients
+ *
+ * @param p Particle.
+ */
__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
struct part* p) {
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_face.h b/src/hydro/Gizmo/hydro_slope_limiters_face.h
index 7ef9319877fbb60309afdf7b92acf296f9e117b4..6745638369fc71dcffc1c480c8d5872f70cfaf5f 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_face.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters_face.h
@@ -17,72 +17,105 @@
*
******************************************************************************/
-__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
- float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
- float r) {
+/**
+ * @brief Slope limit a single quantity at the interface
+ *
+ * @param phi_i Value of the quantity at the particle position.
+ * @param phi_j Value of the quantity at the neighbouring particle position.
+ * @param phi_mid0 Extrapolated value of the quantity at the interface position.
+ * @param xij_norm Distance between the particle position and the interface
+ * position.
+ * @param r Distance between the particle and its neighbour.
+ * @return The slope limited difference between the quantity at the particle
+ * position and the quantity at the interface position.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
+ float xij_norm, float r) {
- float xij_i_norm;
- float phi_i, phi_j;
- float delta1, delta2;
- float phiminus, phiplus;
- float phimin, phimax;
- float phibar;
- /* free parameters, values from Hopkins */
- float psi1 = 0.5, psi2 = 0.25;
- float phi_mid0, phi_mid;
- int k;
+ float delta1, delta2, phimin, phimax, phibar, phiplus, phiminus, phi_mid;
+ const float psi1 = 0.5f;
+ const float psi2 = 0.25f;
- for (k = 0; k < 10; k++) {
- if (k < 5) {
- phi_i = Wi[k];
- phi_j = Wj[k];
- phi_mid0 = Wi[k] + dWi[k];
- xij_i_norm = sqrtf(xij_i[0] * xij_i[0] + xij_i[1] * xij_i[1] +
- xij_i[2] * xij_i[2]);
- } else {
- phi_i = Wj[k - 5];
- phi_j = Wi[k - 5];
- phi_mid0 = Wj[k - 5] + dWj[k - 5];
- xij_i_norm = sqrtf(xij_j[0] * xij_j[0] + xij_j[1] * xij_j[1] +
- xij_j[2] * xij_j[2]);
- }
+ delta1 = psi1 * fabs(phi_i - phi_j);
+ delta2 = psi2 * fabs(phi_i - phi_j);
- delta1 = psi1 * fabs(phi_i - phi_j);
- delta2 = psi2 * fabs(phi_i - phi_j);
+ phimin = fmin(phi_i, phi_j);
+ phimax = fmax(phi_i, phi_j);
- phimin = fmin(phi_i, phi_j);
- phimax = fmax(phi_i, phi_j);
+ phibar = phi_i + xij_norm / r * (phi_j - phi_i);
- phibar = phi_i + xij_i_norm / r * (phi_j - phi_i);
+ /* if sign(phimax+delta1) == sign(phimax) */
+ if ((phimax + delta1) * phimax > 0.0f) {
+ phiplus = phimax + delta1;
+ } else {
+ phiplus = phimax / (1.0f + delta1 / fabs(phimax));
+ }
- /* if sign(phimax+delta1) == sign(phimax) */
- if ((phimax + delta1) * phimax > 0.0f) {
- phiplus = phimax + delta1;
- } else {
- phiplus = phimax / (1.0f + delta1 / fabs(phimax));
- }
+ /* if sign(phimin-delta1) == sign(phimin) */
+ if ((phimin - delta1) * phimin > 0.0f) {
+ phiminus = phimin - delta1;
+ } else {
+ phiminus = phimin / (1.0f + delta1 / fabs(phimin));
+ }
- /* if sign(phimin-delta1) == sign(phimin) */
- if ((phimin - delta1) * phimin > 0.0f) {
- phiminus = phimin - delta1;
+ if (phi_i == phi_j) {
+ phi_mid = phi_i;
+ } else {
+ if (phi_i < phi_j) {
+ phi_mid = fmax(phiminus, fmin(phibar + delta2, phi_mid0));
} else {
- phiminus = phimin / (1.0f + delta1 / fabs(phimin));
+ phi_mid = fmin(phiplus, fmax(phibar - delta2, phi_mid0));
}
+ }
- if (phi_i == phi_j) {
- phi_mid = phi_i;
- } else {
- if (phi_i < phi_j) {
- phi_mid = fmax(phiminus, fmin(phibar + delta2, phi_mid0));
- } else {
- phi_mid = fmin(phiplus, fmax(phibar - delta2, phi_mid0));
- }
- }
+ return phi_mid - phi_i;
+}
- if (k < 5) {
- dWi[k] = phi_mid - phi_i;
- } else {
- dWj[k - 5] = phi_mid - phi_i;
- }
- }
+/**
+ * @brief Slope limit the slopes at the interface between two particles
+ *
+ * @param Wi Hydrodynamic variables of particle i.
+ * @param Wj Hydrodynamic variables of particle j.
+ * @param dWi Difference between the hydrodynamic variables of particle i at the
+ * position of particle i and at the interface position.
+ * @param dWj Difference between the hydrodynamic variables of particle j at the
+ * position of particle j and at the interface position.
+ * @param xij_i Relative position vector of the interface w.r.t. particle i.
+ * @param xij_j Relative position vector of the interface w.r.t. partilce j.
+ * @param r Distance between particle i and particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
+ float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
+ float r) {
+
+ float xij_i_norm, xij_j_norm;
+
+ xij_i_norm =
+ sqrtf(xij_i[0] * xij_i[0] + xij_i[1] * xij_i[1] + xij_i[2] * xij_i[2]);
+
+ xij_j_norm =
+ sqrtf(xij_j[0] * xij_j[0] + xij_j[1] * xij_j[1] + xij_j[2] * xij_j[2]);
+
+ dWi[0] = hydro_slope_limit_face_quantity(Wi[0], Wj[0], Wi[0] + dWi[0],
+ xij_i_norm, r);
+ dWi[1] = hydro_slope_limit_face_quantity(Wi[1], Wj[1], Wi[1] + dWi[1],
+ xij_i_norm, r);
+ dWi[2] = hydro_slope_limit_face_quantity(Wi[2], Wj[2], Wi[2] + dWi[2],
+ xij_i_norm, r);
+ dWi[3] = hydro_slope_limit_face_quantity(Wi[3], Wj[3], Wi[3] + dWi[3],
+ xij_i_norm, r);
+ dWi[4] = hydro_slope_limit_face_quantity(Wi[4], Wj[4], Wi[4] + dWi[4],
+ xij_i_norm, r);
+
+ dWj[0] = hydro_slope_limit_face_quantity(Wj[0], Wi[0], Wj[0] + dWj[0],
+ xij_j_norm, r);
+ dWj[1] = hydro_slope_limit_face_quantity(Wj[1], Wi[1], Wj[1] + dWj[1],
+ xij_j_norm, r);
+ dWj[2] = hydro_slope_limit_face_quantity(Wj[2], Wi[2], Wj[2] + dWj[2],
+ xij_j_norm, r);
+ dWj[3] = hydro_slope_limit_face_quantity(Wj[3], Wi[3], Wj[3] + dWj[3],
+ xij_j_norm, r);
+ dWj[4] = hydro_slope_limit_face_quantity(Wj[4], Wi[4], Wj[4] + dWj[4],
+ xij_j_norm, r);
}
diff --git a/src/runner.c b/src/runner.c
index f4ac9f825df807416c56534f9167a777f796b853..d8c42920675aa662106a23d754dd11e790844914 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -428,7 +428,39 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_init);
}
-void runner_do_extra_ghost(struct runner *r, struct cell *c) {}
+/**
+ * @brief Intermediate task after the gradient loop that does final operations
+ * on the gradient quantities and optionally slope limits the gradients
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ */
+void runner_do_extra_ghost(struct runner *r, struct cell *c) {
+ struct part *restrict parts = c->parts;
+ const int count = c->count;
+ const int ti_current = r->e->ti_current;
+
+ /* Recurse? */
+ if (c->split) {
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) runner_do_extra_ghost(r, c->progeny[k]);
+ return;
+ } else {
+
+ /* Loop over the parts in this cell. */
+ for (int i = 0; i < count; i++) {
+
+ /* Get a direct pointer on the part. */
+ struct part *restrict p = &parts[i];
+
+ if (p->ti_end <= ti_current) {
+
+ /* Get ready for a force calculation */
+ hydro_end_gradient(p);
+ }
+ }
+ }
+}
/**
* @brief Intermediate task after the density to check that the smoothing