From ad691f87126b0afa3c88f9ae90fcbc4b86eb05a3 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Fri, 12 Jul 2019 14:14:39 +0200
Subject: [PATCH 01/16] Add rayleigh taylor

---
 examples/HydroTests/Rayleigh-Taylor/makeIC.py | 171 ++++++++++++++++++
 .../Rayleigh-Taylor/rayleigh_taylor.yml       |  37 ++++
 examples/HydroTests/Rayleigh-Taylor/run.sh    |  15 ++
 3 files changed, 223 insertions(+)
 create mode 100644 examples/HydroTests/Rayleigh-Taylor/makeIC.py
 create mode 100644 examples/HydroTests/Rayleigh-Taylor/rayleigh_taylor.yml
 create mode 100644 examples/HydroTests/Rayleigh-Taylor/run.sh

diff --git a/examples/HydroTests/Rayleigh-Taylor/makeIC.py b/examples/HydroTests/Rayleigh-Taylor/makeIC.py
new file mode 100644
index 000000000..53b9eceea
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor/makeIC.py
@@ -0,0 +1,171 @@
+###############################################################################
+# 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 vortex in a periodic box
+# following the conditions given in Hopkins 2013
+
+# Parameters
+N = [50, 100]  # Particles along one edge
+gamma = 5./3.  # Gas adiabatic index
+delta = 0.025  # interface size
+dv = 0.025  # velocity perturbation
+rho1 = 1    # Lower region density
+rho2 = 2    # Upper region density
+P0 = 0      # integration constant of the hydrostatic equation
+g = -0.5    # gravitational acceleration
+box_size = [0.5, 1.]  # size of the box
+fileOutputName = "rayleigh_taylor.hdf5"
+
+
+def density(y):
+    """
+    Mass density as function of the position y.
+    """
+    tmp = 1 + np.exp(-(y - 0.5 * box_size[1]) / delta)
+    return rho1 + (rho2 - rho1) / tmp
+
+
+def mass(y):
+    """
+    Integrated Mass in x and y.
+    """
+    tmp = np.log(np.exp(y / delta) + np.exp(box_size[1] * 0.5 / delta))
+    tmp -= 0.5 * box_size[1] / delta
+    m = box_size[0] * (rho1 * y + (rho2 - rho1) * delta * tmp)
+    return m
+
+
+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 pressure(y):
+    """
+    Pressure as function of the position y.
+    Here we assume hydrostatic equation.
+    """
+    tmp = np.log(np.exp(y / delta) + np.exp(box_size[1] * 0.5 / delta))
+    tmp -= 0.5 * box_size[1] / delta
+    tmp *= g * (rho2 - rho1) * delta
+    return P0 - rho1 * g * y + tmp
+
+
+# ---------------------------------------------------
+
+if (box_size[0] / N[0] != box_size[1] / N[1]):
+    raise Exception(
+        "Assuming the same number of particle per unit of distance")
+
+
+# Start by generating grids of particles
+numPart = N[0] * N[1]
+
+m_tot = mass(box_size[1])
+m_part = m_tot / numPart
+
+coords = np.zeros((numPart, 3))
+h = np.ones(numPart) * 1.2348 * box_size[0] / N[0]
+m = np.ones(numPart) * m_part
+u = np.zeros(numPart)
+vel = np.zeros((numPart, 3))
+
+
+# generate grid of particles
+y_prev = 0
+# loop over y
+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] - 1e-5 * (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
+        u[index] = pressure(y_j) / (rho1 * (gamma-1.))
+
+ids = np.linspace(1, numPart, numPart)
+
+# Velocity perturbation
+ind = coords[:, 1] < 0.7
+ind = np.logical_and(ind, coords[:, 1] > 0.3)
+vel[ind, 1] = (1 + np.cos(8 * np.pi * (coords[index, 0] + 0.5 * box_size[0])))
+vel[ind, 1] += (1 + np.cos(5 * np.pi * (coords[index, 1] - 0.5 * box_size[1])))
+vel[ind, 1] *= dv
+
+# 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, 0, 0, 0, 0, 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))
+
+fileOutput.close()
diff --git a/examples/HydroTests/Rayleigh-Taylor/rayleigh_taylor.yml b/examples/HydroTests/Rayleigh-Taylor/rayleigh_taylor.yml
new file mode 100644
index 000000000..dc9f69e7a
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor/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:   6.   # 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).
+
+# 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/run.sh b/examples/HydroTests/Rayleigh-Taylor/run.sh
new file mode 100644
index 000000000..1c8ff3eaf
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor/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=1 rayleigh_taylor.yml 2>&1 | tee output.log
+
+python makeMovie.py
-- 
GitLab


From 18d8c44af74dad6ad9f882fcb4952414a7ff13eb Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Tue, 16 Jul 2019 09:49:43 +0200
Subject: [PATCH 02/16] Rayleigh: IC from abel 2011

---
 examples/HydroTests/Rayleigh-Taylor/makeIC.py | 259 +++++++++++-------
 examples/HydroTests/Rayleigh-Taylor/run.sh    |   5 +-
 2 files changed, 156 insertions(+), 108 deletions(-)

diff --git a/examples/HydroTests/Rayleigh-Taylor/makeIC.py b/examples/HydroTests/Rayleigh-Taylor/makeIC.py
index 53b9eceea..317f7d6fd 100644
--- a/examples/HydroTests/Rayleigh-Taylor/makeIC.py
+++ b/examples/HydroTests/Rayleigh-Taylor/makeIC.py
@@ -20,39 +20,52 @@
 import h5py
 import numpy as np
 from scipy.optimize import bisect
+from scipy.integrate import quad
 
-# Generates a swift IC file for the Rayleigh-Taylor vortex in a periodic box
-# following the conditions given in Hopkins 2013
+# Generates a swift IC file for the Rayleigh-Taylor instability in a periodic
+# box following the conditions given in Abel 2010
 
 # Parameters
-N = [50, 100]  # Particles along one edge
-gamma = 5./3.  # Gas adiabatic index
-delta = 0.025  # interface size
-dv = 0.025  # velocity perturbation
-rho1 = 1    # Lower region density
-rho2 = 2    # Upper region density
-P0 = 0      # integration constant of the hydrostatic equation
+N = [100, 200]  # Particles along one edge
+gamma = 7./5.  # Gas adiabatic index
+delta = 0.1  # interface size
+dv = 0.1    # velocity perturbation
+rho1 = 2    # Upper region density
+rho2 = 1    # Lower region density
 g = -0.5    # gravitational acceleration
 box_size = [0.5, 1.]  # size of the box
+
+fixed = [0.1, 0.9]  # y-range of non fixed particles
+perturbation_box = [0.3, 0.7]  # 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.
     """
-    tmp = 1 + np.exp(-(y - 0.5 * box_size[1]) / delta)
-    return rho1 + (rho2 - rho1) / tmp
+    tmp = 1 + np.exp(-2.*(y - 0.5 * box_size[1]) / delta)
+    return rho2 + (rho1 - rho2) / tmp
 
 
 def mass(y):
     """
-    Integrated Mass in x and y.
+    Integral of the density
     """
-    tmp = np.log(np.exp(y / delta) + np.exp(box_size[1] * 0.5 / delta))
-    tmp -= 0.5 * box_size[1] / delta
-    m = box_size[0] * (rho1 * y + (rho2 - rho1) * delta * tmp)
-    return m
+    return box_size[0] * quad(density, 0, y)[0]
+
+
+P0 = rho1 / 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):
@@ -73,99 +86,133 @@ def pressure(y):
     Pressure as function of the position y.
     Here we assume hydrostatic equation.
     """
-    tmp = np.log(np.exp(y / delta) + np.exp(box_size[1] * 0.5 / delta))
-    tmp -= 0.5 * box_size[1] / delta
-    tmp *= g * (rho2 - rho1) * delta
-    return P0 - rho1 * g * y + tmp
+    return P0 + g * density(y) * (y - 0.5 * box_size[1])
 
 
-# ---------------------------------------------------
-
-if (box_size[0] / N[0] != box_size[1] / N[1]):
-    raise Exception(
-        "Assuming the same number of particle per unit of distance")
+def smoothing_length(y, m):
+    """
+    Compute the smoothing length
+    """
+    return 1.23 * np.sqrt(m / density(y))
 
 
-# Start by generating grids of particles
-numPart = N[0] * N[1]
+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))
 
