diff --git a/src/Makefile.am b/src/Makefile.am
index eec459727ba4951777a6b1d9c214bdd13986a1e3..043c084df6ef9a966450f0dd172bb44e79355263 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -137,7 +137,11 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/Gizmo/hydro_io.h hydro/Gizmo/hydro_debug.h \
hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h \
hydro/Gizmo/hydro_part.h \
- hydro/Gizmo/MFV/hydro.h hydro/Gizmo/MFV/hydro_iact.h \
+ hydro/Gizmo/hydro_gradients.h \
+ hydro/Gizmo/hydro_getters.h \
+ hydro/Gizmo/hydro_setters.h \
+ hydro/Gizmo/hydro_flux.h \
+ hydro/Gizmo/MFV/hydro.h \
hydro/Gizmo/MFV/hydro_io.h hydro/Gizmo/MFV/hydro_debug.h \
hydro/Gizmo/MFV/hydro_part.h \
hydro/Gizmo/MFV/hydro_gradients_gizmo.h \
@@ -149,7 +153,10 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/Gizmo/MFV/hydro_unphysical.h \
hydro/Gizmo/MFV/hydro_velocities.h \
hydro/Gizmo/MFV/hydro_parameters.h \
- hydro/Gizmo/MFM/hydro.h hydro/Gizmo/MFM/hydro_iact.h \
+ hydro/Gizmo/MFV/hydro_getters.h \
+ hydro/Gizmo/MFV/hydro_setters.h \
+ hydro/Gizmo/MFV/hydro_flux.h \
+ hydro/Gizmo/MFM/hydro.h \
hydro/Gizmo/MFM/hydro_io.h hydro/Gizmo/MFM/hydro_debug.h \
hydro/Gizmo/MFM/hydro_part.h \
hydro/Gizmo/MFM/hydro_gradients_gizmo.h \
@@ -160,6 +167,9 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/Gizmo/MFM/hydro_slope_limiters.h \
hydro/Gizmo/MFM/hydro_unphysical.h \
hydro/Gizmo/MFM/hydro_parameters.h \
+ hydro/Gizmo/MFM/hydro_getters.h \
+ hydro/Gizmo/MFM/hydro_setters.h \
+ hydro/Gizmo/MFM/hydro_flux.h \
hydro/Shadowswift/hydro_debug.h \
hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
hydro/Shadowswift/hydro_iact.h \
diff --git a/src/hydro/Gizmo/MFM/hydro_flux.h b/src/hydro/Gizmo/MFM/hydro_flux.h
new file mode 100644
index 0000000000000000000000000000000000000000..9cdc2dbeecd5dbc9303d852bb5c39d1de0999605
--- /dev/null
+++ b/src/hydro/Gizmo/MFM/hydro_flux.h
@@ -0,0 +1,47 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFM_HYDRO_FLUX_H
+#define SWIFT_GIZMO_MFM_HYDRO_FLUX_H
+
+#include "riemann.h"
+
+/**
+ * @brief Compute the flux for the Riemann problem with the given left and right
+ * state, and interface normal, surface area and velocity.
+ *
+ * @param WL Left state variables.
+ * @param WR Right state variables.
+ * @param n_unit Unit vector of the interface.
+ * @param vLR Velocity of the interface.
+ * @param Anorm Surface area of the interface.
+ * @param fluxes Array to store the result in (of size 5 or more).
+ */
+__attribute__((always_inline)) INLINE static void hydro_compute_flux(
+ const float *WL, const float *WR, const float *n_unit, const float *vLR,
+ const float Anorm, float *fluxes) {
+
+ riemann_solve_for_middle_state_flux(WL, WR, n_unit, vLR, fluxes);
+
+ fluxes[1] *= Anorm;
+ fluxes[2] *= Anorm;
+ fluxes[3] *= Anorm;
+ fluxes[4] *= Anorm;
+}
+
+#endif /* SWIFT_GIZMO_MFM_HYDRO_FLUX_H */
diff --git a/src/hydro/Gizmo/MFM/hydro_getters.h b/src/hydro/Gizmo/MFM/hydro_getters.h
new file mode 100644
index 0000000000000000000000000000000000000000..bbf87d048f8c40e55956a35f084b44208e129357
--- /dev/null
+++ b/src/hydro/Gizmo/MFM/hydro_getters.h
@@ -0,0 +1,46 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFM_HYDRO_GETTERS_H
+#define SWIFT_GIZMO_MFM_HYDRO_GETTERS_H
+
+/**
+ * @brief Get a 5-element state vector W containing the primitive hydrodynamic
+ * variables.
+ *
+ * @param p Particle.
+ * @param W Pointer to the array in which the result needs to be stored (of size
+ * 5 or more).
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_get_primitive_variables(const struct part *restrict p, float *W) {
+
+ W[0] = p->rho;
+ W[1] = p->v[0];
+ W[2] = p->v[1];
+ W[3] = p->v[2];
+ W[4] = p->P;
+}
+
+/**
+ * @brief Check if the gradient matrix for this particle is well behaved.
+ */
+#define hydro_part_geometry_well_behaved(p) \
+ (p->geometry.wcorr > const_gizmo_min_wcorr)
+
+#endif /* SWIFT_GIZMO_MFM_HYDRO_GETTERS_H */
diff --git a/src/hydro/Gizmo/MFM/hydro_iact.h b/src/hydro/Gizmo/MFM/hydro_iact.h
deleted file mode 100644
index e8b8efaeb71553cfdb776ee0510de9c25ef70bb7..0000000000000000000000000000000000000000
--- a/src/hydro/Gizmo/MFM/hydro_iact.h
+++ /dev/null
@@ -1,491 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_GIZMO_MFM_HYDRO_IACT_H
-#define SWIFT_GIZMO_MFM_HYDRO_IACT_H
-
-#include "adiabatic_index.h"
-#include "hydro_gradients.h"
-#include "riemann.h"
-
-#include "./hydro_parameters.h"
-
-#define GIZMO_VOLUME_CORRECTION
-
-/**
- * @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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_density(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- float wi, wj, wi_dx, wj_dx;
-
- /* Get r and h inverse. */
- const float r = sqrtf(r2);
-
- /* Compute density of pi. */
- const float hi_inv = 1.0f / hi;
- const float xi = r * hi_inv;
- kernel_deval(xi, &wi, &wi_dx);
-
- pi->density.wcount += wi;
- pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
-
- /* these are eqns. (1) and (2) in the summary */
- pi->geometry.volume += wi;
- for (int k = 0; k < 3; k++)
- for (int l = 0; l < 3; l++)
- pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
-
- pi->geometry.centroid[0] -= dx[0] * wi;
- pi->geometry.centroid[1] -= dx[1] * wi;
- pi->geometry.centroid[2] -= dx[2] * wi;
-
- /* Compute density of pj. */
- const float hj_inv = 1.0f / hj;
- const float xj = r * hj_inv;
- kernel_deval(xj, &wj, &wj_dx);
-
- pj->density.wcount += wj;
- pj->density.wcount_dh -= (hydro_dimension * wj + xj * wj_dx);
-
- /* these are eqns. (1) and (2) in the summary */
- pj->geometry.volume += wj;
- for (int k = 0; k < 3; k++)
- for (int l = 0; l < 3; l++)
- pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
-
- pj->geometry.centroid[0] += dx[0] * wj;
- pj->geometry.centroid[1] += dx[1] * wj;
- pj->geometry.centroid[2] += dx[2] * wj;
-}
-
-/**
- * @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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- const struct part *restrict pj, float a, float H) {
-
- float wi, wi_dx;
-
- /* Get r and h inverse. */
- const float r = sqrtf(r2);
-
- const float hi_inv = 1.0f / hi;
- const float xi = r * hi_inv;
- kernel_deval(xi, &wi, &wi_dx);
-
- pi->density.wcount += wi;
- pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
-
- /* these are eqns. (1) and (2) in the summary */
- pi->geometry.volume += wi;
- for (int k = 0; k < 3; k++)
- for (int l = 0; l < 3; l++)
- pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
-
- pi->geometry.centroid[0] -= dx[0] * wi;
- pi->geometry.centroid[1] -= dx[1] * wi;
- pi->geometry.centroid[2] -= dx[2] * wi;
-}
-
-/**
- * @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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_gradient(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- 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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- 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), 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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, int mode, float a, float H) {
-
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
-
- /* Initialize local variables */
- float Bi[3][3];
- float Bj[3][3];
- float vi[3], vj[3];
- for (int k = 0; k < 3; k++) {
- for (int l = 0; l < 3; l++) {
- Bi[k][l] = pi->geometry.matrix_E[k][l];
- Bj[k][l] = pj->geometry.matrix_E[k][l];
- }
- vi[k] = pi->v[k]; /* particle velocities */
- vj[k] = pj->v[k];
- }
- const float Vi = pi->geometry.volume;
- const float Vj = pj->geometry.volume;
- float Wi[5], Wj[5];
- Wi[0] = pi->rho;
- Wi[1] = pi->v[0];
- Wi[2] = pi->v[1];
- Wi[3] = pi->v[2];
- Wi[4] = pi->P;
- Wj[0] = pj->rho;
- Wj[1] = pj->v[0];
- Wj[2] = pj->v[1];
- Wj[3] = pj->v[2];
- Wj[4] = pj->P;
-
- /* calculate the maximal signal velocity */
- float vmax;
- if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
- const float ci = gas_soundspeed_from_pressure(Wi[0], Wi[4]);
- const float cj = gas_soundspeed_from_pressure(Wj[0], Wj[4]);
- vmax = ci + cj;
- } else
- vmax = 0.0f;
-
- /* Velocity on the axis linking the particles */
- float dvdr = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
- (Wi[3] - Wj[3]) * dx[2];
-
- /* We only care about this velocity for particles moving towards each other
- */
- const float dvdotdx = min(dvdr, 0.0f);
-
- /* Get the signal velocity */
- vmax -= const_viscosity_beta * dvdotdx * r_inv;
-
- /* Store the signal velocity */
- pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
- if (mode == 1) pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
-
- /* Compute kernel of pi. */
- float wi, wi_dx;
- const float hi_inv = 1.0f / hi;
- const float hi_inv_dim = pow_dimension(hi_inv);
- const float xi = r * hi_inv;
- kernel_deval(xi, &wi, &wi_dx);
-
- /* Compute kernel of pj. */
- float wj, wj_dx;
- const float hj_inv = 1.0f / hj;
- const float hj_inv_dim = pow_dimension(hj_inv);
- const float xj = r * hj_inv;
- kernel_deval(xj, &wj, &wj_dx);
-
- /* Compute h_dt. We are going to use an SPH-like estimate of div_v for that */
- const float hidp1 = pow_dimension_plus_one(hi_inv);
- const float hjdp1 = pow_dimension_plus_one(hj_inv);
- const float wi_dr = hidp1 * wi_dx;
- const float wj_dr = hjdp1 * wj_dx;
- dvdr *= r_inv;
- if (pj->rho > 0.0f)
- pi->force.h_dt -= pj->conserved.mass * dvdr / pj->rho * wi_dr;
- if (mode == 1 && pi->rho > 0.0f)
- pj->force.h_dt -= pi->conserved.mass * dvdr / pi->rho * wj_dr;
-
- /* Compute (square of) area */
- /* eqn. (7) */
- float Anorm2 = 0.0f;
- float A[3];
- if (pi->geometry.wcorr > const_gizmo_min_wcorr &&
- pj->geometry.wcorr > const_gizmo_min_wcorr) {
- /* in principle, we use Vi and Vj as weights for the left and right
- contributions to the generalized surface vector.
- However, if Vi and Vj are very different (because they have very
- different
- smoothing lengths), then the expressions below are more stable. */
- float Xi = Vi;
- float Xj = Vj;
-#ifdef GIZMO_VOLUME_CORRECTION
- if (fabsf(Vi - Vj) / min(Vi, Vj) > 1.5f * hydro_dimension) {
- Xi = (Vi * hj + Vj * hi) / (hi + hj);
- Xj = Xi;
- }
-#endif
- for (int k = 0; k < 3; k++) {
- /* we add a minus sign since dx is pi->x - pj->x */
- A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
- wi * hi_inv_dim -
- Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
- wj * hj_inv_dim;
- Anorm2 += A[k] * A[k];
- }
- } else {
- /* ill condition gradient matrix: revert to SPH face area */
- const float Anorm =
- -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * r_inv;
- A[0] = -Anorm * dx[0];
- A[1] = -Anorm * dx[1];
- A[2] = -Anorm * dx[2];
- Anorm2 = Anorm * Anorm * r2;
- }
-
- /* if the interface has no area, nothing happens and we return */
- /* continuing results in dividing by zero and NaN's... */
- if (Anorm2 == 0.0f) return;
-
- /* Compute the area */
- const float Anorm_inv = 1.0f / sqrtf(Anorm2);
- const float Anorm = Anorm2 * Anorm_inv;
-
-#ifdef SWIFT_DEBUG_CHECKS
- /* For stability reasons, we do require A and dx to have opposite
- directions (basically meaning that the surface normal for the surface
- always points from particle i to particle j, as it would in a real
- moving-mesh code). If not, our scheme is no longer upwind and hence can
- become unstable. */
- const float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
- /* In GIZMO, Phil Hopkins reverts to an SPH integration scheme if this
- happens. We curently just ignore this case and display a message. */
- const float rdim = pow_dimension(r);
- if (dA_dot_dx > 1.e-6f * rdim) {
- message("Ill conditioned gradient matrix (%g %g %g %g %g)!", dA_dot_dx,
- Anorm, Vi, Vj, r);
- }
-#endif
-
- /* compute the normal vector of the interface */
- const float n_unit[3] = {A[0] * Anorm_inv, A[1] * Anorm_inv,
- A[2] * Anorm_inv};
-
- /* Compute interface position (relative to pi, since we don't need the actual
- * position) eqn. (8) */
- const float xfac = -hi / (hi + hj);
- const float xij_i[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
-
- /* Compute interface velocity */
- /* eqn. (9) */
- const float vij[3] = {vi[0] + (vi[0] - vj[0]) * xfac,
- vi[1] + (vi[1] - vj[1]) * xfac,
- vi[2] + (vi[2] - vj[2]) * xfac};
-
- /* complete calculation of position of interface */
- /* NOTE: dx is not necessarily just pi->x - pj->x but can also contain
- correction terms for periodicity. If we do the interpolation,
- we have to use xij w.r.t. the actual particle.
- => we need a separate xij for pi and pj... */
- /* tldr: we do not need the code below, but we do need the same code as above
- but then with i and j swapped */
- // for ( k = 0 ; k < 3 ; k++ )
- // xij[k] += pi->x[k];
-
- hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
-
- /* Boost the primitive variables to the frame of reference of the interface */
- /* Note that velocities are indices 1-3 in W */
- Wi[1] -= vij[0];
- Wi[2] -= vij[1];
- Wi[3] -= vij[2];
- Wj[1] -= vij[0];
- Wj[2] -= vij[1];
- Wj[3] -= vij[2];
-
- /* we don't need to rotate, we can use the unit vector in the Riemann problem
- * itself (see GIZMO) */
-
- float totflux[5];
- riemann_solve_for_middle_state_flux(Wi, Wj, n_unit, vij, totflux);
-
- /* Multiply with the interface surface area */
- totflux[1] *= Anorm;
- totflux[2] *= Anorm;
- totflux[3] *= Anorm;
- totflux[4] *= Anorm;
-
- /* Update conserved variables */
- /* We shamelessly exploit the fact that the mass flux is zero and omit all
- terms involving it */
- /* eqn. (16) */
- pi->flux.momentum[0] -= totflux[1];
- pi->flux.momentum[1] -= totflux[2];
- pi->flux.momentum[2] -= totflux[3];
- pi->flux.energy -= totflux[4];
-
-#ifndef GIZMO_TOTAL_ENERGY
- pi->flux.energy += totflux[1] * pi->v[0];
- pi->flux.energy += totflux[2] * pi->v[1];
- pi->flux.energy += totflux[3] * pi->v[2];
-#endif
-
- /* Note that this used to be much more complicated in early implementations of
- * the GIZMO scheme, as we wanted manifest conservation of conserved variables
- * and had to do symmetric flux exchanges. Now we don't care about manifest
- * conservation anymore and just assume the current fluxes are representative
- * for the flux over the entire time step. */
- if (mode == 1) {
- pj->flux.momentum[0] += totflux[1];
- pj->flux.momentum[1] += totflux[2];
- pj->flux.momentum[2] += totflux[3];
- pj->flux.energy += totflux[4];
-
-#ifndef GIZMO_TOTAL_ENERGY
- pj->flux.energy -= totflux[1] * pj->v[0];
- pj->flux.energy -= totflux[2] * pj->v[1];
- pj->flux.energy -= totflux[3] * pj->v[2];
-#endif
- }
-}
-
-/**
- * @brief Flux calculation between particle i and particle j
- *
- * This method calls runner_iact_fluxes_common with mode 1.
- *
- * @param r2 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_force(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1, a, H);
-}
-
-/**
- * @brief Flux calculation between particle i and particle j: non-symmetric
- * version
- *
- * This method calls runner_iact_fluxes_common with mode 0.
- *
- * @param r2 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
-}
-
-#endif /* SWIFT_GIZMO_MFM_HYDRO_IACT_H */
diff --git a/src/hydro/Gizmo/MFM/hydro_setters.h b/src/hydro/Gizmo/MFM/hydro_setters.h
new file mode 100644
index 0000000000000000000000000000000000000000..68767745d248c8e4eab6f9c2f9dc63a0c3092b85
--- /dev/null
+++ b/src/hydro/Gizmo/MFM/hydro_setters.h
@@ -0,0 +1,71 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFM_HYDRO_SETTERS_H
+#define SWIFT_GIZMO_MFM_HYDRO_SETTERS_H
+
+#include "const.h"
+
+/**
+ * @brief Update the fluxes for the particle with the given contributions,
+ * assuming the particle is to the left of the interparticle interface.
+ *
+ * @param p Particle.
+ * @param fluxes Fluxes accross the interface.
+ * @param dx Distance between the particles that share the interface.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_update_fluxes_left(
+ struct part *restrict p, const float *fluxes, const float *dx) {
+
+ p->flux.momentum[0] -= fluxes[1];
+ p->flux.momentum[1] -= fluxes[2];
+ p->flux.momentum[2] -= fluxes[3];
+ p->flux.energy -= fluxes[4];
+
+#ifndef GIZMO_TOTAL_ENERGY
+ p->flux.energy += fluxes[1] * p->v[0];
+ p->flux.energy += fluxes[2] * p->v[1];
+ p->flux.energy += fluxes[3] * p->v[2];
+#endif
+}
+
+/**
+ * @brief Update the fluxes for the particle with the given contributions,
+ * assuming the particle is to the right of the interparticle interface.
+ *
+ * @param p Particle.
+ * @param fluxes Fluxes accross the interface.
+ * @param dx Distance between the particles that share the interface.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_update_fluxes_right(struct part *restrict p, const float *fluxes,
+ const float *dx) {
+
+ p->flux.momentum[0] += fluxes[1];
+ p->flux.momentum[1] += fluxes[2];
+ p->flux.momentum[2] += fluxes[3];
+ p->flux.energy += fluxes[4];
+
+#ifndef GIZMO_TOTAL_ENERGY
+ p->flux.energy -= fluxes[1] * p->v[0];
+ p->flux.energy -= fluxes[2] * p->v[1];
+ p->flux.energy -= fluxes[3] * p->v[2];
+#endif
+}
+
+#endif /* SWIFT_GIZMO_MFM_HYDRO_SETTERS_H */
diff --git a/src/hydro/Gizmo/MFV/hydro_flux.h b/src/hydro/Gizmo/MFV/hydro_flux.h
new file mode 100644
index 0000000000000000000000000000000000000000..bc66e6088f09357ff2398c3b15f744029b86111a
--- /dev/null
+++ b/src/hydro/Gizmo/MFV/hydro_flux.h
@@ -0,0 +1,48 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFV_HYDRO_FLUX_H
+#define SWIFT_GIZMO_MFV_HYDRO_FLUX_H
+
+#include "riemann.h"
+
+/**
+ * @brief Compute the flux for the Riemann problem with the given left and right
+ * state, and interface normal, surface area and velocity.
+ *
+ * @param WL Left state variables.
+ * @param WR Right state variables.
+ * @param n_unit Unit vector of the interface.
+ * @param vLR Velocity of the interface.
+ * @param Anorm Surface area of the interface.
+ * @param fluxes Array to store the result in (of size 5 or more).
+ */
+__attribute__((always_inline)) INLINE static void hydro_compute_flux(
+ const float *WL, const float *WR, const float *n_unit, const float *vLR,
+ const float Anorm, float *fluxes) {
+
+ riemann_solve_for_flux(WL, WR, n_unit, vLR, fluxes);
+
+ fluxes[0] *= Anorm;
+ fluxes[1] *= Anorm;
+ fluxes[2] *= Anorm;
+ fluxes[3] *= Anorm;
+ fluxes[4] *= Anorm;
+}
+
+#endif /* SWIFT_GIZMO_MFV_HYDRO_FLUX_H */
diff --git a/src/hydro/Gizmo/MFV/hydro_getters.h b/src/hydro/Gizmo/MFV/hydro_getters.h
new file mode 100644
index 0000000000000000000000000000000000000000..e793380d7ff88624db6634eeb878749ff1371e73
--- /dev/null
+++ b/src/hydro/Gizmo/MFV/hydro_getters.h
@@ -0,0 +1,46 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFV_HYDRO_GETTERS_H
+#define SWIFT_GIZMO_MFV_HYDRO_GETTERS_H
+
+/**
+ * @brief Get a 5-element state vector W containing the primitive hydrodynamic
+ * variables.
+ *
+ * @param p Particle.
+ * @param W Pointer to the array in which the result needs to be stored (of size
+ * 5 or more).
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_get_primitive_variables(const struct part *restrict p, float *W) {
+
+ W[0] = p->primitives.rho;
+ W[1] = p->primitives.v[0];
+ W[2] = p->primitives.v[1];
+ W[3] = p->primitives.v[2];
+ W[4] = p->primitives.P;
+}
+
+/**
+ * @brief Check if the gradient matrix for this particle is well behaved.
+ */
+#define hydro_part_geometry_well_behaved(p) \
+ (p->density.wcorr > const_gizmo_min_wcorr)
+
+#endif /* SWIFT_GIZMO_MFV_HYDRO_GETTERS_H */
diff --git a/src/hydro/Gizmo/MFV/hydro_iact.h b/src/hydro/Gizmo/MFV/hydro_iact.h
deleted file mode 100644
index 1b8aac4ff89e4945a314ff72c44bff8dbfe1502c..0000000000000000000000000000000000000000
--- a/src/hydro/Gizmo/MFV/hydro_iact.h
+++ /dev/null
@@ -1,506 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_GIZMO_MFV_HYDRO_IACT_H
-#define SWIFT_GIZMO_MFV_HYDRO_IACT_H
-
-#include "adiabatic_index.h"
-#include "hydro_gradients.h"
-#include "riemann.h"
-
-#include "./hydro_parameters.h"
-
-#define GIZMO_VOLUME_CORRECTION
-
-/**
- * @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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_density(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- float wi, wj, wi_dx, wj_dx;
-
- /* Get r and h inverse. */
- const float r = sqrtf(r2);
-
- /* Compute density of pi. */
- const float hi_inv = 1.f / hi;
- const float xi = r * hi_inv;
- kernel_deval(xi, &wi, &wi_dx);
-
- pi->density.wcount += wi;
- pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
-
- /* these are eqns. (1) and (2) in the summary */
- pi->geometry.volume += wi;
- for (int k = 0; k < 3; k++)
- for (int l = 0; l < 3; l++)
- pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
-
- pi->geometry.centroid[0] -= dx[0] * wi;
- pi->geometry.centroid[1] -= dx[1] * wi;
- pi->geometry.centroid[2] -= dx[2] * wi;
-
- /* Compute density of pj. */
- const float hj_inv = 1.f / hj;
- const float xj = r * hj_inv;
- kernel_deval(xj, &wj, &wj_dx);
-
- pj->density.wcount += wj;
- pj->density.wcount_dh -= (hydro_dimension * wj + xj * wj_dx);
-
- /* these are eqns. (1) and (2) in the summary */
- pj->geometry.volume += wj;
- for (int k = 0; k < 3; k++)
- for (int l = 0; l < 3; l++)
- pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
-
- pj->geometry.centroid[0] += dx[0] * wj;
- pj->geometry.centroid[1] += dx[1] * wj;
- pj->geometry.centroid[2] += dx[2] * wj;
-}
-
-/**
- * @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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- const struct part *restrict pj, float a, float H) {
-
- float wi, wi_dx;
-
- /* Get r and h inverse. */
- const float r = sqrtf(r2);
-
- const float hi_inv = 1.f / hi;
- const float xi = r * hi_inv;
- kernel_deval(xi, &wi, &wi_dx);
-
- pi->density.wcount += wi;
- pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
-
- /* these are eqns. (1) and (2) in the summary */
- pi->geometry.volume += wi;
- for (int k = 0; k < 3; k++)
- for (int l = 0; l < 3; l++)
- pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
-
- pi->geometry.centroid[0] -= dx[0] * wi;
- pi->geometry.centroid[1] -= dx[1] * wi;
- pi->geometry.centroid[2] -= dx[2] * wi;
-}
-
-/**
- * @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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_gradient(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- 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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- 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), 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 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, int mode, float a, float H) {
-
- const float r_inv = 1.f / sqrtf(r2);
- const float r = r2 * r_inv;
-
- /* Initialize local variables */
- float Bi[3][3];
- float Bj[3][3];
- float vi[3], vj[3];
- for (int k = 0; k < 3; k++) {
- for (int l = 0; l < 3; l++) {
- Bi[k][l] = pi->geometry.matrix_E[k][l];
- Bj[k][l] = pj->geometry.matrix_E[k][l];
- }
- vi[k] = pi->v[k]; /* particle velocities */
- vj[k] = pj->v[k];
- }
- const float Vi = pi->geometry.volume;
- const float Vj = pj->geometry.volume;
- float Wi[5], Wj[5];
- 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;
-
- /* calculate the maximal signal velocity */
- float vmax;
- if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
- const float ci = gas_soundspeed_from_pressure(Wi[0], Wi[4]);
- const float cj = gas_soundspeed_from_pressure(Wj[0], Wj[4]);
- vmax = ci + cj;
- } else
- vmax = 0.0f;
-
- float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2];
-
- /* Velocity on the axis linking the particles */
- float dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
- (Wi[3] - Wj[3]) * dx[2];
- dvdotdx = min(dvdotdx, dvdr);
-
- /* We only care about this velocity for particles moving towards each others
- */
- dvdotdx = min(dvdotdx, 0.f);
-
- /* Get the signal velocity */
- vmax -= const_viscosity_beta * dvdotdx * r_inv;
-
- /* Store the signal velocity */
- pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
- if (mode == 1) pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
-
- /* Compute kernel of pi. */
- float wi, wi_dx;
- const float hi_inv = 1.f / hi;
- const float hi_inv_dim = pow_dimension(hi_inv);
- const float xi = r * hi_inv;
- kernel_deval(xi, &wi, &wi_dx);
-
- /* Compute kernel of pj. */
- float wj, wj_dx;
- const float hj_inv = 1.f / hj;
- const float hj_inv_dim = pow_dimension(hj_inv);
- const float xj = r * hj_inv;
- kernel_deval(xj, &wj, &wj_dx);
-
- /* Compute h_dt. We are going to use an SPH-like estimate of div_v for that */
- const float hidp1 = pow_dimension_plus_one(hi_inv);
- const float hjdp1 = pow_dimension_plus_one(hj_inv);
- const float wi_dr = hidp1 * wi_dx;
- const float wj_dr = hjdp1 * wj_dx;
- dvdr *= r_inv;
- if (pj->primitives.rho > 0.)
- pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
- if (mode == 1 && pi->primitives.rho > 0.)
- pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
-
- /* Compute (square of) area */
- /* eqn. (7) */
- float Anorm2 = 0.0f;
- float A[3];
- if (pi->density.wcorr > const_gizmo_min_wcorr &&
- pj->density.wcorr > const_gizmo_min_wcorr) {
- /* in principle, we use Vi and Vj as weights for the left and right
- contributions to the generalized surface vector.
- However, if Vi and Vj are very different (because they have very
- different
- smoothing lengths), then the expressions below are more stable. */
- float Xi = Vi;
- float Xj = Vj;
-#ifdef GIZMO_VOLUME_CORRECTION
- if (fabsf(Vi - Vj) / min(Vi, Vj) > 1.5f * hydro_dimension) {
- Xi = (Vi * hj + Vj * hi) / (hi + hj);
- Xj = Xi;
- }
-#endif
- for (int k = 0; k < 3; k++) {
- /* we add a minus sign since dx is pi->x - pj->x */
- A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
- wi * hi_inv_dim -
- Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
- wj * hj_inv_dim;
- Anorm2 += A[k] * A[k];
- }
- } else {
- /* ill condition gradient matrix: revert to SPH face area */
- const float Anorm =
- -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * r_inv;
- A[0] = -Anorm * dx[0];
- A[1] = -Anorm * dx[1];
- A[2] = -Anorm * dx[2];
- Anorm2 = Anorm * Anorm * r2;
- }
-
- /* if the interface has no area, nothing happens and we return */
- /* continuing results in dividing by zero and NaN's... */
- if (Anorm2 == 0.f) return;
-
- /* Compute the area */
- const float Anorm_inv = 1. / sqrtf(Anorm2);
- const float Anorm = Anorm2 * Anorm_inv;
-
-#ifdef SWIFT_DEBUG_CHECKS
- /* For stability reasons, we do require A and dx to have opposite
- directions (basically meaning that the surface normal for the surface
- always points from particle i to particle j, as it would in a real
- moving-mesh code). If not, our scheme is no longer upwind and hence can
- become unstable. */
- const float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
- /* In GIZMO, Phil Hopkins reverts to an SPH integration scheme if this
- happens. We curently just ignore this case and display a message. */
- const float rdim = pow_dimension(r);
- if (dA_dot_dx > 1.e-6f * rdim) {
- message("Ill conditioned gradient matrix (%g %g %g %g %g)!", dA_dot_dx,
- Anorm, Vi, Vj, r);
- }
-#endif
-
- /* compute the normal vector of the interface */
- const float n_unit[3] = {A[0] * Anorm_inv, A[1] * Anorm_inv,
- A[2] * Anorm_inv};
-
- /* Compute interface position (relative to pi, since we don't need the actual
- * position) eqn. (8) */
- const float xfac = -hi / (hi + hj);
- const float xij_i[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
-
- /* Compute interface velocity */
- /* eqn. (9) */
- const float vij[3] = {vi[0] + xfac * (vi[0] - vj[0]),
- vi[1] + xfac * (vi[1] - vj[1]),
- vi[2] + xfac * (vi[2] - vj[2])};
-
- hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
-
- /* Boost the primitive variables to the frame of reference of the interface */
- /* Note that velocities are indices 1-3 in W */
- Wi[1] -= vij[0];
- Wi[2] -= vij[1];
- Wi[3] -= vij[2];
- Wj[1] -= vij[0];
- Wj[2] -= vij[1];
- Wj[3] -= vij[2];
-
- /* we don't need to rotate, we can use the unit vector in the Riemann problem
- * itself (see GIZMO) */
-
- float totflux[5];
- riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
-
- /* Multiply with the interface surface area */
- totflux[0] *= Anorm;
- totflux[1] *= Anorm;
- totflux[2] *= Anorm;
- totflux[3] *= Anorm;
- totflux[4] *= Anorm;
-
- /* Store mass flux */
- const float mflux_i = totflux[0];
- pi->gravity.mflux[0] += mflux_i * dx[0];
- pi->gravity.mflux[1] += mflux_i * dx[1];
- pi->gravity.mflux[2] += mflux_i * dx[2];
-
- /* Update conserved variables */
- /* eqn. (16) */
- pi->conserved.flux.mass -= totflux[0];
- pi->conserved.flux.momentum[0] -= totflux[1];
- pi->conserved.flux.momentum[1] -= totflux[2];
- pi->conserved.flux.momentum[2] -= totflux[3];
- pi->conserved.flux.energy -= totflux[4];
-
-#ifndef GIZMO_TOTAL_ENERGY
- const float ekin_i = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
- pi->primitives.v[1] * pi->primitives.v[1] +
- pi->primitives.v[2] * pi->primitives.v[2]);
- pi->conserved.flux.energy += totflux[1] * pi->primitives.v[0];
- pi->conserved.flux.energy += totflux[2] * pi->primitives.v[1];
- pi->conserved.flux.energy += totflux[3] * pi->primitives.v[2];
- pi->conserved.flux.energy -= totflux[0] * ekin_i;
-#endif
-
- /* Note that this used to be much more complicated in early implementations of
- * the GIZMO scheme, as we wanted manifest conservation of conserved variables
- * and had to do symmetric flux exchanges. Now we don't care about manifest
- * conservation anymore and just assume the current fluxes are representative
- * for the flux over the entire time step. */
- if (mode == 1) {
- /* Store mass flux */
- const float mflux_j = totflux[0];
- pj->gravity.mflux[0] += mflux_j * dx[0];
- pj->gravity.mflux[1] += mflux_j * dx[1];
- pj->gravity.mflux[2] += mflux_j * dx[2];
-
- pj->conserved.flux.mass += totflux[0];
- pj->conserved.flux.momentum[0] += totflux[1];
- pj->conserved.flux.momentum[1] += totflux[2];
- pj->conserved.flux.momentum[2] += totflux[3];
- pj->conserved.flux.energy += totflux[4];
-
-#ifndef GIZMO_TOTAL_ENERGY
- const float ekin_j = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
- pj->primitives.v[1] * pj->primitives.v[1] +
- pj->primitives.v[2] * pj->primitives.v[2]);
- pj->conserved.flux.energy -= totflux[1] * pj->primitives.v[0];
- pj->conserved.flux.energy -= totflux[2] * pj->primitives.v[1];
- pj->conserved.flux.energy -= totflux[3] * pj->primitives.v[2];
- pj->conserved.flux.energy += totflux[0] * ekin_j;
-#endif
- }
-}
-
-/**
- * @brief Flux calculation between particle i and particle j
- *
- * This method calls runner_iact_fluxes_common with mode 1.
- *
- * @param r2 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_force(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1, a, H);
-}
-
-/**
- * @brief Flux calculation between particle i and particle j: non-symmetric
- * version
- *
- * This method calls runner_iact_fluxes_common with mode 0.
- *
- * @param r2 Comoving squared distance between particle i and particle j.
- * @param dx Comoving distance vector between the particles (dx = pi->x -
- * pj->x).
- * @param hi Comoving smoothing-length of particle i.
- * @param hj Comoving smoothing-length of particle j.
- * @param pi Particle i.
- * @param pj Particle j.
- * @param a Current scale factor.
- * @param H Current Hubble parameter.
- */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
- float r2, const float *dx, float hi, float hj, struct part *restrict pi,
- struct part *restrict pj, float a, float H) {
-
- runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
-}
-
-#endif /* SWIFT_GIZMO_MFV_HYDRO_IACT_H */
diff --git a/src/hydro/Gizmo/MFV/hydro_setters.h b/src/hydro/Gizmo/MFV/hydro_setters.h
new file mode 100644
index 0000000000000000000000000000000000000000..d3ae4177612cf47b9ca9b21c1aec74361176a515
--- /dev/null
+++ b/src/hydro/Gizmo/MFV/hydro_setters.h
@@ -0,0 +1,89 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFV_HYDRO_SETTERS_H
+#define SWIFT_GIZMO_MFV_HYDRO_SETTERS_H
+
+#include "const.h"
+
+/**
+ * @brief Update the fluxes for the particle with the given contributions,
+ * assuming the particle is to the left of the interparticle interface.
+ *
+ * @param p Particle.
+ * @param fluxes Fluxes accross the interface.
+ * @param dx Distance between the particles that share the interface.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_update_fluxes_left(
+ struct part *restrict p, const float *fluxes, const float *dx) {
+
+ p->gravity.mflux[0] += fluxes[0] * dx[0];
+ p->gravity.mflux[1] += fluxes[0] * dx[1];
+ p->gravity.mflux[2] += fluxes[0] * dx[2];
+
+ p->conserved.flux.mass -= fluxes[0];
+ p->conserved.flux.momentum[0] -= fluxes[1];
+ p->conserved.flux.momentum[1] -= fluxes[2];
+ p->conserved.flux.momentum[2] -= fluxes[3];
+ p->conserved.flux.energy -= fluxes[4];
+
+#ifndef GIZMO_TOTAL_ENERGY
+ const float ekin = 0.5f * (p->primitives.v[0] * p->primitives.v[0] +
+ p->primitives.v[1] * p->primitives.v[1] +
+ p->primitives.v[2] * p->primitives.v[2]);
+ p->conserved.flux.energy += fluxes[1] * p->primitives.v[0];
+ p->conserved.flux.energy += fluxes[2] * p->primitives.v[1];
+ p->conserved.flux.energy += fluxes[3] * p->primitives.v[2];
+ p->conserved.flux.energy -= fluxes[0] * ekin;
+#endif
+}
+
+/**
+ * @brief Update the fluxes for the particle with the given contributions,
+ * assuming the particle is to the right of the interparticle interface.
+ *
+ * @param p Particle.
+ * @param fluxes Fluxes accross the interface.
+ * @param dx Distance between the particles that share the interface.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_part_update_fluxes_right(struct part *restrict p, const float *fluxes,
+ const float *dx) {
+
+ p->gravity.mflux[0] += fluxes[0] * dx[0];
+ p->gravity.mflux[1] += fluxes[0] * dx[1];
+ p->gravity.mflux[2] += fluxes[0] * dx[2];
+
+ p->conserved.flux.mass += fluxes[0];
+ p->conserved.flux.momentum[0] += fluxes[1];
+ p->conserved.flux.momentum[1] += fluxes[2];
+ p->conserved.flux.momentum[2] += fluxes[3];
+ p->conserved.flux.energy += fluxes[4];
+
+#ifndef GIZMO_TOTAL_ENERGY
+ const float ekin = 0.5f * (p->primitives.v[0] * p->primitives.v[0] +
+ p->primitives.v[1] * p->primitives.v[1] +
+ p->primitives.v[2] * p->primitives.v[2]);
+ p->conserved.flux.energy -= fluxes[1] * p->primitives.v[0];
+ p->conserved.flux.energy -= fluxes[2] * p->primitives.v[1];
+ p->conserved.flux.energy -= fluxes[3] * p->primitives.v[2];
+ p->conserved.flux.energy += fluxes[0] * ekin;
+#endif
+}
+
+#endif /* SWIFT_GIZMO_MFV_HYDRO_SETTERS_H */
diff --git a/src/hydro/Gizmo/hydro_flux.h b/src/hydro/Gizmo/hydro_flux.h
new file mode 100644
index 0000000000000000000000000000000000000000..0fba48fe6aafa6073d8af90ac57a4a6cc44508a7
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_flux.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_HYDRO_FLUX_H
+#define SWIFT_GIZMO_HYDRO_FLUX_H
+
+#if defined(GIZMO_MFV_SPH)
+#include "MFV/hydro_flux.h"
+#elif defined(GIZMO_MFM_SPH)
+#include "MFM/hydro_flux.h"
+#endif
+
+#endif /* SWIFT_GIZMO_HYDRO_FLUX_H */
diff --git a/src/hydro/Gizmo/hydro_getters.h b/src/hydro/Gizmo/hydro_getters.h
new file mode 100644
index 0000000000000000000000000000000000000000..59217744441c00680b20ed9ab13f35ba7a22c5cb
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_getters.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_HYDRO_GETTERS_H
+#define SWIFT_GIZMO_HYDRO_GETTERS_H
+
+#if defined(GIZMO_MFV_SPH)
+#include "MFV/hydro_getters.h"
+#elif defined(GIZMO_MFM_SPH)
+#include "MFM/hydro_getters.h"
+#endif
+
+#endif /* SWIFT_GIZMO_HYDRO_GETTERS_H */
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
new file mode 100644
index 0000000000000000000000000000000000000000..da0b096bdda889cc1959fb310402707549edf43b
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_HYDRO_GRADIENTS_H
+#define SWIFT_GIZMO_HYDRO_GRADIENTS_H
+
+#if defined(GIZMO_MFV_SPH)
+#include "MFV/hydro_gradients.h"
+#elif defined(GIZMO_MFM_SPH)
+#include "MFM/hydro_gradients.h"
+#endif
+
+#endif /* SWIFT_GIZMO_HYDRO_GRADIENTS_H */
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 61f307037baf15a03732e4e9106bd794fb6ce049..243a3d37c343ec6ae4157ab25815b08e6d0c6ddc 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -19,10 +19,446 @@
#ifndef SWIFT_GIZMO_HYDRO_IACT_H
#define SWIFT_GIZMO_HYDRO_IACT_H
-#if defined(GIZMO_MFV_SPH)
-#include "MFV/hydro_iact.h"
-#elif defined(GIZMO_MFM_SPH)
-#include "MFM/hydro_iact.h"
+#include "hydro_flux.h"
+#include "hydro_getters.h"
+#include "hydro_gradients.h"
+#include "hydro_setters.h"
+
+#define GIZMO_VOLUME_CORRECTION
+
+/**
+ * @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 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ float wi, wj, wi_dx, wj_dx;
+
+ /* Get r and h inverse. */
+ const float r = sqrtf(r2);
+
+ /* Compute density of pi. */
+ const float hi_inv = 1.0f / hi;
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
+
+ /* these are eqns. (1) and (2) in the summary */
+ pi->geometry.volume += wi;
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+
+ pi->geometry.centroid[0] -= dx[0] * wi;
+ pi->geometry.centroid[1] -= dx[1] * wi;
+ pi->geometry.centroid[2] -= dx[2] * wi;
+
+ /* Compute density of pj. */
+ const float hj_inv = 1.0f / hj;
+ const float xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= (hydro_dimension * wj + xj * wj_dx);
+
+ /* these are eqns. (1) and (2) in the summary */
+ pj->geometry.volume += wj;
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
+
+ pj->geometry.centroid[0] += dx[0] * wj;
+ pj->geometry.centroid[1] += dx[1] * wj;
+ pj->geometry.centroid[2] += dx[2] * wj;
+}
+
+/**
+ * @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 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ const struct part *restrict pj, float a, float H) {
+
+ float wi, wi_dx;
+
+ /* Get r and h inverse. */
+ const float r = sqrtf(r2);
+
+ const float hi_inv = 1.0f / hi;
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
+
+ /* these are eqns. (1) and (2) in the summary */
+ pi->geometry.volume += wi;
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+
+ pi->geometry.centroid[0] -= dx[0] * wi;
+ pi->geometry.centroid[1] -= dx[1] * wi;
+ pi->geometry.centroid[2] -= dx[2] * wi;
+}
+
+/**
+ * @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 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_gradient(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ 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 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ 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), 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 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, int mode, float a, float H) {
+
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Initialize local variables */
+ float Bi[3][3];
+ float Bj[3][3];
+ float vi[3], vj[3];
+ for (int k = 0; k < 3; k++) {
+ for (int l = 0; l < 3; l++) {
+ Bi[k][l] = pi->geometry.matrix_E[k][l];
+ Bj[k][l] = pj->geometry.matrix_E[k][l];
+ }
+ vi[k] = pi->v[k]; /* particle velocities */
+ vj[k] = pj->v[k];
+ }
+ const float Vi = pi->geometry.volume;
+ const float Vj = pj->geometry.volume;
+ float Wi[5], Wj[5];
+ hydro_part_get_primitive_variables(pi, Wi);
+ hydro_part_get_primitive_variables(pj, Wj);
+
+ /* calculate the maximal signal velocity */
+ float vmax;
+ if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
+ const float ci = gas_soundspeed_from_pressure(Wi[0], Wi[4]);
+ const float cj = gas_soundspeed_from_pressure(Wj[0], Wj[4]);
+ vmax = ci + cj;
+ } else {
+ vmax = 0.0f;
+ }
+
+ float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Velocity on the axis linking the particles */
+ /* This velocity will be the same as dvdr for MFM, so hopefully this gets
+ optimised out. */
+ float dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
+ (Wi[3] - Wj[3]) * dx[2];
+
+ /* We only care about this velocity for particles moving towards each others
+ */
+ dvdotdx = min3(dvdr, dvdotdx, 0.f);
+
+ /* Get the signal velocity */
+ vmax -= const_viscosity_beta * dvdotdx * r_inv;
+
+ /* Store the signal velocity */
+ pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
+ if (mode == 1) {
+ pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
+ }
+
+ /* Compute kernel of pi. */
+ float wi, wi_dx;
+ const float hi_inv = 1.0f / hi;
+ const float hi_inv_dim = pow_dimension(hi_inv);
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* Compute kernel of pj. */
+ float wj, wj_dx;
+ const float hj_inv = 1.0f / hj;
+ const float hj_inv_dim = pow_dimension(hj_inv);
+ const float xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* Compute h_dt. We are going to use an SPH-like estimate of div_v for that */
+ const float hidp1 = pow_dimension_plus_one(hi_inv);
+ const float hjdp1 = pow_dimension_plus_one(hj_inv);
+ const float wi_dr = hidp1 * wi_dx;
+ const float wj_dr = hjdp1 * wj_dx;
+ dvdr *= r_inv;
+ if (Wj[0] > 0.0f) {
+ pi->force.h_dt -= pj->conserved.mass * dvdr / Wj[0] * wi_dr;
+ }
+ if (mode == 1 && Wi[0] > 0.0f) {
+ pj->force.h_dt -= pi->conserved.mass * dvdr / Wi[0] * wj_dr;
+ }
+
+ /* Compute (square of) area */
+ /* eqn. (7) */
+ float Anorm2 = 0.0f;
+ float A[3];
+ if (hydro_part_geometry_well_behaved(pi) &&
+ hydro_part_geometry_well_behaved(pj)) {
+ /* in principle, we use Vi and Vj as weights for the left and right
+ contributions to the generalized surface vector.
+ However, if Vi and Vj are very different (because they have very
+ different
+ smoothing lengths), then the expressions below are more stable. */
+ float Xi = Vi;
+ float Xj = Vj;
+#ifdef GIZMO_VOLUME_CORRECTION
+ if (fabsf(Vi - Vj) / min(Vi, Vj) > 1.5f * hydro_dimension) {
+ Xi = (Vi * hj + Vj * hi) / (hi + hj);
+ Xj = Xi;
+ }
+#endif
+ for (int k = 0; k < 3; k++) {
+ /* we add a minus sign since dx is pi->x - pj->x */
+ A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
+ wi * hi_inv_dim -
+ Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
+ wj * hj_inv_dim;
+ Anorm2 += A[k] * A[k];
+ }
+ } else {
+ /* ill condition gradient matrix: revert to SPH face area */
+ const float Anorm =
+ -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * r_inv;
+ A[0] = -Anorm * dx[0];
+ A[1] = -Anorm * dx[1];
+ A[2] = -Anorm * dx[2];
+ Anorm2 = Anorm * Anorm * r2;
+ }
+
+ /* if the interface has no area, nothing happens and we return */
+ /* continuing results in dividing by zero and NaN's... */
+ if (Anorm2 == 0.0f) {
+ return;
+ }
+
+ /* Compute the area */
+ const float Anorm_inv = 1.0f / sqrtf(Anorm2);
+ const float Anorm = Anorm2 * Anorm_inv;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* For stability reasons, we do require A and dx to have opposite
+ directions (basically meaning that the surface normal for the surface
+ always points from particle i to particle j, as it would in a real
+ moving-mesh code). If not, our scheme is no longer upwind and hence can
+ become unstable. */
+ const float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
+ /* In GIZMO, Phil Hopkins reverts to an SPH integration scheme if this
+ happens. We curently just ignore this case and display a message. */
+ const float rdim = pow_dimension(r);
+ if (dA_dot_dx > 1.e-6f * rdim) {
+ message("Ill conditioned gradient matrix (%g %g %g %g %g)!", dA_dot_dx,
+ Anorm, Vi, Vj, r);
+ }
#endif
+ /* compute the normal vector of the interface */
+ const float n_unit[3] = {A[0] * Anorm_inv, A[1] * Anorm_inv,
+ A[2] * Anorm_inv};
+
+ /* Compute interface position (relative to pi, since we don't need the actual
+ * position) eqn. (8) */
+ const float xfac = -hi / (hi + hj);
+ const float xij_i[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
+
+ /* Compute interface velocity */
+ /* eqn. (9) */
+ const float vij[3] = {vi[0] + (vi[0] - vj[0]) * xfac,
+ vi[1] + (vi[1] - vj[1]) * xfac,
+ vi[2] + (vi[2] - vj[2]) * xfac};
+
+ /* complete calculation of position of interface */
+ /* NOTE: dx is not necessarily just pi->x - pj->x but can also contain
+ correction terms for periodicity. If we do the interpolation,
+ we have to use xij w.r.t. the actual particle.
+ => we need a separate xij for pi and pj... */
+ /* tldr: we do not need the code below, but we do need the same code as above
+ but then with i and j swapped */
+ // for ( k = 0 ; k < 3 ; k++ )
+ // xij[k] += pi->x[k];
+
+ hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
+
+ /* Boost the primitive variables to the frame of reference of the interface */
+ /* Note that velocities are indices 1-3 in W */
+ Wi[1] -= vij[0];
+ Wi[2] -= vij[1];
+ Wi[3] -= vij[2];
+ Wj[1] -= vij[0];
+ Wj[2] -= vij[1];
+ Wj[3] -= vij[2];
+
+ /* we don't need to rotate, we can use the unit vector in the Riemann problem
+ * itself (see GIZMO) */
+
+ float totflux[5];
+ hydro_compute_flux(Wi, Wj, n_unit, vij, Anorm, totflux);
+
+ hydro_part_update_fluxes_left(pi, totflux, dx);
+
+ /* Note that this used to be much more complicated in early implementations of
+ * the GIZMO scheme, as we wanted manifest conservation of conserved variables
+ * and had to do symmetric flux exchanges. Now we don't care about manifest
+ * conservation anymore and just assume the current fluxes are representative
+ * for the flux over the entire time step. */
+ if (mode == 1) {
+ hydro_part_update_fluxes_right(pj, totflux, dx);
+ }
+}
+
+/**
+ * @brief Flux calculation between particle i and particle j
+ *
+ * This method calls runner_iact_fluxes_common with mode 1.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1, a, H);
+}
+
+/**
+ * @brief Flux calculation between particle i and particle j: non-symmetric
+ * version
+ *
+ * This method calls runner_iact_fluxes_common with mode 0.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
+}
+
#endif /* SWIFT_GIZMO_HYDRO_IACT_H */
diff --git a/src/hydro/Gizmo/hydro_setters.h b/src/hydro/Gizmo/hydro_setters.h
new file mode 100644
index 0000000000000000000000000000000000000000..351796979f84417647719dd2076b9f32fc1f551b
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_setters.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_HYDRO_SETTERS_H
+#define SWIFT_GIZMO_HYDRO_SETTERS_H
+
+#if defined(GIZMO_MFV_SPH)
+#include "MFV/hydro_setters.h"
+#elif defined(GIZMO_MFM_SPH)
+#include "MFM/hydro_setters.h"
+#endif
+
+#endif /* SWIFT_GIZMO_HYDRO_SETTERS_H */