Skip to content
GitLab
Commit 2aaae570 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added the Keplerian Ring hydro test along with two new relevant softened potentials. See ls.

parent e405189e
Expand all Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ 2aaae570
......@@ -275,9 +275,14 @@ _minted*
*.sympy
sympy-plots-for-*.tex/
# python
*.pyc
# todonotes
*.tdo
# xindy
*.xdy
# macOS
*.DS_Store
\ No newline at end of file
configure.ac
View file @ 2aaae570
......@@ -980,7 +980,7 @@ esac
# External potential
AC_ARG_WITH([ext-potential],
[AS_HELP_STRING([--with-ext-potential=<pot>],
[external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch, sine-wave default: none@:>@]
[external potential @<:@none, point-mass, point-mass-ring, point-mass-softened, isothermal, softened-isothermal, disc-patch, sine-wave, default: none@:>@]
)],
[with_potential="$withval"],
[with_potential="none"]
......@@ -1001,6 +1001,12 @@ case "$with_potential" in
sine-wave)
AC_DEFINE([EXTERNAL_POTENTIAL_SINE_WAVE], [1], [Sine wave external potential in 1D])
;;
point-mass-ring)
AC_DEFINE([EXTERNAL_POTENTIAL_POINTMASS_RING], [1], [Point mass potential for Keplerian Ring (Hopkins 2015).])
;;
point-mass-softened)
AC_DEFINE([EXTERNAL_POTENTIAL_POINTMASS_SOFT], [1], [Softened point-mass potential with form 1/(r^2 + softening^2).])
;;
*)
AC_MSG_ERROR([Unknown external potential: $with_potential])
;;
......
examples/KeplerianRing/README.md 0 → 100644
View file @ 2aaae570
Keplerian Ring Test
===================
This test uses the '2D Keplerian Ring' to assess the accuracy of SPH
calculations. For more information on the test, please see Cullen & Dehnen 2009
(http://arxiv.org/abs/1006.1524) Section 4.3.
It is key that for this test in particular that the initial conditions are
perfectly arranged (here we use the same methodology as described in the paper
referenced above).
The test uses:
+ $M = 10^10$ for a central point mass
+ $r$ up to 5 (particles created to edge of box when relevant)
+ Approximately 16000 particles
+ $c_s = 0.01 \ll v_\phi = 10$, enforced by giving them a mass of 1 unit and an
internal energy of 0.015.
Please note that the initial condition generator **requires python3 rather than
python2**, as well as the plotting script.
Code Setup
----------
You will need to recompile SWIFT with an external potential. This can be done
by using
./configure --with-ext-potential=point-mass
in the root directory of the project. We suggest leaving all other code options
as their defaults.
Please also consider using:
./configure --with-ext-potential=point-mass-softened
if you are running with the initial conditions generated with the script and
using a nonzero softening.
Initial Conditions Generation
-----------------------------
The initial condition generation is handled by ```makeIC.py```, for which
the options are available with
python3 makeIC.py --help
however the defaults are very sensible. If you wish to change the central mass,
you will need to change the external point mass in ```keplerian_ring.yml``` as
well.
There are three main types of generation, handled through `--generationmethod`:
+ `spiral`: using the original spiral method with a density distribution
created using different mass particles.
+ `gaussian`: a guassian density profile with density distribution changed
through the _distribution_ of equal mass particles.
+ `grid`: a grid with a flat density profile (similar to Hopkins 2013) imposed
on it through different particle masses.
All of them have their benefits and drawbacks, so give each a go.
Plotting
--------
Once you have ran swift (we suggest that you use the following)
../swift -g -S -s -t 16 keplerian_ring.yml 2>&1 | tee output.log
there will be around 350 ```.hdf5``` files in your directory. To check out
the results of the example use the plotting script:
python3 plotSolution.py
which will produce a png file in the directory. Check out `plotSolution.py
--help` for more options.
You can also try the `make_movie.py` script which should make a movie.
Theory
------
We use the 'spiral' technique to generate the initial conditions. To do this,
we first calculate
$$
m_i = \frac{i - \frac{1}{2}}{N}
$$
for each particle $i$ with $N$ the total number of required particles. Then the
radius of each particle is given by
$$
r_i = f^{-1}(m_i)
$$
where for us the PDF $f$ is a Gaussian.
To choose the radial positions, we then choose an initial angle $\theta_1$, and
from there
$$
\delta \theta_i = \left[2\pi\left( 1 - \frac{r_{i-1}}{r_i} \right)\right]~.
$$
### Inverse CDF of Gaussian
We need:
$$
f^{-1}(m_i) = r_i = r_{central} + \sqrt{2}\sigma \mathrm{erf}^{-1}(2m_i - 1),
$$
and thankfully scipy has an inverse error function built in.
Implementation Details
----------------------
+ The $m_i$ are generated by ```generate_m_i```.
+ The inverse Gaussian PDF is implemented as ```inverse_gaussian```.
+ The $\theta_i$ are generated by ```generate_theta_i```.
Useful Information
------------------
$G$ in simulation units (by default using this yaml) is given by:
$$
4.300970\times10^{-6}
$$
This is used to convert the mass of the central potential (used in the
parameterfile) to the $GM$ required in the initial conditions generator.
examples/KeplerianRing/keplerian_ring.yml 0 → 100644
View file @ 2aaae570
# 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: 50 # 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-3 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: keplerian_ring # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.2 # 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).
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
InitialConditions:
file_name: initial_conditions.hdf5 # The file to read
shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
shift_z: 0.
# External potential parameters
PointMassPotential:
position_x: 5. # location of external point mass in internal units
position_y: 5. # here just take this to be the centre of the ring
position_z: 5.
mass: 1.498351e7 # mass of external point mass in internal units
timestep_mult: 0.03 # controls time step
softening: 0.05
examples/KeplerianRing/makeIC.py 0 → 100644
View file @ 2aaae570
This diff is collapsed.
examples/KeplerianRing/make_movie.py 0 → 100644
View file @ 2aaae570
"""
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2017
#
# 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/>.
#
#
# -----------------------------------------------------------------------------
#
# This program creates test plots for the initial condition generator provided
# for the Keplerian Ring example.
#
###############################################################################
"""
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.animation as anim
import numpy as np
import h5py as h5
from tqdm import tqdm
def get_colours(numbers, norm=None, cm_name="viridis"):
if norm is None:
norm = matplotlib.colors.Normalize(vmax=numbers.max(), vmin=numbers.min())
cmap = matplotlib.cm.get_cmap(cm_name)
return cmap(norm(numbers)), norm
def load_data(filename, silent=True, extra=None, logextra=True, exclude=None):
if not silent:
print(f"Loading data from {filename}")
with h5.File(filename, "r") as file_handle:
coords = file_handle['PartType0']['Coordinates'][...]
time = float(file_handle['Header'].attrs['Time'])
boxsize = file_handle['Header'].attrs['BoxSize']
# For some old runs we have z=0
if np.sum(coords[:, 2]) == 0:
centre_of_box = np.array(list(boxsize[:2]/2) + [0])
else:
centre_of_box = boxsize/2
if exclude is not None:
distance_from_centre = coords - centre_of_box
r2 = np.sum(distance_from_centre * distance_from_centre, 1)
mask = r2 < (exclude * exclude)
masked = np.ma.array(coords, mask=np.array([mask]*3))
coords = masked.compressed()
if extra is not None:
other = file_handle['PartType0'][extra][...]
if exclude is not None:
masked = np.ma.array(other, mask=mask)
other = masked.compressed()
if logextra:
other = np.log(other)
return coords, time, other
else:
return coords, time
def rms(x):
return np.sqrt(sum(x**2))
def rotation_velocity_at_r(r, params):
"""
Gets the rotation velocity at a given radius r by assuming it is keplerian.
Assumes we are in cgs units, which may one day be our downfall.
"""
unit_length = float(params[r"InternalUnitSystem:UnitLength_in_cgs"])
if unit_length != 1.:
print(f"Your unit length: {unit_length}")
raise InternalUnitSystemError(
"This function is only created to handle CGS units."
)
central_mass = float(params["PointMassPotential:mass"])
G = 6.674e-8
v = np.sqrt( G * central_mass / r)
return v
def get_rotation_period_at_r(r, params):
"""
Gets the rotation period at a given radius r, assuming a keplerian
orbit.
"""
v = rotation_velocity_at_r(r, params)
return 2*np.pi / v
def get_metadata(filename, r=1):
""" The metadata should be extracted from the first snapshot. """
with h5.File(filename, "r") as file_handle:
header = file_handle['Header'].attrs
code = file_handle['Code'].attrs
hydro = file_handle['HydroScheme'].attrs
params = file_handle['Parameters'].attrs
period = get_rotation_period_at_r(r, params)
return_values = {
"header" : dict(header),
"code" : dict(code),
"period" : float(period),
"hydro" : dict(hydro),
"params" : dict(params)
}
return return_values
def plot_single(number, scatter, text, metadata, ax, extra=None, norm=None):
filename = "keplerian_ring_{:04d}.hdf5".format(number)
if extra is not None:
coordinates, time, other = load_data(filename, extra=extra)
else:
coordinates, time = load_data(filename)
text.set_text(
"Time: {:1.2f} | Rotations {:1.2f}".format(
time,
time/metadata['period'],
)
)
data = coordinates[:, 0:2]
scatter.set_offsets(data)
if extra is not None:
colours, _ = get_colours(other, norm)
scatter.set_color(colours)
return scatter,
if __name__ == "__main__":
import os
import argparse as ap
parser = ap.ArgumentParser(
description="""
Plots a movie of the current snapshots in your
directory. Can also colourmap your information.
"""
)
parser.add_argument(
"-e",
"--exclude_central",
help="""
Region from the centre of the ring to exclude from the
vmin/vmax of the colourmap. Note that these particles are
still plotted -- they are just excluded for the purposes
of colourmapping. Default: 1 simulation unit.
""",
default=1.,
required=False
)
parser.add_argument(
"-c",
"--cmap",
help="""
The item from the GADGET hdf5 file to clourmap with.
Examples include Density, InternalEnergy.
Default: don't use a colourmap. (Much faster).
""",
required=False,
default=None
)
parser.add_argument(
"-m",
"--max",
help="""
Maximum radii to plot.
Default: 2.5.
""",
required=False,
default=2.5,
)
args = vars(parser.parse_args())
# Look for the number of files in the directory.
i = 0
while True:
if os.path.isfile("keplerian_ring_{:04d}.hdf5".format(i)):
i += 1
else:
break
if i > 10000:
break
# Now deal with the colourmapping (if present)
if args["cmap"] is not None:
_, _, numbers0 = load_data("keplerian_ring_0000.hdf5", extra=args["cmap"], exclude=float(args["exclude_central"]))
_, _, numbersend = load_data("keplerian_ring_{:04d}.hdf5".format(i-1), extra=args["cmap"], exclude=float(args["exclude_central"]))
vmax = max([numbers0.max(), numbersend.max()])
vmin = min([numbers0.min(), numbersend.min()])
norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
else:
norm = None
# Initial plot setup
metadata = get_metadata("keplerian_ring_0000.hdf5")
n_particle = metadata['header']['NumPart_Total'][0]
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111)
scatter = ax.scatter([0]*n_particle, [0]*n_particle, s=0.5, marker="o")
diff = float(args["max"])
left = metadata['header']['BoxSize'][0]/2 - diff
right = metadata['header']['BoxSize'][0]/2 + diff
ax.set_xlim(left, right)
ax.set_ylim(left, right)
offset = 0.25
time_text = ax.text(
offset + left,
offset + left,
"Time: {:1.2f} | Rotations {:1.2f}".format(
0,
0/metadata['period'],
)
)
ax.text(
offset + left,
right-offset-0.35,
"Code: {} {} | {} {} \nHydro {}\n$\eta$={:1.4f}".format(
metadata['code']['Git Branch'].decode("utf-8"),
metadata['code']['Git Revision'].decode("utf-8"),
metadata['code']['Compiler Name'].decode("utf-8"),
metadata['code']['Compiler Version'].decode("utf-8"),
metadata['hydro']['Scheme'].decode("utf-8"),
metadata['hydro']['Kernel eta'][0],
)
)
ax.set_title("Keplerian Ring Test")
ax.set_xlabel("$x$ position")
ax.set_ylabel("$y$ position")
anim = anim.FuncAnimation(
fig,
plot_single,
tqdm(np.arange(i)),
fargs = [
scatter,
time_text,
metadata,
ax,
args["cmap"],
norm
],
interval=50,
repeat_delay=3000,
blit=True,
)
anim.save("keplerian_ring.mp4", dpi=int(640/8))
examples/KeplerianRing/plotSolution.py 0 → 100644
View file @ 2aaae570
"""
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2017
#
# 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/>.
#
###############################################################################
"""
# Plotting script for the Keplerian Ring example.
# We use yt for the projection rendering of the ring,
# and then our own density as a function of radius calculation.
import matplotlib
matplotlib.use("Agg")
matplotlib.rc("text", usetex=True)
try:
import yt
ytavail = True
except ImportError:
print("yt not found. Falling back on homebrew plots.")
ytavail = False
import h5py
import os
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.gridspec as gridspec
try:
from tqdm import tqdm
except ImportError:
tqdm = lambda x: x
from make_movie import get_metadata
yt.funcs.mylog.setLevel(50)
tqdm.monitor_interval = 0
class InternalUnitSystemError(Exception):
pass
def get_axes_grid(figure):
"""
Grab our axes grid.
"""
gs = gridspec.GridSpec(2, 30)
grid = []
grid.append(figure.add_subplot(gs[0, 0:10]))
grid.append(figure.add_subplot(gs[0, 10:20]))
grid.append(figure.add_subplot(gs[0, 20:]))
grid.append(figure.add_subplot(gs[1, 0:9]))
grid.append(figure.add_subplot(gs[1, 11:20]))
grid.append(figure.add_subplot(gs[1, 21:]))
return grid
def get_yt_actual_data(plot, name="density"):
"""
Extracts the image data and colourmap from a yt plot.
This is used to put on our own grid.
"""
data = plot.plots[name].image.get_array()
cmap = plot.plots[name].image.cmap
return data, cmap