-m_tot = mass(box_size[1])
-m_part = m_tot / numPart
 
-coords = np.zeros((numPart, 3))
-h = np.ones(numPart) * 1.2348 * box_size[0] / N[0]
-m = np.ones(numPart) * m_part
-u = np.zeros(numPart)
-vel = np.zeros((numPart, 3))
-
-
-# generate grid of particles
-y_prev = 0
-# loop over y
-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] - 1e-5 * (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
-        u[index] = pressure(y_j) / (rho1 * (gamma-1.))
-
-ids = np.linspace(1, numPart, numPart)
-
-# Velocity perturbation
-ind = coords[:, 1] < 0.7
-ind = np.logical_and(ind, coords[:, 1] > 0.3)
-vel[ind, 1] = (1 + np.cos(8 * np.pi * (coords[index, 0] + 0.5 * box_size[0])))
-vel[ind, 1] += (1 + np.cos(5 * np.pi * (coords[index, 1] - 0.5 * box_size[1])))
-vel[ind, 1] *= dv
-
-# 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, 0, 0, 0, 0, 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))
-
-fileOutput.close()
+def vy(x, y):
+    """
+    Velocity along the direction y
+    """
+    ind = y < perturbation_box[1]
+    ind = np.logical_and(ind, y > perturbation_box[0])
+
+    v = np.zeros(len(x))
+
+    v[ind] = (1 + np.cos(8 * np.pi * (x[ind] + 0.5 * box_size[0])))
+    v[ind] *= (1 + np.cos(5 * np.pi * (y[ind] - 0.5 * box_size[1])))
+    v[ind] *= 0.25 * 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
+    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] - 1e-5 * (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 configure with --fix-below-id={}".format(uni_id+1))
+    ids[ids == 0] = np.linspace(uni_id, numPart, numPart-uni_id+1)
+
+    u = pressure(coords[:, 1]) / (density(coords[:, 1]) * (gamma-1.))
+    if (u.min() < 0):
+        raise Exception("P0 too small")
+
+    # smoothing length
+    h = smoothing_length(coords[:, 1], m)
+
+    # density
+    rho = density(coords[:, 1])
+
+    # 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, 0, 0, 0, 0, 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()
+
+    print("Time scale: ", 1./growth_rate())
diff --git a/examples/HydroTests/Rayleigh-Taylor/run.sh b/examples/HydroTests/Rayleigh-Taylor/run.sh
index 1c8ff3eaf..ee4dd4619 100644
--- a/examples/HydroTests/Rayleigh-Taylor/run.sh
+++ b/examples/HydroTests/Rayleigh-Taylor/run.sh
@@ -10,6 +10,7 @@ then
 fi
 
 # Run SWIFT
-../../swift --hydro --external-gravity  --threads=1 rayleigh_taylor.yml 2>&1 | tee output.log
+rm rayleigh_taylor_*
+../../swift --hydro --external-gravity  --threads=8 rayleigh_taylor.yml 2>&1 | tee output.log
 
-python makeMovie.py
+# python makeMovie.py
-- 
GitLab


From c0d859289568ecf18eb441acf25b935b7165fee6 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Tue, 16 Jul 2019 15:40:02 +0200
Subject: [PATCH 03/16] Rayleigh: add constant grav acceleration + boundary
 particles

---
 configure.ac                       |  19 +++-
 src/potential.h                    |   2 +
 src/potential/constant/potential.h | 145 +++++++++++++++++++++++++++++
 src/runner.c                       |  12 +++
 4 files changed, 177 insertions(+), 1 deletion(-)
 create mode 100644 src/potential/constant/potential.h

diff --git a/configure.ac b/configure.ac
index c715fb34b..79f0841e6 100644
--- a/configure.ac
+++ b/configure.ac
@@ -372,6 +372,19 @@ 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],
+     [Enable the boundaries particles (id == 0)]
+   )],
+   [boundary_particles="$enableval"],
+   [boundary_particles="no"]
+)
+if test "$no_acceleration_below_id" != "no"; then
+   AC_DEFINE_UNQUOTED([SWIFT_NO_GRAVITY_BELOW_ID], [1] ,[Particles with smaller ID than this will have zero gravity forces])
+   AC_DEFINE_UNQUOTED([SWIFT_BOUNDARY_PARTICLES], [1] ,[Particles with ID == 0 are considered as boundaries.])
+fi
+
 # Check whether we have any of the ARM v8.1 tick timers
 AX_ASM_ARM_PMCCNTR
 AX_ASM_ARM_CNTVCT
