diff --git a/configure.ac b/configure.ac
index 6f47f86d304bf666fc261bd6c941ec081cc17a6f..498820bfd17e9dca0e076c32e9520eb41c24ef85 100644
--- a/configure.ac
+++ b/configure.ac
@@ -362,6 +362,21 @@ elif test "$no_gravity_below_id" != "no"; then
AC_DEFINE_UNQUOTED([SWIFT_NO_GRAVITY_BELOW_ID], [$enableval] ,[Particles with smaller ID than this will have zero gravity forces])
fi
+# Check if we want to use boundary particles.
+AC_ARG_ENABLE([boundary-particles],
+ [AS_HELP_STRING([--enable-boundary-particles],
+ [Set all particles with an ID smaller than @<:@N@:>@ as boundary particles (i.e. receive zero gravity + hydro forces).]
+ )],
+ [boundary_particles="$enableval"],
+ [boundary_particles="no"]
+)
+if test "$boundary_particles" == "yes"; then
+ AC_MSG_ERROR(Need to specify the ID below which particles get zero forces when using --enable-boundary-particles!)
+elif test "$boundary_particles" != "no"; then
+ AC_DEFINE_UNQUOTED([SWIFT_NO_GRAVITY_BELOW_ID], [$enableval] ,[Particles with smaller ID than this will have zero gravity forces])
+ AC_DEFINE_UNQUOTED([SWIFT_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.])
+fi
+
# Check whether we have any of the ARM v8.1 tick timers
AX_ASM_ARM_PMCCNTR
AX_ASM_ARM_CNTVCT
@@ -1832,7 +1847,7 @@ esac
# External potential
AC_ARG_WITH([ext-potential],
[AS_HELP_STRING([--with-ext-potential=<pot>],
- [external potential @<:@none, point-mass, point-mass-ring, point-mass-softened, isothermal, nfw, hernquist, disc-patch, sine-wave, default: none@:>@]
+ [external potential @<:@none, point-mass, point-mass-ring, point-mass-softened, isothermal, nfw, hernquist, disc-patch, sine-wave, constant, default: none@:>@]
)],
[with_potential="$withval"],
[with_potential="none"]
@@ -1865,6 +1880,9 @@ case "$with_potential" in
point-mass-softened)
AC_DEFINE([EXTERNAL_POTENTIAL_POINTMASS_SOFT], [1], [Softened point-mass potential with form 1/(r^2 + softening^2).])
;;
+ constant)
+ AC_DEFINE([EXTERNAL_POTENTIAL_CONSTANT], [1], [Constant gravitational acceleration.])
+ ;;
*)
AC_MSG_ERROR([Unknown external potential: $with_potential])
;;
@@ -2043,5 +2061,6 @@ AC_MSG_RESULT([
Naive stars interactions : $enable_naive_interactions_stars
Gravity checks : $gravity_force_checks
Custom icbrtf : $enable_custom_icbrtf
+ Boundary particles : $boundary_particles
------------------------])
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/README b/examples/HydroTests/Rayleigh-Taylor_2D/README
new file mode 100644
index 0000000000000000000000000000000000000000..b1ac67c0c24d98d6b3e03cec70697114bb2fb434
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/README
@@ -0,0 +1,17 @@
+Rayleigh Taylor
+---------------
+
+This is a common test for hydrodynamics (see Abel 2011, Saitoh and Makino 2013, Hopkins 2013).
+It consists in a low density layer of fluid at the bottom and a high density layer at the top.
+Due to a constant vertical gravitational acceleration, the two fluids mix togerther through the
+Rayleigh Taylor instability (e.g. nuclear mushroom cloud).
+
+In this example, we follow the implementation of Saitoh and Makino 2013, but with a longer box in order
+to avoid an interaction with the wall (see last image in Figure 10).
+
+The code needs to be compiled with the following options: `--with-hydro-dimension=2`,
+`--with-ext-potential=constant`, `--enable-boundary-particles` and `--with-adiabatic-index=7/5`.
+I also recommend to use `--with-kernel=wendland-C2`.
+
+The option `--enable-boundary-particles` requires the ID of the last boundary particle.
+This value is provided by the script makeIC.py and depends on the resolution.
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/generate_movie_from_png.sh b/examples/HydroTests/Rayleigh-Taylor_2D/generate_movie_from_png.sh
new file mode 100644
index 0000000000000000000000000000000000000000..f97e447dc21abeabac6fe98eb23e17efd9051cf4
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/generate_movie_from_png.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+d1=gizmo/png
+d2=anarchy/png
+d3=pu/png
+
+for i in {0..1040}
+do
+ echo $i
+ tmp=`printf "%04i" $i`
+ convert $d1/rayleigh_taylor_$tmp.png $d2/rayleigh_taylor_$tmp.png $d3/rayleigh_taylor_$tmp.png +append movie/rayleigh_taylor_$tmp.png
+done
+
+ffmpeg -framerate 10 -i movie/rayleigh_taylor_%04d.png -c:v libx264 -r 30 -pix_fmt yuv420p rayleigh_taylor.mp4
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/makeIC.py b/examples/HydroTests/Rayleigh-Taylor_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..d75d5f215f8b71d859693b842a1e637b52095f45
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/makeIC.py
@@ -0,0 +1,250 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# 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/>.
+#
+##############################################################################
+
+import h5py
+import numpy as np
+from scipy.optimize import bisect
+
+# Generates a swift IC file for the Rayleigh-Taylor instability in a periodic
+# box following the conditions given in Saitoh and Makino 2013: 1202.4277v3
+
+# Parameters
+N = [128, 192] # Particles along one edge
+gamma = 7./5. # Gas adiabatic index
+dv = 0.025 # velocity perturbation
+rho_h = 2 # high density region
+rho_l = 1 # Low density region
+g = -0.5 # gravitational acceleration
+box_size = [1., 1.5] # size of the box
+
+fixed = [0.1, 1.4] # y-range of non fixed particles
+perturbation_box = [0.3, 1.2] # y-range for the velocity perturbation
+fileOutputName = "rayleigh_taylor.hdf5"
+
+
+# ---------------------------------------------------
+
+if (N[0] / box_size[0] != N[1] / box_size[1]):
+ raise Exception("Suppose the same ratio for each directions")
+
+
+def density(y):
+ """
+ Mass density as function of the position y.
+ """
+ if isinstance(y, float) or isinstance(y, int):
+ y = np.array(y)
+
+ ind = y < 0.5 * box_size[1]
+ rho = np.zeros(y.shape)
+ tmp = (gamma - 1.) * g / (gamma * P0)
+ alpha = 1. / (gamma - 1.)
+
+ rho[ind] = rho_l * (1 + rho_l * tmp * (y[ind] - 0.5 * box_size[1]))**alpha
+ rho[~ind] = rho_h * (1 + rho_h * tmp * (y[~ind] - 0.5 * box_size[1]))**alpha
+
+ return rho
+
+
+def mass(y):
+ """
+ Integral of the density
+ """
+ if isinstance(y, float) or isinstance(y, int):
+ y = np.array(y)
+
+ ind = y < 0.5 * box_size[1]
+ m = np.zeros(y.shape)
+
+ B = (gamma - 1.) * g / (gamma * P0)
+ alpha = 1. / (gamma - 1.)
+
+ m[ind] = (1 + B * rho_l * (y[ind] - 0.5 * box_size[1]))**(alpha + 1)
+
+ m[~ind] = (1 + B * rho_h * (y[~ind] - 0.5 * box_size[1]))**(alpha + 1)
+
+ m -= (1 - 0.5 * B * box_size[1] * rho_l)**(alpha + 1)
+
+ return box_size[0] * m / (B * (alpha + 1))
+
+
+P0 = rho_h / gamma # integration constant of the hydrostatic equation
+numPart = N[0] * N[1]
+
+m_tot = mass(box_size[1])
+m_part = m_tot / numPart
+
+
+def inv_mass(m):
+ """
+ Inverse of the function `mass`.
+ """
+ left = 0
+ right = box_size[1]
+
+ def f(y, x):
+ return x - mass(y)
+
+ return bisect(f, left, right, args=(m))
+
+
+def entropy(y):
+ """
+ Entropy as function of the position y.
+ Here we assume isoentropic medium.
+ """
+ if isinstance(y, float) or isinstance(y, int):
+ y = np.array(y)
+
+ ind = y < 0.5 * box_size[1]
+
+ a = np.zeros(y.shape)
+ a[ind] = P0 * rho_l**(-gamma)
+ a[~ind] = P0 * rho_h**(-gamma)
+ return a
+
+
+def smoothing_length(rho, m):
+ """
+ Compute the smoothing length
+ """
+ return 1.23 * np.sqrt(m / rho)
+
+
+def growth_rate():
+ """
+ Compute the growth rate of the instability.
+ Assumes a wavelength equal to the boxsize.
+ """
+ ymin = 0.
+ ymax = box_size[1]
+ A = density(ymax) - density(ymin)
+ A /= density(ymax) + density(ymin)
+ return np.sqrt(A * np.abs(g) * ymax / (2. * np.pi))
+
+
+def vy(x, y):
+ """
+ Velocity along the direction y
+ """
+ ind = y > perturbation_box[0]
+ ind = np.logical_and(ind, y < perturbation_box[1])
+
+ v = np.zeros(len(x))
+
+ v[ind] = 1 + np.cos(4 * np.pi * x[ind])
+ v[ind] *= 1 + np.cos(5 * np.pi * (y[ind] - 0.5 * box_size[1]))
+ v[ind] *= dv
+ return v
+
+
+if __name__ == "__main__":
+ # Start by generating grids of particles
+
+ coords = np.zeros((numPart, 3))
+ m = np.ones(numPart) * m_part
+ u = np.zeros(numPart)
+ vel = np.zeros((numPart, 3))
+ ids = np.zeros(numPart)
+
+ # generate grid of particles
+ y_prev = 0
+ uni_id = 1
+
+ # loop over y
+ eps = 1e-3 * box_size[1] / N[1]
+ for j in range(N[1]):
+ m_y = m_part * N[0] + mass(y_prev)
+ if m_y > m_tot:
+ y_j = box_size[1] - eps * (box_size[1] - y_prev)
+ else:
+ y_j = inv_mass(m_y)
+ y_prev = y_j
+
+ # loop over x
+ for i in range(N[0]):
+
+ index = j * N[0] + i
+
+ x = i * box_size[0] / float(N[0])
+
+ coords[index, 0] = x
+ coords[index, 1] = y_j
+ if (y_j < fixed[0] or y_j > fixed[1]):
+ ids[index] = uni_id
+ uni_id += 1
+
+ print("You need to compile the code with "
+ "--enable-boundary-particles=%i" % uni_id)
+ ids[ids == 0] = np.linspace(uni_id, numPart, numPart-uni_id+1)
+
+ # density
+ rho = density(coords[:, 1])
+
+ # internal energy
+ a = entropy(coords[:, 1])
+ u = a * rho**(gamma-1) / (gamma - 1.)
+
+ # smoothing length
+ h = smoothing_length(rho, m)
+
+ # Velocity perturbation
+ vel[:, 1] = vy(coords[:, 0], coords[:, 1])
+
+ # File
+ fileOutput = h5py.File(fileOutputName, 'w')
+
+ # Header
+ grp = fileOutput.create_group("/Header")
+ grp.attrs["BoxSize"] = box_size
+ grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+ grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+ grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+ grp.attrs["Time"] = 0.0
+ grp.attrs["NumFileOutputsPerSnapshot"] = 1
+ grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ grp.attrs["Flag_Entropy_ICs"] = 0
+ grp.attrs["Dimension"] = 2
+
+ # Units
+ grp = fileOutput.create_group("/Units")
+ grp.attrs["Unit length in cgs (U_L)"] = 1.
+ grp.attrs["Unit mass in cgs (U_M)"] = 1.
+ grp.attrs["Unit time in cgs (U_t)"] = 1.
+ grp.attrs["Unit current in cgs (U_I)"] = 1.
+ grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+ # Particle group
+ grp = fileOutput.create_group("/PartType0")
+ ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ ds[()] = coords
+ ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ ds[()] = vel
+ ds = grp.create_dataset('Masses', (numPart, 1), 'f')
+ ds[()] = m.reshape((numPart, 1))
+ ds = grp.create_dataset('SmoothingLength', (numPart, 1), 'f')
+ ds[()] = h.reshape((numPart, 1))
+ ds = grp.create_dataset('InternalEnergy', (numPart, 1), 'f')
+ ds[()] = u.reshape((numPart, 1))
+ ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ ds[()] = ids.reshape((numPart, 1))
+ ds = grp.create_dataset('Density', (numPart, 1), 'f')
+ ds[()] = rho.reshape((numPart, 1))
+
+ fileOutput.close()
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
new file mode 100644
index 0000000000000000000000000000000000000000..b00c2946994e0e3d8e93b0824a4bc44b9079967e
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
@@ -0,0 +1,280 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+# 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/>.
+#
+##############################################################################
+
+from swiftsimio import load
+from swiftsimio.visualisation import scatter
+from matplotlib.animation import FuncAnimation
+from matplotlib.colors import LogNorm
+import matplotlib.pyplot as plt
+
+from numpy import max, min
+import numpy as np
+
+try:
+ # Try and load this, otherwise we're stuck with serial
+ from p_tqdm import p_map
+except:
+ print("Try installing the p_tqdm package to make movie frames in parallel")
+ p_map = map
+ pass
+
+
+info_frames = 40
+generate_png = False
+text_args = dict(color="black")
+
+
+class Metadata(object):
+ """
+ Copy the useful data in order to decrease the memory usage
+ """
+ def __init__(self, data):
+ metadata = data.metadata
+ self.t = metadata.t
+ try:
+ self.viscosity_info = metadata.viscosity_info
+ except:
+ self.viscosity_info = "No info"
+ try:
+ self.diffusion_info = metadata.diffusion_info
+ except:
+ self.diffusion_info = "No info"
+
+ self.code_info = metadata.code_info
+ self.compiler_info = metadata.compiler_info
+ self.hydro_info = metadata.hydro_info
+
+
+def project(data, m_res, property, ylim):
+ x, y, _ = data.gas.coordinates.value.T
+
+ mask = np.logical_and(y >= ylim[0], y <= ylim[1])
+
+ x = x[mask]
+ y = y[mask] - np.float64(ylim[0])
+
+ h = data.gas.smoothing_length[mask]
+
+ if property == "density":
+ property = "masses"
+
+ if property is not None:
+ quant = getattr(data.gas, property).value[mask]
+ else:
+ quant = np.ones_like(x)
+
+ image = scatter(x=x, y=y, m=quant, h=h, res=m_res)
+
+ return image.T
+
+
+def load_and_make_image(filename, res, property):
+ image = np.zeros(res, dtype=np.float32)
+ m_res = min(res)
+ border = int(0.2 * m_res)
+
+ # first part of the image
+ ylim = np.array([0., 1.])
+
+ data = load(filename)
+ image[:m_res, :m_res] = project(data, m_res, property, ylim)
+ if property != "density":
+ image[:m_res, :m_res] /= project(data, m_res, None, ylim)
+
+ # second part of the image
+ ylim = np.array([0.5, 1.5])
+
+ left = -m_res + border
+ image[left:, :] = project(data, m_res, property, ylim)[left:, :]
+ if property != "density":
+ image[left:, :] /= project(data, m_res, None, ylim)[left:, :]
+
+ metadata = Metadata(data)
+ return image, metadata
+
+
+def time_formatter(metadata):
+ return f"$t = {metadata.t:2.2f}$"
+
+
+def create_movie(filename, start, stop, resolution, property, output_filename):
+ """
+ Creates the movie with:
+
+ snapshots named {filename}_{start}.hdf5 to {filename}_{stop}.hdf5
+ at {resolution} x {resolution} pixel size and smoothing the given
+ {property} and outputting to {output_filename}.
+ """
+
+ if property != "density":
+ name = property
+ else:
+ name = "Fluid Density $\\rho$"
+
+ def baked_in_load(n):
+ f = filename + "_{:04d}.hdf5".format(n)
+ return load_and_make_image(f, resolution, property)
+
+ # Make frames in parallel (reading also parallel!)
+ frames, metadata = zip(*p_map(baked_in_load, list(range(start, stop))))
+
+ vmax = max(list(p_map(max, frames)))
+ vmin = min(list(p_map(min, frames)))
+
+ fig, ax = plt.subplots(figsize=(8, 1.5 * 8), dpi=resolution[0] // 8)
+ ax.axis("off")
+ fig.subplots_adjust(0, 0, 1, 1)
+
+ norm = LogNorm(vmin=vmin, vmax=vmax, clip="black")
+
+ image = ax.imshow(np.zeros_like(frames[0]), origin="lower",
+ norm=norm)
+
+ description_text = ax.text(
+ 0.5, 0.5,
+ get_simulation_information(metadata[0]),
+ va="center", ha="center",
+ **text_args, transform=ax.transAxes,
+ )
+
+ time_text = ax.text(
+ (1 - 0.025 * 0.25), 0.975,
+ time_formatter(metadata[0]),
+ **text_args,
+ va="top", ha="right",
+ transform=ax.transAxes,
+ )
+
+ ax.text(
+ 0.025 * 0.25, 0.975, name, **text_args, va="top", ha="left",
+ transform=ax.transAxes
+ )
+
+ def frame(n):
+ if n >= 0:
+ image.set_array(frames[n])
+ description_text.set_text("")
+ time_text.set_text(time_formatter(metadata[n]))
+
+ if generate_png:
+ name = filename + "_{:04d}".format(n+info_frames)
+ fig.savefig(name + ".png")
+
+ else:
+ return (image,)
+
+ if generate_png:
+ for i in range(-info_frames, stop-start):
+ frame(i)
+ else:
+ animation = FuncAnimation(fig, frame,
+ range(-info_frames, stop-start),
+ interval=40)
+ animation.save(output_filename)
+
+
+def get_simulation_information(metadata):
+ """
+ Generates a string from the SWIFT metadata
+ """
+ viscosity = metadata.viscosity_info
+ diffusion = metadata.diffusion_info
+
+ output = (
+ "$\\bf{Rayleigh-Taylor}$ $\\bf{Instability}$\n\n"
+ "$\\bf{SWIFT}$\n"
+ + metadata.code_info
+ + "\n\n"
+ + "$\\bf{Compiler}$\n"
+ + metadata.compiler_info
+ + "\n\n"
+ + "$\\bf{Hydrodynamics}$\n"
+ + metadata.hydro_info
+ + "\n\n"
+ + "$\\bf{Viscosity}$\n"
+ + viscosity
+ + "\n\n"
+ + "$\\bf{Diffusion}$\n"
+ + diffusion
+ )
+
+ return output
+
+
+if __name__ == "__main__":
+ from argparse import ArgumentParser
+
+ parser = ArgumentParser(description="Creates a movie of the whole box.")
+
+ parser.add_argument(
+ "-s",
+ "--snapshot",
+ help="Snapshot name. Default: rayleigh_taylor",
+ type=str,
+ default="rayleigh_taylor",
+ )
+
+ parser.add_argument(
+ "-i", "--initial", help="Initial snapshot. Default: 0",
+ type=int, default=0
+ )
+
+ parser.add_argument(
+ "-f", "--final", help="Final snapshot. Default: 40",
+ type=int, default=40
+ )
+
+ parser.add_argument(
+ "-o",
+ "--output",
+ help="Output filename. Default: rayleigh_taylor.mp4",
+ type=str,
+ default="rayleigh_taylor.mp4",
+ )
+
+ parser.add_argument(
+ "-p",
+ "--property",
+ help="(swiftsimio) Property to plot. Default: density",
+ type=str,
+ default="density",
+ )
+
+ parser.add_argument(
+ "-r",
+ "--resolution",
+ help="Resolution to make the movie at. Default: 512",
+ type=int,
+ default=512,
+ )
+
+ vars = parser.parse_args()
+
+ yres = int(1.5 * vars.resolution)
+ vars.resolution = [yres, vars.resolution]
+
+ create_movie(
+ filename=vars.snapshot,
+ start=vars.initial,
+ stop=vars.final,
+ resolution=vars.resolution,
+ property=vars.property,
+ output_filename=vars.output,
+ )
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py b/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py
new file mode 100644
index 0000000000000000000000000000000000000000..e89c4e8525fe4c88e517acbd453b0941f8f573c8
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py
@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+import makeIC
+import matplotlib.pyplot as plt
+from swiftsimio import load
+import numpy as np
+
+
+N = 1000
+filename = "rayleigh_taylor_0000.hdf5"
+
+f = load(filename)
+
+# Get data from snapshot
+x, y, _ = f.gas.coordinates.value.T
+rho = f.gas.density.value
+a = f.gas.entropy.value
+
+# Get analytical solution
+y_an = np.linspace(0, makeIC.box_size[1], N)
+rho_an = makeIC.density(y_an)
+a_an = makeIC.entropy(y_an)
+
+plt.figure()
+plt.plot(a, y, ".r", label="Entropy")
+plt.plot(a_an, y_an, "--r")
+
+plt.plot(rho, y, ".k", label="Density")
+plt.plot(rho_an, y_an, "--k")
+
+plt.ylabel("Position")
+plt.xlabel("Density / Entropy")
+
+plt.xlim(0, 2.5)
+plt.ylim(0, makeIC.box_size[1])
+plt.legend()
+
+plt.show()
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml b/examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml
new file mode 100644
index 0000000000000000000000000000000000000000..9e843c49448c27ff8c1452f7e28adbc6fe9a9cd6
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml
@@ -0,0 +1,37 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 10. # The end time of the simulation (in internal units).
+ dt_min: 1e-8 # 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).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: rayleigh_taylor # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./rayleigh_taylor.hdf5 # The file to read
+ periodic: 1
+
+ConstantPotential:
+ g_cgs: [0, -0.5, 0]
\ No newline at end of file
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/run.sh b/examples/HydroTests/Rayleigh-Taylor_2D/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..0c8d3eaab9a401c4be3a586e3b7b8900406bbdbb
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/run.sh
@@ -0,0 +1,15 @@
+#!/bin/bash
+
+set -e
+
+# Generate the initial conditions if they are not present.
+if [ ! -e rayleigh_taylor.hdf5 ]
+then
+ echo "Generating initial conditions for the Rayleigh Taylor example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../../swift --hydro --external-gravity --threads=8 rayleigh_taylor.yml 2>&1 | tee output.log
+
+python makeMovie.py -i 0 -f 1001
diff --git a/src/potential.h b/src/potential.h
index 59567fe92296068f838c39a3eb5ff55c14005d48..9d1418cadafcabe1d66658c43f2c1ef15778fc2b 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -46,6 +46,8 @@
#include "./potential/point_mass_ring/potential.h"
#elif defined(EXTERNAL_POTENTIAL_POINTMASS_SOFT)
#include "./potential/point_mass_softened/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_CONSTANT)
+#include "./potential/constant/potential.h"
#else
#error "Invalid choice of external potential"
#endif
diff --git a/src/potential/constant/potential.h b/src/potential/constant/potential.h
new file mode 100644
index 0000000000000000000000000000000000000000..f91d2b3fc8a13792901fee680678c9b93af20cbf
--- /dev/null
+++ b/src/potential/constant/potential.h
@@ -0,0 +1,143 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * 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_POTENTIAL_CONSTANT_H
+#define SWIFT_POTENTIAL_CONSTANT_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief External Potential Properties - Constant acceleration.
+ */
+struct external_potential {
+
+ /*! Value of the acceleration */
+ double g[3];
+};
+
+/**
+ * @brief Computes the time-step due to the acceleration from an constant
+ * acceleration.
+ *
+ * @param time The current time.
+ * @param potential The #external_potential used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static float external_gravity_timestep(
+ double time, const struct external_potential* restrict potential,
+ const struct phys_const* restrict phys_const,
+ const struct gpart* restrict g) {
+
+ return FLT_MAX;
+}
+
+/**
+ * @brief Computes the gravitational acceleration from an constant acceleration.
+ *
+ * @param time The current time.
+ * @param potential The #external_potential used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static void external_gravity_acceleration(
+ double time, const struct external_potential* potential,
+ const struct phys_const* const phys_const, struct gpart* g) {
+
+ g->a_grav[0] += potential->g[0];
+ g->a_grav[1] += potential->g[1];
+ g->a_grav[2] += potential->g[2];
+}
+
+/**
+ * @brief Computes the gravitational potential energy of a particle in an
+ * constant acceleration.
+ *
+ * @param time The current time (unused here).
+ * @param potential The #external_potential used in the run.
+ * @param phys_const Physical constants in internal units.
+ * @param g Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+external_gravity_get_potential_energy(
+ double time, const struct external_potential* potential,
+ const struct phys_const* const phys_const, const struct gpart* g) {
+
+ const float gh = g->x[0] * potential->g[0] + g->x[1] * potential->g[1] +
+ g->x[2] * potential->g[2];
+
+ return g->mass * gh;
+}
+
+/**
+ * @brief Initialises the external potential properties in the internal system
+ * of units.
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param potential The external potential properties to initialize
+ */
+static INLINE void potential_init_backend(
+ struct swift_params* parameter_file, const struct phys_const* phys_const,
+ const struct unit_system* us, const struct space* s,
+ struct external_potential* potential) {
+
+ /* Read in the acceleration */
+ parser_get_param_double_array(parameter_file, "ConstantPotential:g_cgs", 3,
+ potential->g);
+
+ /* Change the unit system */
+ const double unit_length = units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
+ const double unit_time = units_cgs_conversion_factor(us, UNIT_CONV_TIME);
+ const double unit_g = unit_length / (unit_time * unit_time);
+
+ for (int i = 0; i < 3; i++) {
+ // Need to divide by G due to gravity_end_force
+ potential->g[i] /= unit_g * phys_const->const_newton_G;
+ }
+}
+
+/**
+ * @brief Prints the properties of the external potential to stdout.
+ *
+ * @param potential The external potential properties.
+ */
+static INLINE void potential_print_backend(
+ const struct external_potential* potential) {
+
+ message(
+ "External potential is 'Constant' with properties are g = (%e, "
+ "%e, %e)",
+ potential->g[0], potential->g[1], potential->g[2]);
+}
+
+#endif /* SWIFT_CONSTANT_H */
diff --git a/src/runner.c b/src/runner.c
index 7f645665eb786ab578fd8e0cfccaa72420c08b1b..fdf371db2fb5b8a9ace982d4c35f7e37c4d2e3a1 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3614,6 +3614,26 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
/* Finish the force loop */
hydro_end_force(p, cosmo);
+
+#ifdef SWIFT_BOUNDARY_PARTICLES
+
+ /* Get the ID of the part */
+ const long long id = p->id;
+
+ /* Cancel hdyro forces of these particles */
+ if (id < SWIFT_BOUNDARY_PARTICLES) {
+
+ /* Don't move ! */
+ hydro_reset_acceleration(p);
+
+#if defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
+
+ /* Some values need to be reset in the Gizmo case. */
+ hydro_prepare_force(p, &c->hydro.xparts[k], cosmo,
+ e->hydro_properties, 0);
+#endif
+ }
+#endif
}
}
}