diff --git a/.gitignore b/.gitignore
index b6f80a5ad0a549843b6261c439bcc2a704b43124..8c3ede8f3125f5024c1fe01a405024d1cf5a7f19 100644
--- a/.gitignore
+++ b/.gitignore
@@ -26,17 +26,20 @@ examples/swift_mpi
examples/*/*.xmf
examples/*/*.h5
examples/*/*.png
+examples/*/*.mp4
examples/*/*.txt
examples/*/*.dot
+examples/*/restart/*
examples/*/used_parameters.yml
examples/*/*/*.xmf
examples/*/*/*.png
+examples/*/*/*.mp4
examples/*/*/*.txt
examples/*/*/*.dot
examples/*/*/*.rst
examples/*/*/*.hdf5
examples/*/snapshots*
-examples/*/restart
+examples/*/restart/*
examples/*/*/used_parameters.yml
examples/*/*.mpg
examples/*/gravity_checks_*.dat
@@ -298,3 +301,6 @@ sympy-plots-for-*.tex/
# macOS
*.DS_Store
+
+#ctags
+*tags
diff --git a/configure.ac b/configure.ac
index c4134686dc5bca99082981a19029063e3a4db1b1..56018272debbb9e9ccb2f425eff45b97c358aec9 100644
--- a/configure.ac
+++ b/configure.ac
@@ -855,7 +855,7 @@ fi
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
- [Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo-mfv, gizmo-mfm, shadowfax, minimal-multi-mat, debug default: gadget2@:>@]
+ [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, minimal-multi-mat, debug default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
@@ -876,9 +876,12 @@ case "$with_hydro" in
minimal)
AC_DEFINE([MINIMAL_SPH], [1], [Minimal SPH])
;;
- hopkins)
+ pressure-entropy)
AC_DEFINE([HOPKINS_PE_SPH], [1], [Pressure-Entropy SPH])
;;
+ pressure-energy)
+ AC_DEFINE([HOPKINS_PU_SPH], [1], [Pressure-Energy SPH])
+ ;;
default)
AC_DEFINE([DEFAULT_SPH], [1], [Default SPH])
;;
diff --git a/examples/KelvinHelmholtz_2D/kelvinHelmholtz.yml b/examples/KelvinHelmholtz_2D/kelvinHelmholtz.yml
index 19f036c364c4e9f40580e906e2840cad6df378a0..ccc7526b391374a4da0883f6615a65c7b93a0948 100644
--- a/examples/KelvinHelmholtz_2D/kelvinHelmholtz.yml
+++ b/examples/KelvinHelmholtz_2D/kelvinHelmholtz.yml
@@ -8,8 +8,8 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
- time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1.5 # The end time of the simulation (in internal units).
+ time_begin: 0.0 # The starting time of the simulation (in internal units).
+ time_end: 4.5 # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
@@ -17,7 +17,7 @@ TimeIntegration:
Snapshots:
basename: kelvinHelmholtz # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.25 # Time difference between consecutive outputs (in internal units)
+ delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/KelvinHelmholtz_2D/makeIC.py b/examples/KelvinHelmholtz_2D/makeIC.py
index bd0f39ed90faf0d67ff4a508bff83067bf748d43..744b39de8260720521ae8e77ed5d0a12161f2b6a 100644
--- a/examples/KelvinHelmholtz_2D/makeIC.py
+++ b/examples/KelvinHelmholtz_2D/makeIC.py
@@ -24,7 +24,7 @@ import sys
# Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
# Parameters
-L2 = 128 # Particles along one edge in the low-density region
+L2 = 256 # Particles along one edge in the low-density region
gamma = 5./3. # Gas adiabatic index
P1 = 2.5 # Central region pressure
P2 = 2.5 # Outskirts pressure
diff --git a/examples/KelvinHelmholtz_2D/makeMovie.py b/examples/KelvinHelmholtz_2D/makeMovie.py
new file mode 100644
index 0000000000000000000000000000000000000000..84fe99106bf607830e89b6aa663135b48b6c0744
--- /dev/null
+++ b/examples/KelvinHelmholtz_2D/makeMovie.py
@@ -0,0 +1,112 @@
+"""
+Makes a movie of the KH 2D data.
+
+You will need to run your movie with far higher time-resolution than usual to
+get a nice movie; around 450 snapshots over 6s is required.
+
+Edit this file near the bottom with the number of snaps you have.
+
+Written by Josh Borrow (joshua.borrow@durham.ac.uk)
+"""
+
+import os
+import h5py as h5
+import numpy as np
+import scipy.interpolate as si
+
+
+def load_and_extract(filename):
+ """
+ Load the data and extract relevant info.
+ """
+
+ with h5.File(filename, "r") as f:
+ x, y, _ = f["PartType0/Coordinates"][...].T
+ density = f["PartType0/Density"][...]
+
+ return x, y, density
+
+
+def make_plot(filename, array, nx, ny, dx, dy):
+ """
+ Load the data and plop it on the grid using nearest
+ neighbour searching for finding the 'correct' value of
+ the density.
+ """
+
+ data_x, data_y, density = load_and_extract(filename)
+
+ # Make the grid
+ x = np.linspace(*dx, nx)
+ y = np.linspace(*dy, ny)
+
+ xv, yv = np.meshgrid(x, y)
+
+ mesh = si.griddata((data_x, data_y), density, (xv, yv), method="nearest")
+
+ array.set_array(mesh)
+
+ return array,
+
+
+def frame(n, *args):
+ """
+ Make a single frame. Requires the global variables plot and dpi.
+ """
+
+ global plot, dpi
+
+ fn = "{}_{:04d}.hdf5".format(filename, n)
+
+ return make_plot(fn, plot, dpi, dpi, (0, 1), (0, 1))
+
+
+if __name__ == "__main__":
+ import matplotlib
+ matplotlib.use("Agg")
+
+ from tqdm import tqdm
+ from matplotlib.animation import FuncAnimation
+ from scipy.stats import gaussian_kde
+
+ import matplotlib.pyplot as plt
+
+ filename = "kelvinhelmholtz"
+ dpi = 512
+
+
+ # Look for the number of files in the directory.
+ i = 0
+ while True:
+ if os.path.isfile("{}_{:04d}.hdf5".format(filename, i)):
+ i += 1
+ else:
+ break
+
+ if i > 10000:
+ raise FileNotFoundError(
+ "Could not find the snapshots in the directory")
+
+ frames = tqdm(np.arange(0, i))
+
+ # Creation of first frame
+ fig, ax = plt.subplots(1, 1, figsize=(1, 1), frameon=False)
+
+ data_x, data_y, density = load_and_extract("kelvinhelmholtz_0000.hdf5")
+
+ x = np.linspace(0, 1, dpi)
+ y = np.linspace(0, 1, dpi)
+ xv, yv = np.meshgrid(x, y)
+
+ mesh = si.griddata((data_x, data_y), density, (xv, yv), method="nearest")
+
+ # Global variable for set_array
+ plot = ax.imshow(mesh, extent=[0, 1, 0, 1], animated=True, interpolation="none")
+
+ anim = FuncAnimation(fig, frame, frames, interval=40, blit=False)
+
+ # Remove all whitespace
+ fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None)
+
+ # Actually make the movie
+ anim.save("khmovie.mp4", dpi=dpi, bitrate=4096)
diff --git a/examples/KelvinHelmholtz_2D/run.sh b/examples/KelvinHelmholtz_2D/run.sh
index 3d83bcc6dd7c226885ed83c3188f3177c4807154..dbb39caf383279dbc71c2baa125499d115538654 100755
--- a/examples/KelvinHelmholtz_2D/run.sh
+++ b/examples/KelvinHelmholtz_2D/run.sh
@@ -8,7 +8,8 @@ then
fi
# Run SWIFT
-../swift -s -t 1 kelvinHelmholtz.yml 2>&1 | tee output.log
+../swift -s -t 4 kelvinHelmholtz.yml 2>&1 | tee output.log
# Plot the solution
python plotSolution.py 6
+python makeMovie.py
diff --git a/src/debug.c b/src/debug.c
index 3e43d05f13b521ea23e76508922d7a3244c50ab1..93d14952f523be5f1d1fa90484e9e7951f8e3f6e 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -46,6 +46,8 @@
#include "./hydro/Gadget2/hydro_debug.h"
#elif defined(HOPKINS_PE_SPH)
#include "./hydro/PressureEntropy/hydro_debug.h"
+#elif defined(HOPKINS_PU_SPH)
+#include "./hydro/PressureEnergy/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 f7921db56f4063f8326780f411eb5f410471faa0..950f63526a1590fa0fdcf2bfb5e650a2dfe14431 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -41,6 +41,10 @@
#include "./hydro/PressureEntropy/hydro.h"
#include "./hydro/PressureEntropy/hydro_iact.h"
#define SPH_IMPLEMENTATION "Pressure-Entropy SPH (Hopkins 2013)"
+#elif defined(HOPKINS_PU_SPH)
+#include "./hydro/PressureEnergy/hydro.h"
+#include "./hydro/PressureEnergy/hydro_iact.h"
+#define SPH_IMPLEMENTATION "Pressure-Energy SPH (Hopkins 2013)"
#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
new file mode 100644
index 0000000000000000000000000000000000000000..1e0d8208e82f7f48691d1df7603a6b02d1471c12
--- /dev/null
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -0,0 +1,657 @@
+/*******************************************************************************
+ * 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_HYDRO_H
+#define SWIFT_PRESSURE_ENERGY_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 constant
+ * viscosity term 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 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;
+}
+
+/**
+ * @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 * kernel_norm * 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.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ const float fac_mu = cosmo->a_factor_mu;
+
+ /* 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 / p->h);
+
+ /* 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);
+
+ /* Update variables. */
+ 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 = 0.0f;
+}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__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).
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+ const struct cosmology *cosmo,
+ const struct hydro_props *restrict hydro_properties) {
+
+ /* 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;
+
+ /* 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
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {}
+
+/**
+ * @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_MINIMAL_HYDRO_H */
diff --git a/src/hydro/PressureEnergy/hydro_debug.h b/src/hydro/PressureEnergy/hydro_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..6324167f12726e155eeaa3359be9741aca3a1e42
--- /dev/null
+++ b/src/hydro/PressureEnergy/hydro_debug.h
@@ -0,0 +1,43 @@
+/*******************************************************************************
+ * 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_HYDRO_DEBUG_H
+#define SWIFT_PRESSURE_ENERGY_HYDRO_DEBUG_H
+/**
+ * @file PressureEnergy/hydro_debug.h
+ * @brief P-U conservative implementation of SPH (Debugging routines)
+ */
+
+__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\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);
+}
+
+#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..65c46a55554d4a8f09b32bb6eb1deb1fdcfc932a
--- /dev/null
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -0,0 +1,418 @@
+/*******************************************************************************
+ * 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_MINIMAL_HYDRO_IACT_H
+#define SWIFT_MINIMAL_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 constant
+ * viscosity term 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] + a2_Hubble * r2;
+
+ /* Are the part*icles moving towards each others ? */
+ const float omega_ij = min(dvdr, 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 visc = -0.25f * const_viscosity_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;
+
+ /* 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] + a2_Hubble * r2;
+
+ /* Are the part*icles moving towards each others ? */
+ const float omega_ij = min(dvdr, 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 visc = -0.25f * const_viscosity_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;
+
+ /* Assemble the energy equation term */
+ const float du_dt_i = sph_du_term_i + visc_du_term;
+
+ /* Internal energy time derivatibe */
+ 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_MINIMAL_HYDRO_IACT_H */
diff --git a/src/hydro/PressureEnergy/hydro_io.h b/src/hydro/PressureEnergy/hydro_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..776e7653ac3152e1594f25a33796a470dfcf69d3
--- /dev/null
+++ b/src/hydro/PressureEnergy/hydro_io.h
@@ -0,0 +1,196 @@
+/*******************************************************************************
+ * 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_HYDRO_IACT_H
+#define SWIFT_PRESSURE_ENERGY_HYDRO_IACT_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 constant
+ * viscosity term 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.
+ */
+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);
+}
+
+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);
+}
+
+void convert_S(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_entropy(p);
+}
+
+void convert_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_pressure(p);
+}
+
+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];
+ }
+}
+
+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->a2_inv;
+ ret[1] *= cosmo->a2_inv;
+ ret[2] *= cosmo->a2_inv;
+}
+
+/**
+ * @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.
+ */
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
+
+ *num_fields = 9;
+
+ /* 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);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+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",
+ "Minimal treatment as in Monaghan (1992)");
+
+ /* 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.
+ */
+int writeEntropyFlag() { return 0; }
+
+#endif /* SWIFT_MINIMAL_HYDRO_IO_H */
diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..bc7d14b612556dc722ecca67dd6ce823192e00f0
--- /dev/null
+++ b/src/hydro/PressureEnergy/hydro_part.h
@@ -0,0 +1,183 @@
+/*******************************************************************************
+ * 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_HYDRO_PART_H
+#define SWIFT_PRESSURE_ENERGY_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 constant
+ * viscosity term 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;
+
+ /* 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_MINIMAL_HYDRO_PART_H */
diff --git a/src/hydro_io.h b/src/hydro_io.h
index 9b5641deda4bedc7c47f6007021e5a2563898ebf..d752bb8bc03f619fe759fc8f5de32a01b3a61abe 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -29,6 +29,8 @@
#include "./hydro/Gadget2/hydro_io.h"
#elif defined(HOPKINS_PE_SPH)
#include "./hydro/PressureEntropy/hydro_io.h"
+#elif defined(HOPKINS_PU_SPH)
+#include "./hydro/PressureEnergy/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 a9a79d3189b48f436f4cd58c7707565f2f348398..e6750ea864bf3785df0b4ebe011e0ad741d7b5c7 100644
--- a/src/part.h
+++ b/src/part.h
@@ -51,6 +51,9 @@
#elif defined(HOPKINS_PE_SPH)
#include "./hydro/PressureEntropy/hydro_part.h"
#define hydro_need_extra_init_loop 1
+#elif defined(HOPKINS_PU_SPH)
+#include "./hydro/PressureEnergy/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/Makefile.am b/tests/Makefile.am
index 77aad7beb58c8ab80da04c28fc68b64afba3d99a..891eef3f518f83c17b66623e3dac1832512d31f3 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -51,6 +51,9 @@ testReading_SOURCES = testReading.c
testSymmetry_SOURCES = testSymmetry.c
+# Added because of issues using memcmp on clang 4.x
+testSymmetry_CFLAGS = $(AM_CFLAGS) -fno-builtin-memcmp
+
testTimeIntegration_SOURCES = testTimeIntegration.c
testSPHStep_SOURCES = testSPHStep.c
diff --git a/tests/difffloat.py b/tests/difffloat.py
index ddcf7bcb29758afa3429dea8bcf50e1c5c0477dc..55295309e1309ab35542a50f7dd5b15c86c1508f 100644
--- a/tests/difffloat.py
+++ b/tests/difffloat.py
@@ -57,7 +57,7 @@ if fileTol != "":
if shape(data1) != shape(data2):
- print "Non-matching array sizes in the files", file1, "and", file2, "."
+ print("Non-matching array sizes in the files", file1, "and", file2, ".")
sys.exit(1)
n_lines = shape(data1)[0]
@@ -65,18 +65,18 @@ n_columns = shape(data1)[1]
if fileTol != "":
if n_linesTol != 3:
- print "Incorrect number of lines in tolerance file '%s'."%fileTol
+ print("Incorrect number of lines in tolerance file '%s'."%fileTol)
if n_columnsTol != n_columns:
- print "Incorrect number of columns in tolerance file '%s'."%fileTol
+ print("Incorrect number of columns in tolerance file '%s'."%fileTol)
if fileTol == "":
- print "Absolute difference tolerance:", abs_tol
- print "Relative difference tolerance:", rel_tol
+ print("Absolute difference tolerance:", abs_tol)
+ print("Relative difference tolerance:", rel_tol)
absTol = ones(n_columns) * abs_tol
relTol = ones(n_columns) * rel_tol
limTol = zeros(n_columns)
else:
- print "Tolerances read from file"
+ print("Tolerances read from file")
absTol = dataTol[0,:]
relTol = dataTol[1,:]
limTol = dataTol[2,:]
@@ -85,10 +85,10 @@ n_lines_to_check = 0
if number_to_check > 0:
n_lines_to_check = number_to_check**3
n_lines_to_check = min(n_lines_to_check, n_lines)
- print "Checking the first %d particles."%n_lines_to_check
+ print("Checking the first %d particles."%n_lines_to_check)
else:
n_lines_to_check = n_lines
- print "Checking all particles in the file."
+ print("Checking all particles in the file.")
error = False
for i in range(n_lines_to_check):
@@ -103,26 +103,26 @@ for i in range(n_lines_to_check):
rel_diff = 0.
if( abs_diff > 1.1*absTol[j]):
- print "Absolute difference larger than tolerance (%e) for particle %d, column %s:"%(absTol[j], data1[i,0], part_props[j])
- print "%10s: a = %e"%("File 1", data1[i,j])
- print "%10s: b = %e"%("File 2", data2[i,j])
- print "%10s: |a-b| = %e"%("Difference", abs_diff)
- print ""
+ print("Absolute difference larger than tolerance (%e) for particle %d, column %s:"%(absTol[j], data1[i,0], part_props[j]))
+ print("%10s: a = %e"%("File 1", data1[i,j]))
+ print("%10s: b = %e"%("File 2", data2[i,j]))
+ print("%10s: |a-b| = %e"%("Difference", abs_diff))
+ print("")
error = True
if abs(data1[i,j]) + abs(data2[i,j]) < limTol[j] : continue
if( rel_diff > 1.1*relTol[j]):
- print "Relative difference larger than tolerance (%e) for particle %d, column %s:"%(relTol[j], data1[i,0], part_props[j])
- print "%10s: a = %e"%("File 1", data1[i,j])
- print "%10s: b = %e"%("File 2", data2[i,j])
- print "%10s: |a-b|/|a+b| = %e"%("Difference", rel_diff)
- print ""
+ print("Relative difference larger than tolerance (%e) for particle %d, column %s:"%(relTol[j], data1[i,0], part_props[j]))
+ print("%10s: a = %e"%("File 1", data1[i,j]))
+ print("%10s: b = %e"%("File 2", data2[i,j]))
+ print("%10s: |a-b|/|a+b| = %e"%("Difference", rel_diff))
+ print("")
error = True
if error:
exit(1)
else:
- print "No differences found"
+ print("No differences found")
exit(0)
diff --git a/tests/test125cells.c b/tests/test125cells.c
index af962a7b1cbe5e8abc3f8c7cdc9e3aeab602c0a8..a50b847308422cbf10f58c737811934979d21899 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -118,7 +118,7 @@ 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)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(MINIMAL_MULTI_MAT_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
@@ -406,7 +406,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].v[2], main_cell->parts[pid].h,
hydro_get_comoving_density(&main_cell->parts[pid]),
#if defined(MINIMAL_SPH) || defined(MINIMAL_MULTI_MAT_SPH) || \
- defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
+ defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) || \
+ defined(HOPKINS_PU_SPH)
0.f,
#else
main_cell->parts[pid].density.div_v,
@@ -423,7 +424,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
#elif defined(DEFAULT_SPH)
main_cell->parts[pid].force.v_sig, 0.f,
main_cell->parts[pid].force.u_dt
-#elif defined(MINIMAL_SPH)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH)
main_cell->parts[pid].force.v_sig, 0.f, main_cell->parts[pid].u_dt
#else
0.f, 0.f, 0.f
@@ -663,6 +664,7 @@ int main(int argc, char *argv[]) {
runner_do_sort(&runner, cells[j], 0x1FFF, 0, 0);
/* Do the density calculation */
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Initialise the particle cache. */
#ifdef WITH_VECTORIZATION
@@ -706,10 +708,13 @@ int main(int argc, char *argv[]) {
for (int j = 0; j < 27; ++j)
runner_doself1_density(&runner, inner_cells[j]);
+#endif
+
/* Ghost to finish everything on the central cells */
for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j], 0);
/* Do the force calculation */
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
#ifdef WITH_VECTORIZATION
/* Initialise the cache. */
@@ -745,6 +750,7 @@ int main(int argc, char *argv[]) {
DOSELF2(&runner, main_cell);
timings[26] += getticks() - self_tic;
+#endif
/* Finally, give a gentle kick */
runner_do_end_force(&runner, main_cell, 0);
@@ -790,17 +796,18 @@ int main(int argc, char *argv[]) {
const ticks tic = getticks();
- /* Kick the central cell */
- // runner_do_kick1(&runner, main_cell, 0);
+/* Kick the central cell */
+// runner_do_kick1(&runner, main_cell, 0);
- /* And drift it */
- // runner_do_drift_particles(&runner, main_cell, 0);
+/* And drift it */
+// runner_do_drift_particles(&runner, main_cell, 0);
- /* Initialise the particles */
- // for (int j = 0; j < 125; ++j) runner_do_drift_particles(&runner, cells[j],
- // 0);
+/* Initialise the particles */
+// for (int j = 0; j < 125; ++j) runner_do_drift_particles(&runner, cells[j],
+// 0);
- /* Do the density calculation */
+/* Do the density calculation */
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Run all the pairs (only once !)*/
for (int i = 0; i < 5; i++) {
@@ -835,10 +842,13 @@ int main(int argc, char *argv[]) {
/* And now the self-interaction for the central cells*/
for (int j = 0; j < 27; ++j) self_all_density(&runner, inner_cells[j]);
+#endif
+
/* Ghost to finish everything on the central cells */
for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j], 0);
- /* Do the force calculation */
+/* Do the force calculation */
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Do the pairs (for the central 27 cells) */
for (int i = 1; i < 4; i++) {
@@ -855,6 +865,8 @@ int main(int argc, char *argv[]) {
/* And now the self-interaction for the main cell */
self_all_force(&runner, main_cell);
+#endif
+
/* Finally, give a gentle kick */
runner_do_end_force(&runner, main_cell, 0);
// runner_do_kick2(&runner, main_cell, 0);
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 11affc56a0aebb750bd7f562df9001e60bd2084a..e60262df71f6dc455e944b82a337261a57bcc9bc 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -263,12 +263,15 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
hydro_get_comoving_density(&main_cell->parts[pid]),
#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
0.f,
+#elif defined(HOPKINS_PU_SPH)
+ main_cell->parts[pid].density.pressure_bar_dh,
#else
main_cell->parts[pid].density.rho_dh,
#endif
main_cell->parts[pid].density.wcount,
main_cell->parts[pid].density.wcount_dh,
-#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH) || \
+ defined(HOPKINS_PU_SPH)
main_cell->parts[pid].density.div_v,
main_cell->parts[pid].density.rot_v[0],
main_cell->parts[pid].density.rot_v[1],
diff --git a/tests/testActivePair.c b/tests/testActivePair.c
index cede31a876f214b286eeecc34c81b9522332e947..0453f6d5896eaa53b0f44a567d353d7d8e8fb7df 100644
--- a/tests/testActivePair.c
+++ b/tests/testActivePair.c
@@ -110,7 +110,7 @@ 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)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH)
part->u = 1.f;
#elif defined(HOPKINS_PE_SPH)
part->entropy = 1.f;
@@ -194,13 +194,13 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo) {
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
p->density.div_v = 0.f;
-#endif
+#endif /* GADGET-2 */
#ifdef MINIMAL_SPH
p->rho = 1.f;
p->density.rho_dh = 0.f;
p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
p->density.wcount_dh = 0.f;
-#endif
+#endif /* MINIMAL */
#ifdef HOPKINS_PE_SPH
p->rho = 1.f;
p->rho_bar = 1.f;
@@ -208,7 +208,15 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo) {
p->density.pressure_dh = 0.f;
p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
p->density.wcount_dh = 0.f;
-#endif
+#endif /* PRESSURE-ENTROPY */
+#ifdef HOPKINS_PU_SPH
+ p->rho = 1.f;
+ p->pressure_bar = 0.6666666;
+ p->density.rho_dh = 0.f;
+ p->density.pressure_bar_dh = 0.f;
+ p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
+ p->density.wcount_dh = 0.f;
+#endif /* PRESSURE-ENERGY */
/* And prepare for a round of force tasks. */
hydro_prepare_force(p, xp, cosmo);
diff --git a/tests/testDump.c b/tests/testDump.c
index d4a3b3c1bacdc8071a32fce6d5f1f746530e589c..fa68ef9869f2f3ac2ee790b9815e42f73976ac9f 100644
--- a/tests/testDump.c
+++ b/tests/testDump.c
@@ -99,6 +99,8 @@ int main(int argc, char *argv[]) {
#else
+#include <stdio.h>
+
int main(int argc, char *argv[]) {
printf("No posix_fallocate, not testing anything.\n");
return 0;
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index f32fcfdc166bc0344ea08893341b8947838859de..5473dc2588d66e0df2e3e3caddfc04ba3e6f7a2c 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -24,7 +24,10 @@
#include <unistd.h>
#include "swift.h"
-#ifdef WITH_VECTORIZATION
+/* Other schemes need to be added here if they are not vectorized, otherwise
+ * this test will simply not compile. */
+
+#if defined(GADGET2_SPH) && defined(WITH_VECTORIZATION)
#define array_align sizeof(float) * VEC_SIZE
#define ACC_THRESHOLD 1e-5
@@ -92,7 +95,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
-#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH)
+#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH)
p->mass = 1.0f;
#endif
}
@@ -104,8 +107,9 @@ struct part *make_particles(size_t count, double *offset, double spacing,
*/
void prepare_force(struct part *parts, size_t count) {
-#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) && \
- !defined(MINIMAL_SPH) && !defined(MINIMAL_MULTI_MAT_SPH)
+#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) && \
+ !defined(MINIMAL_SPH) && !defined(MINIMAL_MULTI_MAT_SPH) && \
+ !defined(HOPKINS_PU_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
@@ -145,7 +149,7 @@ 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)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH)
p->force.v_sig, 0.f, p->u_dt
#else
0.f, 0.f, 0.f
@@ -549,7 +553,9 @@ 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)
pOrhoi2q[i] = pi_vec.force.P_over_rho2;
+#endif
balsaraiq[i] = pi_vec.force.balsara;
ciq[i] = pi_vec.force.soundspeed;
@@ -560,7 +566,9 @@ 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)
pOrhoj2q[i] = pj_vec[i].force.P_over_rho2;
+#endif
balsarajq[i] = pj_vec[i].force.balsara;
cjq[i] = pj_vec[i].force.soundspeed;
}
@@ -640,7 +648,9 @@ 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)
VEC_HADD(entropy_dtSum, piq[0]->entropy_dt);
+#endif
vec_time += getticks() - vec_tic;
}
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index 917652f195d5838dd2f6659cc748adb3211543de..1ab493a7c149070dc667a2377ab205df7f873856 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -165,14 +165,14 @@ void test() {
printParticle_single(&pi2, &xpi);
print_bytes(&pj, sizeof(struct part));
print_bytes(&pj2, sizeof(struct part));
- error("Particles 'pi' do not match after force (byte = %d)", i_not_ok);
+ error("Particles 'pi' do not match after density (byte = %d)", i_not_ok);
}
if (j_not_ok) {
printParticle_single(&pj, &xpj);
printParticle_single(&pj2, &xpj);
print_bytes(&pj, sizeof(struct part));
print_bytes(&pj2, sizeof(struct part));
- error("Particles 'pj' do not match after force (byte = %d)", j_not_ok);
+ error("Particles 'pj' do not match after density (byte = %d)", j_not_ok);
}
/* --- Test the force loop --- */
diff --git a/theory/SPH/Flavours/plotSoundspeed.py b/theory/SPH/Flavours/plotSoundspeed.py
new file mode 100644
index 0000000000000000000000000000000000000000..5a08ea08f4ad4cbecc1b6b22608de32691eac67b
--- /dev/null
+++ b/theory/SPH/Flavours/plotSoundspeed.py
@@ -0,0 +1,91 @@
+"""
+Makes a movie of the Sedov 2D data. Adapted from
+
+ KelvinHelmholtz_2D/makeMovie.py
+
+You will need to run your movie with far higher time-resolution than usual to
+get a nice movie; around 450 snapshots over 6s is required.
+
+Edit this file near the bottom with the number of snaps you have.
+
+Written by Josh Borrow (joshua.borrow@durham.ac.uk)
+"""
+
+import numpy as np
+
+if __name__ == "__main__":
+ import matplotlib
+ matplotlib.use("Agg")
+ params = {'axes.labelsize': 9,
+ 'axes.titlesize': 10,
+ 'font.size': 12,
+ 'legend.fontsize': 12,
+ 'xtick.labelsize': 9,
+ 'ytick.labelsize': 9,
+ 'text.usetex': True,
+ 'figure.figsize' : (3.15,2.60),
+ 'figure.subplot.left' : 0.17,
+ 'figure.subplot.right' : 0.99 ,
+ 'figure.subplot.bottom' : 0.08 ,
+ 'figure.subplot.top' : 0.99 ,
+ 'figure.subplot.wspace' : 0. ,
+ 'figure.subplot.hspace' : 0. ,
+ 'lines.markersize' : 6,
+ 'lines.linewidth' : 3.,
+ 'text.latex.unicode': True}
+ matplotlib.rcParams.update(params)
+ matplotlib.rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+ from matplotlib.colors import LogNorm
+
+ import matplotlib.pyplot as plt
+
+ filename = "sedov"
+ dpi = 1024
+
+
+ # Creation of first frame
+ fig, ax = plt.subplots(1, 1, frameon=False)
+
+ with np.load("sedov_soundspeed_ratio_data.npz") as file:
+ mesh = file.items()[0][1]
+
+ # Global variable for set_array
+ img = ax.imshow(mesh,
+ extent=[0,1,0,1],
+ animated=True,
+ interpolation="none",
+ norm=LogNorm())
+
+ circle = matplotlib.patches.Circle([0.5, 0.5],
+ radius=0.18242863869665918,
+ animated=True,
+ lw=1,
+ fill=False,
+ ec="red")
+
+ ax.add_artist(circle)
+
+ fig.colorbar(img,
+ label=r"$c_{s, {\rm smoothed}}$ / $c_{s, {\rm gas}}$",
+ pad=0.0)
+
+
+ plt.tick_params(axis='x',
+ which='both',
+ bottom=False,
+ top=False,
+ labelbottom=False)
+ plt.tick_params(axis='y',
+ which='both',
+ left=False,
+ right=False,
+ labelleft=False)
+
+ plt.xlim(0.2, 0.8)
+ plt.ylim(0.2, 0.8)
+
+ # Actually make the movie
+ plt.tight_layout()
+
+ plt.savefig("sedov_blast_soundspeed.pdf", dpi=300)
diff --git a/theory/SPH/Flavours/run.sh b/theory/SPH/Flavours/run.sh
index 6d0791823d93f7feb8f469747f81b24032612523..39677c961cfa6557ae0db52d04ff15e39949c80c 100755
--- a/theory/SPH/Flavours/run.sh
+++ b/theory/SPH/Flavours/run.sh
@@ -1,4 +1,6 @@
#!/bin/bash
+python plotSoundspeed.py
+
pdflatex -jobname=sph_flavours sph_flavours_standalone.tex
bibtex sph_flavours.aux
pdflatex -jobname=sph_flavours sph_flavours_standalone.tex
diff --git a/theory/SPH/Flavours/sedov_soundspeed_ratio_data.npz b/theory/SPH/Flavours/sedov_soundspeed_ratio_data.npz
new file mode 100644
index 0000000000000000000000000000000000000000..1111a2696b76d6ea5a733b29eec47648d9c15a55
Binary files /dev/null and b/theory/SPH/Flavours/sedov_soundspeed_ratio_data.npz differ
diff --git a/theory/SPH/Flavours/sph_flavours.tex b/theory/SPH/Flavours/sph_flavours.tex
index 3c80fefb4989505b76cfaf8b38676ac4276b8da8..5d62af3aab777e66f0b33b89e861d2b21e10b38c 100644
--- a/theory/SPH/Flavours/sph_flavours.tex
+++ b/theory/SPH/Flavours/sph_flavours.tex
@@ -6,7 +6,8 @@ vector quantity $\Wij \equiv \frac{1}{2}\left(W(\vec{x}_{ij}, h_i) +
\nabla_x W(\vec{x}_{ij},h_j)\right)$.
-%#######################################################################################################
+%##############################################################################
+
\subsection{\MinimalSPH}
\label{sec:sph:minimal}
@@ -33,7 +34,8 @@ and the derivative of its density with respect to $h$:
\begin{equation}
\label{eq:sph:minimal:rho_dh}
- \rho_{\partial h_i} \equiv \dd{\rho}{h}(\vec{x}_i) = \sum_j m_j \dd{W}{h}(\vec{x}_{ij}
+ \rho_{\partial h_i} \equiv \dd{\rho}{h}(\vec{x}_i) = \sum_j m_j
+\dd{W}{h}(\vec{x}_{ij}
, h_i).
\end{equation}
This corresponds to $x_i = \tilde{x}_i = m_i$, and $y_i =\tilde{y}_i = \rho_i$
@@ -56,8 +58,8 @@ $\rho_i$ and $u_i$:
\subsubsection{Hydrodynamical accelerations (\nth{2} neighbour loop)}
-We can then proceed with the second loop over
-neighbours. The signal velocity $v_{{\rm sig},ij}$ between two particles is given by
+We can then proceed with the second loop over neighbours. The signal velocity
+$v_{{\rm sig},ij}$ between two particles is given by
\begin{align}
\mu_{ij} &=
@@ -91,7 +93,7 @@ and the change in internal energy,
&+\left. \frac{1}{2}\nu_{ij}\vec{v}_{ij}\cdot\Big. \Wij\right], \nonumber
\end{align}
where in both cases the first line corresponds to the standard SPH
-term and the second line to the viscuous accelerations.
+term and the second line to the viscous accelerations.
We also compute an estimator of the change in smoothing length to be
used in the prediction step. This is an estimate of the local
@@ -119,13 +121,17 @@ For each particle, we compute a time-step given by the CFL condition:
\Delta t = 2 C_{\rm CFL} \frac{H_i}{v_{{\rm sig},i}},
\label{eq:sph:minimal:dt}
\end{equation}
-where $C_{\rm CFL}$ is a free dimensionless parameter and $H_i = \gamma h_i$ is the
-kernel support size. Particles can then be ``kicked'' forward in time:
+where $C_{\rm CFL}$ is a free dimensionless parameter and $H_i = \gamma h_i$ is
+the kernel support size. Particles can then be ``kicked'' forward in time:
\begin{align}
- \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t \label{eq:sph:minimal:kick_v}\\
- u_i &\rightarrow u_i + \frac{du_i}{dt} \Delta t \label{eq:sph:minimal:kick_u}\\
- P_i &\rightarrow P_{\rm eos}\left(\rho_i, u_i\right) \label{eq:sph:minimal:kick_P}, \\
- c_i &\rightarrow c_{\rm eos}\left(\rho_i, u_i\right) \label{eq:sph:minimal:kick_c},
+ \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t
+\label{eq:sph:minimal:kick_v}\\
+ u_i &\rightarrow u_i + \frac{du_i}{dt} \Delta t
+\label{eq:sph:minimal:kick_u}\\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, u_i\right)
+\label{eq:sph:minimal:kick_P}, \\
+ c_i &\rightarrow c_{\rm eos}\left(\rho_i, u_i\right)
+\label{eq:sph:minimal:kick_c},
\end{align}
where we used the pre-defined equation of state to compute the new
value of the pressure and sound-speed.
@@ -136,12 +142,14 @@ Inactive particles need to have their quantities predicted forward in
time in the ``drift'' operation. We update them as follows:
\begin{align}
- \vec{x}_i &\rightarrow \vec{x}_i + \vec{v}_i \Delta t \label{eq:sph:minimal:drift_x} \\
+ \vec{x}_i &\rightarrow \vec{x}_i + \vec{v}_i \Delta t
+\label{eq:sph:minimal:drift_x} \\
h_i &\rightarrow h_i \exp\left(\frac{1}{h_i} \frac{dh_i}{dt}
\Delta t\right), \label{eq:sph:minimal:drift_h}\\
\rho_i &\rightarrow \rho_i \exp\left(-\frac{3}{h_i} \frac{dh_i}{dt}
\Delta t\right), \label{eq:sph:minimal:drift_rho} \\
- P_i &\rightarrow P_{\rm eos}\left(\rho_i, u_i + \frac{du_i}{dt} \Delta t\right), \label{eq:sph:minimal:drift_P}\\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, u_i + \frac{du_i}{dt} \Delta
+t\right), \label{eq:sph:minimal:drift_P}\\
c_i &\rightarrow c_{\rm eos}\left(\rho_i, u_i + \frac{du_i}{dt}
\Delta t\right) \label{eq:sph:minimal:drift_c},
\end{align}
@@ -149,7 +157,7 @@ where, as above, the last two updated quantities are obtained using
the pre-defined equation of state. Note that the thermal energy $u_i$
itself is \emph{not} updated.
-%#######################################################################################################
+%##############################################################################
\subsection{Gadget-2 SPH}
\label{sec:sph:gadget2}
@@ -158,15 +166,16 @@ This flavour of SPH is the one implemented in the \gadget-2 code
\citep{Springel2005}. The basic equations were derived by
\cite{Springel2002} and also includes a \cite{Balsara1995} switch for
the suppression of viscosity. The implementation here follows closely the
-presentation of \cite{Springel2005}. Specifically, we use their equations (5), (7),
-(8), (9), (10), (13), (14) and (17). We summarize them here for completeness.
+presentation of \cite{Springel2005}. Specifically, we use their equations (5),
+(7), (8), (9), (10), (13), (14) and (17). We summarise them here for
+completeness.
\subsubsection{Density and other fluid properties (\nth{1} neighbour loop)}
For a set of particles $i$ with positions $\vec{x}_i$ with velocities
$\vec{v}_i$, masses $m_i$, entropic function per unit mass $A_i$ and
-smoothing length $h_i$, we compute the density, derivative of the density with respect
-to $h$ and the ``h-terms'' in a similar way to the minimal-SPH case
+smoothing length $h_i$, we compute the density, derivative of the density with
+respect to $h$ and the ``h-terms'' in a similar way to the minimal-SPH case
(Eq. \ref{eq:sph:minimal:rho}, \ref{eq:sph:minimal:rho_dh} and
\ref{eq:sph:minimal:f_i}). From these the pressure and sound-speed can
be computed using the pre-defined equation of state:
@@ -175,13 +184,16 @@ be computed using the pre-defined equation of state:
P_i &= P_{\rm eos}(\rho_i, A_i), \label{eq:sph:gadget2:P}\\
c_i &= c_{\rm eos}(\rho_i, A_i). \label{eq:sph:gadget2:c}
\end{align}
-We additionally compute the divergence and
-curl of the velocity field using standard SPH expressions:
+We additionally compute the divergence and curl of the velocity field using
+standard SPH expressions:
\begin{align}
- \nabla\cdot\vec{v}_i &\equiv\nabla\cdot \vec{v}(\vec{x}_i) = \frac{1}{\rho_i} \sum_j m_j
- \vec{v}_{ij}\cdot\nabla_x W(\vec{x}_{ij}, h_i) \label{eq:sph:gadget2:div_v},\\
- \nabla\times\vec{v}_i &\equiv \nabla\times \vec{v}(\vec{x}_i) = \frac{1}{\rho_i} \sum_j m_j
+ \nabla\cdot\vec{v}_i &\equiv\nabla\cdot \vec{v}(\vec{x}_i) = \frac{1}{\rho_i}
+\sum_j m_j
+ \vec{v}_{ij}\cdot\nabla_x W(\vec{x}_{ij}, h_i)
+\label{eq:sph:gadget2:div_v},\\
+ \nabla\times\vec{v}_i &\equiv \nabla\times \vec{v}(\vec{x}_i) =
+\frac{1}{\rho_i} \sum_j m_j
\vec{v}_{ij}\times\nabla_x W(\vec{x}_{ij}, h_i) \label{eq:sph:gadget2:rot_v}.
\end{align}
These are used to construct the \cite{Balsara1995} switch $B_i$:
@@ -190,7 +202,8 @@ These are used to construct the \cite{Balsara1995} switch $B_i$:
B_i = \frac{|\nabla\cdot\vec{v}_i|}{|\nabla\cdot\vec{v}_i| +
|\nabla\times\vec{v}_i| + 10^{-4}c_i / h_i}, \label{eq:sph:gadget2:balsara}
\end{equation}
-where the last term in the denominator is added to prevent numerical instabilities.
+where the last term in the denominator is added to prevent numerical
+instabilities.
\subsubsection{Hydrodynamical accelerations (\nth{2} neighbour loop)}
@@ -199,7 +212,8 @@ case with the exception of the viscosity term which is now modified to
include the switch. Instead of Eq. \ref{eq:sph:minimal:nu_ij}, we get:
\begin{equation}
-\nu_{ij} = -\frac{1}{2}\frac{\alpha \bar B_{ij} \mu_{ij} v_{{\rm sig},ij}}{\bar\rho_{ij}},
+\nu_{ij} = -\frac{1}{2}\frac{\alpha \bar B_{ij} \mu_{ij} v_{{\rm
+sig},ij}}{\bar\rho_{ij}},
\label{eq:sph:gadget2:nu_ij}
\end{equation}
whilst equations \ref{eq:sph:minimal:v_sig},
@@ -207,8 +221,8 @@ whilst equations \ref{eq:sph:minimal:v_sig},
\ref{eq:sph:minimal:v_sig_update} remain unchanged. The only other
change is the equation of motion for the thermodynamic variable which
now has to be describing the evolution of the entropic function and
-not the evolution of the thermal energy. Instead of
-eq. \ref{eq:sph:minimal:du_dt}, we have
+not the evolution of the thermal energy. Instead of eq.
+\ref{eq:sph:minimal:du_dt}, we have
\begin{equation}
\frac{dA_i}{dt} = \frac{1}{2} A_{\rm eos}\left(\rho_i, \sum_j
@@ -228,10 +242,14 @@ by an integration for the the entropy:
\begin{align}
- \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t \label{eq:sph:gadget2:kick_v}\\
- A_i &\rightarrow A_i + \frac{dA_i}{dt} \Delta t \label{eq:sph:gadget2:kick_A}\\
- P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i\right) \label{eq:sph:gadget2:kick_P}, \\
- c_i &\rightarrow c_{\rm eos}\left(\rho_i, A_i\right) \label{eq:sph:gadget2:kick_c},
+ \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t
+\label{eq:sph:gadget2:kick_v}\\
+ A_i &\rightarrow A_i + \frac{dA_i}{dt} \Delta t
+\label{eq:sph:gadget2:kick_A}\\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i\right)
+\label{eq:sph:gadget2:kick_P}, \\
+ c_i &\rightarrow c_{\rm eos}\left(\rho_i, A_i\right)
+\label{eq:sph:gadget2:kick_c},
\end{align}
where, once again, we made use of the equation of state relating
thermodynamical quantities.
@@ -242,12 +260,14 @@ The prediction step is also identical to the \MinimalSPH case with the
entropic function replacing the thermal energy.
\begin{align}
- \vec{x}_i &\rightarrow \vec{x}_i + \vec{v}_i \Delta t \label{eq:sph:gadget2:drift_x} \\
+ \vec{x}_i &\rightarrow \vec{x}_i + \vec{v}_i \Delta t
+\label{eq:sph:gadget2:drift_x} \\
h_i &\rightarrow h_i \exp\left(\frac{1}{h_i} \frac{dh_i}{dt}
\Delta t\right), \label{eq:sph:gadget2:drift_h}\\
\rho_i &\rightarrow \rho_i \exp\left(-\frac{3}{h_i} \frac{dh_i}{dt}
\Delta t\right), \label{eq:sph:gadget2:drift_rho} \\
- P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt} \Delta t\right), \label{eq:sph:gadget2:drift_P}\\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt} \Delta
+t\right), \label{eq:sph:gadget2:drift_P}\\
c_i &\rightarrow c_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt}
\Delta t\right) \label{eq:sph:gadget2:drift_c},
\end{align}
@@ -255,8 +275,39 @@ where, as above, the last two updated quantities are obtained using
the pre-defined equation of state. Note that the entropic function $A_i$
itself is \emph{not} updated.
+\subsection{Weighted-Pressure SPH Validity}
+
+A new class of Lagrangian SPH methods were introduced to the astrophysical
+community by \citet{Hopkins2013} and \citet{Saitoh2013}. Two of these methods,
+Pressure-Entropy (used in the original ANARCHY implementation in EAGLE) and
+Pressure-Energy, are implemented for use in \swift{}. Before considering the
+use of these methods, though, it is important to pause for a moment and
+consider where it is valid to use them in a cosmological context. These methods
+(as implemented in \swift{}) are only valid for cases that use an \emph{ideal
+gas equation of state}, i.e. one in which
+\begin{equation}
+ P = (\gamma - 1) u \rho \propto \rho.
+ \nonumber
+\end{equation}
+Implementations that differ from this, such as the original ANARCHY scheme in
+EAGLE, may have some problems with energy conservation \cite[see][]{Hosono2013}
+and other properties as at their core they assume that pressure is proportional
+to the local energy density, i.e. $P \propto \rho u$. This is most easily shown
+in Pressure-Energy SPH where the weighted pressure $\bar{P}$ is written as
+\begin{equation}
+ \bar{P} = \sum_j \frac{P_j}{\rho_j} W_{ij} = \sum_j m_j (\gamma - 1) u_j
+W_{ij},
+\end{equation}
+and the right-hand side is what is actually calculated using the scheme. It is
+clear that this does not give a valid weighted pressure for any scheme using a
+non-ideal equation of state. Fortunately, there is a general prescription for
+including non-ideal equations of state in the P-X formalisms, but this is yet
+to be implemented in \swift{} and requires an extra density loop. Attempting to
+use these weighted schemes with a non-ideal equation of state will lead to an
+incorrect calculation of both the pressure and the equation fo motion. How
+incorrect this estimate is, however, remains to be seen.
-%#######################################################################################################
+%##############################################################################
\subsection{Pressure-Entropy SPH}
\label{sec:sph:pe}
@@ -284,7 +335,8 @@ which enters many equations. This allows us to construct the
entropy-weighted density $\bar\rho_i$:
\begin{equation}
- \bar\rho_i = \frac{1}{\tilde{A_i}} \sum_j m_j \tilde{A_j} W(\vec{x}_{ij}, h_i),
+ \bar\rho_i = \frac{1}{\tilde{A_i}} \sum_j m_j \tilde{A_j} W(\vec{x}_{ij},
+h_i),
\label{eq:sph:pe:rho_bar}
\end{equation}
which can then be used to construct an entropy-weighted sound-speed
@@ -318,8 +370,11 @@ P_{\partial h_i}$ and $\rho_{\partial h_i}$ (eq. \ref{eq:sph:minimal:rho_dh}):
The accelerations are given by the following term:
\begin{align}
- \frac{d\vec{v}_i}{dt} = -\sum_j m_j &\left[\frac{\bar P_i}{\bar\rho_i^2} \left(\frac{\tilde A_j}{\tilde A_i} - \frac{f_i}{\tilde A_i}\right)\nabla_x W(\vec{x}_{ij}, h_i) \right. \nonumber \\
- &+\frac{P_j}{\rho_j^2} \left(\frac{\tilde A_i}{\tilde A_j} - \frac{f_j}{\tilde A_j}\right)\nabla_x W(\vec{x}_{ij},h_j) \\
+ \frac{d\vec{v}_i}{dt} = -\sum_j m_j &\left[\frac{\bar P_i}{\bar\rho_i^2}
+\left(\frac{\tilde A_j}{\tilde A_i} - \frac{f_i}{\tilde A_i}\right)\nabla_x
+W(\vec{x}_{ij}, h_i) \right. \nonumber \\
+ &+\frac{P_j}{\rho_j^2} \left(\frac{\tilde A_i}{\tilde A_j} -
+\frac{f_j}{\tilde A_j}\right)\nabla_x W(\vec{x}_{ij},h_j) \\
&+ \left. \bigg.\nu_{ij} \Wij \right], \label{eq:sph:pe:dv_dt}
\end{align}
where the viscosity term $\nu_{ij}$ has been computed as in
@@ -342,9 +397,11 @@ internal energy (Eq. \ref{eq:sph:minimal:kick_u}) which gets replaced
by an integration for the the entropy:
\begin{align}
- \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t \label{eq:sph:pe:kick_v}\\
+ \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t
+\label{eq:sph:pe:kick_v}\\
A_i &\rightarrow A_i + \frac{dA_i}{dt} \Delta t \label{eq:sph:pe:kick_A}\\
- P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i\right) \label{eq:sph:pe:kick_P}, \\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i\right)
+\label{eq:sph:pe:kick_P}, \\
c_i &\rightarrow c_{\rm eos}\left(\rho_i,
A_i\right) \label{eq:sph:pe:kick_c}, \\
\tilde A_i &= A_i^{1/\gamma}
@@ -359,14 +416,16 @@ The prediction step is also identical to the \MinimalSPH case with the
entropic function replacing the thermal energy.
\begin{align}
- \vec{x}_i &\rightarrow \vec{x}_i + \vec{v}_i \Delta t \label{eq:sph:pe:drift_x} \\
+ \vec{x}_i &\rightarrow \vec{x}_i + \vec{v}_i \Delta t
+\label{eq:sph:pe:drift_x} \\
h_i &\rightarrow h_i \exp\left(\frac{1}{h_i} \frac{dh_i}{dt}
\Delta t\right), \label{eq:sph:pe:drift_h}\\
\rho_i &\rightarrow \rho_i \exp\left(-\frac{3}{h_i} \frac{dh_i}{dt}
\Delta t\right), \label{eq:sph:pe:drift_rho} \\
\tilde A_i &\rightarrow \left(A_i + \frac{dA_i}{dt}
\Delta t\right)^{1/\gamma} \label{eq:sph:pe:drift_A_tilde}, \\
- P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt} \Delta t\right), \label{eq:sph:pe:drift_P}\\
+ P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt} \Delta
+t\right), \label{eq:sph:pe:drift_P}\\
c_i &\rightarrow c_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt}
\Delta t\right) \label{eq:sph:pe:drift_c},
\end{align}
@@ -377,7 +436,160 @@ itself is \emph{not} updated.
\subsection{Pressure-Energy SPH}
\label{sec:sph:pu}
-Section 2.2.2 of \cite{Hopkins2013}.\\ \tbd
+Section 2.2.2 of \cite{Hopkins2013} describes the equations for Pressure-Energy
+(P-U) SPH; they are reproduced here with some more details.
+
+P-U SPH depends on the calculation of a smoothed pressure, and follows the
+evolution of the internal energy, as opposed to the entropy.
+
+For P-U, the following choice of parameters in the formalism of \S
+\ref{sec:derivation} provides convenient properties:
+\begin{align}
+ x_i =&~ (\gamma - 1) m_i u_i, \\
+ \tilde{x}_i =&~ 1,
+ \label{eq:sph:pu:xichoice}
+\end{align}
+leading to the following requirements to ensure correct volume elements:
+\begin{align}
+ y_i =& \sum_{j} (\gamma - 1) m_j u_j W_{ij} = \bar{P}_i,\\
+ \tilde{y}_i =& \sum_{j} W_{ij} = \bar{n}_i,
+ \label{eq:sph:pu:yichoice}
+\end{align}
+with the resulting variables representing a smoothed pressure and particle
+number density. This choice of variables leads to the following equation of
+motion:
+\begin{align}
+ \frac{\mathrm{d} \mathbf{v}_i}{\mathrm{d} t} = -\sum_j (\gamma - 1)^2 m_j u_j
+u_i
+ &\left[\frac{f_{ij}}{\bar{P}_i} \nabla_i W_{ij}(h_i) ~+ \right. \nonumber \\
+ &\frac{f_{ji}}{\bar{P}_j} \nabla_i W_{ji}(h_j) ~+ \nonumber \\
+ & \left.\nu_{ij}\bar{\nabla_i W_{ij}}\right]~.
+ \label{eq:sph:pu:eom}
+\end{align}
+which includes the Monaghan artificial viscosity term and Balsara switch in
+the final term.
+
+The $h$-terms are given as
+\begin{align}
+ f_{ij} = 1 - \left[\frac{h_i}{n_d (\gamma - 1) \bar{n}_i \left\{m_j
+u_j\right\}}
+ \frac{\partial \bar{P}_i}{\partial h_i} \right]
+ \left( 1 + \frac{h_i}{n_d \bar{n}_i}
+ \frac{\partial \bar{n}_i}{\partial h_i} \right)^{-1}
+ \label{eq:sph:pu:fij}
+\end{align}
+with $n_d$ the number of dimensions. In practice, the majority of $f_{ij}$ is
+precomputed in {\tt hydro\_prepare\_force} as only the curly-bracketed term
+depends on the $j$ particle. This cuts out on the majority of operations,
+including expensive divisions.
+
+In a similar fashion to \MinimalSPH, the internal energy must also be
+evolved. Following \cite{Hopkins2013}, this is calculated as
+\begin{align}
+ \frac{\mathrm{d}u_i}{\mathrm{d}t} = \sum_j (\gamma - 1)^2 m_j u_j u_i
+ \frac{f_{ij}}{\bar{P}_i}(\mathbf{v}_i - \mathbf{v}_j) \cdot
+ \nabla_i W_{ij}(h_i)~.
+ \label{eq:sph:pu:dudt}
+\end{align}
+The sound-speed in P-U requires some consideration. To see what the `correct'
+sound-speed
+is, it is worth looking at the equation of motion (Equation
+\ref{eq:sph:pu:eom}) in
+contrast with that of the EoM for Density-Energy SPH (Equation
+\ref{eq:sph:minimal:dv_dt})
+to see what terms are applicable. For Density-Energy SPH, we see that
+\begin{align}
+ \frac{\mathrm{d}\mathbf{v}_i}{\mathrm{d} t} \sim \frac{c_{s, i}}{\rho_i}
+\nabla_i W_{ij},
+ \nonumber
+\end{align}
+and for Pressure-Energy SPH
+\begin{align}
+ \frac{\mathrm{d}\mathbf{v}_i}{\mathrm{d} t} \sim (\gamma - 1)^2
+ \frac{u_i u_j}{\bar{P}_i} \nabla_i W_{ij}.
+ \nonumber
+\end{align}
+From this it is reasonable to assume that the sound-speed, i.e. the speed at
+which information propagates in the system through pressure waves, is given by
+the expression
+\begin{align}
+ c_{s, i} = (\gamma - 1) u_i \sqrt{\gamma \frac{\rho_i}{\bar{P_i}}}.
+ \label{eq:sph:pu:soundspeedfromeom}
+\end{align}
+This expression is dimensionally consistent with a sound-speed, and includes
+the gas density information (through $\rho$), traditionally used for
+sound-speeds, as well as including the extra information from the smoothed
+pressure $\bar{P}$. However, this scheme causes some problems, which can be
+illustrated using the Sedov Blast test. Such a sound-speed leads to a
+considerably \emph{higher} time-step in front of the shock wave (where the
+smoothed pressure is higher, but the SPH density is relatively constant),
+leading to integration problems. An alternative to this is to use the smoothed
+pressure in the place of the ``real" pressure. Whilst it is well understood
+that $\bar{P}$ should not be used to replace the real pressure in general, here
+(in the sound-speed) it is only used as part of the time-stepping condition.
+Using
+\begin{align}
+ c_{s, i} = \sqrt{\gamma \frac{\bar{P}_i}{\rho_i}}
+ \label{eq:sph:pu:soundspeed}
+\end{align}
+instead of Equation \ref{eq:sph:pu:soundspeedfromeom} leads to a much improved
+time-stepping condition that actually allows particles to be woken up before
+being hit by a shock (see Figure \ref{fig:sph:soundspeed}).
+
+\begin{figure}
+ \centering
+ \includegraphics[width=\columnwidth]{sedov_blast_soundspeed.pdf}
+ \caption{The ratio of the sound-speed calculated in Equation
+ \ref{eq:sph:pu:soundspeed} to the gas sound-speed,
+ $c_s = \sqrt{\gamma(\gamma - 1) u}$ with $u$ the internal energy. Note how
+ the sound-speed increases ahead of the shock, leading to a much smaller
+ time-step for these particles ($\Delta t \propto c_s^{-1}$), waking them up
+ just before they are hit by a shock. The physical reasoning behind using
+ this particular metric for the sound-speed is weak, but as a time-step
+ criterion it appears to be useful. The smoothed pressure calculation
+ ``catches" the hot particles from the shock that is incoming and is
+ increased for those in front of the shock. Thankfully, particles that are
+ behind the shock seem to be relatively unaffected. The red line shows the
+ shock front.}
+ \label{fig:sph:soundspeed}
+\end{figure}
+
+
+\subsubsection{Time integration}
+
+Time integration follows exactly the same scheme as \MinimalSPH, as the
+fundamental equations are exactly the same (just slightly different quantities
+enter the equations of motion). The CFL condition is used, along with the
+sound-speed that is discussed above, such that
+\begin{align}
+ \Delta t = 2 C_{\rm CFL} \frac{H_i}{v_{{\rm sig},i}},
+ \label{eq:sph:pu:dt}
+\end{align}
+where $C_{\rm CFL}$ is a free dimensionless parameter and $H_i = \gamma h_i$ is
+the kernel support size. There is an additional requirement placed on the
+maximal change in $u$ such that
+\begin{align}
+ \Delta t = C_{u} \frac{u}{du/dt},
+ \label{eq:sph:pu:dt_du}
+\end{align}
+with $C_{u}$ a free dimensionless parameter that describes the maximal change
+in $u$ that is allowed.
+
+\subsubsection{Particle properties prediction}
+
+The prediction of particle properties follows exactly the same scheme as
+\MinimalSPH, with the exception of course of the additional smoothed quantity
+$\bar{P}$. This is drifted in the same way as the density; they should both
+follow from the same differential equation with an additional $u$ factor on
+both sides, such that
+\begin{align}
+ \bar{P}_i \rightarrow \bar{P}_i
+ \exp\left(-\frac{3}{h_i}\frac{dh_i}{dt} \Delta t\right).
+ \label{eq:sph:pu:drift}
+\end{align}
+
+%##############################################################################
+
\subsection{Anarchy SPH}
Dalla Vecchia (\textit{in prep.}), also described in section 2.2.2 of
\cite{Schaller2015}.\\
diff --git a/theory/SPH/bibliography.bib b/theory/SPH/bibliography.bib
index 4bcb0e1939a257d54c4c0aa99495d7568838b4f8..2c34da68ec1ed8e4b0321d9312773a308f629b7f 100644
--- a/theory/SPH/bibliography.bib
+++ b/theory/SPH/bibliography.bib
@@ -34,12 +34,14 @@ archivePrefix = "arXiv",
@ARTICLE{Hopkins2013,
author = {{Hopkins}, P.~F.},
- title = "{A general class of Lagrangian smoothed particle hydrodynamics methods and implications for fluid mixing problems}",
+title = "{A general class of Lagrangian smoothed particle hydrodynamics methods
+and implications for fluid mixing problems}",
journal = {\mnras},
archivePrefix = "arXiv",
eprint = {1206.5006},
primaryClass = "astro-ph.IM",
- keywords = {hydrodynamics, instabilities, turbulence, methods: numerical, cosmology: theory},
+keywords = {hydrodynamics, instabilities, turbulence, methods: numerical,
+cosmology: theory},
year = 2013,
month = feb,
volume = 428,
@@ -51,10 +53,12 @@ archivePrefix = "arXiv",
@ARTICLE{Springel2002,
author = {{Springel}, V. and {Hernquist}, L.},
- title = "{Cosmological smoothed particle hydrodynamics simulations: the entropy equation}",
+title = "{Cosmological smoothed particle hydrodynamics simulations: the entropy
+equation}",
journal = {\mnras},
eprint = {astro-ph/0111016},
- keywords = {methods: numerical, galaxies: evolution, galaxies: starburst, methods: numerical, galaxies: evolution, galaxies: starburst},
+keywords = {methods: numerical, galaxies: evolution, galaxies: starburst,
+methods: numerical, galaxies: evolution, galaxies: starburst},
year = 2002,
month = jul,
volume = 333,
@@ -66,7 +70,8 @@ archivePrefix = "arXiv",
@ARTICLE{Balsara1995,
author = {{Balsara}, D.~S.},
- title = "{von Neumann stability analysis of smooth particle hydrodynamics--suggestions for optimal algorithms}",
+title = "{von Neumann stability analysis of smooth particle
+hydrodynamics--suggestions for optimal algorithms}",
journal = {Journal of Computational Physics},
year = 1995,
volume = 121,
@@ -81,11 +86,13 @@ archivePrefix = "arXiv",
author = {{Schaller}, M. and {Dalla Vecchia}, C. and {Schaye}, J. and
{Bower}, R.~G. and {Theuns}, T. and {Crain}, R.~A. and {Furlong}, M. and
{McCarthy}, I.~G.},
- title = "{The EAGLE simulations of galaxy formation: the importance of the hydrodynamics scheme}",
+title = "{The EAGLE simulations of galaxy formation: the importance of the
+hydrodynamics scheme}",
journal = {\mnras},
archivePrefix = "arXiv",
eprint = {1509.05056},
- keywords = {methods: numerical, galaxies: clusters: intracluster medium, galaxies: formation, cosmology: theory},
+keywords = {methods: numerical, galaxies: clusters: intracluster medium,
+galaxies: formation, cosmology: theory},
year = 2015,
month = dec,
volume = 454,
@@ -100,7 +107,8 @@ archivePrefix = "arXiv",
@ARTICLE{Dehnen2012,
author = {{Dehnen}, W. and {Aly}, H.},
- title = "{Improving convergence in smoothed particle hydrodynamics simulations without pairing instability}",
+title = "{Improving convergence in smoothed particle hydrodynamics simulations
+without pairing instability}",
journal = {\mnras},
archivePrefix = "arXiv",
eprint = {1204.2471},
@@ -114,3 +122,44 @@ archivePrefix = "arXiv",
adsurl = {http://adsabs.harvard.edu/abs/2012MNRAS.425.1068D},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}
+
+
+
+@article{Saitoh2013,
+archivePrefix = {arXiv},
+arxivId = {1202.4277},
+author = {Saitoh, Takayuki R. and Makino, Junichiro},
+doi = {10.1088/0004-637X/768/1/44},
+eprint = {1202.4277},
+file = {:Users/josh/Downloads/Saitoh{\_}2013{\_}ApJ{\_}768{\_}44.pdf:pdf},
+isbn = {0004-637X$\backslash$r1538-4357},
+issn = {15384357},
+journal = {Astrophysical Journal},
+keywords = {galaxies: ISM,galaxies: evolution,methods: numerical},
+number = {1},
+title = {{A density-independent formulation of smoothed particle
+hydrodynamics}},
+volume = {768},
+year = {2013}
+}
+
+
+
+@article{Hosono2013,
+archivePrefix = {arXiv},
+arxivId = {1307.0916},
+author = {Hosono, Natsuki and Saitoh, Takayuki R. and Makino, Junichiro},
+doi = {10.1093/pasj/65.5.108},
+eprint = {1307.0916},
+file = {:Users/josh/Downloads/pasj65-0108.pdf:pdf},
+issn = {0004-6264},
+keywords = {hydrodynamics,methods,numerical},
+number = {May},
+pages = {1--11},
+title = {{Density Independent Smoothed Particle Hydrodynamics for Non-Ideal
+Equation of State}},
+url =
+{http://arxiv.org/abs/1307.0916{\%}0Ahttp://dx.doi.org/10.1093/pasj/65.5.108},
+year = {2013}
+}
+
diff --git a/theory/SPH/run.sh b/theory/SPH/run.sh
index 651aadad79c2471f58221e2b4fcad7a09ab12256..8d33be12825a31b2906e1259e31185fee2cc74bf 100755
--- a/theory/SPH/run.sh
+++ b/theory/SPH/run.sh
@@ -4,6 +4,10 @@ python kernels.py
cp kernels.pdf ..
cp kernel_derivatives.pdf ..
cd ..
+cd Flavours
+python plotSoundspeed.py
+cp sedov_blast_soundspeed.pdf ..
+cd ..
pdflatex swift_sph.tex
bibtex swift_sph.aux
pdflatex swift_sph.tex
diff --git a/theory/SPH/swift_sph.tex b/theory/SPH/swift_sph.tex
index 4be8b965b2888bb04b5c150b41ccc38ef2ec3f95..e9c185c3cd0b845bff75be2092846bffbdcfd1a9 100644
--- a/theory/SPH/swift_sph.tex
+++ b/theory/SPH/swift_sph.tex
@@ -29,7 +29,7 @@
\section{Equation of state}
\input{EoS/eos}
-\section{Derivation of the Equation of Motion}
+\section{Derivation of the Equation of Motion}\label{sec:derivation}
\input{Derivation/sph_derivation.tex}
\section{SPH flavours}