@@ -1813,7 +1826,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"]
@@ -1846,6 +1859,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])
    ;;
@@ -2024,5 +2040,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/src/potential.h b/src/potential.h
index 59567fe92..9d1418cad 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 000000000..bfa047d16
--- /dev/null
+++ b/src/potential/constant/potential.h
@@ -0,0 +1,145 @@
+/*******************************************************************************
+ * 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 <math.h>
+#include <float.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 d42ebb4e9..d1c15cabc 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3614,6 +3614,18 @@ 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 gpart */
+        const long long id = p->id;
+
+        /* Cancel gravity forces of these particles */
+        if (id == 0) {
+
+          /* Don't move ! */
+	  hydro_reset_acceleration(p);
+        }
+#endif
       }
     }
   }
-- 
GitLab


From 973b1aec57a03dac5f1da4e5f5c500b7988b2c39 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Wed, 17 Jul 2019 14:47:27 +0200
Subject: [PATCH 04/16] Rayleigh: code is working

---
 examples/HydroTests/Rayleigh-Taylor_2D/README |  14 ++
 .../makeIC.py                                 | 111 ++++++----
 .../Rayleigh-Taylor_2D/makeMovie.py           | 202 ++++++++++++++++++
 .../Rayleigh-Taylor_2D/plotInitialProfile.py  |  37 ++++
 .../rayleigh_taylor.yml                       |   6 +-
 .../run.sh                                    |   3 +-
 6 files changed, 329 insertions(+), 44 deletions(-)
 create mode 100644 examples/HydroTests/Rayleigh-Taylor_2D/README
 rename examples/HydroTests/{Rayleigh-Taylor => Rayleigh-Taylor_2D}/makeIC.py (67%)
 create mode 100644 examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
 create mode 100644 examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py
 rename examples/HydroTests/{Rayleigh-Taylor => Rayleigh-Taylor_2D}/rayleigh_taylor.yml (86%)
 rename examples/HydroTests/{Rayleigh-Taylor => Rayleigh-Taylor_2D}/run.sh (88%)

diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/README b/examples/HydroTests/Rayleigh-Taylor_2D/README
new file mode 100644
index 000000000..1b3345f84
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/README
@@ -0,0 +1,14 @@
+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`.
\ No newline at end of file
diff --git a/examples/HydroTests/Rayleigh-Taylor/makeIC.py b/examples/HydroTests/Rayleigh-Taylor_2D/makeIC.py
similarity index 67%
rename from examples/HydroTests/Rayleigh-Taylor/makeIC.py
rename to examples/HydroTests/Rayleigh-Taylor_2D/makeIC.py
index 317f7d6fd..3d2bec517 100644
--- a/examples/HydroTests/Rayleigh-Taylor/makeIC.py
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/makeIC.py
@@ -20,23 +20,21 @@
 import h5py
 import numpy as np
 from scipy.optimize import bisect
-from scipy.integrate import quad
 
 # Generates a swift IC file for the Rayleigh-Taylor instability in a periodic
-# box following the conditions given in Abel 2010
+# box following the conditions given in Saitoh and Makino 2013: 1202.4277v3
 
 # Parameters
-N = [100, 200]  # Particles along one edge
+N = [128, 192]  # Particles along one edge
 gamma = 7./5.  # Gas adiabatic index
-delta = 0.1  # interface size
-dv = 0.1    # velocity perturbation
-rho1 = 2    # Upper region density
-rho2 = 1    # Lower region density
+dv = 0.025    # velocity perturbation
+rho_h = 2    # high density region
+rho_l = 1    # Low density region
 g = -0.5    # gravitational acceleration
-box_size = [0.5, 1.]  # size of the box
+box_size = [1., 1.5]  # size of the box
 
-fixed = [0.1, 0.9]  # y-range of non fixed particles
-perturbation_box = [0.3, 0.7]  # y-range for the velocity perturbation
+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"
 
 
@@ -50,18 +48,43 @@ def density(y):
     """
     Mass density as function of the position y.
     """
-    tmp = 1 + np.exp(-2.*(y - 0.5 * box_size[1]) / delta)
-    return rho2 + (rho1 - rho2) / tmp
+    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
     """
-    return box_size[0] * quad(density, 0, y)[0]
+    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 = rho1 / gamma  # integration constant of the hydrostatic equation
+P0 = rho_h / gamma  # integration constant of the hydrostatic equation
 numPart = N[0] * N[1]
 
 m_tot = mass(box_size[1])
@@ -81,19 +104,27 @@ def inv_mass(m):
     return bisect(f, left, right, args=(m))
 
 
-def pressure(y):
+def entropy(y):
     """
-    Pressure as function of the position y.
-    Here we assume hydrostatic equation.
+    Entropy as function of the position y.
+    Here we assume isoentropic medium.
     """
-    return P0 + g * density(y) * (y - 0.5 * box_size[1])
+    if isinstance(y, float) or isinstance(y, int):
+        y = np.array(y)
 
+    ind = y < 0.5 * box_size[1]
 
-def smoothing_length(y, m):
+    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 / density(y))
+    return 1.23 * np.sqrt(m / rho)
 
 
 def growth_rate():
@@ -101,7 +132,7 @@ def growth_rate():
     Compute the growth rate of the instability.
     Assumes a wavelength equal to the boxsize.
     """
-    ymin = 0
+    ymin = 0.
     ymax = box_size[1]
     A = density(ymax) - density(ymin)
     A /= density(ymax) + density(ymin)
@@ -112,14 +143,14 @@ def vy(x, y):
     """
     Velocity along the direction y
     """
-    ind = y < perturbation_box[1]
-    ind = np.logical_and(ind, y > perturbation_box[0])
+    ind = y > perturbation_box[0]
+    ind = np.logical_and(ind, y < perturbation_box[1])
 
     v = np.zeros(len(x))
 
-    v[ind] = (1 + np.cos(8 * np.pi * (x[ind] + 0.5 * box_size[0])))
-    v[ind] *= (1 + np.cos(5 * np.pi * (y[ind] - 0.5 * box_size[1])))
-    v[ind] *= 0.25 * dv
+    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
 
 
@@ -130,16 +161,18 @@ if __name__ == "__main__":
     m = np.ones(numPart) * m_part
     u = np.zeros(numPart)
     vel = np.zeros((numPart, 3))
-    ids = np.zeros(numPart)
+    ids = np.ones(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] - 1e-5 * (box_size[1] - y_prev)
+            y_j = box_size[1] - eps * (box_size[1] - y_prev)
         else:
             y_j = inv_mass(m_y)
         y_prev = y_j
@@ -154,22 +187,22 @@ if __name__ == "__main__":
             coords[index, 0] = x
             coords[index, 1] = y_j
             if (y_j < fixed[0] or y_j > fixed[1]):
-                ids[index] = uni_id
+                ids[index] = 0
                 uni_id += 1
 
-    print("You need to configure with --fix-below-id={}".format(uni_id+1))
-    ids[ids == 0] = np.linspace(uni_id, numPart, numPart-uni_id+1)
-
-    u = pressure(coords[:, 1]) / (density(coords[:, 1]) * (gamma-1.))
-    if (u.min() < 0):
-        raise Exception("P0 too small")
-
-    # smoothing length
-    h = smoothing_length(coords[:, 1], m)
+    N = numPart-uni_id+1
+    ids[ids != 0] = np.linspace(1, N+1, N)
 
     # 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])
 
@@ -185,7 +218,7 @@ if __name__ == "__main__":
     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, 0, 0, 0, 0, 0]
+    grp.attrs["Flag_Entropy_ICs"] = 0
     grp.attrs["Dimension"] = 2
 
     # Units
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
new file mode 100644
index 000000000..ac3e29754
--- /dev/null
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
@@ -0,0 +1,202 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2019 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, mask
+from swiftsimio.visualisation import project_gas_pixel_grid
+from matplotlib.animation import FuncAnimation
+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
+
+    # map = p_map
+except:
+    print("Try installing the p_tqdm package to make movie frames in parallel")
+    pass
+
+reverse_axis = False
+
+
+def project(data, m_res, property):
+    tmp = project_gas_pixel_grid(data, m_res, property)
+
+    if reverse_axis:
+        return tmp
+    else:
+        return tmp.T
+
+
+def load_and_make_image(filename, res, property):
+    image = np.zeros(res, dtype=np.float32)
+    image_none = np.zeros(res, dtype=np.float32)
+    m_res = min(res)
+
+    # first part of the image
+    m = mask(filename)
+    ylim = np.array([0., 2./3.]) * m.metadata.boxsize[1]
+    m.constrain_spatial([None, ylim, None])
+
+    # data = load(filename, mask=m)
+    data = load(filename)
+    y = data.gas.coordinates.value[:, 1]
+    ind = y > ylim[0]
+    ind = np.logical_and(ind, y < ylim[1])
+    data.gas.particle_ids.values = data.gas.particle_ids[ind]
+    data.gas.density.values = data.gas.density[ind]
+    data.gas.coordinates.values = data.gas.coordinates[ind, :]
+    data.gas.smoothing_length.values = data.gas.smoothing_length[ind]
+    image[:m_res, :m_res] = project(data, m_res, property)
+    tmp = project(data, m_res, None)
+    image_none[:m_res, :m_res] = tmp
+
+    # second part of the image
+    m = mask(filename)
+    ylim = np.array([1./3., 1.]) * m.metadata.boxsize[1]
+    m.constrain_spatial([None, ylim, None])
+
+    # data = load(filename, mask=m)
+    data = load(filename)
+    y = data.gas.coordinates.value[:, 1]
+    ind = y > ylim[0]
+    ind = np.logical_and(ind, y < ylim[1])
+    data.gas.particle_ids = data.gas.particle_ids[ind]
+    data.gas.density = data.gas.density[ind]
+    data.gas.coordinates = data.gas.coordinates[ind, :]
+    data.gas.smoothing_length = data.gas.smoothing_length[ind]
+    image[:m_res, :m_res] = project(data, m_res, property)
+    tmp2 = project(data, m_res, None)
+    image_none[:m_res, :m_res] = tmp2
+
+    if reverse_axis:
+        image[:, -m_res:] = project(data, m_res, property)
+        image_none[:, -m_res:] = project(data, m_res, None)
+    else:
+        image[-m_res:, :] = project(data, m_res, property)
+        image_none[-m_res:, :] = project(data, m_res, None)
+
+    if np.sum(image_none == 0):
+        plt.imshow(tmp)
+        plt.figure()
+        plt.imshow(tmp2)
+        plt.colorbar()
+        plt.show()
+    return image / image_none
+
+
+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}.
+    """
+
+    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 = map(baked_in_load, list(range(start, stop)))
+
+    vmax = max(list(map(max, frames)))
+    vmin = min(list(map(min, frames)))
+
+    fig, ax = plt.subplots(figsize=(1, 1))
+    ax.axis("off")
+    fig.subplots_adjust(0, 0, 1, 1)
+
+    image = ax.imshow(frames[0], origin="lower", vmax=vmax, vmin=vmin)
+
+    def frame(n):
+        image.set_array(frames[n])
+
+        return (image,)
+
+    animation = FuncAnimation(fig, frame, range(0, stop-start), interval=40)
+    animation.save(output_filename, dpi=resolution[0])
+
+
+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)
+    if reverse_axis:
+        vars.resolution = [vars.resolution, yres]
+    else:
+        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 000000000..e89c4e852
--- /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/rayleigh_taylor.yml b/examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml
similarity index 86%
rename from examples/HydroTests/Rayleigh-Taylor/rayleigh_taylor.yml
rename to examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml
index dc9f69e7a..568567756 100644
--- a/examples/HydroTests/Rayleigh-Taylor/rayleigh_taylor.yml
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml
@@ -9,8 +9,8 @@ InternalUnitSystem:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   6.   # The end time of the simulation (in internal units).
-  dt_min:     1e-6  # The minimal time-step size 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
@@ -26,7 +26,7 @@ Statistics:
 # 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.
+  CFL_condition:         0.3      # Courant-Friedrich-Levy condition for time integration.
   
 # Parameters related to the initial conditions
 InitialConditions:
diff --git a/examples/HydroTests/Rayleigh-Taylor/run.sh b/examples/HydroTests/Rayleigh-Taylor_2D/run.sh
similarity index 88%
rename from examples/HydroTests/Rayleigh-Taylor/run.sh
rename to examples/HydroTests/Rayleigh-Taylor_2D/run.sh
index ee4dd4619..3044c1cc1 100644
--- a/examples/HydroTests/Rayleigh-Taylor/run.sh
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/run.sh
@@ -10,7 +10,6 @@ then
 fi
 
 # Run SWIFT
-rm rayleigh_taylor_*
 ../../swift --hydro --external-gravity  --threads=8 rayleigh_taylor.yml 2>&1 | tee output.log
 
-# python makeMovie.py
+python makeMovie.py -i 0 -t 1001
-- 
GitLab


From 6f4caf6fb770e67990deb05a0ee5bd672b4723b4 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Wed, 17 Jul 2019 16:49:24 +0200
Subject: [PATCH 05/16] RayleighTaylor: improve movie

---
 .../Rayleigh-Taylor_2D/makeMovie.py           | 95 +++++++++----------
 examples/HydroTests/Rayleigh-Taylor_2D/run.sh |  2 +-
 2 files changed, 45 insertions(+), 52 deletions(-)

diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
index ac3e29754..5e1d0d7c4 100644
--- a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
@@ -1,6 +1,7 @@
 ###############################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+# 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
@@ -17,8 +18,8 @@
 #
 ##############################################################################
 
-from swiftsimio import load, mask
-from swiftsimio.visualisation import project_gas_pixel_grid
+from swiftsimio import load
+from swiftsimio.visualisation import scatter
 from matplotlib.animation import FuncAnimation
 import matplotlib.pyplot as plt
 
@@ -29,7 +30,7 @@ try:
     # Try and load this, otherwise we're stuck with serial
     from p_tqdm import p_map
 
-    # map = p_map
+    map = p_map
 except:
     print("Try installing the p_tqdm package to make movie frames in parallel")
     pass
@@ -37,70 +38,62 @@ except:
 reverse_axis = False
 
 
-def project(data, m_res, property):
-    tmp = project_gas_pixel_grid(data, m_res, property)
+def project(data, m_res, property, ylim):
+    x, y, _ = data.gas.coordinates.value.T
 
     if reverse_axis:
-        return tmp
+        x, y = y, x
+
+    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)
+
+    if reverse_axis:
+        return image
     else:
-        return tmp.T
+        return image.T
 
 
 def load_and_make_image(filename, res, property):
     image = np.zeros(res, dtype=np.float32)
-    image_none = np.zeros(res, dtype=np.float32)
     m_res = min(res)
+    border = int(0.2 * m_res)
 
     # first part of the image
-    m = mask(filename)
-    ylim = np.array([0., 2./3.]) * m.metadata.boxsize[1]
-    m.constrain_spatial([None, ylim, None])
+    ylim = np.array([0., 1.])
 
-    # data = load(filename, mask=m)
     data = load(filename)
-    y = data.gas.coordinates.value[:, 1]
-    ind = y > ylim[0]
-    ind = np.logical_and(ind, y < ylim[1])
-    data.gas.particle_ids.values = data.gas.particle_ids[ind]
-    data.gas.density.values = data.gas.density[ind]
-    data.gas.coordinates.values = data.gas.coordinates[ind, :]
-    data.gas.smoothing_length.values = data.gas.smoothing_length[ind]
-    image[:m_res, :m_res] = project(data, m_res, property)
-    tmp = project(data, m_res, None)
-    image_none[:m_res, :m_res] = tmp
+    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
-    m = mask(filename)
-    ylim = np.array([1./3., 1.]) * m.metadata.boxsize[1]
-    m.constrain_spatial([None, ylim, None])
-
-    # data = load(filename, mask=m)
-    data = load(filename)
-    y = data.gas.coordinates.value[:, 1]
-    ind = y > ylim[0]
-    ind = np.logical_and(ind, y < ylim[1])
-    data.gas.particle_ids = data.gas.particle_ids[ind]
-    data.gas.density = data.gas.density[ind]
-    data.gas.coordinates = data.gas.coordinates[ind, :]
-    data.gas.smoothing_length = data.gas.smoothing_length[ind]
-    image[:m_res, :m_res] = project(data, m_res, property)
-    tmp2 = project(data, m_res, None)
-    image_none[:m_res, :m_res] = tmp2
+    ylim = np.array([0.5, 1.5])
 
+    left = -m_res + border
     if reverse_axis:
-        image[:, -m_res:] = project(data, m_res, property)
-        image_none[:, -m_res:] = project(data, m_res, None)
+        image[:, left:] = project(data, m_res, property, ylim)[:, left:]
+        if property != "density":
+            image[:, left:] /= project(data, m_res, None, ylim)[:, left:]
     else:
-        image[-m_res:, :] = project(data, m_res, property)
-        image_none[-m_res:, :] = project(data, m_res, None)
-
-    if np.sum(image_none == 0):
-        plt.imshow(tmp)
-        plt.figure()
-        plt.imshow(tmp2)
-        plt.colorbar()
-        plt.show()
-    return image / image_none
+        image[left:, :] = project(data, m_res, property, ylim)[left:, :]
+        if property != "density":
+            image[left:, :] /= project(data, m_res, None, ylim)[left:, :]
+
+    return image
 
 
 def create_movie(filename, start, stop, resolution, property, output_filename):
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/run.sh b/examples/HydroTests/Rayleigh-Taylor_2D/run.sh
index 3044c1cc1..0c8d3eaab 100644
--- a/examples/HydroTests/Rayleigh-Taylor_2D/run.sh
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/run.sh
@@ -12,4 +12,4 @@ fi
 # Run SWIFT
 ../../swift --hydro --external-gravity  --threads=8 rayleigh_taylor.yml 2>&1 | tee output.log
 
-python makeMovie.py -i 0 -t 1001
+python makeMovie.py -i 0 -f 1001
-- 
GitLab


From 24bf5c20a22ef31583c51b282ea8cca9a3e0364e Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Thu, 18 Jul 2019 10:29:17 +0200
Subject: [PATCH 06/16] Improve rayleigh movie

---
 .../Rayleigh-Taylor_2D/makeMovie.py           | 129 ++++++++++++++----
 1 file changed, 102 insertions(+), 27 deletions(-)

diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
index 5e1d0d7c4..10e53f226 100644
--- a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
@@ -21,6 +21,7 @@
 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
@@ -35,15 +36,35 @@ except:
     print("Try installing the p_tqdm package to make movie frames in parallel")
     pass
 
-reverse_axis = False
+
+info_frames = 40
+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
 
-    if reverse_axis:
-        x, y = y, x
-
     mask = np.logical_and(y >= ylim[0], y <= ylim[1])
 
     x = x[mask]
@@ -61,10 +82,7 @@ def project(data, m_res, property, ylim):
 
     image = scatter(x=x, y=y, m=quant, h=h, res=m_res)
 
-    if reverse_axis:
-        return image
-    else:
-        return image.T
+    return image.T
 
 
 def load_and_make_image(filename, res, property):
@@ -84,16 +102,16 @@ def load_and_make_image(filename, res, property):
     ylim = np.array([0.5, 1.5])
 
     left = -m_res + border
-    if reverse_axis:
-        image[:, left:] = project(data, m_res, property, ylim)[:, left:]
-        if property != "density":
-            image[:, left:] /= project(data, m_res, None, ylim)[:, left:]
-    else:
-        image[left:, :] = project(data, m_res, property, ylim)[left:, :]
-        if property != "density":
-            image[left:, :] /= project(data, m_res, None, ylim)[left:, :]
+    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
+
 
-    return image
+def time_formatter(metadata):
+    return f"$t = {metadata.t:2.2f}$"
 
 
 def create_movie(filename, start, stop, resolution, property, output_filename):
@@ -105,29 +123,89 @@ def create_movie(filename, start, stop, resolution, property, output_filename):
     {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 = map(baked_in_load, list(range(start, stop)))
+    frames, metadata = zip(*map(baked_in_load, list(range(start, stop))))
 
     vmax = max(list(map(max, frames)))
     vmin = min(list(map(min, frames)))
 
-    fig, ax = plt.subplots(figsize=(1, 1))
+    fig, ax = plt.subplots(figsize=(8, 1.5 * 8), dpi=resolution[0] // 8)
     ax.axis("off")
     fig.subplots_adjust(0, 0, 1, 1)
 
-    image = ax.imshow(frames[0], origin="lower", vmax=vmax, vmin=vmin)
+    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,
+    )
+
+    info_text = ax.text(
+        0.025 * 0.25, 0.975, name, **text_args, va="top", ha="left",
+        transform=ax.transAxes
+    )
 
     def frame(n):
-        image.set_array(frames[n])
+        if n >= 0:
+            image.set_array(frames[n])
+            description_text.set_text("")
+            time_text.set_text(time_formatter(metadata[n]))
 
         return (image,)
 
-    animation = FuncAnimation(fig, frame, range(0, stop-start), interval=40)
-    animation.save(output_filename, dpi=resolution[0])
+    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__":
@@ -180,10 +258,7 @@ if __name__ == "__main__":
     vars = parser.parse_args()
 
     yres = int(1.5 * vars.resolution)
-    if reverse_axis:
-        vars.resolution = [vars.resolution, yres]
-    else:
-        vars.resolution = [yres, vars.resolution]
+    vars.resolution = [yres, vars.resolution]
 
     create_movie(
         filename=vars.snapshot,
-- 
GitLab


From 4805968a34277b72a1fd18086eacff568f7113e5 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Thu, 18 Jul 2019 11:10:24 +0200
Subject: [PATCH 07/16] Rayleigh: update the CFL condition

---
 examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml b/examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml
index 568567756..9e843c494 100644
--- a/examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/rayleigh_taylor.yml
@@ -26,7 +26,7 @@ Statistics:
 # 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.3      # Courant-Friedrich-Levy condition for time integration.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
   
 # Parameters related to the initial conditions
 InitialConditions:
-- 
GitLab


From 7415e08fde053490f6c43aebc63104ba33357d70 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Fri, 19 Jul 2019 08:36:51 +0200
Subject: [PATCH 08/16] Rayleigh: add prepare force for gizmo

---
 src/runner.c | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/runner.c b/src/runner.c
index d1c15cabc..6d7ee331a 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3624,6 +3624,8 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
 
           /* Don't move ! */
 	  hydro_reset_acceleration(p);
+	  hydro_prepare_force(p, &c->hydro.xparts[k], cosmo,
+			      e->hydro_properties, 0);
         }
 #endif
       }
-- 
GitLab


From cf2c14d99dc56c82a5719c9182241602049a92dc Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Fri, 19 Jul 2019 09:00:33 +0200
Subject: [PATCH 09/16] Rayleigh: add comment

---
 src/runner.c | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/runner.c b/src/runner.c
index 6d7ee331a..fb31c57ac 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3624,6 +3624,8 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
 
           /* Don't move ! */
 	  hydro_reset_acceleration(p);
+
+	  /* Required for GIZMO */
 	  hydro_prepare_force(p, &c->hydro.xparts[k], cosmo,
 			      e->hydro_properties, 0);
         }
-- 
GitLab


From 4add1e0dda57c9706ec51e0adf7b68c01126afa9 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Mon, 22 Jul 2019 13:47:39 +0200
Subject: [PATCH 10/16] Rayleigh: Add argument to --enable-boundary-particles

---
 configure.ac                                  | 10 +++---
 .../generate_movie_from_png.sh                | 14 ++++++++
 .../HydroTests/Rayleigh-Taylor_2D/makeIC.py   | 11 +++----
 .../Rayleigh-Taylor_2D/makeMovie.py           | 32 ++++++++++++-------
 src/runner.c                                  |  2 +-
 5 files changed, 47 insertions(+), 22 deletions(-)
 create mode 100644 examples/HydroTests/Rayleigh-Taylor_2D/generate_movie_from_png.sh

diff --git a/configure.ac b/configure.ac
index 79f0841e6..8e717e17c 100644
--- a/configure.ac
+++ b/configure.ac
@@ -375,14 +375,16 @@ fi
 # Check if we want to use boundary particles.
 AC_ARG_ENABLE([boundary-particles],
    [AS_HELP_STRING([--enable-boundary-particles],
-     [Enable the boundaries particles (id == 0)]
+     [Set all particles with an ID smaller than @<:@N@:>@ as boundary particles.]
    )],
    [boundary_particles="$enableval"],
    [boundary_particles="no"]
 )
-if test "$no_acceleration_below_id" != "no"; then
-   AC_DEFINE_UNQUOTED([SWIFT_NO_GRAVITY_BELOW_ID], [1] ,[Particles with smaller ID than this will have zero gravity forces])
-   AC_DEFINE_UNQUOTED([SWIFT_BOUNDARY_PARTICLES], [1] ,[Particles with ID == 0 are considered as boundaries.])
+if test "$no_gravity_below_id" == "yes"; then
+   AC_MSG_ERROR(Need to specify the ID below which particles get zero forces when using --enable-boundary-particles!)
+elif test "$no_acceleration_below_id" != "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 ID == 0 are considered as boundaries.])
 fi
 
 # Check whether we have any of the ARM v8.1 tick timers
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 000000000..f97e447dc
--- /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
index 3d2bec517..d75d5f215 100644
--- a/examples/HydroTests/Rayleigh-Taylor_2D/makeIC.py
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/makeIC.py
@@ -161,7 +161,7 @@ if __name__ == "__main__":
     m = np.ones(numPart) * m_part
     u = np.zeros(numPart)
     vel = np.zeros((numPart, 3))
-    ids = np.ones(numPart)
+    ids = np.zeros(numPart)
 
     # generate grid of particles
     y_prev = 0
@@ -187,11 +187,12 @@ if __name__ == "__main__":
             coords[index, 0] = x
             coords[index, 1] = y_j
             if (y_j < fixed[0] or y_j > fixed[1]):
-                ids[index] = 0
+                ids[index] = uni_id
                 uni_id += 1
 
-    N = numPart-uni_id+1
-    ids[ids != 0] = np.linspace(1, N+1, N)
+    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])
@@ -247,5 +248,3 @@ if __name__ == "__main__":
     ds[()] = rho.reshape((numPart, 1))
 
     fileOutput.close()
-
-    print("Time scale: ", 1./growth_rate())
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
index 10e53f226..b00c29469 100644
--- a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py
@@ -30,14 +30,14 @@ import numpy as np
 try:
     # Try and load this, otherwise we're stuck with serial
     from p_tqdm import p_map
-
-    map = 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")
 
 
@@ -133,10 +133,10 @@ def create_movie(filename, start, stop, resolution, property, output_filename):
         return load_and_make_image(f, resolution, property)
 
     # Make frames in parallel (reading also parallel!)
-    frames, metadata = zip(*map(baked_in_load, list(range(start, stop))))
+    frames, metadata = zip(*p_map(baked_in_load, list(range(start, stop))))
 
-    vmax = max(list(map(max, frames)))
-    vmin = min(list(map(min, frames)))
+    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")
@@ -162,7 +162,7 @@ def create_movie(filename, start, stop, resolution, property, output_filename):
         transform=ax.transAxes,
     )
 
-    info_text = ax.text(
+    ax.text(
         0.025 * 0.25, 0.975, name, **text_args, va="top", ha="left",
         transform=ax.transAxes
     )
@@ -173,11 +173,21 @@ def create_movie(filename, start, stop, resolution, property, output_filename):
             description_text.set_text("")
             time_text.set_text(time_formatter(metadata[n]))
 
-        return (image,)
+        if generate_png:
+            name = filename + "_{:04d}".format(n+info_frames)
+            fig.savefig(name + ".png")
+
+        else:
+            return (image,)
 
-    animation = FuncAnimation(fig, frame, range(-info_frames, stop-start),
-                              interval=40)
-    animation.save(output_filename)
+    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):
@@ -186,7 +196,7 @@ def get_simulation_information(metadata):
     """
     viscosity = metadata.viscosity_info
     diffusion = metadata.diffusion_info
-    
+
     output = (
         "$\\bf{Rayleigh-Taylor}$ $\\bf{Instability}$\n\n"
         "$\\bf{SWIFT}$\n"
diff --git a/src/runner.c b/src/runner.c
index fb31c57ac..7123eb281 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3620,7 +3620,7 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
         const long long id = p->id;
 
         /* Cancel gravity forces of these particles */
-        if (id == 0) {
+        if (id < SWIFT_BOUNDARY_PARTICLES) {
 
           /* Don't move ! */
 	  hydro_reset_acceleration(p);
-- 
GitLab


From 9460633c8b06c898eaae99ac837ad736cdd64583 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Mon, 22 Jul 2019 13:50:27 +0200
Subject: [PATCH 11/16] Update comment

---
 configure.ac | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configure.ac b/configure.ac
index 8e717e17c..c0bbd97cd 100644
--- a/configure.ac
+++ b/configure.ac
@@ -384,7 +384,7 @@ if test "$no_gravity_below_id" == "yes"; then
    AC_MSG_ERROR(Need to specify the ID below which particles get zero forces when using --enable-boundary-particles!)
 elif test "$no_acceleration_below_id" != "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 ID == 0 are considered as boundaries.])
+   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
-- 
GitLab


From ff8627ae6c573f93bd2ca4bc0431fe223c54e4e1 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Mon, 22 Jul 2019 14:06:24 +0200
Subject: [PATCH 12/16] Fix configure

---
 configure.ac | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/configure.ac b/configure.ac
index c0bbd97cd..9bee8288c 100644
--- a/configure.ac
+++ b/configure.ac
@@ -380,9 +380,9 @@ AC_ARG_ENABLE([boundary-particles],
    [boundary_particles="$enableval"],
    [boundary_particles="no"]
 )
-if test "$no_gravity_below_id" == "yes"; then
+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 "$no_acceleration_below_id" != "no"; then
+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
-- 
GitLab


From 6e544c7653afbb66eb5c75206b7481ced533ce58 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Mon, 22 Jul 2019 15:55:03 +0200
Subject: [PATCH 13/16] Rayleigh: Make required modifications

---
 src/potential/constant/potential.h | 2 +-
 src/runner.c                       | 6 ++++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/potential/constant/potential.h b/src/potential/constant/potential.h
index bfa047d16..5ec10439a 100644
--- a/src/potential/constant/potential.h
+++ b/src/potential/constant/potential.h
@@ -120,7 +120,7 @@ static INLINE void potential_init_backend(
   /* 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;
+  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
diff --git a/src/runner.c b/src/runner.c
index 7123eb281..fd7749db6 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3616,18 +3616,20 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
         hydro_end_force(p, cosmo);
 #ifdef SWIFT_BOUNDARY_PARTICLES
 
-        /* Get the ID of the gpart */
+        /* Get the ID of the part */
         const long long id = p->id;
 
-        /* Cancel gravity forces of these particles */
+        /* Cancel hdyro forces of these particles */
         if (id < SWIFT_BOUNDARY_PARTICLES) {
 
           /* Don't move ! */
 	  hydro_reset_acceleration(p);
 
+#ifdef EXTRA_HYDRO_LOOP
 	  /* Required for GIZMO */
 	  hydro_prepare_force(p, &c->hydro.xparts[k], cosmo,
 			      e->hydro_properties, 0);
+#endif
         }
 #endif
       }
-- 
GitLab


From 5cda99ff986bee083a95df94358ad7c7df393b2c Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Tue, 23 Jul 2019 07:56:51 +0200
Subject: [PATCH 14/16] Rayleigh: add info boundary particle

---
 examples/HydroTests/Rayleigh-Taylor_2D/README | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/README b/examples/HydroTests/Rayleigh-Taylor_2D/README
index 1b3345f84..b1ac67c0c 100644
--- a/examples/HydroTests/Rayleigh-Taylor_2D/README
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/README
@@ -11,4 +11,7 @@ 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`.
\ No newline at end of file
+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.
-- 
GitLab


