From 4d69cea150cba01d93042612aec2d1a8db4eeecb Mon Sep 17 00:00:00 2001
From: Josh Borrow <joshua.borrow@durham.ac.uk>
Date: Wed, 7 Nov 2018 18:02:19 +0000
Subject: [PATCH] Add P-U with M&M Switch
---
configure.ac | 5 +-
examples/KelvinHelmholtz_2D/makeMovie.py | 1 +
examples/SodShock_1D/plotSolution.py | 31 +-
src/debug.c | 2 +
src/hydro.h | 6 +
src/hydro/PressureEnergy/hydro.h | 6 +-
src/hydro/PressureEnergy/hydro_iact.h | 4 +-
.../PressureEnergyMorrisMonaghanAV/hydro.h | 766 ++++++++++++++++++
.../hydro_debug.h | 47 ++
.../hydro_iact.h | 427 ++++++++++
.../PressureEnergyMorrisMonaghanAV/hydro_io.h | 216 +++++
.../hydro_part.h | 187 +++++
src/hydro_io.h | 2 +
src/part.h | 3 +
tests/test125cells.c | 8 +-
tests/test27cells.c | 4 +-
tests/testActivePair.c | 5 +-
tests/testInteractions.c | 11 +-
theory/Cosmology/artificialvisc.tex | 75 ++
theory/Cosmology/bibliography.bib | 27 +-
theory/Cosmology/cosmology_standalone.tex | 2 +
theory/SPH/Flavours/bibliography.bib | 13 +
theory/SPH/Flavours/sph_flavours.tex | 36 +
23 files changed, 1858 insertions(+), 26 deletions(-)
create mode 100644 src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
create mode 100644 src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
create mode 100644 src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
create mode 100644 src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
create mode 100644 src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
create mode 100644 theory/Cosmology/artificialvisc.tex
diff --git a/configure.ac b/configure.ac
index ab2ed97305..b91303ef1e 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 84fe99106b..a52784891a 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 e001a8d87a..0149d66a0c 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 a980ef9457..809d7048c4 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 b3716996cc..15c45c1dcf 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 3274c1b1b5..4c3cc5c1c5 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 2ed7fe8cb8..4146e61a53 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 0000000000..060694a6af
--- /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 0000000000..ead5fcc0c8
--- /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 0000000000..747fca714c
--- /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 0000000000..1600679bc2
--- /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 0000000000..da63912368
--- /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 b6e0c36cc7..1a2d6319b7 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 625b928106..64babf4a37 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 93a85bea87..cfe7e18b07 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 ba98d91a42..6930638efe 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 dd7b4c6277..a4eb04330f 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 19eb44fa79..0241a46634 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 0000000000..55cbe27567
--- /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 6979bf7dd2..94ac688f70 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 31a96d3a00..56f31297e5 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 2bc11dacca..02ebed25a4 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 5d62af3aab..81ac2153ed 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}.\\
--
GitLab