Commit cf616d4b authored by Loic Hausammann's avatar Loic Hausammann Committed by Matthieu Schaller
Browse files

Rayleigh

parent 6052ad4a
......@@ -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
------------------------])
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.
#!/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
###############################################################################
# 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()
###############################################################################
# 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,
)