From 49ed4d59ec40d9334643481f5f3be6153a44b76c Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Tue, 23 Jul 2019 09:11:05 +0100
Subject: [PATCH 15/16] Applied code formatting tool.

---
 src/potential/constant/potential.h | 16 +++++++---------
 src/runner.c                       |  8 ++++----
 2 files changed, 11 insertions(+), 13 deletions(-)

diff --git a/src/potential/constant/potential.h b/src/potential/constant/potential.h
index 5ec10439a..f91d2b3fc 100644
--- a/src/potential/constant/potential.h
+++ b/src/potential/constant/potential.h
@@ -23,8 +23,8 @@
 #include "../config.h"
 
 /* Some standard headers. */
-#include <math.h>
 #include <float.h>
+#include <math.h>
 
 /* Local includes. */
 #include "error.h"
@@ -91,11 +91,9 @@ 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];
-    
+  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;
 }
 
@@ -114,15 +112,15 @@ static INLINE void potential_init_backend(
     struct external_potential* potential) {
 
   /* Read in the acceleration */
-  parser_get_param_double_array(parameter_file, "ConstantPotential:g_cgs",
-                                3, potential->g);
+  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++) {
+  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;
   }
diff --git a/src/runner.c b/src/runner.c
index fd7749db6..e3fde24dc 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3623,12 +3623,12 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
         if (id < SWIFT_BOUNDARY_PARTICLES) {
 
           /* Don't move ! */
-	  hydro_reset_acceleration(p);
+          hydro_reset_acceleration(p);
 
 #ifdef EXTRA_HYDRO_LOOP
-	  /* Required for GIZMO */
-	  hydro_prepare_force(p, &c->hydro.xparts[k], cosmo,
-			      e->hydro_properties, 0);
+          /* Required for GIZMO */
+          hydro_prepare_force(p, &c->hydro.xparts[k], cosmo,
+                              e->hydro_properties, 0);
 #endif
         }
 #endif
-- 
GitLab


From 3c5ce231077bb4f797e5512f02016e8250f1b900 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Tue, 23 Jul 2019 09:13:34 +0100
Subject: [PATCH 16/16] Applied code formatting tool.

---
 configure.ac | 2 +-
 src/runner.c | 6 ++++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/configure.ac b/configure.ac
index 9bee8288c..d94544ffe 100644
--- a/configure.ac
+++ b/configure.ac
@@ -375,7 +375,7 @@ 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.]
+     [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"]
diff --git a/src/runner.c b/src/runner.c
index e3fde24dc..c8235ae33 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3614,6 +3614,7 @@ 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 */
@@ -3625,8 +3626,9 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
           /* Don't move ! */
           hydro_reset_acceleration(p);
 
-#ifdef EXTRA_HYDRO_LOOP
-          /* Required for GIZMO */
+#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
-- 
GitLab