diff --git a/configure.ac b/configure.ac
index ab2ed97305041f1a150d25a7d85488fcab6eed67..b91303ef1ec191a8bf0ddeb3fdf3bbafd483a6a5 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1150,7 +1150,7 @@ esac
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
- [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, debug default: gadget2@:>@]
+ [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, debug default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
@@ -1177,6 +1177,9 @@ case "$with_hydro" in
pressure-energy)
AC_DEFINE([HOPKINS_PU_SPH], [1], [Pressure-Energy SPH])
;;
+ pressure-energy-monaghan)
+ AC_DEFINE([HOPKINS_PU_SPH_MONAGHAN], [1], [Pressure-Energy SPH with M&M Variable A.V.])
+ ;;
default)
AC_DEFINE([DEFAULT_SPH], [1], [Default SPH])
;;
diff --git a/examples/KelvinHelmholtz_2D/makeMovie.py b/examples/KelvinHelmholtz_2D/makeMovie.py
index 84fe99106bf607830e89b6aa663135b48b6c0744..a52784891ab4689dcd59dc27945e573e602785f3 100644
--- a/examples/KelvinHelmholtz_2D/makeMovie.py
+++ b/examples/KelvinHelmholtz_2D/makeMovie.py
@@ -91,6 +91,7 @@ if __name__ == "__main__":
# Creation of first frame
fig, ax = plt.subplots(1, 1, figsize=(1, 1), frameon=False)
+ ax.axis("off") # Remove annoying black frame.
data_x, data_y, density = load_and_extract("kelvinhelmholtz_0000.hdf5")
diff --git a/examples/SodShock_1D/plotSolution.py b/examples/SodShock_1D/plotSolution.py
index e001a8d87a03cb246be63ab10d245f95eb1a7ce7..0149d66a0c28c777a4265da10d86ed160086995d 100644
--- a/examples/SodShock_1D/plotSolution.py
+++ b/examples/SodShock_1D/plotSolution.py
@@ -70,11 +70,11 @@ snap = int(sys.argv[1])
sim = h5py.File("sodShock_%04d.hdf5"%snap, "r")
boxSize = sim["/Header"].attrs["BoxSize"][0]
time = sim["/Header"].attrs["Time"][0]
-scheme = sim["/HydroScheme"].attrs["Scheme"]
-kernel = sim["/HydroScheme"].attrs["Kernel function"]
+scheme = str(sim["/HydroScheme"].attrs["Scheme"])
+kernel = str(sim["/HydroScheme"].attrs["Kernel function"])
neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta = sim["/HydroScheme"].attrs["Kernel eta"]
-git = sim["Code"].attrs["Git Revision"]
+git = str(sim["Code"].attrs["Git Revision"])
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0]
@@ -82,6 +82,11 @@ u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
rho = sim["/PartType0/Density"][:]
+try:
+ alpha = sim["/PartType0/Viscosity"][:]
+ plot_alpha = True
+except:
+ plot_alpha = False
N = 1000 # Number of points
x_min = -1.
@@ -259,14 +264,23 @@ ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
xlim(-0.5, 0.5)
ylim(0.8, 2.2)
-# Entropy profile ---------------------------------
+# Entropy/alpha profile ---------------------------------
subplot(235)
-plot(x, S, '.', color='r', ms=4.0)
-plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+
+if plot_alpha:
+ plot(x, alpha, '.', color='r', ms=4.0)
+ ylabel(r"${\rm{Viscosity}}~\alpha$", labelpad=0)
+ # Show location of shock
+ plot([x_56, x_56], [-100, 100], color="k", alpha=0.5, ls="dashed", lw=1.2)
+ ylim(0, 1)
+else:
+ plot(x, S, '.', color='r', ms=4.0)
+ plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+ ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+ ylim(0.8, 3.8)
+
xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=0)
xlim(-0.5, 0.5)
-ylim(0.8, 3.8)
# Information -------------------------------------
subplot(236, frameon=False)
@@ -284,5 +298,6 @@ ylim(0, 1)
xticks([])
yticks([])
+tight_layout()
savefig("SodShock.png", dpi=200)
diff --git a/src/debug.c b/src/debug.c
index a980ef9457a31d4b36f99926a24d7bc55822670a..809d7048c45888eb41bd277bdb4971d469a98dc3 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -48,6 +48,8 @@
#include "./hydro/PressureEntropy/hydro_debug.h"
#elif defined(HOPKINS_PU_SPH)
#include "./hydro/PressureEnergy/hydro_debug.h"
+#elif defined(HOPKINS_PU_SPH_MONAGHAN)
+#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_debug.h"
#elif defined(GIZMO_MFV_SPH)
diff --git a/src/hydro.h b/src/hydro.h
index b3716996cc4da68f9445adccd12315b32d81a34c..15c45c1dcfa1217d842904dfc1303acea607e3ab 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -45,6 +45,12 @@
#include "./hydro/PressureEnergy/hydro.h"
#include "./hydro/PressureEnergy/hydro_iact.h"
#define SPH_IMPLEMENTATION "Pressure-Energy SPH (Hopkins 2013)"
+#elif defined(HOPKINS_PU_SPH_MONAGHAN)
+#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro.h"
+#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h"
+#define SPH_IMPLEMENTATION \
+ "Pressure-Energy SPH (Hopkins 2013) with a Morris and Monaghan (1997) " \
+ "variable artificial viscosity."
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro.h"
#include "./hydro/Default/hydro_iact.h"
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index 3274c1b1b58adf0fa6a4c94132fa5a87186e59ce..4c3cc5c1c588e19de0d4833fc867ae9c0aed1209 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -472,8 +472,10 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
- /* Finish calculation of the velocity divergence */
- p->density.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2;
+ /* Finish calculation of the velocity divergence, including hubble flow term
+ */
+ p->density.div_v *=
+ h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension;
}
/**
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index 2ed7fe8cb8112c42de933a2ca315966e67108f0a..4146e61a53dd7ece57e263cb90308e2579aa3930 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -284,7 +284,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
wj_dr * dvdr * r_inv;
/* Viscosity term */
- const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term;
@@ -404,7 +404,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
wi_dr * dvdr * r_inv;
/* Viscosity term */
- const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term;
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..060694a6afa850c4d4815899fde1450316da81f5
--- /dev/null
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
@@ -0,0 +1,766 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ * Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ * 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_PRESSURE_ENERGY_MORRIS_HYDRO_H
+#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H
+
+/**
+ * @file PressureEnergy/hydro.h
+ * @brief P-U conservative implementation of SPH (Non-neighbour loop
+ * equations)
+ *
+ * The thermal variable is the internal energy (u). A simple variable
+ * viscosity term (Morris & Monaghan 1997) with a Balsara switch is
+ * implemented.
+ *
+ * No thermal conduction term is implemented.
+ *
+ * This implementation corresponds to the one presented in the SWIFT
+ * documentation and in Hopkins, "A general class of Lagrangian smoothed
+ * particle hydrodynamics methods and implications for fluid mixing problems",
+ * MNRAS, 2013.
+ */
+
+#include "adiabatic_index.h"
+#include "approx_math.h"
+#include "cosmology.h"
+#include "dimension.h"
+#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "hydro_space.h"
+#include "kernel_hydro.h"
+#include "minmax.h"
+
+#include <float.h>
+
+/**
+ * @brief Returns the comoving internal energy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part *restrict p) {
+
+ return p->u;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return p->u * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
+ const struct part *restrict p) {
+
+ return p->pressure_bar;
+}
+
+/**
+ * @brief Returns the physical pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties and
+ * convert it to physical coordinates.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_pressure * p->pressure_bar;
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
+ const struct part *restrict p) {
+
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the physical entropy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ /* Note: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the comoving sound speed of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part *restrict p) {
+
+ /* Compute the sound speed -- see theory section for justification */
+ /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
+ const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho);
+
+ return square_rooted;
+}
+
+/**
+ * @brief Returns the physical sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_soundspeed(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_sound_speed * p->force.soundspeed;
+}
+
+/**
+ * @brief Returns the comoving internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_internal_energy(const struct part *restrict p) {
+
+ return p->u;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return p->u * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_entropy(const struct part *restrict p) {
+
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the physical entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_entropy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ /* Note: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the comoving density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
+ const struct part *restrict p) {
+
+ return p->rho;
+}
+
+/**
+ * @brief Returns the comoving density of a particle.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return cosmo->a3_inv * p->rho;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+ const struct part *restrict p) {
+
+ return p->mass;
+}
+
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
+ * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+ const struct part *restrict p, const struct xpart *xp, float dt_kick_hydro,
+ float dt_kick_grav, float v[3]) {
+
+ v[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro +
+ xp->a_grav[0] * dt_kick_grav;
+ v[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro +
+ xp->a_grav[1] * dt_kick_grav;
+ v[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro +
+ xp->a_grav[2] * dt_kick_grav;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
+ const struct part *restrict p) {
+
+ return p->u_dt;
+}
+
+/**
+ * @brief Sets the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
+ struct part *restrict p, float du_dt) {
+
+ p->u_dt = du_dt;
+}
+
+/**
+ * @brief Computes the hydro time-step of a given particle
+ *
+ * This function returns the time-step of a particle given its hydro-dynamical
+ * state. A typical time-step calculation would be the use of the CFL condition.
+ *
+ * @param p Pointer to the particle data
+ * @param xp Pointer to the extended particle data
+ * @param hydro_properties The SPH parameters
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct hydro_props *restrict hydro_properties,
+ const struct cosmology *restrict cosmo) {
+
+ const float CFL_condition = hydro_properties->CFL_condition;
+
+ /* CFL condition */
+ const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
+ (cosmo->a_factor_sound_speed * p->force.v_sig);
+
+ const float dt_u_change =
+ (p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt) : FLT_MAX;
+
+ return fminf(dt_cfl, dt_u_change);
+}
+
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
+/**
+ * @brief Prepares a particle for the density calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the various density loop over neighbours. Typically, all fields of the
+ * density sub-structure of a particle get zeroed in here.
+ *
+ * @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
+ */
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part *restrict p, const struct hydro_space *hs) {
+
+ p->density.wcount = 0.f;
+ p->density.wcount_dh = 0.f;
+ p->rho = 0.f;
+ p->density.rho_dh = 0.f;
+ p->pressure_bar = 0.f;
+ p->density.pressure_bar_dh = 0.f;
+
+ p->density.div_v = 0.f;
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and add the self-contribution term.
+ * Additional quantities such as velocity gradients will also get the final
+ * terms added to them here.
+ *
+ * Also adds/multiplies the cosmological terms if need be.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part *restrict p, const struct cosmology *cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+ /* Final operation on the density (add self-contribution). */
+ p->rho += p->mass * kernel_root;
+ p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
+ p->pressure_bar += p->mass * p->u * kernel_root;
+ p->density.pressure_bar_dh -= hydro_dimension * p->mass * p->u * kernel_root;
+ p->density.wcount += kernel_root;
+ p->density.wcount_dh -= hydro_dimension * kernel_root;
+
+ /* Finish the calculation by inserting the missing h-factors */
+ p->rho *= h_inv_dim;
+ p->density.rho_dh *= h_inv_dim_plus_one;
+ p->pressure_bar *= (h_inv_dim * hydro_gamma_minus_one);
+ p->density.pressure_bar_dh *= (h_inv_dim_plus_one * hydro_gamma_minus_one);
+ p->density.wcount *= h_inv_dim;
+ p->density.wcount_dh *= h_inv_dim_plus_one;
+
+ const float rho_inv = 1.f / p->rho;
+ const float a_inv2 = cosmo->a2_inv;
+
+ /* Finish calculation of the velocity curl components */
+ p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+ p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+ p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+
+ /* Finish calculation of the velocity divergence */
+ p->density.div_v *=
+ h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * In the desperate case where a particle has no neighbours (likely because
+ * of the h_max ceiling), set the particle fields to something sensible to avoid
+ * NaNs in the next calculations.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+ /* Re-set problematic values */
+ p->rho = p->mass * kernel_root * h_inv_dim;
+ p->pressure_bar =
+ p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim;
+ p->density.wcount = kernel_root * h_inv_dim;
+ p->density.rho_dh = 0.f;
+ p->density.wcount_dh = 0.f;
+ p->density.pressure_bar_dh = 0.f;
+
+ p->density.div_v = 0.f;
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * This function is called in the ghost task to convert some quantities coming
+ * from the density loop over neighbours into quantities ready to be used in the
+ * force loop over neighbours. Quantities are typically read from the density
+ * sub-structure and written to the force sub-structure.
+ * Examples of calculations done here include the calculation of viscosity term
+ * constants, thermal conduction terms, hydro conversions, etc.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties.
+ * @param dt_alpha The time-step used to evolve non-cosmological quantities such
+ * as the artificial viscosity.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+ const float dt_alpha) {
+
+ const float fac_mu = cosmo->a_factor_mu;
+
+ const float h_inv = 1.f / p->h;
+
+ /* Compute the norm of the curl */
+ const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
+ p->density.rot_v[1] * p->density.rot_v[1] +
+ p->density.rot_v[2] * p->density.rot_v[2]);
+
+ /* Compute the norm of div v */
+ const float abs_div_v = fabsf(p->density.div_v);
+
+ /* Compute the sound speed -- see theory section for justification */
+ const float soundspeed = hydro_get_comoving_soundspeed(p);
+
+ /* Compute the Balsara switch */
+ const float balsara =
+ abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu * h_inv);
+
+ /* Compute the "grad h" term */
+ const float common_factor = p->h / (hydro_dimension * p->density.wcount);
+ const float grad_h_term = (p->density.pressure_bar_dh * common_factor *
+ hydro_one_over_gamma_minus_one) /
+ (1.f + common_factor * p->density.wcount_dh);
+
+ /* Artificial viscosity updates */
+
+ const float inverse_tau = hydro_props->viscosity.length * soundspeed * h_inv;
+ const float source_term = -1.f * min(p->density.div_v, 0.f);
+
+ /* Compute da/dt */
+ const float alpha_time_differential =
+ source_term + (hydro_props->viscosity.alpha_min - p->alpha) * inverse_tau;
+
+ /* Update variables. */
+ p->alpha += alpha_time_differential * dt_alpha;
+ p->force.f = grad_h_term;
+ p->force.soundspeed = soundspeed;
+ p->force.balsara = balsara;
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * Resets all hydro acceleration and time derivative fields in preparation
+ * for the sums taking place in the various force tasks.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part *restrict p) {
+
+ /* Reset the acceleration. */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ /* Reset the time derivatives. */
+ p->u_dt = 0.0f;
+ p->force.h_dt = 0.0f;
+ p->force.v_sig = p->force.soundspeed;
+}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+ struct part *restrict p, const struct xpart *restrict xp) {
+
+ /* Re-set the predicted velocities */
+ p->v[0] = xp->v_full[0];
+ p->v[1] = xp->v_full[1];
+ p->v[2] = xp->v_full[2];
+
+ /* Re-set the entropy */
+ p->u = xp->u_full;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * Additional hydrodynamic quantites are drifted forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * Note the different time-step sizes used for the different quantities as they
+ * include cosmological factors.
+ *
+ * @param p The particle.
+ * @param xp The extended data of the particle.
+ * @param dt_drift The drift time-step for positions.
+ * @param dt_therm The drift time-step for thermal quantities.
+ */
+__attribute__((always_inline)) INLINE static void hydro_predict_extra(
+ struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
+ float dt_therm) {
+
+ const float h_inv = 1.f / p->h;
+
+ /* Predict smoothing length */
+ const float w1 = p->force.h_dt * h_inv * dt_drift;
+ if (fabsf(w1) < 0.2f)
+ p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+ else
+ p->h *= expf(w1);
+
+ /* Predict density and weighted pressure */
+ const float w2 = -hydro_dimension * w1;
+ if (fabsf(w2) < 0.2f) {
+ const float expf_approx =
+ approx_expf(w2); /* 4th order expansion of exp(w) */
+ p->rho *= expf_approx;
+ p->pressure_bar *= expf_approx;
+ } else {
+ const float expf_exact = expf(w2);
+ p->rho *= expf_exact;
+ p->pressure_bar *= expf_exact;
+ }
+
+ /* Predict the internal energy */
+ p->u += p->u_dt * dt_therm;
+
+ /* Compute the new sound speed */
+ const float soundspeed = hydro_get_comoving_soundspeed(p);
+
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Finishes the force calculation.
+ *
+ * Multiplies the force and accelerations by the appropiate constants
+ * and add the self-contribution term. In most cases, there is little
+ * to do here.
+ *
+ * Cosmological terms are also added/multiplied here.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part *restrict p, const struct cosmology *cosmo) {
+
+ p->force.h_dt *= p->h * hydro_dimension_inv;
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * Additional hydrodynamic quantites are kicked forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * @param p The particle to act upon.
+ * @param xp The particle extended data to act upon.
+ * @param dt_therm The time-step for this kick (for thermodynamic quantities).
+ * @param dt_grav The time-step for this kick (for gravity quantities).
+ * @param dt_hydro The time-step for this kick (for hydro quantities).
+ * @param dt_kick_corr The time-step for this kick (for gravity corrections).
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+ float dt_grav, float dt_hydro, float dt_kick_corr,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
+
+ /* Do not decrease the energy by more than a factor of 2*/
+ if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
+ p->u_dt = -0.5f * xp->u_full / dt_therm;
+ }
+ xp->u_full += p->u_dt * dt_therm;
+
+ /* Apply the minimal energy limit */
+ const float min_energy =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+ if (xp->u_full < min_energy) {
+ xp->u_full = min_energy;
+ p->u_dt = 0.f;
+ }
+
+ /* Compute the sound speed */
+ const float soundspeed = hydro_get_comoving_soundspeed(p);
+
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Converts hydro quantity of a particle at the start of a run
+ *
+ * This function is called once at the end of the engine_init_particle()
+ * routine (at the start of a calculation) after the densities of
+ * particles have been computed.
+ * This can be used to convert internal energy into entropy for instance.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle to act upon
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme.
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
+
+ /* Convert the physcial internal energy to the comoving one. */
+ /* u' = a^(3(g-1)) u */
+ const float factor = 1.f / cosmo->a_factor_internal_energy;
+ p->u *= factor;
+ xp->u_full = p->u;
+
+ /* Apply the minimal energy limit */
+ const float min_energy =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+ if (xp->u_full < min_energy) {
+ xp->u_full = min_energy;
+ p->u = min_energy;
+ p->u_dt = 0.f;
+ }
+
+ /* Start out with a 'regular' AV for comparison to other schemes */
+ p->alpha = hydro_props->viscosity.alpha;
+
+ /* Note that unlike Minimal the pressure and sound speed cannot be calculated
+ * here because they are smoothed properties in this scheme. */
+}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions or assignments between the particle
+ * and extended particle fields.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part *restrict p, struct xpart *restrict xp) {
+
+ p->time_bin = 0;
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+ xp->a_grav[0] = 0.f;
+ xp->a_grav[1] = 0.f;
+ xp->a_grav[2] = 0.f;
+ xp->u_full = p->u;
+
+ hydro_reset_acceleration(p);
+ hydro_init_part(p, NULL);
+}
+
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->u = u_init;
+}
+
+#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H */
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..ead5fcc0c842d8018f784a1084941bdb9ebcb6ca
--- /dev/null
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
@@ -0,0 +1,47 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ * Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ * 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_PRESSURE_ENERGY_MORRIS_HYDRO_DEBUG_H
+#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_DEBUG_H
+/**
+ * @file PressureEnergy/hydro_debug.h
+ * @brief P-U conservative implementation of SPH (Debugging routines)
+ *
+ * The thermal variable is the internal energy (u). A simple variable
+ * viscosity term (Morris & Monaghan 1997) with a Balsara switch is
+ * implemented.
+ */
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ const struct part* p, const struct xpart* xp) {
+ printf(
+ "x=[%.3e,%.3e,%.3e], "
+ "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
+ "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+ "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, \n"
+ "p_dh=%.3e, p_bar=%.3e \n"
+ "time_bin=%d, alpha=%.3e\n",
+ p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
+ xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+ p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
+ p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+ p->density.pressure_bar_dh, p->pressure_bar, p->time_bin, p->alpha);
+}
+
+#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_DEBUG_H */
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..747fca714ce20d9c2b018e14ac24a6492c51a75f
--- /dev/null
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
@@ -0,0 +1,427 @@
+/*******************************************************************************
+ * This file is part* of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ * Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ * 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_PRESSURE_ENERGY_MORRIS_HYDRO_IACT_H
+#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IACT_H
+
+/**
+ * @file PressureEnergy/hydro_iact.h
+ * @brief P-U implementation of SPH (Neighbour loop equations)
+ *
+ * The thermal variable is the internal energy (u). A simple variable
+ * viscosity term (Morris & Monaghan 1997) with a Balsara switch is
+ * implemented.
+ *
+ * No thermal conduction term is implemented.
+ *
+ * See PressureEnergy/hydro.h for references.
+ */
+
+#include "adiabatic_index.h"
+#include "minmax.h"
+
+/**
+ * @brief Density interaction between two part*icles.
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle.
+ * @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* pi,
+ struct part* pj, float a, float H) {
+
+ float wi, wj, wi_dx, wj_dx;
+ float dv[3], curlvr[3];
+
+ const float r = sqrtf(r2);
+
+ /* Get the masses. */
+ const float mi = pi->mass;
+ const float mj = pj->mass;
+
+ /* Compute density of pi. */
+ const float hi_inv = 1.f / hi;
+ const float ui = r * hi_inv;
+
+ kernel_deval(ui, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+
+ pi->pressure_bar += mj * wi * pj->u;
+ pi->density.pressure_bar_dh -=
+ mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ /* Compute density of pj. */
+ const float hj_inv = 1.f / hj;
+ const float uj = r * hj_inv;
+ kernel_deval(uj, &wj, &wj_dx);
+
+ pj->rho += mi * wj;
+ pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+ pj->pressure_bar += mi * wj * pi->u;
+ pj->density.pressure_bar_dh -=
+ mi * pi->u * (hydro_dimension * wj + uj * wj_dx);
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+
+ /* Now we need to compute the div terms */
+ const float r_inv = 1.f / r;
+ const float faci = mj * wi_dx * r_inv;
+ const float facj = mi * wj_dx * r_inv;
+
+ /* Compute dv dot r */
+ dv[0] = pi->v[0] - pj->v[0];
+ dv[1] = pi->v[1] - pj->v[1];
+ dv[2] = pi->v[2] - pj->v[2];
+ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+ pi->density.div_v -= faci * dvdr;
+ pj->density.div_v -= facj * dvdr;
+
+ /* Compute dv cross r */
+ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+ pi->density.rot_v[0] += faci * curlvr[0];
+ pi->density.rot_v[1] += faci * curlvr[1];
+ pi->density.rot_v[2] += faci * curlvr[2];
+
+ /* Negative because of the change in sign of dx & dv. */
+ pj->density.rot_v[0] += facj * curlvr[0];
+ pj->density.rot_v[1] += facj * curlvr[1];
+ pj->density.rot_v[2] += facj * curlvr[2];
+}
+
+/**
+ * @brief Density interaction between two part*icles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @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* pi,
+ const struct part* pj, float a, float H) {
+
+ float wi, wi_dx;
+ float dv[3], curlvr[3];
+
+ /* Get the masses. */
+ const float mj = pj->mass;
+
+ /* Get r and r inverse. */
+ const float r = sqrtf(r2);
+
+ const float h_inv = 1.f / hi;
+ const float ui = r * h_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+
+ pi->pressure_bar += mj * wi * pj->u;
+
+ pi->density.pressure_bar_dh -=
+ mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ const float r_inv = 1.f / r;
+ const float faci = mj * wi_dx * r_inv;
+
+ /* Compute dv dot r */
+ dv[0] = pi->v[0] - pj->v[0];
+ dv[1] = pi->v[1] - pj->v[1];
+ dv[2] = pi->v[2] - pj->v[2];
+ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+ pi->density.div_v -= faci * dvdr;
+
+ /* Compute dv cross r */
+ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+ pi->density.rot_v[0] += faci * curlvr[0];
+ pi->density.rot_v[1] += faci * curlvr[1];
+ pi->density.rot_v[2] += faci * curlvr[2];
+}
+
+/**
+ * @brief Force interaction between two part*icles.
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle.
+ * @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* pi,
+ struct part* pj, float a, float H) {
+
+ /* Cosmological factors entering the EoMs */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+
+ /* Recover some data */
+ const float mj = pj->mass;
+ const float mi = pi->mass;
+
+ const float miui = mi * pi->u;
+ const float mjuj = mj * pj->u;
+
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+ /* Compute gradient terms */
+ const float f_ij = 1.f - (pi->force.f / mjuj);
+ const float f_ji = 1.f - (pj->force.f / miui);
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ const float wi_dr = hid_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hjd_inv * wj_dx;
+
+ /* Compute dv dot r. */
+ const 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];
+
+ /* Includes the hubble flow term; not used for du/dt */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+
+ /* Are the part*icles moving towards each others ? */
+ const float omega_ij = min(dvdr_Hubble, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Compute sound speeds and signal velocity */
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Balsara term */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
+
+ /* Construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float alpha = 0.5f * (pi->alpha + pj->alpha);
+ const float visc =
+ -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
+
+ /* Convolve with the kernel */
+ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+
+ /* SPH acceleration term */
+ const float sph_acc_term =
+ pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
+ ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
+ r_inv;
+
+ /* Assemble the acceleration */
+ const float acc = sph_acc_term + visc_acc_term;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ pj->a_hydro[0] += mi * acc * dx[0];
+ pj->a_hydro[1] += mi * acc * dx[1];
+ pj->a_hydro[2] += mi * acc * dx[2];
+
+ /* Get the time derivative for u. */
+ const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
+ pj->u * pi->u * (f_ij / pi->pressure_bar) *
+ wi_dr * dvdr * r_inv;
+ const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one *
+ pi->u * pj->u * (f_ji / pj->pressure_bar) *
+ wj_dr * dvdr * r_inv;
+
+ /* Viscosity term */
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
+
+ /* Assemble the energy equation term */
+ const float du_dt_i = sph_du_term_i + visc_du_term;
+ const float du_dt_j = sph_du_term_j + visc_du_term;
+
+ /* Internal energy time derivative */
+ pi->u_dt += du_dt_i * mj;
+ pj->u_dt += du_dt_j * mi;
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+ pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = max(pi->force.v_sig, v_sig);
+ pj->force.v_sig = max(pj->force.v_sig, v_sig);
+}
+
+/**
+ * @brief Force interaction between two part*icles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @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* pi,
+ const struct part* pj, float a, float H) {
+
+ /* Cosmological factors entering the EoMs */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+
+ /* Recover some data */
+ // const float mi = pi->mass;
+ const float mj = pj->mass;
+ const float mi = pi->mass;
+
+ const float miui = mi * pi->u;
+ const float mjuj = mj * pj->u;
+
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+ /* Compute gradient terms */
+ const float f_ij = 1.f - (pi->force.f / mjuj);
+ const float f_ji = 1.f - (pj->force.f / miui);
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ const float wi_dr = hid_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hjd_inv * wj_dx;
+
+ /* Compute dv dot r. */
+ const 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];
+
+ /* Includes the hubble flow term; not used for du/dt */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+
+ /* Are the part*icles moving towards each others ? */
+ const float omega_ij = min(dvdr_Hubble, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Compute sound speeds and signal velocity */
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Balsara term */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
+
+ /* Construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float alpha = 0.5f * (pi->alpha + pj->alpha);
+ const float visc =
+ -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
+
+ /* Convolve with the kernel */
+ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+
+ /* SPH acceleration term */
+ const float sph_acc_term =
+ pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
+ ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
+ r_inv;
+
+ /* Assemble the acceleration */
+ const float acc = sph_acc_term + visc_acc_term;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ /* Get the time derivative for u. */
+ const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
+ pj->u * pi->u * (f_ij / pi->pressure_bar) *
+ wi_dr * dvdr * r_inv;
+
+ /* Viscosity term */
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
+
+ /* Assemble the energy equation term */
+ const float du_dt_i = sph_du_term_i + visc_du_term;
+
+ /* Internal energy time derivative */
+ pi->u_dt += du_dt_i * mj;
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = max(pi->force.v_sig, v_sig);
+}
+
+#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IACT_H */
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..1600679bc2e840d0b3b958531c279f5f29293b48
--- /dev/null
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
@@ -0,0 +1,216 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ * Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ * 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_PRESSURE_ENERGY_MORRIS_HYDRO_IO_H
+#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IO_H
+/**
+ * @file PressureEnergy/hydro_io.h
+ * @brief P-U implementation of SPH (i/o routines)
+ *
+ * The thermal variable is the internal energy (u). A simple variable
+ * viscosity term (Morris & Monaghan 1997) with a Balsara switch is
+ * implemented.
+ *
+ * No thermal conduction term is implemented.
+ *
+ * See PressureEnergy/hydro.h for references.
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
+#include "io_properties.h"
+#include "kernel_hydro.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void hydro_read_particles(struct part* parts,
+ struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 8;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, parts, v);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ parts, mass);
+ list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, h);
+ list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+ list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
+ UNIT_CONV_ACCELERATION, parts, a_hydro);
+ list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_DENSITY, parts, rho);
+}
+
+INLINE static void convert_u(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_internal_energy(p);
+}
+
+INLINE static void convert_S(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_entropy(p);
+}
+
+INLINE static void convert_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_pressure(p);
+}
+
+INLINE static void convert_part_pos(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, double* ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = p->x[0];
+ ret[1] = p->x[1];
+ ret[2] = p->x[2];
+ }
+}
+
+INLINE static void convert_part_vel(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a_inv;
+ ret[1] *= cosmo->a_inv;
+ ret[2] *= cosmo->a_inv;
+}
+
+INLINE static void convert_part_potential(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, float* ret) {
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+INLINE static void hydro_write_particles(const struct part* parts,
+ const struct xpart* xparts,
+ struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 11;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
+ list[2] =
+ io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
+ list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ parts, h);
+ list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS,
+ parts, xparts, convert_u);
+ list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] =
+ io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+ list[7] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
+ parts, pressure_bar);
+ list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+ UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+ parts, xparts, convert_S);
+ list[9] = io_make_output_field("Viscosity", FLOAT, 1, UNIT_CONV_NO_UNITS,
+ parts, alpha);
+ list[10] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+INLINE static void hydro_write_flavour(hid_t h_grpsph) {
+
+ /* Viscosity and thermal conduction */
+ /* Nothing in this minimal model... */
+ io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
+ io_write_attribute_s(h_grpsph, "Viscosity Model",
+ "Variable viscosity as in Morris and Monaghan (1997)");
+
+ /* Time integration properties */
+ io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
+ const_max_u_change);
+}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+INLINE static int writeEntropyFlag(void) { return 0; }
+
+#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IO_H */
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..da6391236811e2a907281c3db05462bb57602fe0
--- /dev/null
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
@@ -0,0 +1,187 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
+ * Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ * 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_PRESSURE_ENERGY_MORRIS_HYDRO_PART_H
+#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_PART_H
+/**
+ * @file PressureEnergy/hydro_part.h
+ * @brief P-U implementation of SPH (Particle definition)
+ *
+ * The thermal variable is the internal energy (u). A simple variable
+ * viscosity term (Morris & Monaghan 1997) with a Balsara switch is
+ * implemented.
+ *
+ * No thermal conduction term is implemented.
+ *
+ * See PressureEnergy/hydro.h for references.
+ */
+
+#include "chemistry_struct.h"
+#include "cooling_struct.h"
+
+/**
+ * @brief Particle fields not needed during the SPH loops over neighbours.
+ *
+ * This structure contains the particle fields that are not used in the
+ * density or force loops. Quantities should be used in the kick, drift and
+ * potentially ghost tasks only.
+ */
+struct xpart {
+
+ /*! Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /*! Offset between the current position and position at the last sort. */
+ float x_diff_sort[3];
+
+ /*! Velocity at the last full step. */
+ float v_full[3];
+
+ /*! Gravitational acceleration at the last full step. */
+ float a_grav[3];
+
+ /*! Internal energy at the last full step. */
+ float u_full;
+
+ /*! Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
+
+} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Particle fields for the SPH particles
+ *
+ * The density and force substructures are used to contain variables only used
+ * within the density and force loops over neighbours. All more permanent
+ * variables should be declared in the main part of the part structure,
+ */
+struct part {
+
+ /*! Particle unique ID. */
+ long long id;
+
+ /*! Pointer to corresponding gravity part. */
+ struct gpart* gpart;
+
+ /*! Particle position. */
+ double x[3];
+
+ /*! Particle predicted velocity. */
+ float v[3];
+
+ /*! Particle acceleration. */
+ float a_hydro[3];
+
+ /*! Particle mass. */
+ float mass;
+
+ /*! Particle smoothing length. */
+ float h;
+
+ /*! Particle internal energy. */
+ float u;
+
+ /*! Time derivative of the internal energy. */
+ float u_dt;
+
+ /*! Particle density. */
+ float rho;
+
+ /*! Particle pressure (weighted) */
+ float pressure_bar;
+
+ /*! Artificial viscosity */
+ float alpha;
+
+ /* Store density/force specific stuff. */
+ union {
+
+ /**
+ * @brief Structure for the variables only used in the density loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the density
+ * loop over neighbours and the ghost task.
+ */
+ struct {
+
+ /*! Neighbour number count. */
+ float wcount;
+
+ /*! Derivative of the neighbour number with respect to h. */
+ float wcount_dh;
+
+ /*! Derivative of density with respect to h */
+ float rho_dh;
+
+ /*! Derivative of the weighted pressure with respect to h */
+ float pressure_bar_dh;
+
+ /*! Particle velocity curl. */
+ float rot_v[3];
+
+ /*! Particle velocity divergence. */
+ float div_v;
+ } density;
+
+ /**
+ * @brief Structure for the variables only used in the force loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the force
+ * loop over neighbours and the ghost, drift and kick tasks.
+ */
+ struct {
+
+ /*! "Grad h" term -- only partial in P-U */
+ float f;
+
+ /*! Particle soundspeed. */
+ float soundspeed;
+
+ /*! Particle signal velocity */
+ float v_sig;
+
+ /*! Time derivative of smoothing length */
+ float h_dt;
+
+ /*! Balsara switch */
+ float balsara;
+ } force;
+ };
+
+ /* Chemistry information */
+ struct chemistry_part_data chemistry_data;
+
+ /*! Time-step length */
+ timebin_t time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_PART_H */
diff --git a/src/hydro_io.h b/src/hydro_io.h
index b6e0c36cc7415a1f628a109795aa98b4f583036c..1a2d6319b7caf6c09b9af406cbdd323f27607791 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -31,6 +31,8 @@
#include "./hydro/PressureEntropy/hydro_io.h"
#elif defined(HOPKINS_PU_SPH)
#include "./hydro/PressureEnergy/hydro_io.h"
+#elif defined(HOPKINS_PU_SPH_MONAGHAN)
+#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_io.h"
#elif defined(GIZMO_MFV_SPH)
diff --git a/src/part.h b/src/part.h
index 625b928106cc18c219a250bf3ddac9bd36ddd3d2..64babf4a37696d7cb49b4804ee77773b4e1981fc 100644
--- a/src/part.h
+++ b/src/part.h
@@ -54,6 +54,9 @@
#elif defined(HOPKINS_PU_SPH)
#include "./hydro/PressureEnergy/hydro_part.h"
#define hydro_need_extra_init_loop 0
+#elif defined(HOPKINS_PU_SPH_MONAGHAN)
+#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h"
+#define hydro_need_extra_init_loop 0
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_part.h"
#define hydro_need_extra_init_loop 0
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 93a85bea87eb1b7c204f0cf9e6ea37ecea6d18f4..cfe7e18b07b30c309531642a769241841de3a593 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -111,7 +111,8 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
part->entropy = pressure / pow_gamma(density);
#elif defined(DEFAULT_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
-#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
+ defined(HOPKINS_PU_SPH_MONAGHAN)
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(PLANETARY_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
@@ -402,7 +403,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
hydro_get_comoving_density(&main_cell->hydro.parts[pid]),
#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) || \
defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) || \
- defined(HOPKINS_PU_SPH)
+ defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
0.f,
#else
main_cell->hydro.parts[pid].density.div_v,
@@ -422,7 +423,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
#elif defined(DEFAULT_SPH)
main_cell->hydro.parts[pid].force.v_sig, 0.f,
main_cell->hydro.parts[pid].force.u_dt
-#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
+ defined(HOPKINS_PU_SPH_MONAGHAN)
main_cell->hydro.parts[pid].force.v_sig, 0.f,
main_cell->hydro.parts[pid].u_dt
#else
diff --git a/tests/test27cells.c b/tests/test27cells.c
index ba98d91a4250865f301bd96594b1364d95034bbb..6930638efe28aafab1c614e8ce71b353eb0c5b8f 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -275,7 +275,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
hydro_get_comoving_density(&main_cell->hydro.parts[pid]),
#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
0.f,
-#elif defined(HOPKINS_PU_SPH)
+#elif defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
main_cell->hydro.parts[pid].density.pressure_bar_dh,
#else
main_cell->hydro.parts[pid].density.rho_dh,
@@ -283,7 +283,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->hydro.parts[pid].density.wcount,
main_cell->hydro.parts[pid].density.wcount_dh,
#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH) || \
- defined(HOPKINS_PU_SPH)
+ defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
main_cell->hydro.parts[pid].density.div_v,
main_cell->hydro.parts[pid].density.rot_v[0],
main_cell->hydro.parts[pid].density.rot_v[1],
diff --git a/tests/testActivePair.c b/tests/testActivePair.c
index dd7b4c6277142f4e074fc154aac455c355200816..a4eb04330fba7c984333fc7c705235223810ec4d 100644
--- a/tests/testActivePair.c
+++ b/tests/testActivePair.c
@@ -111,7 +111,8 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
/* Set the thermodynamic variable */
#if defined(GADGET2_SPH)
part->entropy = 1.f;
-#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
+ defined(HOPKINS_PU_SPH_MONAGHAN)
part->u = 1.f;
#elif defined(HOPKINS_PE_SPH)
part->entropy = 1.f;
@@ -211,7 +212,7 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo,
p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
p->density.wcount_dh = 0.f;
#endif /* PRESSURE-ENTROPY */
-#ifdef HOPKINS_PU_SPH
+#if defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
p->rho = 1.f;
p->pressure_bar = 0.6666666;
p->density.rho_dh = 0.f;
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 19eb44fa79759f50da2611cd764a1b59bd1dc579..0241a46634bdd12fb8a89bc3912c3b160198c389 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -113,7 +113,7 @@ void prepare_force(struct part *parts, size_t count) {
#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) && \
!defined(MINIMAL_SPH) && !defined(PLANETARY_SPH) && \
- !defined(HOPKINS_PU_SPH)
+ !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
@@ -153,7 +153,8 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
p->force.v_sig, p->entropy_dt, 0.f
#elif defined(DEFAULT_SPH)
p->force.v_sig, 0.f, p->force.u_dt
-#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
+ defined(HOPKINS_PU_SPH_MONAGHAN)
p->force.v_sig, 0.f, p->u_dt
#else
0.f, 0.f, 0.f
@@ -557,7 +558,7 @@ void test_force_interactions(struct part test_part, struct part *parts,
vizq[i] = pi_vec.v[2];
rhoiq[i] = pi_vec.rho;
grad_hiq[i] = pi_vec.force.f;
-#if !defined(HOPKINS_PU_SPH)
+#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN)
pOrhoi2q[i] = pi_vec.force.P_over_rho2;
#endif
balsaraiq[i] = pi_vec.force.balsara;
@@ -570,7 +571,7 @@ void test_force_interactions(struct part test_part, struct part *parts,
vjzq[i] = pj_vec[i].v[2];
rhojq[i] = pj_vec[i].rho;
grad_hjq[i] = pj_vec[i].force.f;
-#if !defined(HOPKINS_PU_SPH)
+#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN)
pOrhoj2q[i] = pj_vec[i].force.P_over_rho2;
#endif
balsarajq[i] = pj_vec[i].force.balsara;
@@ -652,7 +653,7 @@ void test_force_interactions(struct part test_part, struct part *parts,
VEC_HADD(a_hydro_zSum, piq[0]->a_hydro[2]);
VEC_HADD(h_dtSum, piq[0]->force.h_dt);
VEC_HMAX(v_sigSum, piq[0]->force.v_sig);
-#if !defined(HOPKINS_PU_SPH)
+#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN)
VEC_HADD(entropy_dtSum, piq[0]->entropy_dt);
#endif
diff --git a/theory/Cosmology/artificialvisc.tex b/theory/Cosmology/artificialvisc.tex
new file mode 100644
index 0000000000000000000000000000000000000000..55cbe2756714f875c2a9c52d7afae49499b4328b
--- /dev/null
+++ b/theory/Cosmology/artificialvisc.tex
@@ -0,0 +1,75 @@
+\subsection{Cosmological factors for properties entering the artificial viscosity}
+\label{ssec:artificialvisc}
+
+There are multiple properties that enter into the more complex artificial
+viscosity schemes, such as those by \citet{Morris1997} (henceforth M\&M) and
+\citet{Cullen2010} (henceforth C\&D).
+
+\subsubsection{M\&M basic scheme}
+\label{sssec:mandm}
+
+This relies on the velocity divergence as a shock indicator, i.e. the property
+$\nabla \cdot \mathbf{v}$. The interpretation of this is the velocity divergence of
+the fluid overall, i.e. the physical velocity divergence. Starting with
+\begin{equation}
+\mathbf{v}_p = a \dot{\mathbf{r}}' + \dot{a} \mathbf{r}', \nonumber
+\end{equation}
+with the divergence,
+\begin{equation}
+\nabla \cdot \mathbf{v}_p =
+ \nabla \cdot \left(a \dot{\mathbf{r}}'\right) +
+ \nabla \cdot \left(\dot{a} \mathbf{r}'\right). \nonumber
+\end{equation}
+The quantity on the left is the one that we want to enter the source term for the
+artificial viscosity. Transforming to the co-moving derivative on the right hand side
+to enable it to be calculated in the code,
+\begin{equation}
+\nabla \cdot \mathbf{v}_p =
+ \nabla' \cdot \dot{\mathbf{r}}' + n_d H(a),
+\label{eqn:divvwithcomovingcoordinates}
+\end{equation}
+with $n_d$ the number of spatial dimensions, and the final transformation
+being the one to internal code velocity units,
+\begin{equation}
+\nabla \cdot \mathbf{v}_p =
+ \frac{1}{a^2} \nabla' \cdot \mathbf{v}' + n_d H(a).
+\label{eqn:divvcodeunits}
+\end{equation}
+We note that there is no similar hubble flow term in the expression for
+$\nabla \times \mathbf{v}_p$.
+
+In some more complex schemes, such as the one presented by \cite{Cullen2010},
+the time differential of the velocity divergence is used as a way to differentiate
+the pre- and post-shock region.
+
+Building on the above, we take the time differential of both sides,
+\begin{equation}
+ \frac{{\mathrm d}}{{\mathrm d} t} \nabla \cdot \mathbf{v}_p =
+ \frac{{\mathrm d}}{{\mathrm d} t} \left(
+ \frac{1}{a^2} \nabla' \cdot \mathbf{v}' + n_d H(a)
+ \right).
+ \nonumber
+\end{equation}
+Collecting the factors, we see
+\begin{align}
+ \frac{{\mathrm d}}{{\mathrm d} t} \nabla \cdot \mathbf{v}_p =
+ \frac{1}{a^2} &\left(
+ \frac{{\mathrm d}}{{\mathrm d} t} \nabla ' \cdot \mathbf{v}' -
+ 2H(a) \nabla' \cdot \mathbf{v}'
+ \right) \\
+ + n_d &\left(
+ \frac{\ddot{a}}{a} - \frac{\dot{a}}{a^2}
+ \right).
+ \label{eqn:divvdtcodeunits}
+\end{align}
+This looks like quite a mess, but in most cases we calculate this implicitly
+from the velocity divergence itself, and so we do not actually need to take
+into account these factors; i.e. we actually calculate
+\begin{equation}
+ \frac{\mathrm d}{{\mathrm d} t} \nabla \cdot \mathbf{v}_p =
+ \frac{
+ \nabla \cdot \mathbf{v}_p (t + {\mathrm d}t) - \nabla \cdot \mathbf{v}_p (t)
+ }{dt},
+ \label{eqn:divvdtcodeunitsimplicit}
+\end{equation}
+meaning that the above is taken into account self-consistently.
\ No newline at end of file
diff --git a/theory/Cosmology/bibliography.bib b/theory/Cosmology/bibliography.bib
index 6979bf7dd23bdb8543ac8752c12432837480d4ed..94ac688f7043e1e911561e460b1fc607ce1e1421 100644
--- a/theory/Cosmology/bibliography.bib
+++ b/theory/Cosmology/bibliography.bib
@@ -137,4 +137,29 @@ issn = "0021-9991",
doi = "https://doi.org/10.1006/jcph.1997.5732",
url = "http://www.sciencedirect.com/science/article/pii/S0021999197957326",
author = "J.J. Monaghan"
-}
\ No newline at end of file
+}
+
+@article{Cullen2010,
+author = {Cullen, Lee and Dehnen, Walter},
+title = {{Inviscid smoothed particle hydrodynamics}},
+journal = {Monthly Notices of the Royal Astronomical Society},
+year = {2010},
+volume = {408},
+number = {2},
+pages = {669--683},
+month = oct,
+annote = {14 pages (15 in arXiv), 15 figures, accepted for publication in MNRAS}
+}
+
+@article{Morris1997,
+author = {Morris, J P and Monaghan, J J},
+title = {{A Switch to Reduce SPH Viscosity}},
+journal = {Journal of Computational Physics},
+year = {1997},
+volume = {136},
+number = {1},
+pages = {41--50},
+month = sep
+}
+
+
diff --git a/theory/Cosmology/cosmology_standalone.tex b/theory/Cosmology/cosmology_standalone.tex
index 31a96d3a002aae423b2f8e16ef3044e357fdea6a..56f31297e55cef298b20204db0ee35867ad1789b 100644
--- a/theory/Cosmology/cosmology_standalone.tex
+++ b/theory/Cosmology/cosmology_standalone.tex
@@ -42,6 +42,8 @@ Making cosmology great again.
\input{timesteps}
+\input{artificialvisc}
+
\bibliographystyle{mnras}
\bibliography{./bibliography.bib}
diff --git a/theory/SPH/Flavours/bibliography.bib b/theory/SPH/Flavours/bibliography.bib
index 2bc11dacca90fe03d05c2e847503105d80eb1317..02ebed25a407ae5adba87d9f46d3f004bf9fbae2 100644
--- a/theory/SPH/Flavours/bibliography.bib
+++ b/theory/SPH/Flavours/bibliography.bib
@@ -97,4 +97,17 @@ archivePrefix = "arXiv",
+@article{Morris1997,
+abstract = {Smoothed particle hydrodynamics is a Lagrangian particle method for fluid dynamics which simulates shocks by using an artificial viscosity. Unlike Eulerian methods it is not convenient to reduce the effects of viscosity by means of switches based on spatial gradients. In this paper we introduce the idea of time-varying coefficients which fits more naturally with a particle formulation. Each particle has a viscosity parameter which evolves according to a simple source and decay equation. The source causes the parameter to grow when the particle enters a shock and the decay term causes it to decay to a small value beyond the shock. Tests on one-dimensional shocks and a two-dimensional shock-bubble interaction confirm that the method gives good results. {\textcopyright} 1997 Academic Press.},
+author = {Morris, J. P. and Monaghan, J. J.},
+doi = {10.1006/jcph.1997.5690},
+isbn = {0021-9991},
+issn = {00219991},
+journal = {Journal of Computational Physics},
+number = {1},
+pages = {41--50},
+title = {{A switch to reduce SPH viscosity}},
+volume = {136},
+year = {1997}
+}
diff --git a/theory/SPH/Flavours/sph_flavours.tex b/theory/SPH/Flavours/sph_flavours.tex
index 5d62af3aab777e66f0b33b89e861d2b21e10b38c..81ac2153ed23f29f14345e0377774548420c84c9 100644
--- a/theory/SPH/Flavours/sph_flavours.tex
+++ b/theory/SPH/Flavours/sph_flavours.tex
@@ -590,6 +590,42 @@ both sides, such that
%##############################################################################
+\subsection{Variable artificial viscosity}
+
+Here we consider a modified version of the Pressure-Energy scheme described
+above but one that uses a variable artificial viscosity. The prescription used
+in this scheme was originally introduced by \citet{Morris1997} and is almost
+identical to the above equations, but tracks an individual viscosity paramaeter
+$\alpha_i$ for each particle. This viscosity is then updated each time-step to
+a more appropraite value. The hope is that the artificial viscosity will be
+high in regions that contain shocks, but as low as possible in regions where it
+is uneccesary such as shear flows. This is already accomplished somewhat with
+the inclusion of a \citet{Balsara1995} switch, but a fixed $\alpha$ still leads
+to spurious transport of angular momentum and vorticity.
+
+The equation governing the growth of the viscosity is
+\begin{align}
+ \frac{\mathrm{d} \alpha_i}
+ {\mathrm{d} t} =
+ - (\alpha_i - \alpha_{\rm min}) \ell \frac{c_{s, i}}{h},
+ \label{eq:sph:pu:alphadt}
+\end{align}
+with $\alpha_{\rm min}=0.1$ the minimal artificial viscosity parameter, and
+$\ell=0.1$ the viscosity ``length" that governs how quickly the viscosity
+decays. This equation is solved implicitly in a similar way to
+$\mathrm{d}\mathbf{v}/ \mathrm{d}t$ and $\mathrm{d}u/\mathrm{d}t$ - i.e.
+$\alpha_{i} (t+\Delta t_i) = \alpha_{i}(t) + \dot{\alpha}_i \Delta t_i$.
+
+To ensure that the scheme is conservative, the viscosity coefficients must be
+combined in a fully conservative way; this is performed by taking the mean
+viscosity parameter of the two particles that are being interacted, such that
+\begin{align}
+ \alpha_{ij} = \frac{\alpha_i + \alpha_j}{2}.
+\end{align}
+The rest of the artificial viscosity implementation, including the
+\citet{Balsara1995} switch, is the same - just with $\alpha \rightarrow
+\alpha_{ij}$.
+
\subsection{Anarchy SPH}
Dalla Vecchia (\textit{in prep.}), also described in section 2.2.2 of
\cite{Schaller2015}.\\