diff --git a/.gitignore b/.gitignore
index f2726b05eb3bb1d782333eb5bf4bd98fde02ff27..434d88e2771adfceb347838560812ce418d3a22e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -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
diff --git a/configure.ac b/configure.ac
index edda4e608d97227633dc55c726a429ffb67d8aad..c4e9f9d971154ea3889cc3fbb696dc20e5d02fc3 100644
--- a/configure.ac
+++ b/configure.ac
@@ -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])
;;
diff --git a/examples/KeplerianRing/README.md b/examples/KeplerianRing/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..d486cbfe412c7ff9ca121c87bbc548d0ece078dd
--- /dev/null
+++ b/examples/KeplerianRing/README.md
@@ -0,0 +1,130 @@
+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.
diff --git a/examples/KeplerianRing/keplerian_ring.yml b/examples/KeplerianRing/keplerian_ring.yml
new file mode 100644
index 0000000000000000000000000000000000000000..421d1e89255195cd05a66975d838c59a4ad77c72
--- /dev/null
+++ b/examples/KeplerianRing/keplerian_ring.yml
@@ -0,0 +1,46 @@
+# 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
diff --git a/examples/KeplerianRing/makeIC.py b/examples/KeplerianRing/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..e99b7f3b5176f0c7c857e7407634caa78ed6c431
--- /dev/null
+++ b/examples/KeplerianRing/makeIC.py
@@ -0,0 +1,740 @@
+"""
+###############################################################################
+# 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/>.
+#
+###############################################################################
+"""
+
+
+import numpy as np
+import write_gadget as wg
+import h5py as h5
+
+from scipy.special import erfinv
+
+
+class Particles(object):
+ """
+ Holder class for particle properties. Also contains some methods to e.g.
+ set their keplerian velocities based on their positions. These properties
+ are set using the 'generationmethod' functions below.
+ """
+ def __init__(self, meta):
+ self.gravitymass = meta["gravitymass"]
+ self.nparts = meta["nparts"]**2
+ self.particlemass = meta["particlemass"]
+ self.softening = meta["softening"]
+ self.boxsize = meta["boxsize"]
+
+ self.smoothing = np.zeros(self.nparts) + meta["smoothing"]
+ self.internalenergy = np.zeros(self.nparts) + meta["internalenergy"]
+
+ self.positions = np.array([])
+ self.radii = np.array([])
+ self.theta = np.array([])
+ self.phi = np.array([])
+ self.velocities = np.array([])
+ self.ids = np.array([])
+ self.densities = np.array([])
+ self.masses = np.array([])
+
+ return
+
+
+ def calculate_masses(self):
+ """
+ Calculate the individual masses for the particles such that we have
+ a uniform thickness disk, given their densities.
+ """
+ mass_factor = self.particlemass / self.densities.max()
+ self.masses = self.densities * mass_factor
+
+ return self.masses
+
+
+ def calculate_velocities(self, angle=0):
+ """
+ Calculates keplerian orbit velocities for the, well, keplerian ring.
+ Requires that radius and phi are set already.
+
+ Please give angle in radians.
+ """
+ force_modifier = np.sqrt(self.gravitymass / (self.radii**2 + self.softening**2)**(3/2)) * self.radii
+ try:
+ v_x = np.sin(self.theta) * np.sin(self.phi) * np.cos(angle)
+ except ValueError:
+ # Phi are not yet set. Make them and then move on to v_x
+ if angle:
+ raise ValueError("Unable to find phi.")
+
+ # If we have angle of inclination 0 we can set phi.
+ self.phi = np.zeros_like(self.theta) + np.pi / 2
+ v_x = np.sin(self.theta) * np.sin(self.phi) * np.cos(angle)
+
+ v_y = np.cos(self.phi) * np.sin(angle) - np.cos(self.theta) * np.sin(self.phi) * np.cos(angle)
+ v_z = np.sin(self.theta) * np.sin(self.phi) * np.sin(angle)
+
+ self.velocities = (force_modifier * np.array([v_x, v_y, v_z])).T
+
+ return self.velocities
+
+
+ def generate_ids(self):
+ """
+ Generate consecutive IDs from 0 based on the number of particles
+ currently in the object.
+ """
+ self.ids = np.arange(self.nparts)
+
+ return self.ids
+
+
+ def convert_polar_to_cartesian(self, centre_of_ring=(5, 5), boxsize=None):
+ """
+ Converts self.radii, self.theta to self.positions.
+
+ If self.phi is set, it also uses that to conver to z.
+ """
+ if len(self.phi) == 0:
+ x = self.radii * np.cos(self.theta) + centre_of_ring[0]
+ y = self.radii * np.sin(self.theta) + centre_of_ring[1]
+ z = np.zeros_like(x) + boxsize/2
+ else:
+ x = self.radii * np.cos(self.theta) * np.sin(self.phi) + centre_of_ring[0]
+ y = self.radii * np.sin(self.theta) * np.sin(self.phi) + centre_of_ring[1]
+ try:
+ z = self.radii * np.cos(self.phi) + centre_of_ring[2]
+ except AttributeError:
+ # Just set the central z value to the middle of the box
+ z = self.radii * np.cos(self.phi) + boxsize/2
+
+ self.positions = np.array([x, y, z]).T
+
+ return self.positions
+
+
+ def convert_cartesian_to_polar(self, centre_of_ring=(5, 5)):
+ """
+ Converts self.positions to self.radii, self.theta, and self.phi.
+ """
+ x, y, z = self.positions.T
+
+ try:
+ x = x - centre_of_ring[0]
+ y = y - centre_of_ring[1]
+ z = z - centre_of_ring[2]
+ except AttributeError:
+ raise AttributeError("Invalid array for centre_of_ring provided.")
+
+ xsquare = x*x
+ ysquare = y*y
+ zsquare = z*z
+
+ r = np.sqrt(xsquare + ysquare + zsquare)
+
+ # We need to mask over the x = 0 cases (theta = 0).
+ mask = np.isclose(x, 0, 1e-6)
+
+ masked_x = np.ma.array(x, mask=mask)
+ theta_for_unmasked = np.arctan2(y, masked_x)
+
+ theta = theta_for_unmasked + mask * 0
+
+ # Thankfully we have already ensured that r != 0.
+ phi = np.arccos(z/r)
+
+
+ self.radii = r
+ self.theta = theta
+ self.phi = phi
+
+ return r, theta, phi
+
+
+ def wiggle_positions(self, tol=1e-3):
+ """
+ 'Wiggle' the positions to avoid precision issues.
+
+ Note that this does not touch r, theta, phi.
+ """
+ self.positions += np.random.random(self.positions.shape) * tol
+
+ return
+
+
+ def exclude_particles(self, range):
+ """
+ Exclude all particles that are *not* within range (of radii).
+ """
+ mask = np.logical_or(
+ self.radii < range[0],
+ self.radii > range[1],
+ )
+
+ x = np.ma.array(self.positions[:, 0], mask=mask).compressed()
+ y = np.ma.array(self.positions[:, 1], mask=mask).compressed()
+ z = np.ma.array(self.positions[:, 2], mask=mask).compressed()
+
+ self.positions = np.array([x, y, z]).T
+
+ try:
+ v_x = np.ma.array(self.velocities[:, 0], mask=mask).compressed()
+ v_y = np.ma.array(self.velocities[:, 1], mask=mask).compressed()
+ v_z = np.ma.array(self.velocities[:, 2], mask=mask).compressed()
+
+ self.velocities = np.array([v_x, v_y, v_z])
+ except IndexError:
+ # We haven't filled them yet anyway...
+ pass
+
+ try:
+ self.ids = np.ma.array(self.ids, mask=mask).compressed()
+ except np.ma.core.MaskError:
+ # Not filled yet.
+ pass
+
+ try:
+ self.densities = np.ma.array(self.densities, mask=mask).compressed()
+ except np.ma.core.MaskError:
+ # Not filled yet.
+ pass
+
+ # QSP Fix has modified us, so first we need to chop off extras.
+ # Then, as they are all the same, we don't need to remove 'specific'
+ # values, and so we can just chop off the ends.
+ self.smoothing = self.smoothing[:len(x)]
+ self.internalenergy = self.internalenergy[:len(x)]
+
+ try:
+ self.masses = np.ma.array(self.masses, mask=mask).compressed()
+ except np.ma.core.MaskError:
+ # Not filled yet.
+ pass
+
+ self.radii = np.ma.array(self.radii, mask=mask).compressed()
+ self.theta = np.ma.array(self.theta, mask=mask).compressed()
+
+ try:
+ self.phi = np.ma.array(self.phi, mask=mask).compressed()
+ except np.ma.core.MaskError:
+ # We have not allocated our phis,
+ pass
+
+ self.nparts = len(self.radii)
+
+ return
+
+
+ def tilt_particles(self, angle, center=(5, 5, 5)):
+ """
+ Tilts the particles around the x-axis around the point given.
+
+ Assumes that the particles already have set r, theta, phi.
+ """
+ angle_radians = angle * np.pi / 180
+ rotation_matrix = np.array(
+ [
+ [ np.cos(angle_radians), 0, np.sin(angle_radians) ],
+ [ 0, 1, 0 ],
+ [ -np.sin(angle_radians), 0, np.cos(angle_radians) ]
+ ]
+ )
+
+ self.positions = (
+ (
+ np.matmul(
+ rotation_matrix,
+ (self.positions - np.array(center)).T
+ )
+ ).T
+ ) + np.array(center)
+
+ self.velocities = (
+ (
+ np.matmul(
+ rotation_matrix,
+ (self.velocities).T
+ )
+ ).T
+ )
+
+ self.convert_cartesian_to_polar(center)
+
+ return
+
+
+ def save_to_gadget(self, filename, boxsize=10.):
+ """
+ Save the particle data to a GADGET .hdf5 file.
+
+ Uses the internal options, but you must specify a filename.
+ """
+ with h5.File(filename, "w") as handle:
+ wg.write_header(
+ handle,
+ boxsize=boxsize,
+ flag_entropy=0,
+ np_total=np.array([self.nparts, 0, 0, 0, 0, 0]),
+ np_total_hw=np.array([0, 0, 0, 0, 0, 0]),
+ other={
+ "MassTable" : np.array([self.particlemass, 0, 0, 0, 0, 0]),
+ "Time" : 0,
+ }
+ )
+
+ wg.write_runtime_pars(
+ handle,
+ periodic_boundary=1,
+ )
+
+ wg.write_units(
+ handle,
+ current=1.,
+ length=1.,
+ mass=1,
+ temperature=1.,
+ time=1.,
+ )
+
+ wg.write_block(
+ handle,
+ 0, # gas
+ self.positions,
+ self.velocities,
+ self.ids,
+ mass=self.masses,
+ int_energy=self.internalenergy,
+ smoothing=self.smoothing,
+ other={"Density": self.densities},
+ )
+
+ return
+
+
+def __sigma(r):
+ """
+ Density distribution of the ring, this comes directly from Hopkins 2015.
+ """
+ if r < 0.5:
+ return 0.01 + (r/0.5)**3
+ elif r <= 2:
+ return 1.01
+ else:
+ return 0.01 + (1 + (r-2)/0.1)**(-3)
+
+
+# This is required because of the if, else statement above that does not
+# play nicely with our numpy arrays.
+sigma = np.vectorize(__sigma)
+
+
+def generate_theta(r, theta_initial=0.):
+ """
+ Generate the theta associated with the particles based on their radii.
+ This uses the method from The effect of Poisson noise on SPH calculations,
+ Annabel Cartwright, Dimitris Stamatellos and Anthony P. Whitworth 2009.
+
+ @param: r | float / array-like
+ - the radii of the particles.
+
+ @param: theta_initial | float | optional
+ - initial radius to start the generation at.
+
+ ---------------------------------------------------------------------------
+
+ @return: theta_i | numpy.array
+ - angles associated with the particles in the plane.
+ """
+ radii_fraction = r[:-1] / r[1:]
+ d_theta_i = np.sqrt(2 * np.pi * (1 - radii_fraction))
+
+ theta_i = np.empty_like(r)
+ theta_i[0] = theta_initial
+
+ # Need to do this awful loop because each entry relies on the one before it.
+ # Unless someone knows how to do this in numpy.
+ for i in range(len(d_theta_i)): # first is theta_initial
+ theta_i[i+1] = theta_i[i] + d_theta_i[i]
+
+ return theta_i
+
+
+def QSP_fix(r_i, theta_i):
+ """
+ The start and end of the disk will have the end of the spiral there. That
+ is no good and introduces some shear forces, so we need to move them to
+ concentric circles. We'll also add some extra particles on the final
+ 'layer' of this giant onion to ensure that all of the circles are complete.
+
+ @param: r_i | numpy.array
+ - the original r_i generated from inverse_gaussian (and perhaps
+ afterwards masked).
+
+ @param: theta_i | numpy.array
+ - the original theta_i generated from generate_theta_i.
+
+ ---------------------------------------------------------------------------
+
+ @return r_i_fixed | numpy.array
+ - the fixed, concentric circle-like r_i. Note that these arrays do not
+ necessarily have the same length as the r_i, theta_i that are
+ passed in and you will have to re-calculate the number of particles
+ in the system.
+
+ @return theta_i_fixed | numpy.array
+ - the fixed, concentric circle-like theta_i
+ """
+
+ # Thankfully, the delta_thetas are not wrapped (i.e. they keep on going
+ # from 2 pi -- we need to split the arrays into 'circles' by splitting
+ # the theta_i array every 2 pi.
+
+ rotations = 1
+ circles = []
+
+ these_r_i = []
+ these_theta_i = []
+
+ for radius, theta in zip(r_i, theta_i):
+ if theta > rotations * 2 * np.pi:
+ circles.append([these_r_i, these_theta_i])
+ these_r_i = []
+ these_theta_i = []
+ rotations += 1
+
+ these_r_i.append(radius)
+ these_theta_i.append(theta)
+
+ # Now we need to do the averaging.
+ # We want to have all particles in each circle a fixed radius away from the
+ # centre, as well as having even spacing between each particle. The final
+ # ring may be a bit dodgy still, but we will see.
+
+ r_i_fixed = []
+ theta_i_fixed = []
+
+ for circle in circles:
+ n_particles = len(circle[0])
+ radius = sum(circle[0]) / n_particles
+ theta_initial = circle[1][0]
+
+ theta_sep = 2 * np.pi / n_particles
+
+ theta = [t * theta_sep for t in range(n_particles)]
+ radii = [radius] * n_particles
+
+ r_i_fixed += radii
+ theta_i_fixed += theta
+
+ return np.array(r_i_fixed), np.array(theta_i_fixed)
+
+
+def gen_particles_grid(meta):
+ """
+ Generates particles on a grid and returns a filled Particles object.
+ """
+ particles = Particles(meta)
+ range = (0, meta["boxsize"])
+ centre_of_ring = [meta["boxsize"]/2.] * 3
+
+ # Because we are using a uniform grid we actually use the same x and y
+ # range for the initial particle setup.
+ step = (range[1] - range[0])/meta["nparts"]
+
+ x_values = np.arange(0, range[1] - range[0], step)
+
+ # These are 2d arrays which isn't actually that helpful.
+ x, y = np.meshgrid(x_values, x_values)
+ x = x.flatten() + centre_of_ring[0] - (range[1] - range[0])/2
+ y = y.flatten() + centre_of_ring[1] - (range[1] - range[0])/2
+ z = np.zeros_like(x) + meta["boxsize"]/2
+
+ particles.positions = np.array([x, y, z]).T
+ particles.radii = np.sqrt((x - centre_of_ring[0])**2 + (y - centre_of_ring[1])**2)
+ particles.theta = np.arctan2(y - centre_of_ring[1], x - centre_of_ring[0])
+ particles.exclude_particles((particles.softening, 100.))
+
+ particles.densities = sigma(particles.radii)
+ particles.calculate_velocities()
+ particles.calculate_masses()
+
+ particles.generate_ids()
+
+ if meta["angle"]:
+ particles.tilt_particles(meta["angle"], centre_of_ring)
+
+ return particles
+
+
+def gen_particles_spiral(meta):
+ """
+ Generates particles on concentric circles and returns a filled Particles
+ object. Based on Cartwright, Stamatellos & Whitworth (2009).
+ """
+ particles = Particles(meta)
+ centre_of_ring = (meta["boxsize"]/2., meta["boxsize"]/2., meta["boxsize"]/2.)
+ max_r = meta["boxsize"]/2.
+
+ m = (np.arange(particles.nparts) + 0.5)/particles.nparts
+ r = max_r * m
+ theta = generate_theta(r)
+
+ particles.radii, particles.theta = QSP_fix(r, theta)
+ # We must do this afterwards as QSP does not always return the same number of parts.
+ phi = np.zeros_like(particles.theta) + np.pi/2
+ particles.phi = phi
+
+ particles.convert_polar_to_cartesian(centre_of_ring, meta["boxsize"])
+ particles.nparts = len(particles.radii)
+
+ particles.exclude_particles((particles.softening, 100.))
+
+ # This way of doing densities does not work for different sized patches.
+ # Therefore we need to weight by effecitve area.
+ dtheta = (particles.theta[1:] - particles.theta[:-1]) / 2
+ areas = np.square(4 * particles.radii[1:] * np.sin(dtheta))
+ areas = np.hstack((np.array([1]), areas))
+
+ # Now do the 'actual' density calculation.
+ particles.densities = areas * sigma(particles.radii)
+ particles.calculate_velocities()
+ particles.calculate_masses()
+
+ particles.generate_ids()
+
+ if meta["angle"]:
+ particles.tilt_particles(meta["angle"], centre_of_ring)
+
+
+ return particles
+
+
+def gen_particles_gaussian(meta):
+ """
+ Generates particles on concentric circles and returns a filled Particles
+ object. Based on Cartwright, Stamatellos & Whitworth (2009).
+
+ This generation function uses a Gaussian PDF.
+ """
+ particles = Particles(meta)
+ centre_of_ring = (meta["boxsize"]/2., meta["boxsize"]/2., meta["boxsize"]/2.)
+ max_r = meta["boxsize"]/2.
+
+ m = (np.arange(particles.nparts) + 0.5)/particles.nparts
+ error_function = erfinv(2*m - 1)
+ r = 2 + 0.5 * error_function
+
+ theta = generate_theta(r)
+
+ particles.radii, particles.theta = QSP_fix(r, theta)
+ # We must do this afterwards as QSP does not always return the same number of parts.
+ phi = np.zeros_like(particles.theta) + np.pi/2
+ particles.phi = phi
+
+ particles.convert_polar_to_cartesian(centre_of_ring, meta["boxsize"])
+ particles.nparts = len(particles.radii)
+
+ particles.exclude_particles((particles.softening, 100.))
+
+ # This way of doing densities does not work for different sized patches.
+ particles.densities = np.zeros_like(particles.radii)
+ particles.calculate_velocities()
+ particles.masses = np.ones_like(particles.radii)
+
+ particles.generate_ids()
+
+ if meta["angle"]:
+ particles.tilt_particles(meta["angle"], centre_of_ring)
+
+
+ return particles
+
+
+if __name__ == "__main__":
+ import argparse as ap
+
+ PARSER = ap.ArgumentParser(
+ description="""
+ Initial conditions generator for the Keplerian Ring
+ example. It has sensible defaults for GIZMO, but if you
+ wish to run the example with SPH you sould use
+ --generationmethod spiral.
+ """
+ )
+
+ PARSER.add_argument(
+ "-m",
+ "--gravitymass",
+ help="""
+ GM for the central point mass. Default: 1.
+ """,
+ required=False,
+ default=1.,
+ )
+
+ PARSER.add_argument(
+ "-f",
+ "--filename",
+ help="""
+ Filename for your initial conditions.
+ Default: initial_conditions.hdf5.
+ """,
+ required=False,
+ default="initial_conditions.hdf5",
+ )
+
+ PARSER.add_argument(
+ "-n",
+ "--nparts",
+ help="""
+ Square-root of the number of particles, i.e. the default
+ nparts=128 leads to a ring with 128^2 particles in it.
+ Note that for the spiral generation method this number is
+ somewhat approximate.
+ """,
+ required=False,
+ default=128,
+ )
+
+ PARSER.add_argument(
+ "-p",
+ "--particlemass",
+ help="""
+ Maximum mass of the gas particles. Default: 10.
+ """,
+ required=False,
+ default=10.,
+ )
+
+ PARSER.add_argument(
+ "-e",
+ "--softening",
+ help="""
+ Gravitational softening for the centre of the ring. We also
+ exclude all particles within this radius from the centre of the
+ ring. Default: 0.05.
+ """,
+ required=False,
+ default=0.05,
+ )
+
+ PARSER.add_argument(
+ "-s",
+ "--smoothing",
+ help="""
+ Initial smoothing length for all of the particles.
+ Default: Boxsize/N
+ """,
+ required=False,
+ default=-1,
+ )
+
+ PARSER.add_argument(
+ "-i",
+ "--internalenergy",
+ help="""
+ Initial internal energy for all of the particles. Note that this
+ should be low to ensure that the ring is very cold (see Inviscid
+ SPH for details). Default: 0.015.
+ """,
+ required=False,
+ default=0.015
+ )
+
+ PARSER.add_argument(
+ "-g",
+ "--generationmethod",
+ help="""
+ Generation method for the particles. Choose between grid, spiral,
+ and gaussian, where spiral ensures that the particles are generated
+ in a way that minimises the energy in SPH. For more details on
+ this method see Cartwright, Stamatellos & Whitworth (2009).
+ Default: grid.
+ """,
+ required=False,
+ default="grid",
+ )
+
+ PARSER.add_argument(
+ "-b",
+ "--boxsize",
+ help="""
+ The box size. Default: 10.
+ """,
+ required=False,
+ default=10.
+ )
+
+ PARSER.add_argument(
+ "-a",
+ "--angle",
+ help="""
+ Angle of incline on the disk (degrees). Default: 0.
+ Note that tilting the ring may cause some issues at the
+ disk boundary as the 'box' is no longer filled with
+ particles.
+ """,
+ required=False,
+ default=0.
+ )
+
+
+ ### --- ### --- Argument Parsing --- ### --- ###
+
+ ARGS = vars(PARSER.parse_args())
+
+ if ARGS["generationmethod"] == "grid":
+ gen_particles = gen_particles_grid
+ elif ARGS["generationmethod"] == "spiral":
+ gen_particles = gen_particles_spiral
+ elif ARGS["generationmethod"] == "gaussian":
+ gen_particles = gen_particles_gaussian
+ else:
+ print(
+ "ERROR: {} is an invalid generation method. Exiting.".format(
+ ARGS["generationmethod"]
+ )
+ )
+ exit(1)
+
+
+ if ARGS["smoothing"] == -1:
+ smoothing = float(ARGS["boxsize"]) / int(ARGS["nparts"])
+ else:
+ smoothing = float(ARGS["smoothing"])
+
+
+ META = {
+ "gravitymass": float(ARGS["gravitymass"]),
+ "nparts": int(ARGS["nparts"]),
+ "particlemass": float(ARGS["particlemass"]),
+ "smoothing": smoothing,
+ "softening": float(ARGS["softening"]),
+ "internalenergy": float(ARGS["internalenergy"]),
+ "boxsize": float(ARGS["boxsize"]),
+ "angle" : float(ARGS["angle"])
+ }
+
+ PARTICLES = gen_particles(META)
+
+ PARTICLES.save_to_gadget(
+ filename=ARGS["filename"],
+ boxsize=ARGS["boxsize"],
+ )
+
diff --git a/examples/KeplerianRing/make_movie.py b/examples/KeplerianRing/make_movie.py
new file mode 100644
index 0000000000000000000000000000000000000000..83e783924086377a0b2aa384d2c4634bfd40443f
--- /dev/null
+++ b/examples/KeplerianRing/make_movie.py
@@ -0,0 +1,310 @@
+"""
+###############################################################################
+# 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))
diff --git a/examples/KeplerianRing/plotSolution.py b/examples/KeplerianRing/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..4d6411106f7f5ddd047148927a03fbb4800e0874
--- /dev/null
+++ b/examples/KeplerianRing/plotSolution.py
@@ -0,0 +1,627 @@
+"""
+###############################################################################
+# 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
+
+
+def chi_square(observed, expected):
+ """
+ The chi squared statistic.
+
+ This also looks for where expected == 0 and masks over those particles to
+ avoid divide by zero errors and unrealistically high chi squared.
+ """
+
+ mask = np.array(expected) != 0
+
+ masked_expected = np.array(expected)[mask]
+ masked_observed = np.array(observed)[mask]
+
+ return sum(((masked_observed - masked_expected)**2)/masked_expected**2)
+
+
+def load_data(filename):
+ """
+ Loads the data and extracts the relevant information for
+ calculating the chi squared statistic and density(r) profiles.
+ """
+
+ with h5py.File(filename, "r") as file:
+ boxsize = np.array(file["Header"].attrs["BoxSize"])
+
+ # Check if z = 0 for all particles. If so we need to set the cetnre
+ # to have z = 0 also.
+ if np.sum(file["PartType0"]["Coordinates"][:, 2]) == 0:
+ centre = [boxsize[0] / 2., boxsize[0] / 2., 0.]
+ else:
+ centre = boxsize / 2.
+
+ radii = np.sqrt(np.sum(((file["PartType0"]["Coordinates"][...] - centre).T)**2, 0))
+ masses = file["PartType0"]["Masses"][...]
+
+ return radii, masses
+
+
+def bin_density_r(radii, density, binrange, binnumber):
+ """
+ Bins the density as a funciton of radius.
+ """
+
+ bins = np.linspace(*binrange, binnumber)
+ indicies = np.digitize(radii, bins)
+
+ binned_masses = np.zeros(len(bins) - 1)
+
+ for index, bin in enumerate(indicies):
+ if bin >= len(bins) - 1:
+ continue
+
+ binned_masses[bin] += density[index]
+
+ areas = [np.pi * (a**2 - b**2) for a, b in zip(bins[1:], bins[:-1])]
+ binned_densities = binned_masses/areas
+
+ return bins, binned_densities
+
+
+def get_density_r(snapshot, filename="keplerian_ring", binrange=(0, 5), binnumber=50):
+ """
+ Gets the binned density as a function of radius.
+ """
+ snap = "{:04d}".format(snapshot)
+ filename = f"{filename}_{snap}.hdf5"
+
+ data = load_data(filename)
+
+ return bin_density_r(*data, binrange, binnumber)
+
+
+def get_mass_outside_inside(snap, radius=2, filename="keplerian_ring"):
+ """
+ Finds the mass outside and inside a radius. This is used to look at mass
+ flow inside and outside of the ring with a postprocessing routine in
+ get_derived_data.
+ """
+ snapshot = "{:04d}".format(snap)
+ filename = f"{filename}_{snapshot}.hdf5"
+
+ radii, masses = load_data(filename)
+
+ below = sum(masses[radii < radius])
+ above = sum(masses[radii > radius])
+
+ return below, above
+
+
+def get_derived_data(minsnap, maxsnap, filename="keplerian_ring", binnumber=50, radius=2):
+ """
+ Gets the derived data from our snapshots, i.e. the
+ density(r) profile and the chi squared (based on the
+ difference between the minsnap and the current snapshot).
+ """
+
+ initial = get_density_r(minsnap, filename, binnumber=binnumber)
+ other_densities = [
+ get_density_r(snap, binnumber=binnumber)[1] for snap in tqdm(
+ range(minsnap+1, maxsnap+1), desc="Densities"
+ )
+ ]
+ densities = [initial[1]] + other_densities
+
+ masses_inside_outside = [
+ get_mass_outside_inside(snap, radius=radius, filename=filename)[0] for snap in tqdm(
+ range(minsnap, maxsnap+1), desc="Mass Flow"
+ )
+ ]
+
+ # Between the initial conditions and the first snapshot we hope that there
+ # has been no mass flow, hence the [0.] +
+ mass_flows = [0.] + [
+ y - x for x, y in zip(
+ masses_inside_outside[1:],
+ masses_inside_outside[:-1]
+ )
+ ]
+
+ cumulative_mass_flows = [
+ sum(mass_flows[:x])/masses_inside_outside[0] for x in range(len(mass_flows))
+ ]
+
+ chisq = [
+ chi_square(
+ dens[int(0.1*binnumber):int(0.4*binnumber)],
+ initial[1][int(0.1*binnumber):int(0.4*binnumber)]
+ ) for dens in tqdm(densities, desc="Chi Squared")
+ ]
+
+ return initial[0], densities, chisq, cumulative_mass_flows
+
+
+def plot_chisq(ax, minsnap, maxsnap, chisq, filename="keplerian_ring"):
+ """
+ Plot the chisq(rotation).
+ """
+ snapshots = np.arange(minsnap, maxsnap + 1)
+ rotations = [convert_snapshot_number_to_rotations_at(1, snap, filename) for snap in snapshots]
+ ax.plot(rotations, np.array(chisq)/max(chisq))
+
+ ax.set_xlabel("Number of rotations")
+ ax.set_ylabel("$\chi^2 / \chi^2_{{max}}$ = {:3.5f}".format(max(chisq)))
+
+ return
+
+
+def plot_mass_flow(ax, minsnap, maxsnap, mass_flow, filename="keplerian_ring"):
+ """
+ Plot the mass_flow(rotation).
+ """
+ snapshots = np.arange(minsnap, maxsnap + 1)
+ rotations = [convert_snapshot_number_to_rotations_at(1, snap, filename) for snap in snapshots]
+
+ ax.plot(rotations, mass_flow)
+
+ ax.set_xlabel("Number of rotations")
+ ax.set_ylabel(r"Mass flow out of ring ($M_{\rm ring}$)")
+
+ return
+
+
+def plot_density_r(ax, bins, densities, snaplist, filename="keplerian_ring"):
+ """
+ Make the density(r) plots.
+
+ Densities is the _full_ list of density profiles, and
+ snaplist is the ones that you wish to plot.
+ """
+ radii = [(x + y)/2 for x, y in zip(bins[1:], bins[:-1])]
+
+ for snap in snaplist:
+ index = snap - snaplist[0]
+ rotations = convert_snapshot_number_to_rotations_at(1, snap, filename)
+ ax.plot(radii, densities[index], label="{:2.2f} Rotations".format(rotations))
+
+ ax.legend()
+ ax.set_xlabel("Radius")
+ ax.set_ylabel("Azimuthally Averaged Surface Density")
+
+ return
+
+
+def plot_extra_info(ax, filename):
+ """
+ Plots all of the extra information on the final axis.
+
+ Give it a filename of any of the snapshots.
+ """
+
+ metadata = get_metadata(filename)
+
+ git = metadata['code']['Git Revision'].decode("utf-8")
+ compiler_name = metadata['code']['Compiler Name'].decode("utf-8")
+ compiler_version = metadata['code']['Compiler Version'].decode("utf-8")
+ scheme = metadata['hydro']['Scheme'].decode("utf-8")
+ kernel = metadata['hydro']['Kernel function'].decode("utf-8")
+ gas_gamma = metadata["hydro"]["Adiabatic index"][0]
+ neighbors = metadata["hydro"]["Kernel target N_ngb"][0]
+ eta = metadata["hydro"]["Kernel eta"][0]
+
+
+ ax.text(-0.49, 0.9, "Keplerian Ring with $\\gamma={:4.4f}$ in 2/3D".format(gas_gamma), fontsize=11)
+ ax.text(-0.49, 0.8, f"Compiler: {compiler_name} {compiler_version}", fontsize=10)
+ ax.text(-0.49, 0.7, "Rotations are quoted at $r=1$", fontsize=10)
+ ax.plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+ ax.text(-0.49, 0.5, f"$\\textsc{{Swift}}$ {git}", fontsize=10)
+ ax.text(-0.49, 0.4, scheme, fontsize=10)
+ ax.text(-0.49, 0.3, kernel, fontsize=10)
+ ax.text(-0.49, 0.2, "${:2.2f}$ neighbours ($\\eta={:3.3f}$)".format(neighbors, eta), fontsize=10)
+
+ ax.set_axis_off()
+ ax.set_xlim(-0.5, 0.5)
+ ax.set_ylim(0, 1)
+
+ return
+
+
+def surface_density_plot_no_yt(ax, snapnum, filename="keplerian_ring", density_limits=None, vlim=None):
+ """
+ Make the surface density plot (sans yt).
+
+ Also returns the max and minimum values for the density so these can
+ be passed to the next call, as well as vlim which are the colourmap
+ max/min.
+ """
+
+ with h5py.File("{}_{:04d}.hdf5".format(filename, snapnum)) as filehandle:
+ density = filehandle["PartType0"]["Density"][...]
+ x, y = filehandle["PartType0"]["Coordinates"][:, 0:2].T
+
+ new_vlim = (density.min(), density.max())
+
+ if vlim is None:
+ vlim = new_vlim
+
+ ax.scatter(x, y, c=density, vmin=vlim[0], vmax=vlim[1], s=0.1)
+
+
+ metadata = get_metadata("{}_{:04d}.hdf5".format(filename, snapnum))
+ period = metadata["period"]
+
+ t = ax.text(
+ 2.5,
+ 7.5,
+ "Snapshot = {:04d}\nRotations = {:1.2f}".format(
+ snapnum,
+ float(metadata["header"]["Time"])/period
+ ),
+ color='black'
+ )
+
+ t.set_bbox(dict(alpha=0.5, color="white"))
+
+ ax.axis("equal")
+
+ ax.set_xlim(2, 8)
+ ax.set_ylim(2, 8)
+
+ # We now want to remove all of the ticklabels.
+
+ for axis in ['x', 'y']:
+ ax.tick_params(
+ axis=axis,
+ which='both',
+ bottom='off',
+ top='off',
+ left='off',
+ right='off',
+ labelleft='off',
+ labelbottom='off'
+ )
+
+ return density_limits, vlim
+
+
+
+def surface_density_plot(ax, snapnum, filename="keplerian_ring", density_limits=None, vlim=None):
+ """
+ Make the surface density plot (via yt).
+
+ Also returns the max and minimum values for the density so these can
+ be passed to the next call, as well as vlim which are the colourmap
+ max/min.
+ """
+
+ unit_base = {
+ 'length': (1.0, 'cm'),
+ 'velocity': (1.0, 'cm/s'),
+ 'mass': (1.0, 'g')
+ }
+
+ filename = "{}_{:04d}.hdf5".format(filename, snapnum)
+
+ try:
+ snap = yt.load(filename, unit_base=unit_base, over_refine_factor=2)
+ except yt.utilities.exceptions.YTOutputNotIdentified:
+ # Probably the file isn't here because we supplied a too high snapshot
+ # number. Just return what we're given.
+ return density_limits, vlim
+
+ projection_plot = yt.ProjectionPlot(
+ snap,
+ "z",
+ ("gas", "cell_mass"),
+ width=5.5
+ )
+
+ max_density = snap.all_data()[("gas", "cell_mass")].max()
+ min_density = snap.all_data()[("gas", "cell_mass")].min()
+
+ new_density_limits = (min_density, max_density)
+
+ if density_limits is None:
+ density_limits = new_density_limits
+
+ projection_plot.set_zlim("cell_mass", *density_limits)
+
+ data = get_yt_actual_data(projection_plot, ("gas", "cell_mass"))
+
+ # Becuase of the way plotting works, we also need a max/min for the colourmap.
+
+ new_vlim = (data[0].min(), data[0].max())
+
+ if vlim is None:
+ vlim = new_vlim
+
+ ax.imshow(
+ data[0],
+ cmap=data[1],
+ vmin=vlim[0],
+ vmax=vlim[1]
+ )
+
+ metadata = get_metadata(filename)
+ period = metadata["period"]
+
+ ax.text(
+ 20,
+ 80,
+ "Snapshot = {:04d}\nRotations = {:1.2f}".format(
+ snapnum,
+ float(snap.current_time)/period
+ ),
+ color='white'
+ )
+
+ # We now want to remove all of the ticklabels.
+
+ for axis in ['x', 'y']:
+ ax.tick_params(
+ axis=axis,
+ which='both',
+ bottom='off',
+ top='off',
+ left='off',
+ right='off',
+ labelleft='off',
+ labelbottom='off'
+ )
+
+ return density_limits, vlim
+
+
+def convert_snapshot_number_to_rotations_at(r, snapnum, filename):
+ """
+ Opens the file and extracts metadata to find the number of rotations.
+ """
+
+ metadata = get_metadata("{}_{:04d}.hdf5".format(filename, snapnum))
+
+ t = metadata["period"]
+ current_time = float(metadata["header"]["Time"])
+
+ return current_time / t
+
+
+if __name__ == "__main__":
+ import argparse as ap
+
+ parser = ap.ArgumentParser(
+ description="""
+ Plotting code for the Keplerian Ring test. Uses yt to make
+ surface density plots.
+ """
+ )
+
+ parser.add_argument(
+ "-b",
+ "--beginning",
+ help="""
+ Initial snapshot to start analysis at.
+ Default: First snapshot in the folder.
+ """,
+ default=-1,
+ required=False
+ )
+
+ parser.add_argument(
+ "-m",
+ "--middle",
+ help="""
+ Middle snapshot in the top row of three surface density plots.
+ Default: (end - beginning) // 2.
+ """,
+ default=-1,
+ required=False
+ )
+
+ parser.add_argument(
+ "-e",
+ "--end",
+ help="""
+ Snapshot to end the analysis at.
+ Default: Last snapshot in the folder.
+ """,
+ default=-1,
+ required=False
+ )
+
+ parser.add_argument(
+ "-f",
+ "--filename",
+ help="""
+ Filename of the output.
+ Default: plot.png
+ """,
+ default="plot.png",
+ required=False
+ )
+
+ parser.add_argument(
+ "-n",
+ "--nbins",
+ help="""
+ Number of bins in histogramming (used for the density(r) plots).
+ Default: 100.
+ """,
+ default=100,
+ required=False
+ )
+
+ parser.add_argument(
+ "-p",
+ "--plotmassflow",
+ help="""
+ Plot the mass flow instead of several density(r) plots.
+ Set this to a nonzero number to do this plot.
+ Default: 0.
+ """,
+ default=0,
+ required=False
+ )
+
+ parser.add_argument(
+ "-s",
+ "--slug",
+ help="""
+ The first part of the output filename. For example, for snapshots
+ with the naming scheme keplerian_ring_xxxx.hdf5, this would be the
+ default value of keplerian_ring.
+ """,
+ default="keplerian_ring",
+ required=False
+ )
+
+ parser.add_argument(
+ "-y",
+ "--yt",
+ help="""
+ Use yt to do the plots at the top of the page. If set to anything
+ other than a 'truthy' value, we will use a homebrew plotting
+ setup. Default: False
+ """,
+ default=False,
+ required=False
+ )
+
+ args = vars(parser.parse_args())
+
+ filename = args["slug"]
+
+ if args["beginning"] == -1:
+ # We look for the maximum number of snapshots.
+ numbers = [
+ int(x[len(filename)+1:-5])
+ for x in os.listdir() if
+ x[:len(filename)] == filename and x[-1] == "5"
+ ]
+
+ snapshots = [min(numbers), (max(numbers) - min(numbers)) // 2, max(numbers)]
+ else:
+ snapshots = [args["beginning"], args["middle"], args["end"]]
+
+ figure = plt.figure(figsize=(12, 10))
+ axes = get_axes_grid(figure)
+
+ density_limits = None
+ vlim = None
+
+ for snap, ax in zip(snapshots, tqdm(axes[0:3], desc="Images")):
+ if args["yt"] and ytavail:
+ density_limits, vlim = surface_density_plot(
+ ax,
+ snap,
+ density_limits=density_limits,
+ vlim=vlim
+ )
+
+ figure.subplots_adjust(hspace=0, wspace=0)
+ else:
+ density_limits, vlim = surface_density_plot_no_yt(
+ ax,
+ snap,
+ density_limits=density_limits,
+ vlim=vlim
+ )
+
+ # Now we need to do the density(r) plot.
+
+ # Derived data includes density profiles and chi squared
+ derived_data = get_derived_data(snapshots[0], snapshots[2], binnumber=int(args["nbins"]))
+
+ if args["plotmassflow"]:
+ plot_mass_flow(axes[3], snapshots[0], snapshots[2], derived_data[3])
+ else:
+ plot_density_r(axes[3], derived_data[0], derived_data[1], snapshots)
+
+ plot_chisq(axes[4], snapshots[0], snapshots[2], derived_data[2])
+
+ plot_extra_info(axes[5], "keplerian_ring_0000.hdf5")
+
+ figure.savefig(args["filename"], dpi=300)
+
+
diff --git a/examples/KeplerianRing/run.sh b/examples/KeplerianRing/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ee8f7b247f38b74a204f9caf2bd543b72c4f94fa
--- /dev/null
+++ b/examples/KeplerianRing/run.sh
@@ -0,0 +1,12 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e initial_conditions.hdf5 ]
+then
+ echo "Generating initial conditions for the keplerian ring example..."
+ echo "Please consider choosing your own options before continuing..."
+ python3 makeIC.py
+fi
+
+rm -rf keplerian_ring_*.hdf5
+../swift -g -s -t 1 -v 1 keplerian_ring.yml 2>&1 | tee output.log
diff --git a/examples/KeplerianRing/testplots.py b/examples/KeplerianRing/testplots.py
new file mode 100644
index 0000000000000000000000000000000000000000..1de4964a5618636b73657eec229cbf6e83ef0ff9
--- /dev/null
+++ b/examples/KeplerianRing/testplots.py
@@ -0,0 +1,45 @@
+"""
+###############################################################################
+# 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.
+#
+###############################################################################
+"""
+
+
+if __name__ == "__main__":
+ import yt
+
+ data = yt.load("initial_conditions.hdf5")
+
+ plot = yt.ParticlePlot(
+ data,
+ "particle_position_x",
+ "particle_position_y",
+ "density"
+ )
+
+ plot.save("test.png")
+
+
diff --git a/examples/KeplerianRing/write_gadget.py b/examples/KeplerianRing/write_gadget.py
new file mode 100644
index 0000000000000000000000000000000000000000..d2519cdaf4d78d9c7327419283072b8001e70fcb
--- /dev/null
+++ b/examples/KeplerianRing/write_gadget.py
@@ -0,0 +1,309 @@
+"""
+###############################################################################
+# 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 is a set of helper functions that write particle data to a
+# GADGET .hdf5 file *handle*. It should also function as a piece of
+# documentation on the required attributes for SWIFT to function when running
+# in GADGET compatibility mode.
+#
+# Example Usage:
+#
+# import write_gadget as wg
+# import h5py as h5
+$ import numpy as np
+#
+#
+# N_PARTICLES = 1000
+#
+#
+# with h5.File("test.hdf5", "w") as f:
+# wg.write_header(
+# f,
+# boxsize=100.,
+# flag_entropy=1,
+# np_total=[N_PARTICLES, 0, 0, 0, 0, 0],
+# np_total_hw=[0, 0, 0, 0, 0, 0]
+# )
+#
+# wg.write_runtime_pars(
+# f,
+# periodic_boundary=1,
+# )
+#
+# wg.write_units(
+# f,
+# current=1.,
+# length=1.,
+# mass=1.,
+# temperature=1.,
+# time=1.
+# )
+#
+# wg.write_block(
+# f,
+# part_type=0,
+# pos=np.array([np.arange(N_PARTICLES)] * 3).T,
+# vel=np.array([np.zeros(N_PARTICLES)] * 3).T,
+# ids=np.arange(N_PARTICLES),
+# mass=np.ones(N_PARTICLES,
+# int_energy=np.zeros(N_PARTICLES),
+# smoothing=np.ones(NP_ARTICLES),
+#
+###############################################################################
+"""
+
+
+def write_header(f, boxsize, flag_entropy, np_total, np_total_hw, other=False):
+ """ Writes the "Header" section of the hdf5 file. The parameters in this
+ function that are required are the ones that are required for SWIFT
+ to function; note that the MassTable is **ignored** and that all
+ particle masses should be placed into the particle data arrays.
+
+ @param: f | file handle
+ - the file handle of the hdf5 object (use h5py.File(filename, "w")
+ to open a file handle of the correct type).
+
+ @param boxsize | float / list (2D / 3D)
+ - the boxsize. If a float is given it is assumed that the box is
+ the same size in all directions.
+
+ @param flag_entropy | int (0/1)
+ - sets Flag_Entropy_ICs. This is a historical variable for cross
+ compatibility with Gadget-3
+
+ @param np_total | list (6D)
+ - the total number of particles required of each type.
+
+ Type/Index | Symbolic Type Name
+ ------------|--------------------
+ 0 | Gas
+ 1 | Halo
+ 2 | Disk
+ 3 | Bulge
+ 4 | Stars
+ 5 | Bndry
+
+ @param np_total_hw | list (6D)
+ - the number of high-word particles in the file.
+
+
+ @param other | dictionary | optional
+ - a dictionary with any other parameters that you wish to pass into
+ the file header. They will be passed such that the key is the
+ name of the attribute in the hdf5 file.
+
+ """
+
+ # We'll first build a dictionary to iterate through.
+
+ default_attributes = {
+ "BoxSize" : boxsize,
+ "Flag_Entropy_ICs" : flag_entropy,
+ "NumPart_Total" : np_total,
+ "NumPart_Total_HighWord" : np_total_hw,
+ "NumFilesPerSnapshot" : 1, # required for backwards compatibility
+ "NumPart_ThisFile" : np_total, # Also required for bw compatibility
+ }
+
+ if other:
+ attributes = dict(default_attributes, **other)
+ else:
+ attributes = default_attributes
+
+ header = f.create_group("Header")
+
+ # For some reason there isn't a direct dictionary interface (at least one
+ # that is documented, so we are stuck doing this loop...
+
+ for name, value in attributes.items():
+ header.attrs[name] = value
+
+ return
+
+
+def write_runtime_pars(f, periodic_boundary, other=False):
+ """ Writes the "RuntimeParams" section in the hdf5 file. The parameters in
+ this function that are required are also required for SWIFT to function.
+ If you wish to pass extra arguments into the runtime parameters you
+ may do that by providing a dictionary to other.
+
+ @param: f | file handle
+ - the file handle of the hdf5 object (use h5py.File(filename, "w")
+ to open a file handle of the correct type).
+
+ @param: periodic_boundary | int (0/1)
+ - the condition for periodic boundary conditions -- they are 'on'
+ if this variable is 1, and off if it is 0. Note that SWIFT
+ currently requires periodic boundary conditions to run (as of
+ August 2017).
+
+ @param other | dictionary | optional
+ - a dictionary with any other parameters that you wish to pass into
+ the RuntimePars. They will be passed such that the key is the
+ name of the attribute in the hdf5 file.
+ """
+
+ # First build the dictionary
+
+ default_attributes = {
+ "PeriodicBoundariesOn" : periodic_boundary,
+ }
+
+ if other:
+ attributes = dict(default_attributes, **other)
+ else:
+ attributes = default_attributes
+
+ runtime = f.create_group("RuntimePars")
+
+ for name, value in attributes.items():
+ runtime.attrs[name] = value
+
+ return
+
+
+def write_units(f, current, length, mass, temperature, time, other=False):
+ """ Writes the "RuntimeParams" section in the hdf5 file. The parameters in
+ this function that are required are also required for SWIFT to function.
+ If you wish to pass extra arguments into the runtime parameters you
+ may do that by providing a dictionary to other.
+
+ @param: f | file handle
+ - the file handle of the hdf5 object (use h5py.File(filename, "w")
+ to open a file handle of the correct type).
+
+ @param: current | float
+ - the current conversion factor in cgs units.
+
+ @param: length | float
+ - the length conversion factor in cgs units.
+
+ @param: mass | float
+ - the mass conversion factor in cgs units.
+
+ @param: temperature | float
+ - the temperature conversion factor in cgs units.
+
+ @param: time | float
+ - the time conversion factor in cgs units.
+
+ @param: other | dictionary | optional
+ - a dictionary with any other parameters that you wish to pass into
+ the units attributes. They will be passed such that the key is
+ the name of the attribute in the hdf5 file.
+ """
+
+ # First build the dictionary
+
+ default_attributes = {
+ "Unit current in cgs (U_I)": current,
+ "Unit length in cgs (U_L)": length,
+ "Unit mass in cgs (U_M)": mass,
+ "Unit temperature in cgs (U_T)": temperature,
+ "Unit time in cgs (U_t)": time,
+ }
+
+ if other:
+ attributes = dict(default_attributes, **other)
+ else:
+ attributes = default_attributes
+
+ units = f.create_group("Units")
+
+ for name, value in attributes.items():
+ units.attrs[name] = value
+
+ return
+
+
+def write_block(f, part_type, pos, vel, ids, mass, int_energy, smoothing, other=False):
+ """ Writes a given block of data to PartType{part_type}. The above
+ required parameters are required for SWIFT to run.
+
+ @param: f | file handle
+ - the file handle of the hdf5 object.
+
+ @param part_type | int (0-5):
+ - the identifiying number of the particle type.
+
+ Type/Index | Symbolic Type Name
+ ------------|--------------------
+ 0 | Gas
+ 1 | Halo
+ 2 | Disk
+ 3 | Bulge
+ 4 | Stars
+ 5 | Bndry
+
+ @param pos | numpy.array
+ - the array of particle positions with shape (n_particles, 3).
+
+ @param vel | numpy.array
+ - the array of particle velocities with shape (n_particles, 3).
+
+ @param ids | numpy.array
+ - the ids of the particles. Please note that particle IDs in
+ SWIFT must be strictly positive. Shape (n_particles, 1)
+
+ @param mass | numpy.array
+ - the masses of the particles. In SWIFT MassTable in the header
+ is ignored and so particle masses must be provided here. Shape
+ (n_particles, 1)
+
+ @param int_energy | numpy.array
+ - the internal energies of the particles. Shape (n_particles, 1)
+
+ @param smoothing | numpy.array
+ - the smoothing lenghts of the individual particles. Please note
+ that this cannot be supplied only in the parameterfile and must
+ be provided on a particle-by-particle basis in SWIFT. Shape
+ (n_particles, 1)
+
+ @param: other | dictionary | optional
+ - a dictionary with any other parameters that you wish to pass into
+ the particle data. They will be passed such that the key is
+ the name of the dataset in the hdf5 file.
+ """
+
+ # Build the dictionary
+
+ default_data = {
+ "Coordinates" : pos,
+ "Velocities" : vel,
+ "ParticleIDs" : ids,
+ "Masses" : mass,
+ "InternalEnergy" : int_energy,
+ "SmoothingLength" : smoothing,
+ }
+
+ if other:
+ data = dict(default_data, **other)
+ else:
+ data = default_data
+
+ particles = f.create_group(f"PartType{part_type}")
+
+ for name, value in data.items():
+ particles.create_dataset(name, data=value)
+
+ return
+
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index b7b87e2224b78d18db751023db2083599d82b901..b7fa7dd4aa00249d009037961888c89abc0d0313 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -659,9 +659,9 @@ runner_iact_nonsym_1_vec_force(
vector viz, vector pirho, vector grad_hi, vector piPOrho2, vector balsara_i,
vector ci, float *Vjx, float *Vjy, float *Vjz, float *Pjrho, float *Grad_hj,
float *PjPOrho2, float *Balsara_j, float *Cj, float *Mj, vector hi_inv,
- vector hj_inv, const float a, const float H, vector *a_hydro_xSum, vector *a_hydro_ySum,
- vector *a_hydro_zSum, vector *h_dtSum, vector *v_sigSum,
- vector *entropy_dtSum, mask_t mask) {
+ vector hj_inv, const float a, const float H, vector *a_hydro_xSum,
+ vector *a_hydro_ySum, vector *a_hydro_zSum, vector *h_dtSum,
+ vector *v_sigSum, vector *entropy_dtSum, mask_t mask) {
#ifdef WITH_VECTORIZATION
@@ -721,13 +721,16 @@ runner_iact_nonsym_1_vec_force(
dvz.v = vec_sub(viz.v, vjz.v);
/* Compute dv dot r. */
- dvdr.v = vec_fma(dvx.v, dx->v, vec_fma(dvy.v, dy->v, vec_fma(dvz.v, dz->v, vec_mul(v_a2_Hubble.v, r2->v))));
+ dvdr.v =
+ vec_fma(dvx.v, dx->v,
+ vec_fma(dvy.v, dy->v,
+ vec_fma(dvz.v, dz->v, vec_mul(v_a2_Hubble.v, r2->v))));
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
omega_ij.v = vec_fmin(dvdr.v, vec_setzero());
- mu_ij.v =
- vec_mul(v_fac_mu.v, vec_mul(ri.v, omega_ij.v)); /* This is 0 or negative */
+ mu_ij.v = vec_mul(v_fac_mu.v,
+ vec_mul(ri.v, omega_ij.v)); /* This is 0 or negative */
/* Compute signal velocity */
v_sig.v = vec_fnma(vec_set1(3.f), mu_ij.v, vec_add(ci.v, cj.v));
@@ -790,9 +793,10 @@ runner_iact_nonsym_2_vec_force(
vector viz, vector pirho, vector grad_hi, vector piPOrho2, vector balsara_i,
vector ci, float *Vjx, float *Vjy, float *Vjz, float *Pjrho, float *Grad_hj,
float *PjPOrho2, float *Balsara_j, float *Cj, float *Mj, vector hi_inv,
- float *Hj_inv, const float a, const float H, vector *a_hydro_xSum, vector *a_hydro_ySum,
- vector *a_hydro_zSum, vector *h_dtSum, vector *v_sigSum,
- vector *entropy_dtSum, mask_t mask, mask_t mask_2, short mask_cond) {
+ float *Hj_inv, const float a, const float H, vector *a_hydro_xSum,
+ vector *a_hydro_ySum, vector *a_hydro_zSum, vector *h_dtSum,
+ vector *v_sigSum, vector *entropy_dtSum, mask_t mask, mask_t mask_2,
+ short mask_cond) {
#ifdef WITH_VECTORIZATION
@@ -897,16 +901,20 @@ runner_iact_nonsym_2_vec_force(
dvz_2.v = vec_sub(viz.v, vjz_2.v);
/* Compute dv dot r. */
- dvdr.v = vec_fma(dvx.v, dx.v, vec_fma(dvy.v, dy.v, vec_fma(dvz.v, dz.v, vec_mul(v_a2_Hubble.v, r2.v))));
- dvdr_2.v = vec_fma(dvx_2.v, dx_2.v,
- vec_fma(dvy_2.v, dy_2.v, vec_fma(dvz_2.v, dz_2.v, vec_mul(v_a2_Hubble.v, r2_2.v))));
+ dvdr.v = vec_fma(
+ dvx.v, dx.v,
+ vec_fma(dvy.v, dy.v, vec_fma(dvz.v, dz.v, vec_mul(v_a2_Hubble.v, r2.v))));
+ dvdr_2.v = vec_fma(
+ dvx_2.v, dx_2.v,
+ vec_fma(dvy_2.v, dy_2.v,
+ vec_fma(dvz_2.v, dz_2.v, vec_mul(v_a2_Hubble.v, r2_2.v))));
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
omega_ij.v = vec_fmin(dvdr.v, vec_setzero());
omega_ij_2.v = vec_fmin(dvdr_2.v, vec_setzero());
- mu_ij.v =
- vec_mul(v_fac_mu.v, vec_mul(ri.v, omega_ij.v)); /* This is 0 or negative */
+ mu_ij.v = vec_mul(v_fac_mu.v,
+ vec_mul(ri.v, omega_ij.v)); /* This is 0 or negative */
mu_ij_2.v = vec_mul(
v_fac_mu.v, vec_mul(ri_2.v, omega_ij_2.v)); /* This is 0 or negative */
diff --git a/src/potential.h b/src/potential.h
index bcb3fd284021fba339ba494c90b81c91bd2ce72f..680d4e235fdf7a7666901f34a82f62feda4ae9bb 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -38,6 +38,10 @@
#include "./potential/disc_patch/potential.h"
#elif defined(EXTERNAL_POTENTIAL_SINE_WAVE)
#include "./potential/sine_wave/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_POINTMASS_RING)
+#include "./potential/point_mass_ring/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_POINTMASS_SOFT)
+#include "./potential/point_mass_softened/potential.h"
#else
#error "Invalid choice of external potential"
#endif
diff --git a/src/potential/point_mass_ring/potential.h b/src/potential/point_mass_ring/potential.h
new file mode 100644
index 0000000000000000000000000000000000000000..ebf047ea7c1f946536300f976713893e66295c59
--- /dev/null
+++ b/src/potential/point_mass_ring/potential.h
@@ -0,0 +1,226 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_POTENTIAL_POINT_MASS_H
+#define SWIFT_POTENTIAL_POINT_MASS_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief External Potential Properties - Point mass case
+ */
+struct external_potential {
+
+ /*! Position of the point mass */
+ double x, y, z;
+
+ /*! Mass */
+ double mass;
+
+ /*! Time-step condition pre-factor */
+ float timestep_mult;
+};
+
+/**
+ * @brief Computes the time-step due to the acceleration from a point mass
+ * based on Hopkins' central potential stuff (i.e. using a 'softened'
+ * gravitational edge).
+ *
+ * We pass in the time for simulations where the potential evolves with time.
+ *
+ * @param time The current time.
+ * @param potential The properties of the external potential.
+ * @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) {
+
+ const float G_newton = phys_const->const_newton_G;
+ const float dx = g->x[0] - potential->x;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+ const float r = sqrtf(dx * dx + dy * dy + dz * dz);
+ const float rinv = 1.f / r;
+ const float rinv2 = rinv * rinv;
+ const float rinv3 = rinv2 * rinv;
+ const float drdv = (g->x[0] - potential->x) * (g->v_full[0]) +
+ (g->x[1] - potential->y) * (g->v_full[1]) +
+ (g->x[2] - potential->z) * (g->v_full[2]);
+ float factor;
+
+ if (r < 0.175) {
+ // We need to flatten gravity as in SWIFT the timestepping is different
+ // than in GIZMO. This causes particles to become trapped close to the
+ // central point mass!
+ factor = 0.;
+ } else if (r < 0.35) {
+ factor = (8.16 * r * r) + 1.f - r / 0.35;
+ } else if (r > 2.1) {
+ factor = 1.f + (r - 2.1) / 0.1;
+ } else {
+ // 0.35 > r > 2.1
+ factor = 1.f;
+ }
+
+ const float dota_x = G_newton * potential->mass * rinv3 *
+ (-g->v_full[0] + 3.f * rinv2 * drdv * dx) * factor;
+ const float dota_y = G_newton * potential->mass * rinv3 *
+ (-g->v_full[1] + 3.f * rinv2 * drdv * dy) * factor;
+ const float dota_z = G_newton * potential->mass * rinv3 *
+ (-g->v_full[2] + 3.f * rinv2 * drdv * dz) * factor;
+
+ const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
+ const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
+ g->a_grav[2] * g->a_grav[2];
+
+ if (fabsf(dota_2) > 0.f)
+ return potential->timestep_mult * sqrtf(a_2 / dota_2);
+ else
+ return FLT_MAX;
+}
+
+/**
+ * @brief Computes the gravitational acceleration of a particle due to a
+ * point mass
+ *
+ * Note that the accelerations are multiplied by Newton's G constant later
+ * on.
+ *
+ * We pass in the time for simulations where the potential evolves with time.
+ *
+ * @param time The current time.
+ * @param potential The proerties of the external potential.
+ * @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* restrict potential,
+ const struct phys_const* restrict phys_const, struct gpart* restrict g) {
+
+ const float dx = g->x[0] - potential->x;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+ const float r = sqrtf(dx * dx + dy * dy + dz * dz);
+ const float rinv = 1.f / r;
+ const float rinv2 = rinv * rinv;
+ const float rinv3 = rinv * rinv2;
+ float factor = 1.f;
+
+ if (r < 0.175) {
+ // We need to flatten gravity as in SWIFT the timestepping is different
+ // than in GIZMO. This causes particles to become trapped close to the
+ // central point mass!
+ factor = 0.;
+ printf("Help me, I'm trapped! (r = %f id = %lld)\n", r,
+ g->id_or_neg_offset);
+ } else if (r < 0.35) {
+ factor = (8.16 * r * r) + 1.f - r / 0.35;
+ printf("Factor is %f, r is %f \n", factor, r);
+ } else if (r > 2.1) {
+ factor = 1.f + (r - 2.1) / 0.1;
+ } else {
+ // 0.35 > r > 2.1
+ factor = 1.f;
+ }
+
+ g->a_grav[0] += -potential->mass * dx * rinv3 * factor;
+ g->a_grav[1] += -potential->mass * dy * rinv3 * factor;
+ g->a_grav[2] += -potential->mass * dz * rinv3 * factor;
+}
+
+/**
+ * @brief Computes the gravitational potential energy of a particle in a point
+ * mass potential.
+ *
+ * @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 dx = g->x[0] - potential->x;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+ const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz);
+ return -phys_const->const_newton_G * potential->mass * rinv;
+}
+
+/**
+ * @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 s The #space we run in.
+ * @param potential The external potential properties to initialize
+ */
+static INLINE void potential_init_backend(
+ const struct swift_params* parameter_file,
+ const struct phys_const* phys_const, const struct unit_system* us,
+ const struct space* s, struct external_potential* potential) {
+
+ potential->x =
+ parser_get_param_double(parameter_file, "PointMassPotential:position_x");
+ potential->y =
+ parser_get_param_double(parameter_file, "PointMassPotential:position_y");
+ potential->z =
+ parser_get_param_double(parameter_file, "PointMassPotential:position_z");
+ potential->mass =
+ parser_get_param_double(parameter_file, "PointMassPotential:mass");
+ potential->timestep_mult = parser_get_param_float(
+ parameter_file, "PointMassPotential:timestep_mult");
+}
+
+/**
+ * @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 'Point mass' with properties (x,y,z) = (%e, %e, "
+ "%e), M = %e timestep multiplier = %e.",
+ potential->x, potential->y, potential->z, potential->mass,
+ potential->timestep_mult);
+}
+
+#endif /* SWIFT_POTENTIAL_POINT_MASS_H */
diff --git a/src/potential/point_mass_softened/potential.h b/src/potential/point_mass_softened/potential.h
new file mode 100644
index 0000000000000000000000000000000000000000..83a79ea3cddbff37fdd70d58d70afcaf46f7bc0e
--- /dev/null
+++ b/src/potential/point_mass_softened/potential.h
@@ -0,0 +1,215 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_POTENTIAL_POINT_MASS_H
+#define SWIFT_POTENTIAL_POINT_MASS_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief External Potential Properties - Point mass case
+ *
+ * This file contains the softened central potential. This has a 'softening
+ * factor' that prevents the potential from becoming singular at the centre
+ * of the potential, i.e. it has the form
+ *
+ * $$ \phi \propto 1/(r^2 + \eta^2) $$
+ *
+ * where $\eta$ is some small value.
+ */
+struct external_potential {
+
+ /*! Position of the point mass */
+ double x, y, z;
+
+ /*! Mass */
+ double mass;
+
+ /*! Time-step condition pre-factor */
+ float timestep_mult;
+
+ /*! Potential softening */
+ float softening;
+};
+
+/**
+ * @brief Computes the time-step due to the acceleration from a point mass
+ *
+ * We pass in the time for simulations where the potential evolves with time.
+ *
+ * @param time The current time.
+ * @param potential The properties of the externa potential.
+ * @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) {
+
+ const float dx = g->x[0] - potential->x;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+
+ const float softening2 = potential->softening * potential->softening;
+ const float r2 = dx * dx + dy * dy + dz * dz;
+ const float rinv = 1.f / sqrtf(r2);
+ const float rinv2_softened = 1.f / (r2 + softening2);
+
+ /* G * M / (r^2 + eta^2)^{3/2} */
+ const float GMr32 = phys_const->const_newton_G * potential->mass *
+ sqrtf(rinv2_softened * rinv2_softened * rinv2_softened);
+
+ /* This is actually dr dot v */
+ const float drdv =
+ (dx) * (g->v_full[0]) + (dy) * (g->v_full[1]) + (dz) * (g->v_full[2]);
+
+ /* We want da/dt */
+ /* da/dt = -GM/(r^2 + \eta^2)^{3/2} *
+ * (\vec{v} - 3 \vec{r} \cdot \vec{v} \hat{r} /
+ * (r^2 + \eta^2)) */
+ const float dota_x =
+ GMr32 * (3.f * drdv * dx * rinv2_softened * rinv - g->v_full[0]);
+ const float dota_y =
+ GMr32 * (3.f * drdv * dy * rinv2_softened * rinv - g->v_full[0]);
+ const float dota_z =
+ GMr32 * (3.f * drdv * dz * rinv2_softened * rinv - g->v_full[0]);
+
+ const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
+ const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
+ g->a_grav[2] * g->a_grav[2];
+
+ if (fabsf(dota_2) > 0.f)
+ return potential->timestep_mult * sqrtf(a_2 / dota_2);
+ else
+ return FLT_MAX;
+}
+
+/**
+ * @brief Computes the gravitational acceleration of a particle due to a
+ * point mass
+ *
+ * Note that the accelerations are multiplied by Newton's G constant later
+ * on.
+ *
+ * We pass in the time for simulations where the potential evolves with time.
+ *
+ * @param time The current time.
+ * @param potential The proerties of the external potential.
+ * @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* restrict potential,
+ const struct phys_const* restrict phys_const, struct gpart* restrict g) {
+
+ const float dx = g->x[0] - potential->x;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+ const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz +
+ potential->softening * potential->softening);
+ const float rinv3 = rinv * rinv * rinv;
+
+ g->a_grav[0] += -potential->mass * dx * rinv3;
+ g->a_grav[1] += -potential->mass * dy * rinv3;
+ g->a_grav[2] += -potential->mass * dz * rinv3;
+}
+
+/**
+ * @brief Computes the gravitational potential energy of a particle in a point
+ * mass potential.
+ *
+ * @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 dx = g->x[0] - potential->x;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+ const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz +
+ potential->softening * potential->softening);
+
+ return -phys_const->const_newton_G * potential->mass * rinv;
+}
+
+/**
+ * @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 s The #space we run in.
+ * @param potential The external potential properties to initialize
+ */
+static INLINE void potential_init_backend(
+ const struct swift_params* parameter_file,
+ const struct phys_const* phys_const, const struct unit_system* us,
+ const struct space* s, struct external_potential* potential) {
+
+ potential->x =
+ parser_get_param_double(parameter_file, "PointMassPotential:position_x");
+ potential->y =
+ parser_get_param_double(parameter_file, "PointMassPotential:position_y");
+ potential->z =
+ parser_get_param_double(parameter_file, "PointMassPotential:position_z");
+ potential->mass =
+ parser_get_param_double(parameter_file, "PointMassPotential:mass");
+ potential->timestep_mult = parser_get_param_float(
+ parameter_file, "PointMassPotential:timestep_mult");
+ potential->softening =
+ parser_get_param_float(parameter_file, "PointMassPotential:softening");
+}
+
+/**
+ * @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 'Point mass' with properties (x,y,z) = (%e, %e, "
+ "%e), M = %e timestep multiplier = %e, softening = %e.",
+ potential->x, potential->y, potential->z, potential->mass,
+ potential->timestep_mult, potential->softening);
+}
+
+#endif /* SWIFT_POTENTIAL_POINT_MASS_H */
diff --git a/src/runner_doiact_vec.c b/src/runner_doiact_vec.c
index c2b807e7061f4b8ac7c2066bb03f2a7c93d9610b..0daea0cf8a8d868ccb14f7d9b27f04aca5a677c8 100644
--- a/src/runner_doiact_vec.c
+++ b/src/runner_doiact_vec.c
@@ -1143,7 +1143,7 @@ void runner_doself2_force_vec(struct runner *r, struct cell *restrict c) {
/* Cosmological terms */
const float a = cosmo->a;
const float H = cosmo->H;
-
+
/* Loop over the particles in the cell. */
for (int pid = 0; pid < count; pid++) {
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 2ded11b162e88380a4ad5c8ad23bb76b22478ec3..9c5dd36415970ff2a53220aa56cecba6fc5fe193 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -610,8 +610,9 @@ void test_force_interactions(struct part test_part, struct part *parts,
(viz_vec), rhoi_vec, grad_hi_vec, pOrhoi2_vec, balsara_i_vec,
ci_vec, &(vjxq[i]), &(vjyq[i]), &(vjzq[i]), &(rhojq[i]),
&(grad_hjq[i]), &(pOrhoj2q[i]), &(balsarajq[i]), &(cjq[i]),
- &(mjq[i]), hi_inv_vec, &(hj_invq[i]), a, H, &a_hydro_xSum, &a_hydro_ySum,
- &a_hydro_zSum, &h_dtSum, &v_sigSum, &entropy_dtSum, mask, mask2, 0);
+ &(mjq[i]), hi_inv_vec, &(hj_invq[i]), a, H, &a_hydro_xSum,
+ &a_hydro_ySum, &a_hydro_zSum, &h_dtSum, &v_sigSum, &entropy_dtSum,
+ mask, mask2, 0);
} else { /* Only use one vector for interaction. */
vector my_r2, my_dx, my_dy, my_dz, hj, hj_inv;
@@ -626,7 +627,7 @@ void test_force_interactions(struct part test_part, struct part *parts,
&my_r2, &my_dx, &my_dy, &my_dz, vix_vec, viy_vec, viz_vec, rhoi_vec,
grad_hi_vec, pOrhoi2_vec, balsara_i_vec, ci_vec, &(vjxq[i]),
&(vjyq[i]), &(vjzq[i]), &(rhojq[i]), &(grad_hjq[i]), &(pOrhoj2q[i]),
- &(balsarajq[i]), &(cjq[i]), &(mjq[i]), hi_inv_vec, hj_inv, a, H,
+ &(balsarajq[i]), &(cjq[i]), &(mjq[i]), hi_inv_vec, hj_inv, a, H,
&a_hydro_xSum, &a_hydro_ySum, &a_hydro_zSum, &h_dtSum, &v_sigSum,
&entropy_dtSum, mask);
}