diff --git a/.gitignore b/.gitignore
index cc936154b298d48a7661070b35dc7d2f326d56ce..83f22370c4ec8421a9b012a49416569e7bd7a58f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
*~
+*.hdf5
Makefile
Makefile.in
@@ -23,14 +24,12 @@ doc/Doxyfile
examples/swift
examples/swift_mpi
examples/*/*.xmf
-examples/*/*.hdf5
examples/*/*.h5
examples/*/*.png
examples/*/*.txt
examples/*/*.dot
examples/*/used_parameters.yml
examples/*/*/*.xmf
-examples/*/*/*.hdf5
examples/*/*/*.png
examples/*/*/*.txt
examples/*/*/*.dot
@@ -73,7 +72,6 @@ tests/test_nonsym_force_1_vec.dat
tests/test_nonsym_force_2_vec.dat
tests/testGreetings
tests/testReading
-tests/input.hdf5
tests/testSingle
tests/testTimeIntegration
tests/testSPHStep
@@ -107,6 +105,8 @@ tests/testDump
tests/testLogger
tests/benchmarkInteractions
tests/testGravityDerivatives
+tests/testPotentialSelf
+tests/testPotentialPair
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
@@ -277,9 +277,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 3826441730116b6d54f0c5da5416fa7d5b596e3d..8cdce3e0c72842abdf075213e4a0e6610ac1f441 100644
--- a/configure.ac
+++ b/configure.ac
@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Init the project.
-AC_INIT([SWIFT],[0.6.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
+AC_INIT([SWIFT],[0.7.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
swift_config_flags="$*"
# Need to define this, instead of using fifth argument of AC_INIT, until 2.64.
@@ -407,13 +407,47 @@ AC_HEADER_STDC
# Check for the libraries we will need.
AC_CHECK_LIB(m,sqrt,,AC_MSG_ERROR(something is wrong with the math library!))
-# Check for GSL
+# Check for GSL. We test for this in the standard directories by default,
+# and only disable if using --with-gsl=no or --without-gsl. When a value
+# is given GSL must be found.
have_gsl="no"
-AC_CHECK_LIB([gslcblas], [cblas_dgemm])
-AC_CHECK_LIB([gsl], [gsl_integration_qag])
-if test "x$ac_cv_lib_gslcblas_cblas_dgemm" != "x"; then
- have_gsl="yes"
+AC_ARG_WITH([gsl],
+ [AS_HELP_STRING([--with-gsl=PATH],
+ [root directory where GSL is installed @<:@yes/no@:>@]
+ )],
+ [with_gsl="$withval"],
+ [with_gsl="test"]
+)
+if test "x$with_gsl" != "xno"; then
+ if test "x$with_gsl" != "xyes" -a "x$with_gsl" != "xtest" -a "x$with_gsl" != "x"; then
+ GSL_LIBS="-L$with_gsl/lib -lgsl -lgslcblas"
+ GSL_INCS="-I$with_gsl/include"
+ else
+ GSL_LIBS="-lgsl -lgslcblas"
+ GSL_INCS=""
+ fi
+ # GSL is not specified, so just check if we have it.
+ if test "x$with_gsl" = "xtest"; then
+ AC_CHECK_LIB([gslcblas],[cblas_dgemm],[have_gsl="yes"],[have_gsl="no"],$GSL_LIBS)
+ if test "x$have_gsl" != "xno"; then
+ AC_DEFINE([HAVE_LIBGSLCBLAS],1,[The GSL CBLAS library appears to be present.])
+ AC_CHECK_LIB([gsl],[gsl_integration_qag],
+ AC_DEFINE([HAVE_LIBGSL],1,[The GSL library appears to be present.]),
+ [have_gsl="no"],$GSL_LIBS)
+ fi
+ else
+ AC_CHECK_LIB([gslcblas],[cblas_dgemm],
+ AC_DEFINE([HAVE_LIBGSLCBLAS],1,[The GSL CBLAS library appears to be present.]),
+ AC_MSG_ERROR(something is wrong with the GSL CBLAS library!), $GSL_LIBS)
+ AC_CHECK_LIB([gsl],[gsl_integration_qag],
+ AC_DEFINE([HAVE_LIBGSL],1,[The GSL library appears to be present.]),
+ AC_MSG_ERROR(something is wrong with the GSL library!), $GSL_LIBS)
+ have_gsl="yes"
+ fi
fi
+AC_SUBST([GSL_LIBS])
+AC_SUBST([GSL_INCS])
+AM_CONDITIONAL([HAVEGSL],[test -n "$GSL_LIBS"])
# Check for pthreads.
AX_PTHREAD([LIBS="$PTHREAD_LIBS $LIBS" CFLAGS="$CFLAGS $PTHREAD_CFLAGS"
@@ -443,7 +477,7 @@ AC_ARG_WITH([parmetis],
[AS_HELP_STRING([--with-parmetis=PATH],
[root directory where ParMETIS is installed @<:@yes/no@:>@]
)],
- [],
+ [with_parmetis="$withval"],
[with_parmetis="no"]
)
@@ -506,13 +540,13 @@ AH_VERBATIM([__STDC_FORMAT_MACROS],
#define __STDC_FORMAT_MACROS 1
#endif])
-# Check for grackle.
+# Check for grackle.
have_grackle="no"
AC_ARG_WITH([grackle],
[AS_HELP_STRING([--with-grackle=PATH],
[root directory where grackle is installed @<:@yes/no@:>@]
)],
- [],
+ [with_grackle="$withval"],
[with_grackle="no"]
)
if test "x$with_grackle" != "xno"; then
@@ -1013,7 +1047,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"]
@@ -1034,6 +1068,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/EAGLE_100/eagle_100.yml b/examples/EAGLE_100/eagle_100.yml
index 4c3caa08f694ba9d40d9278dc7dcc1a339eed2a4..9cf1c9432ed6c79e4b91262cfdd30bad14833558 100644
--- a/examples/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_100/eagle_100.yml
@@ -6,6 +6,15 @@ InternalUnitSystem:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
+# Cosmological parameters
+Cosmology:
+ h: 0.6777 # Reduced Hubble constant
+ a_begin: 0.9090909 # Initial scale-factor of the simulation
+ a_end: 1.0 # Final scale factor of the simulation
+ Omega_m: 0.307 # Matter density parameter
+ Omega_lambda: 0.693 # Dark-energy density parameter
+ Omega_b: 0.0455 # Baryon density parameter
+
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
@@ -19,7 +28,7 @@ Scheduler:
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
+ time_first: 1. # Time of the first output (in internal units)
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
@@ -29,8 +38,8 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
- epsilon: 0.0001 # Softening length (in internal units).
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ epsilon: 0.001 # Softening length (in internal units).
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index af55da76ff513548bee52e51993fa8b78cfdc9e3..d7aedf591a6988f89ad700f0a6729b3682bb6881 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -6,13 +6,6 @@ InternalUnitSystem:
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: 1e-2 # The end time of the simulation (in internal units).
- dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
-
# Cosmological parameters
Cosmology:
h: 0.6777 # Reduced Hubble constant
@@ -21,6 +14,13 @@ Cosmology:
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
Omega_b: 0.0455 # Baryon density parameter
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 1e-2 # The end time of the simulation (in internal units).
+ dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
Scheduler:
max_top_level_cells: 15
@@ -30,6 +30,7 @@ Snapshots:
basename: eagle # Common part of the name of output files
time_first: 1. # Time of the first output (in internal units)
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
+ compression: 4
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index c01abfe94b98e8366ac05f1e3bce50a3931738c5..24daa13e06c2058d8a93243e9838b6f6b2cfe9bb 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -6,6 +6,15 @@ InternalUnitSystem:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
+# Cosmological parameters
+Cosmology:
+ h: 0.6777 # Reduced Hubble constant
+ a_begin: 0.9090909 # Initial scale-factor of the simulation
+ a_end: 1.0 # Final scale factor of the simulation
+ Omega_m: 0.307 # Matter density parameter
+ Omega_lambda: 0.693 # Dark-energy density parameter
+ Omega_b: 0.0455 # Baryon density parameter
+
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
@@ -19,8 +28,9 @@ Scheduler:
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
+ time_first: 1. # Time of the first output (in internal units)
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
+ compression: 4
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml
index 0fe1007414e7190ae9cc6f31aa5a16cec75290db..adf64f6963e2a0bc545637757f5826efd0c027a9 100644
--- a/examples/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_50/eagle_50.yml
@@ -6,6 +6,15 @@ InternalUnitSystem:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
+# Cosmological parameters
+Cosmology:
+ h: 0.6777 # Reduced Hubble constant
+ a_begin: 0.9090909 # Initial scale-factor of the simulation
+ a_end: 1.0 # Final scale factor of the simulation
+ Omega_m: 0.307 # Matter density parameter
+ Omega_lambda: 0.693 # Dark-energy density parameter
+ Omega_b: 0.0455 # Baryon density parameter
+
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
@@ -19,7 +28,7 @@ Scheduler:
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
+ time_first: 1. # Time of the first output (in internal units)
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index ab33b79e80e69bcdac6a10f793930eca203cb26f..d0055607d3cd7c1453e3b93ed8b21cd5ff2673d5 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -6,6 +6,15 @@ InternalUnitSystem:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
+# Cosmological parameters
+Cosmology:
+ h: 0.6777 # Reduced Hubble constant
+ a_begin: 0.9090909 # Initial scale-factor of the simulation
+ a_end: 1.0 # Final scale factor of the simulation
+ Omega_m: 0.307 # Matter density parameter
+ Omega_lambda: 0.693 # Dark-energy density parameter
+ Omega_b: 0.0455 # Baryon density parameter
+
Scheduler:
max_top_level_cells: 8
@@ -15,12 +24,13 @@ TimeIntegration:
time_end: 1e-2 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
-
+
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
+ compression: 4
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/EvrardCollapse_3D/evrard.yml b/examples/EvrardCollapse_3D/evrard.yml
index 006a22e65d3f674f124ce6c4994e752ba39cd1e1..94cd6e7a95d2671f7021ae3cd2ba40477cd8e7fd 100644
--- a/examples/EvrardCollapse_3D/evrard.yml
+++ b/examples/EvrardCollapse_3D/evrard.yml
@@ -18,7 +18,8 @@ Snapshots:
basename: evrard # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
-
+ compression: 4
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
diff --git a/examples/GreshoVortex_3D/gresho.yml b/examples/GreshoVortex_3D/gresho.yml
index df941450196a7de6cd1471e1d258756ca8c36fb1..f42fcbfd00941e7b9c5c09c0d2e3118f5cc1f57d 100644
--- a/examples/GreshoVortex_3D/gresho.yml
+++ b/examples/GreshoVortex_3D/gresho.yml
@@ -21,7 +21,8 @@ Snapshots:
basename: gresho # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-1 # Time difference between consecutive outputs (in internal units)
-
+ compression: 4
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
diff --git a/examples/KelvinHelmholtzGrowthRate_3D/kelvinHelmholtzGrowthRate.yml b/examples/KelvinHelmholtzGrowthRate_3D/kelvinHelmholtzGrowthRate.yml
index 380dc2ab3a530e89b952aa41f425e50709d73ee9..3133e2769e81b80c18760e8258665fc1a6eee6ca 100644
--- a/examples/KelvinHelmholtzGrowthRate_3D/kelvinHelmholtzGrowthRate.yml
+++ b/examples/KelvinHelmholtzGrowthRate_3D/kelvinHelmholtzGrowthRate.yml
@@ -18,7 +18,8 @@ Snapshots:
basename: kelvinHelmholtzGrowthRate # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.04 # Time difference between consecutive outputs (in internal units)
-
+ compression: 4
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
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/examples/Makefile.am b/examples/Makefile.am
index c5537441d5792f00cac9f71a3c9b99a71b13b8e6..23c7c4a44f1221358066df56ee9f8de8404ad6a3 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -19,12 +19,12 @@
MYFLAGS =
# Add the source directory and debug to CFLAGS
-AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS)
+AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS) $(GSL_INCS)
AM_LDFLAGS = $(HDF5_LDFLAGS)
# Extra libraries.
-EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(GRACKLE_LIBS)
+EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS)
# MPI libraries.
MPI_LIBS = $(PARMETIS_LIBS) $(MPI_THREAD_LIBS)
diff --git a/examples/Noh_3D/noh.yml b/examples/Noh_3D/noh.yml
index 1d126f19babd0c9fe28afff907b3fe8259467a24..88119827501e17cc26742e8cad92a3611e83faa7 100644
--- a/examples/Noh_3D/noh.yml
+++ b/examples/Noh_3D/noh.yml
@@ -18,7 +18,8 @@ Snapshots:
basename: noh # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 5e-2 # Time difference between consecutive outputs (in internal units)
-
+ compression: 4
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 3e-3 # Time between statistics output
diff --git a/examples/SedovBlast_3D/sedov.yml b/examples/SedovBlast_3D/sedov.yml
index b80be1954491b19234ce7385baed83931d82433f..057e5a18acba2a1e8211d40f15ffde35b2019a9b 100644
--- a/examples/SedovBlast_3D/sedov.yml
+++ b/examples/SedovBlast_3D/sedov.yml
@@ -18,7 +18,8 @@ Snapshots:
basename: sedov # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
-
+ compression: 4
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-3 # Time between statistics output
diff --git a/examples/SodShockSpherical_3D/sodShock.yml b/examples/SodShockSpherical_3D/sodShock.yml
index a26ab95b21c782ce83310038432ac08df0e024c3..52c06dd65ec0b37a0ec0707315ccd15356a7b2d6 100644
--- a/examples/SodShockSpherical_3D/sodShock.yml
+++ b/examples/SodShockSpherical_3D/sodShock.yml
@@ -18,7 +18,8 @@ Snapshots:
basename: sodShock # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
-
+ compression: 4
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
diff --git a/examples/SodShock_3D/sodShock.yml b/examples/SodShock_3D/sodShock.yml
index 51a188b6d4537d490cb837a03dab15f74c3b083c..26fbfe4faacf2bbbfd9b077bf9f9c075ce93ef6d 100644
--- a/examples/SodShock_3D/sodShock.yml
+++ b/examples/SodShock_3D/sodShock.yml
@@ -18,7 +18,8 @@ Snapshots:
basename: sodShock # 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)
-
+ compression: 4
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
diff --git a/examples/VacuumSpherical_3D/vacuum.yml b/examples/VacuumSpherical_3D/vacuum.yml
index 881b155b62c7f1f2af12a1d013ff5c05f1c16a88..92164a18a46404ad7730f3411dc14953139501bb 100644
--- a/examples/VacuumSpherical_3D/vacuum.yml
+++ b/examples/VacuumSpherical_3D/vacuum.yml
@@ -18,7 +18,8 @@ Snapshots:
basename: vacuum # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.05 # Time difference between consecutive outputs (in internal units)
-
+ compression: 4
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
diff --git a/examples/Vacuum_3D/vacuum.yml b/examples/Vacuum_3D/vacuum.yml
index 5ef5ce3da68febb086a14ad1a2207711f680d9ff..45a6b73d54de69d71e194ec074ebdd00cd6a57c0 100644
--- a/examples/Vacuum_3D/vacuum.yml
+++ b/examples/Vacuum_3D/vacuum.yml
@@ -18,7 +18,8 @@ Snapshots:
basename: vacuum # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
-
+ compression: 4
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
diff --git a/examples/main.c b/examples/main.c
index 828bf871bb034ce76a1a2454ed7e7dbb43e01443..80e4dfb479172d9e2bca6fb19ee676371d130afe 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -606,7 +606,7 @@ int main(int argc, char *argv[]) {
cosmology_init(params, &us, &prog_const, &cosmo);
else
cosmology_init_no_cosmo(&cosmo);
- if (with_cosmology) cosmology_print(&cosmo);
+ if (myrank == 0 && with_cosmology) cosmology_print(&cosmo);
/* Initialise the hydro properties */
if (with_hydro) hydro_props_init(&hydro_properties, params);
diff --git a/examples/plot_gravity_checks.py b/examples/plot_gravity_checks.py
index 41d9d629899f5c34cbb0264b38547a3439c04bf3..de4f37af32cf0d051afb4c5090075654e6fcd65c 100644
--- a/examples/plot_gravity_checks.py
+++ b/examples/plot_gravity_checks.py
@@ -13,7 +13,7 @@ params = {'axes.labelsize': 14,
'xtick.labelsize': 14,
'ytick.labelsize': 14,
'text.usetex': True,
-'figure.figsize': (10, 10),
+'figure.figsize': (12, 10),
'figure.subplot.left' : 0.06,
'figure.subplot.right' : 0.99 ,
'figure.subplot.bottom' : 0.06 ,
@@ -60,11 +60,13 @@ data = np.loadtxt('gravity_checks_exact_step%d.dat'%step)
exact_ids = data[:,0]
exact_pos = data[:,1:4]
exact_a = data[:,4:7]
+exact_pot = data[:,7]
# Sort stuff
sort_index = np.argsort(exact_ids)
exact_ids = exact_ids[sort_index]
exact_pos = exact_pos[sort_index, :]
exact_a = exact_a[sort_index, :]
+exact_pot = exact_pot[sort_index]
exact_a_norm = np.sqrt(exact_a[:,0]**2 + exact_a[:,1]**2 + exact_a[:,2]**2)
# Start the plot
@@ -166,12 +168,14 @@ for i in range(num_order):
ids = data[:,0]
pos = data[:,1:4]
a_grav = data[:, 4:7]
+ pot = data[:, 7]
# Sort stuff
sort_index = np.argsort(ids)
ids = ids[sort_index]
pos = pos[sort_index, :]
a_grav = a_grav[sort_index, :]
+ pot = pot[sort_index]
# Cross-checks
if not np.array_equal(exact_ids, ids):
@@ -185,6 +189,7 @@ for i in range(num_order):
# Compute the error norm
diff = exact_a - a_grav
+ diff_pot = exact_pot - pot
norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
@@ -192,73 +197,87 @@ for i in range(num_order):
error_x = abs(diff[:,0]) / exact_a_norm
error_y = abs(diff[:,1]) / exact_a_norm
error_z = abs(diff[:,2]) / exact_a_norm
+ error_pot = abs(diff_pot) / abs(exact_pot)
# Bin the error
norm_error_hist,_ = np.histogram(norm_error, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
error_x_hist,_ = np.histogram(error_x, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
error_y_hist,_ = np.histogram(error_y, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
error_z_hist,_ = np.histogram(error_z, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
+ error_pot_hist,_ = np.histogram(error_pot, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
norm_median = np.median(norm_error)
median_x = np.median(error_x)
median_y = np.median(error_y)
median_z = np.median(error_z)
+ median_pot = np.median(error_pot)
norm_per99 = np.percentile(norm_error,99)
per99_x = np.percentile(error_x,99)
per99_y = np.percentile(error_y,99)
per99_z = np.percentile(error_z,99)
+ per99_pot = np.percentile(error_pot, 99)
norm_max = np.max(norm_error)
max_x = np.max(error_x)
max_y = np.max(error_y)
max_z = np.max(error_z)
+ max_pot = np.max(error_pot)
print "Order %d ---- "%order[i]
print "Norm: median= %f 99%%= %f max= %f"%(norm_median, norm_per99, norm_max)
print "X : median= %f 99%%= %f max= %f"%(median_x, per99_x, max_x)
print "Y : median= %f 99%%= %f max= %f"%(median_y, per99_y, max_y)
print "Z : median= %f 99%%= %f max= %f"%(median_z, per99_z, max_z)
+ print "Pot : median= %f 99%%= %f max= %f"%(median_pot, per99_pot, max_pot)
print ""
- plt.subplot(221)
- plt.semilogx(bins, norm_error_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
- plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(norm_median, norm_per99), ha="left", va="top", color=cols[i])
- plt.subplot(222)
+ plt.subplot(231)
plt.semilogx(bins, error_x_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_x, per99_x), ha="left", va="top", color=cols[i])
- plt.subplot(223)
+ plt.subplot(232)
plt.semilogx(bins, error_y_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_y, per99_y), ha="left", va="top", color=cols[i])
- plt.subplot(224)
+ plt.subplot(233)
plt.semilogx(bins, error_z_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_z, per99_z), ha="left", va="top", color=cols[i])
+ plt.subplot(234)
+ plt.semilogx(bins, norm_error_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+ plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(norm_median, norm_per99), ha="left", va="top", color=cols[i])
+ plt.subplot(235)
+ plt.semilogx(bins, error_pot_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
+ plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_pot, per99_pot), ha="left", va="top", color=cols[i])
count += 1
-plt.subplot(221)
-plt.xlabel("$|\delta \overrightarrow{a}|/|\overrightarrow{a}_{exact}|$")
-#plt.ylabel("Density")
-plt.xlim(min_error, max_error)
-plt.ylim(0,2.5)
-plt.legend(loc="upper left")
-plt.subplot(222)
+plt.subplot(231)
plt.xlabel("$\delta a_x/|\overrightarrow{a}_{exact}|$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
#plt.legend(loc="center left")
-plt.subplot(223)
+plt.subplot(232)
plt.xlabel("$\delta a_y/|\overrightarrow{a}_{exact}|$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
#plt.legend(loc="center left")
-plt.subplot(224)
+plt.subplot(233)
plt.xlabel("$\delta a_z/|\overrightarrow{a}_{exact}|$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
+plt.subplot(234)
+plt.xlabel("$|\delta \overrightarrow{a}|/|\overrightarrow{a}_{exact}|$")
+#plt.ylabel("Density")
+plt.xlim(min_error, max_error)
+plt.ylim(0,2.5)
+plt.legend(loc="upper left")
+plt.subplot(235)
+plt.xlabel("$\delta \phi/\phi_{exact}$")
+#plt.ylabel("Density")
+plt.xlim(min_error, max_error)
+plt.ylim(0,1.75)
#plt.legend(loc="center left")
diff --git a/src/Makefile.am b/src/Makefile.am
index 6674dfcfaa11dcaa8b039426c4af66f825835adc..885a019ae4990466b924c03737df76b9472de02c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Add the debug flag to the whole thing
-AM_CFLAGS = $(HDF5_CPPFLAGS)
+AM_CFLAGS = $(HDF5_CPPFLAGS) $(GSL_INCS)
# Assign a "safe" version number
AM_LDFLAGS = $(HDF5_LDFLAGS) $(FFTW_LIBS) -version-info 0:0:0
@@ -25,7 +25,7 @@ AM_LDFLAGS = $(HDF5_LDFLAGS) $(FFTW_LIBS) -version-info 0:0:0
GIT_CMD = @GIT_CMD@
# Additional dependencies for shared libraries.
-EXTRA_LIBS = $(HDF5_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(GRACKLE_LIB)
+EXTRA_LIBS = $(HDF5_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(GRACKLE_LIB) $(GSL_LIBS)
# MPI libraries.
MPI_LIBS = $(PARMETIS_LIBS) $(MPI_THREAD_LIBS)
diff --git a/src/cosmology.c b/src/cosmology.c
index 9b301388e08ca0f4aa2a7d5d8f351ab88682c8ba..b1bf71ad25451d2185d911320fb757193e990945 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -667,13 +667,15 @@ void cosmology_struct_dump(const struct cosmology *cosmology, FILE *stream) {
* @brief Restore a cosmology struct from the given FILE as a stream of
* bytes.
*
+ * @param enabled whether cosmology is enabled.
* @param cosmology the struct
* @param stream the file stream
*/
-void cosmology_struct_restore(struct cosmology *cosmology, FILE *stream) {
+void cosmology_struct_restore(int enabled, struct cosmology *cosmology,
+ FILE *stream) {
restart_read_blocks((void *)cosmology, sizeof(struct cosmology), 1, stream,
NULL, "cosmology function");
- /* Re-initialise the tables */
- cosmology_init_tables(cosmology);
+ /* Re-initialise the tables if using a cosmology. */
+ if (enabled) cosmology_init_tables(cosmology);
}
diff --git a/src/cosmology.h b/src/cosmology.h
index b1b33930bad7a49c0f2abe3a2201d71c9be57305..b541e20fa1dd4f55422f9ff1d0952783d8de0ad7 100644
--- a/src/cosmology.h
+++ b/src/cosmology.h
@@ -192,6 +192,7 @@ void cosmology_write_model(hid_t h_grp, const struct cosmology *c);
/* Dump/restore. */
void cosmology_struct_dump(const struct cosmology *cosmology, FILE *stream);
-void cosmology_struct_restore(struct cosmology *cosmology, FILE *stream);
+void cosmology_struct_restore(int enabled, struct cosmology *cosmology,
+ FILE *stream);
#endif /* SWIFT_COSMOLOGY_H */
diff --git a/src/engine.c b/src/engine.c
index 01c1b54ba12885b21458e403e717c51d80e7e3be..e9badb8643cc6f728feac3b2a0bb493f7ce8bca5 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5821,11 +5821,13 @@ void engine_compute_next_snapshot_time(struct engine *e) {
if (e->policy & engine_policy_cosmology) {
const float next_snapshot_time =
exp(e->ti_nextSnapshot * e->time_base) * e->cosmology->a_begin;
- message("Next output time set to a=%e.", next_snapshot_time);
+ if (e->verbose)
+ message("Next output time set to a=%e.", next_snapshot_time);
} else {
const float next_snapshot_time =
e->ti_nextSnapshot * e->time_base + e->time_begin;
- message("Next output time set to t=%e.", next_snapshot_time);
+ if (e->verbose)
+ message("Next output time set to t=%e.", next_snapshot_time);
}
}
}
@@ -5924,7 +5926,7 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
struct cosmology *cosmo =
(struct cosmology *)malloc(sizeof(struct cosmology));
- cosmology_struct_restore(cosmo, stream);
+ cosmology_struct_restore(e->policy & engine_policy_cosmology, cosmo, stream);
e->cosmology = cosmo;
#ifdef WITH_MPI
diff --git a/src/gravity.c b/src/gravity.c
index 05f4f3724414287e5aeaa6e932ff4df7810914d9..0d951b4a9f8babb68d413e1183b1e2b52129abaf 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -54,6 +54,7 @@ struct exact_force_data {
static float fewald_x[Newald + 1][Newald + 1][Newald + 1];
static float fewald_y[Newald + 1][Newald + 1][Newald + 1];
static float fewald_z[Newald + 1][Newald + 1][Newald + 1];
+static float potewald[Newald + 1][Newald + 1][Newald + 1];
/* Factor used to normalize the access to the Ewald table */
float ewald_fac;
@@ -83,6 +84,8 @@ void gravity_exact_force_ewald_init(double boxSize) {
const float factor_exp1 = 2.f * alpha / sqrt(M_PI);
const float factor_exp2 = -M_PI * M_PI / alpha2;
const float factor_sin = 2.f * M_PI;
+ const float factor_cos = 2.f * M_PI;
+ const float factor_pot = M_PI / alpha2;
const float boxSize_inv2 = 1.f / (boxSize * boxSize);
/* Ewald factor to access the table */
@@ -92,6 +95,9 @@ void gravity_exact_force_ewald_init(double boxSize) {
bzero(fewald_x, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
bzero(fewald_y, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
bzero(fewald_z, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
+ bzero(potewald, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
+
+ potewald[0][0][0] = 2.8372975f;
/* Compute the values in one of the octants */
for (int i = 0; i <= Newald; ++i) {
@@ -114,6 +120,7 @@ void gravity_exact_force_ewald_init(double boxSize) {
float f_x = r_x * r_inv3;
float f_y = r_y * r_inv3;
float f_z = r_z * r_inv3;
+ float pot = r_inv + factor_pot;
for (int n_i = -4; n_i <= 4; ++n_i) {
for (int n_j = -4; n_j <= 4; ++n_j) {
@@ -130,15 +137,17 @@ void gravity_exact_force_ewald_init(double boxSize) {
const float r_tilde_inv3 =
r_tilde_inv * r_tilde_inv * r_tilde_inv;
- const float val =
- erfcf(alpha * r_tilde) +
- factor_exp1 * r_tilde * expf(-alpha2 * r_tilde2);
+ const float val_pot = erfcf(alpha * r_tilde);
+
+ const float val_f =
+ val_pot + factor_exp1 * r_tilde * expf(-alpha2 * r_tilde2);
/* First correction term */
- const float f = val * r_tilde_inv3;
+ const float f = val_f * r_tilde_inv3;
f_x -= f * d_x;
f_y -= f * d_y;
f_z -= f * d_z;
+ pot -= val_pot * r_tilde_inv;
}
}
}
@@ -149,16 +158,23 @@ void gravity_exact_force_ewald_init(double boxSize) {
const float h2 = h_i * h_i + h_j * h_j + h_k * h_k;
+ if (h2 == 0.f) continue;
+
const float h2_inv = 1.f / (h2 + FLT_MIN);
const float h_dot_x = h_i * r_x + h_j * r_y + h_k * r_z;
- const float val = 2.f * h2_inv * expf(h2 * factor_exp2) *
- sinf(factor_sin * h_dot_x);
+ const float common = h2_inv * expf(h2 * factor_exp2);
+
+ const float val_pot =
+ (float)M_1_PI * common * cosf(factor_cos * h_dot_x);
+
+ const float val_f = 2.f * common * sinf(factor_sin * h_dot_x);
/* Second correction term */
- f_x -= val * h_i;
- f_y -= val * h_j;
- f_z -= val * h_k;
+ f_x -= val_f * h_i;
+ f_y -= val_f * h_j;
+ f_z -= val_f * h_k;
+ pot -= val_pot;
}
}
}
@@ -167,6 +183,7 @@ void gravity_exact_force_ewald_init(double boxSize) {
fewald_x[i][j][k] = f_x;
fewald_y[i][j][k] = f_y;
fewald_z[i][j][k] = f_z;
+ potewald[i][j][k] = pot;
}
}
}
@@ -196,6 +213,11 @@ void gravity_exact_force_ewald_init(double boxSize) {
H5Dwrite(h_data, H5T_NATIVE_FLOAT, h_space, H5S_ALL, H5P_DEFAULT,
&(fewald_z[0][0][0]));
H5Dclose(h_data);
+ h_data = H5Dcreate(h_file, "Ewald_pot", H5T_NATIVE_FLOAT, h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ H5Dwrite(h_data, H5T_NATIVE_FLOAT, h_space, H5S_ALL, H5P_DEFAULT,
+ &(potewald[0][0][0]));
+ H5Dclose(h_data);
H5Sclose(h_space);
H5Fclose(h_file);
#endif
@@ -207,6 +229,7 @@ void gravity_exact_force_ewald_init(double boxSize) {
fewald_x[i][j][k] *= boxSize_inv2;
fewald_y[i][j][k] *= boxSize_inv2;
fewald_z[i][j][k] *= boxSize_inv2;
+ potewald[i][j][k] *= boxSize_inv2;
}
}
}
@@ -229,11 +252,12 @@ void gravity_exact_force_ewald_init(double boxSize) {
* @param rx x-coordinate of distance vector.
* @param ry y-coordinate of distance vector.
* @param rz z-coordinate of distance vector.
- * @param corr (return) The Ewald correction.
+ * @param corr_f (return) The Ewald correction for the force.
+ * @param corr_p (return) The Ewald correction for the potential.
*/
__attribute__((always_inline)) INLINE static void
gravity_exact_force_ewald_evaluate(double rx, double ry, double rz,
- double corr[3]) {
+ double corr_f[3], double *corr_p) {
const double s_x = (rx < 0.) ? 1. : -1.;
const double s_y = (ry < 0.) ? 1. : -1.;
@@ -258,40 +282,51 @@ gravity_exact_force_ewald_evaluate(double rx, double ry, double rz,
const double tz = 1. - dz;
/* Interpolation in X */
- corr[0] = 0.;
- corr[0] += fewald_x[i + 0][j + 0][k + 0] * tx * ty * tz;
- corr[0] += fewald_x[i + 0][j + 0][k + 1] * tx * ty * dz;
- corr[0] += fewald_x[i + 0][j + 1][k + 0] * tx * dy * tz;
- corr[0] += fewald_x[i + 0][j + 1][k + 1] * tx * dy * dz;
- corr[0] += fewald_x[i + 1][j + 0][k + 0] * dx * ty * tz;
- corr[0] += fewald_x[i + 1][j + 0][k + 1] * dx * ty * dz;
- corr[0] += fewald_x[i + 1][j + 1][k + 0] * dx * dy * tz;
- corr[0] += fewald_x[i + 1][j + 1][k + 1] * dx * dy * dz;
- corr[0] *= s_x;
+ corr_f[0] = 0.;
+ corr_f[0] += fewald_x[i + 0][j + 0][k + 0] * tx * ty * tz;
+ corr_f[0] += fewald_x[i + 0][j + 0][k + 1] * tx * ty * dz;
+ corr_f[0] += fewald_x[i + 0][j + 1][k + 0] * tx * dy * tz;
+ corr_f[0] += fewald_x[i + 0][j + 1][k + 1] * tx * dy * dz;
+ corr_f[0] += fewald_x[i + 1][j + 0][k + 0] * dx * ty * tz;
+ corr_f[0] += fewald_x[i + 1][j + 0][k + 1] * dx * ty * dz;
+ corr_f[0] += fewald_x[i + 1][j + 1][k + 0] * dx * dy * tz;
+ corr_f[0] += fewald_x[i + 1][j + 1][k + 1] * dx * dy * dz;
+ corr_f[0] *= s_x;
/* Interpolation in Y */
- corr[1] = 0.;
- corr[1] += fewald_y[i + 0][j + 0][k + 0] * tx * ty * tz;
- corr[1] += fewald_y[i + 0][j + 0][k + 1] * tx * ty * dz;
- corr[1] += fewald_y[i + 0][j + 1][k + 0] * tx * dy * tz;
- corr[1] += fewald_y[i + 0][j + 1][k + 1] * tx * dy * dz;
- corr[1] += fewald_y[i + 1][j + 0][k + 0] * dx * ty * tz;
- corr[1] += fewald_y[i + 1][j + 0][k + 1] * dx * ty * dz;
- corr[1] += fewald_y[i + 1][j + 1][k + 0] * dx * dy * tz;
- corr[1] += fewald_y[i + 1][j + 1][k + 1] * dx * dy * dz;
- corr[1] *= s_y;
+ corr_f[1] = 0.;
+ corr_f[1] += fewald_y[i + 0][j + 0][k + 0] * tx * ty * tz;
+ corr_f[1] += fewald_y[i + 0][j + 0][k + 1] * tx * ty * dz;
+ corr_f[1] += fewald_y[i + 0][j + 1][k + 0] * tx * dy * tz;
+ corr_f[1] += fewald_y[i + 0][j + 1][k + 1] * tx * dy * dz;
+ corr_f[1] += fewald_y[i + 1][j + 0][k + 0] * dx * ty * tz;
+ corr_f[1] += fewald_y[i + 1][j + 0][k + 1] * dx * ty * dz;
+ corr_f[1] += fewald_y[i + 1][j + 1][k + 0] * dx * dy * tz;
+ corr_f[1] += fewald_y[i + 1][j + 1][k + 1] * dx * dy * dz;
+ corr_f[1] *= s_y;
/* Interpolation in Z */
- corr[2] = 0.;
- corr[2] += fewald_z[i + 0][j + 0][k + 0] * tx * ty * tz;
- corr[2] += fewald_z[i + 0][j + 0][k + 1] * tx * ty * dz;
- corr[2] += fewald_z[i + 0][j + 1][k + 0] * tx * dy * tz;
- corr[2] += fewald_z[i + 0][j + 1][k + 1] * tx * dy * dz;
- corr[2] += fewald_z[i + 1][j + 0][k + 0] * dx * ty * tz;
- corr[2] += fewald_z[i + 1][j + 0][k + 1] * dx * ty * dz;
- corr[2] += fewald_z[i + 1][j + 1][k + 0] * dx * dy * tz;
- corr[2] += fewald_z[i + 1][j + 1][k + 1] * dx * dy * dz;
- corr[2] *= s_z;
+ corr_f[2] = 0.;
+ corr_f[2] += fewald_z[i + 0][j + 0][k + 0] * tx * ty * tz;
+ corr_f[2] += fewald_z[i + 0][j + 0][k + 1] * tx * ty * dz;
+ corr_f[2] += fewald_z[i + 0][j + 1][k + 0] * tx * dy * tz;
+ corr_f[2] += fewald_z[i + 0][j + 1][k + 1] * tx * dy * dz;
+ corr_f[2] += fewald_z[i + 1][j + 0][k + 0] * dx * ty * tz;
+ corr_f[2] += fewald_z[i + 1][j + 0][k + 1] * dx * ty * dz;
+ corr_f[2] += fewald_z[i + 1][j + 1][k + 0] * dx * dy * tz;
+ corr_f[2] += fewald_z[i + 1][j + 1][k + 1] * dx * dy * dz;
+ corr_f[2] *= s_z;
+
+ /* Interpolation for potential */
+ *corr_p = 0.;
+ *corr_p += potewald[i + 0][j + 0][k + 0] * tx * ty * tz;
+ *corr_p += potewald[i + 0][j + 0][k + 1] * tx * ty * dz;
+ *corr_p += potewald[i + 0][j + 1][k + 0] * tx * dy * tz;
+ *corr_p += potewald[i + 0][j + 1][k + 1] * tx * dy * dz;
+ *corr_p += potewald[i + 1][j + 0][k + 0] * dx * ty * tz;
+ *corr_p += potewald[i + 1][j + 0][k + 1] * dx * ty * dz;
+ *corr_p += potewald[i + 1][j + 1][k + 0] * dx * dy * tz;
+ *corr_p += potewald[i + 1][j + 1][k + 1] * dx * dy * dz;
}
#endif
@@ -383,6 +418,7 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
/* Be ready for the calculation */
double a_grav[3] = {0., 0., 0.};
+ double pot = 0.;
/* Interact it with all other particles in the space.*/
for (int j = 0; j < (int)s->nr_gparts; ++j) {
@@ -408,37 +444,42 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
const double r_inv = 1. / sqrt(r2);
const double r = r2 * r_inv;
const double mj = gpj->mass;
- double f;
+ double f, phi;
if (r >= hi) {
/* Get Newtonian gravity */
f = mj * r_inv * r_inv * r_inv;
+ phi = -mj * r_inv;
} else {
const double ui = r * hi_inv;
- double W;
+ double Wf, Wp;
- kernel_grav_eval_double(ui, &W);
+ kernel_grav_eval_force_double(ui, &Wf);
+ kernel_grav_eval_pot_double(ui, &Wp);
/* Get softened gravity */
- f = mj * hi_inv3 * W;
+ f = mj * hi_inv3 * Wf;
+ phi = mj * hi_inv * Wp;
}
a_grav[0] += f * dx;
a_grav[1] += f * dy;
a_grav[2] += f * dz;
+ pot += phi;
/* Apply Ewald correction for periodic BC */
if (periodic && r > 1e-5 * hi) {
- double corr[3];
- gravity_exact_force_ewald_evaluate(dx, dy, dz, corr);
+ double corr_f[3], corr_pot;
+ gravity_exact_force_ewald_evaluate(dx, dy, dz, corr_f, &corr_pot);
- a_grav[0] += mj * corr[0];
- a_grav[1] += mj * corr[1];
- a_grav[2] += mj * corr[2];
+ a_grav[0] += mj * corr_f[0];
+ a_grav[1] += mj * corr_f[1];
+ a_grav[2] += mj * corr_f[2];
+ pot += mj * corr_pot;
}
}
@@ -446,6 +487,7 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
gpi->a_grav_exact[0] = a_grav[0] * const_G;
gpi->a_grav_exact[1] = a_grav[1] * const_G;
gpi->a_grav_exact[2] = a_grav[2] * const_G;
+ gpi->potential_exact = pot * const_G;
counter++;
}
@@ -529,8 +571,9 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
fprintf(file_swift, "# theta= %16.8e\n", e->gravity_properties->theta_crit);
fprintf(file_swift, "# Git Branch: %s\n", git_branch());
fprintf(file_swift, "# Git Revision: %s\n", git_revision());
- fprintf(file_swift, "# %16s %16s %16s %16s %16s %16s %16s\n", "id", "pos[0]",
- "pos[1]", "pos[2]", "a_swift[0]", "a_swift[1]", "a_swift[2]");
+ fprintf(file_swift, "# %16s %16s %16s %16s %16s %16s %16s %16s\n", "id",
+ "pos[0]", "pos[1]", "pos[2]", "a_swift[0]", "a_swift[1]",
+ "a_swift[2]", "potential");
/* Output particle SWIFT accelerations */
for (size_t i = 0; i < s->nr_gparts; ++i) {
@@ -541,9 +584,10 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
gpart_is_starting(gpi, e)) {
- fprintf(file_swift, "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e \n",
+ fprintf(file_swift,
+ "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
- gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2]);
+ gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2], gpi->potential);
}
}
@@ -569,9 +613,9 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
fprintf(file_exact, "# periodic= %d\n", s->periodic);
fprintf(file_exact, "# Git Branch: %s\n", git_branch());
fprintf(file_exact, "# Git Revision: %s\n", git_revision());
- fprintf(file_exact, "# %16s %16s %16s %16s %16s %16s %16s\n", "id",
+ fprintf(file_exact, "# %16s %16s %16s %16s %16s %16s %16s %16s\n", "id",
"pos[0]", "pos[1]", "pos[2]", "a_exact[0]", "a_exact[1]",
- "a_exact[2]");
+ "a_exact[2]", "potential");
/* Output particle exact accelerations */
for (size_t i = 0; i < s->nr_gparts; ++i) {
@@ -582,10 +626,11 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
gpart_is_starting(gpi, e)) {
- fprintf(
- file_exact, "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e \n",
- gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
- gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2]);
+ fprintf(file_exact,
+ "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
+ gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
+ gpi->a_grav_exact[0], gpi->a_grav_exact[1],
+ gpi->a_grav_exact[2], gpi->potential_exact);
}
}
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index bc06b17188f99785404a1162d46384ee26c884ea..099cbcd158436b414b9bc6a912727116b763cff9 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -48,16 +48,29 @@ __attribute__((always_inline)) INLINE static float gravity_get_softening(
}
/**
- * @brief Returns the potential of a particle
+ * @brief Returns the comoving potential of a particle
*
* @param gp The particle of interest
*/
-__attribute__((always_inline)) INLINE static float gravity_get_potential(
- const struct gpart* restrict gp) {
+__attribute__((always_inline)) INLINE static float
+gravity_get_comoving_potential(const struct gpart* restrict gp) {
return gp->potential;
}
+/**
+ * @brief Returns the physical potential of a particle
+ *
+ * @param gp The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_get_physical_potential(const struct gpart* restrict gp,
+ const struct cosmology* cosmo) {
+
+ return gp->potential * cosmo->a_inv;
+}
+
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
@@ -136,6 +149,7 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
gp->a_grav[0] *= const_G;
gp->a_grav[1] *= const_G;
gp->a_grav[2] *= const_G;
+ gp->potential *= const_G;
}
/**
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 1f4310eea49a1b1b298fe00f00e1ae0802d9b688..3951f77d59e71439d2f621dd2a69fd62d39ba4f0 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -38,9 +38,11 @@
* @param h_inv3 Cube of the inverse of the softening length.
* @param mass Mass of the point-mass.
* @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
- float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij) {
+ float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij,
+ float *pot_ij) {
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2);
@@ -50,17 +52,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
/* Get Newtonian gravity */
*f_ij = mass * r_inv * r_inv * r_inv;
+ *pot_ij = -mass * r_inv;
} else {
const float r = r2 * r_inv;
const float ui = r * h_inv;
- float W_ij;
- kernel_grav_eval(ui, &W_ij);
+ float W_f_ij, W_pot_ij;
+ kernel_grav_force_eval(ui, &W_f_ij);
+ kernel_grav_pot_eval(ui, &W_pot_ij);
/* Get softened gravity */
- *f_ij = mass * h_inv3 * W_ij;
+ *f_ij = mass * h_inv3 * W_f_ij;
+ *pot_ij = mass * h_inv * W_pot_ij;
}
}
@@ -78,10 +83,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
* @param mass Mass of the point-mass.
* @param rlr_inv Inverse of the mesh smoothing scale.
* @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
float r2, float h2, float h_inv, float h_inv3, float mass, float rlr_inv,
- float *f_ij) {
+ float *f_ij, float *pot_ij) {
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2);
@@ -92,23 +98,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
/* Get Newtonian gravity */
*f_ij = mass * r_inv * r_inv * r_inv;
+ *pot_ij = -mass * r_inv;
} else {
const float ui = r * h_inv;
- float W_ij;
+ float W_f_ij, W_pot_ij;
- kernel_grav_eval(ui, &W_ij);
+ kernel_grav_force_eval(ui, &W_f_ij);
+ kernel_grav_pot_eval(ui, &W_pot_ij);
/* Get softened gravity */
- *f_ij = mass * h_inv3 * W_ij;
+ *f_ij = mass * h_inv3 * W_f_ij;
+ *pot_ij = mass * h_inv * W_pot_ij;
}
/* Get long-range correction */
const float u_lr = r * rlr_inv;
- float corr_lr;
- kernel_long_grav_eval(u_lr, &corr_lr);
- *f_ij *= corr_lr;
+ float corr_f_lr, corr_pot_lr;
+ kernel_long_grav_force_eval(u_lr, &corr_f_lr);
+ kernel_long_grav_pot_eval(u_lr, &corr_pot_lr);
+ *f_ij *= corr_f_lr;
+ *pot_ij *= corr_pot_lr;
}
/**
@@ -127,33 +138,43 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
* @param f_x (return) The x-component of the acceleration.
* @param f_y (return) The y-component of the acceleration.
* @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
float r_x, float r_y, float r_z, float r2, float h, float h_inv,
- const struct multipole *m, float *f_x, float *f_y, float *f_z) {
+ const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
- runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv3, m->M_000, f_ij);
+ float f_ij;
+ runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
+ &f_ij, pot);
+ *f_x = -f_ij * r_x;
+ *f_y = -f_ij * r_y;
+ *f_z = -f_ij * r_z;
#else
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2);
- struct potential_derivatives_M2P pot;
- compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, &pot);
+ struct potential_derivatives_M2P d;
+ compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, &d);
/* 1st order terms (monopole) */
- *f_x = m->M_000 * pot.D_100;
- *f_y = m->M_000 * pot.D_010;
- *f_z = m->M_000 * pot.D_001;
+ *f_x = m->M_000 * d.D_100;
+ *f_y = m->M_000 * d.D_010;
+ *f_z = m->M_000 * d.D_001;
+ *pot = -m->M_000 * d.D_000;
/* 3rd order terms (quadrupole) */
- *f_x += m->M_200 * pot.D_300 + m->M_020 * pot.D_120 + m->M_002 * pot.D_102;
- *f_y += m->M_200 * pot.D_210 + m->M_020 * pot.D_030 + m->M_002 * pot.D_012;
- *f_z += m->M_200 * pot.D_201 + m->M_020 * pot.D_021 + m->M_002 * pot.D_003;
- *f_x += m->M_110 * pot.D_210 + m->M_101 * pot.D_201 + m->M_011 * pot.D_111;
- *f_y += m->M_110 * pot.D_120 + m->M_101 * pot.D_111 + m->M_011 * pot.D_021;
- *f_z += m->M_110 * pot.D_111 + m->M_101 * pot.D_102 + m->M_011 * pot.D_012;
+ *f_x += m->M_200 * d.D_300 + m->M_020 * d.D_120 + m->M_002 * d.D_102;
+ *f_y += m->M_200 * d.D_210 + m->M_020 * d.D_030 + m->M_002 * d.D_012;
+ *f_z += m->M_200 * d.D_201 + m->M_020 * d.D_021 + m->M_002 * d.D_003;
+ *pot -= m->M_200 * d.D_100 + m->M_020 * d.D_020 + m->M_002 * d.D_002;
+
+ *f_x += m->M_110 * d.D_210 + m->M_101 * d.D_201 + m->M_011 * d.D_111;
+ *f_y += m->M_110 * d.D_120 + m->M_101 * d.D_111 + m->M_011 * d.D_021;
+ *f_z += m->M_110 * d.D_111 + m->M_101 * d.D_102 + m->M_011 * d.D_012;
+ *pot -= m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
#endif
}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 33c444e24e4bc9681dc3138298b5fd9301d084d3..4c7a7b5d4af1c26ee96ac797963f5588d3303d77 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -71,6 +71,8 @@ struct gpart {
/* Brute-force particle acceleration. */
double a_grav_exact[3];
+ /* Brute-force particle potential. */
+ double potential_exact;
#endif
} SWIFT_STRUCT_ALIGN;
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index d46dc3cb1f307aee0b5c7352ea3c5b27b640cfdc..be349cfd4a4583a493219eb020757341c98f540c 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -59,6 +59,9 @@ struct gravity_cache {
/*! #gpart z acceleration. */
float *restrict a_z SWIFT_CACHE_ALIGN;
+ /*! #gpart potential. */
+ float *restrict pot SWIFT_CACHE_ALIGN;
+
/*! Is this #gpart active ? */
int *restrict active SWIFT_CACHE_ALIGN;
@@ -85,6 +88,7 @@ static INLINE void gravity_cache_clean(struct gravity_cache *c) {
free(c->a_x);
free(c->a_y);
free(c->a_z);
+ free(c->pot);
free(c->active);
free(c->use_mpole);
}
@@ -120,6 +124,7 @@ static INLINE void gravity_cache_init(struct gravity_cache *c, int count) {
e += posix_memalign((void **)&c->a_x, SWIFT_CACHE_ALIGNMENT, sizeBytesF);
e += posix_memalign((void **)&c->a_y, SWIFT_CACHE_ALIGNMENT, sizeBytesF);
e += posix_memalign((void **)&c->a_z, SWIFT_CACHE_ALIGNMENT, sizeBytesF);
+ e += posix_memalign((void **)&c->pot, SWIFT_CACHE_ALIGNMENT, sizeBytesF);
e += posix_memalign((void **)&c->active, SWIFT_CACHE_ALIGNMENT, sizeBytesI);
e +=
posix_memalign((void **)&c->use_mpole, SWIFT_CACHE_ALIGNMENT, sizeBytesI);
@@ -278,6 +283,7 @@ __attribute__((always_inline)) INLINE void gravity_cache_write_back(
swift_declare_aligned_ptr(float, a_x, c->a_x, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, a_y, c->a_y, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, a_z, c->a_z, SWIFT_CACHE_ALIGNMENT);
+ swift_declare_aligned_ptr(float, pot, c->pot, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(int, active, c->active, SWIFT_CACHE_ALIGNMENT);
/* Write stuff back to the particles */
@@ -286,6 +292,7 @@ __attribute__((always_inline)) INLINE void gravity_cache_write_back(
gparts[i].a_grav[0] += a_x[i];
gparts[i].a_grav[1] += a_y[i];
gparts[i].a_grav[2] += a_z[i];
+ gparts[i].potential += pot[i];
}
}
}
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index cf8aa54338b2e87e8bf5f2cc453ad7417eea5804..d698c66e418a120e5e7ebbe1cba0a5a9af8cd1f1 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -95,9 +95,16 @@ struct potential_derivatives_M2L {
*/
struct potential_derivatives_M2P {
+ /* 0th order terms */
+ float D_000;
+
/* 1st order terms */
float D_100, D_010, D_001;
+ /* 2nd order terms */
+ float D_200, D_020, D_002;
+ float D_110, D_101, D_011;
+
/* 3rd order terms */
float D_300, D_030, D_003;
float D_210, D_201;
@@ -368,11 +375,22 @@ compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2,
const float r_y3 = r_y2 * r_y;
const float r_z3 = r_z2 * r_z;
+ /* 0th order derivative */
+ pot->D_000 = Dt_1;
+
/* 1st order derivatives */
pot->D_100 = r_x * Dt_3;
pot->D_010 = r_y * Dt_3;
pot->D_001 = r_z * Dt_3;
+ /* 2nd order derivatives */
+ pot->D_200 = r_x2 * Dt_5 + Dt_3;
+ pot->D_020 = r_y2 * Dt_5 + Dt_3;
+ pot->D_002 = r_z2 * Dt_5 + Dt_3;
+ pot->D_110 = r_x * r_y * Dt_5;
+ pot->D_101 = r_x * r_z * Dt_5;
+ pot->D_011 = r_y * r_z * Dt_5;
+
/* 3rd order derivatives */
pot->D_300 = r_x3 * Dt_7 + 3.f * r_x * Dt_5;
pot->D_030 = r_y3 * Dt_7 + 3.f * r_y * Dt_5;
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index 2567545b25f55e03b104844a4142da9cbe397fb5..bfff08dcdf061d86986ffdeabee6e7a84cba1c7b 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -77,7 +77,7 @@ void convert_part_pos(const struct engine* e, const struct part* p,
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 8;
+ *num_fields = 7;
/* List what we want to write */
list[0] = io_make_output_field_convert_part(
@@ -92,9 +92,7 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
- list[6] = io_make_output_field("Acceleration", FLOAT, 3,
- UNIT_CONV_ACCELERATION, parts, a_hydro);
- list[7] =
+ list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
}
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 90484c9cee1b384b0540f56e3c845cc4120cdfe0..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, 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
@@ -687,8 +687,11 @@ runner_iact_nonsym_1_vec_force(
const vector balsara_j = vector_load(Balsara_j);
const vector cj = vector_load(Cj);
- const vector fac_mu =
- vector_set1(1.f); /* Will change with cosmological integration */
+ /* Cosmological terms */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+ const vector v_fac_mu = vector_set1(fac_mu);
+ const vector v_a2_Hubble = vector_set1(a2_Hubble);
/* Load stuff. */
balsara.v = vec_add(balsara_i.v, balsara_j.v);
@@ -718,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_mul(dvz.v, dz->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(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));
@@ -787,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, 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
@@ -853,8 +860,11 @@ runner_iact_nonsym_2_vec_force(
const vector hj_inv = vector_load(Hj_inv);
const vector hj_inv_2 = vector_load(&Hj_inv[VEC_SIZE]);
- const vector fac_mu =
- vector_set1(1.f); /* Will change with cosmological integration */
+ /* Cosmological terms */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+ const vector v_fac_mu = vector_set1(fac_mu);
+ const vector v_a2_Hubble = vector_set1(a2_Hubble);
/* Find the balsara switch. */
balsara.v = vec_add(balsara_i.v, balsara_j.v);
@@ -891,18 +901,22 @@ 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_mul(dvz.v, dz.v)));
- dvdr_2.v = vec_fma(dvx_2.v, dx_2.v,
- vec_fma(dvy_2.v, dy_2.v, vec_mul(dvz_2.v, dz_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(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(
- fac_mu.v, vec_mul(ri_2.v, omega_ij_2.v)); /* This is 0 or negative */
+ v_fac_mu.v, vec_mul(ri_2.v, omega_ij_2.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));
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 1d1f5ed453f249191b6853a3a699de847ed501f2..5c329e9d7297109c378e04f0e4931a0717e07278 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -89,7 +89,7 @@ void convert_part_pos(const struct engine* e, const struct part* p,
void hydro_write_particles(const struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 10;
+ *num_fields = 9;
/* List what we want to write */
list[0] = io_make_output_field_convert_part(
@@ -104,14 +104,12 @@ void hydro_write_particles(const struct part* parts, struct io_props* list,
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
- list[6] = io_make_output_field("Acceleration", FLOAT, 3,
- UNIT_CONV_ACCELERATION, parts, a_hydro);
- list[7] =
+ list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
- list[8] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
+ list[7] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
parts, convert_u);
- list[9] = io_make_output_field_convert_part(
+ list[8] = io_make_output_field_convert_part(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
#ifdef DEBUG_INTERACTIONS_SPH
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 5922834e12f6dd74e951b46c2c93f975d7e8b82e..2ff70fd2f7e95344c17ba0f3237c3b2ee131752c 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -103,7 +103,7 @@ void convert_part_pos(const struct engine* e, const struct part* p,
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 10;
+ *num_fields = 9;
/* List what we want to write */
list[0] = io_make_output_field_convert_part(
@@ -118,13 +118,11 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
- list[6] = io_make_output_field("Acceleration", FLOAT, 3,
- UNIT_CONV_ACCELERATION, parts, a_hydro);
- list[7] =
+ list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
- list[8] = io_make_output_field_convert_part(
+ list[7] = io_make_output_field_convert_part(
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, convert_S);
- list[9] = io_make_output_field_convert_part(
+ list[8] = io_make_output_field_convert_part(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
}
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index 454d999b6dcd604e25da4c513d902854c43c781b..ed8e0e45d240ea37ebfb7df6339afd454dc0e4c4 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -101,7 +101,7 @@ void convert_part_pos(const struct engine* e, const struct part* p,
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 11;
+ *num_fields = 10;
/* List what we want to write */
list[0] = io_make_output_field_convert_part(
@@ -116,17 +116,15 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
- list[6] = io_make_output_field("Acceleration", FLOAT, 3,
- UNIT_CONV_ACCELERATION, parts, a_hydro);
- list[7] =
+ list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
- list[8] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
+ list[7] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
parts, convert_u);
- list[9] = io_make_output_field_convert_part(
+ list[8] = io_make_output_field_convert_part(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
- list[10] = io_make_output_field("WeightedDensity", FLOAT, 1,
- UNIT_CONV_DENSITY, parts, rho_bar);
+ list[9] = io_make_output_field("WeightedDensity", FLOAT, 1, UNIT_CONV_DENSITY,
+ parts, rho_bar);
}
/**
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index 799bda85b0c69dd2757f47fb0225006adb6d1432..24d30136c8213dab6db424f2f459766c45ec6b1e 100644
--- a/src/kernel_gravity.h
+++ b/src/kernel_gravity.h
@@ -27,14 +27,35 @@
#include "minmax.h"
/**
- * @brief Computes the gravity softening function.
+ * @brief Computes the gravity softening function for potential.
*
* This functions assumes 0 < u < 1.
*
* @param u The ratio of the distance to the softening length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
*/
-__attribute__((always_inline)) INLINE static void kernel_grav_eval(
+__attribute__((always_inline)) INLINE static void kernel_grav_pot_eval(
+ float u, float *const W) {
+
+ /* W(u) = 3u^7 - 15u^6 + 28u^5 - 21u^4 + 7u^2 - 3 */
+ *W = 3.f * u - 15.f;
+ *W = *W * u + 28.f;
+ *W = *W * u - 21.f;
+ *W = *W * u;
+ *W = *W * u + 7.f;
+ *W = *W * u;
+ *W = *W * u - 3.f;
+}
+
+/**
+ * @brief Computes the gravity softening function for forces.
+ *
+ * This functions assumes 0 < u < 1.
+ *
+ * @param u The ratio of the distance to the softening length $u = x/h$.
+ * @param W (return) The value of the kernel function $W(x,h)$.
+ */
+__attribute__((always_inline)) INLINE static void kernel_grav_force_eval(
float u, float *const W) {
/* W(u) = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
@@ -48,14 +69,37 @@ __attribute__((always_inline)) INLINE static void kernel_grav_eval(
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/**
- * @brief Computes the gravity softening function in double precision.
+ * @brief Computes the gravity softening function for potential in double
+ * precision.
+ *
+ * This functions assumes 0 < u < 1.
+ *
+ * @param u The ratio of the distance to the softening length $u = x/h$.
+ * @param W (return) The value of the kernel function $W(x,h)$.
+ */
+__attribute__((always_inline)) INLINE static void kernel_grav_eval_pot_double(
+ double u, double *const W) {
+
+ /* W(u) = 3u^7 - 15u^6 + 28u^5 - 21u^4 + 7u^2 - 3 */
+ *W = 3. * u - 15.;
+ *W = *W * u + 28.;
+ *W = *W * u - 21.;
+ *W = *W * u;
+ *W = *W * u + 7.;
+ *W = *W * u;
+ *W = *W * u - 3;
+}
+
+/**
+ * @brief Computes the gravity softening function for forces in double
+ * precision.
*
* This functions assumes 0 < u < 1.
*
* @param u The ratio of the distance to the softening length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
*/
-__attribute__((always_inline)) INLINE static void kernel_grav_eval_double(
+__attribute__((always_inline)) INLINE static void kernel_grav_eval_force_double(
double u, double *const W) {
/* W(u) = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index 2e0097ded217769af832af049b6c00b1305ac339..77bc7c17e2ab8038a4860c74daebd19a2d785070 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -31,12 +31,41 @@
#include <math.h>
/**
- * @brief Computes the long-range correction term for the FFT calculation.
+ * @brief Computes the long-range correction term for the potential calculation
+ * coming from FFT.
*
* @param u The ratio of the distance to the FFT cell scale \f$u = r/r_s\f$.
* @param W (return) The value of the kernel function.
*/
-__attribute__((always_inline)) INLINE static void kernel_long_grav_eval(
+__attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval(
+ const float u, float *const W) {
+
+#ifdef GADGET2_LONG_RANGE_CORRECTION
+
+ const float arg1 = u * 0.5f;
+ const float term1 = erfcf(arg1);
+
+ *W = term1;
+#else
+
+ const float x = 2.f * u;
+ const float exp_x = expf(x);
+ const float alpha = 1.f / (1.f + exp_x);
+
+ /* We want 2 - 2 exp(x) * alpha */
+ *W = 1.f - alpha * exp_x;
+ *W *= 2.f;
+#endif
+}
+
+/**
+ * @brief Computes the long-range correction term for the force calculation
+ * coming from FFT.
+ *
+ * @param u The ratio of the distance to the FFT cell scale \f$u = r/r_s\f$.
+ * @param W (return) The value of the kernel function.
+ */
+__attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
float u, float *const W) {
#ifdef GADGET2_LONG_RANGE_CORRECTION
@@ -54,7 +83,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_eval(
#else
const float arg = 2.f * u;
- const float exp_arg = good_approx_expf(arg);
+ const float exp_arg = expf(arg);
const float term = 1.f / (1.f + exp_arg);
*W = arg * exp_arg * term * term - exp_arg * term + 1.f;
diff --git a/src/multipole.h b/src/multipole.h
index d842081814b6ef7b6f490a1ba47a0152dd84d5a1..852976ad552b3ccf16f3e084612166c9212df2a6 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -2263,30 +2263,38 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
/* Local accumulator */
double a_grav[3] = {0., 0., 0.};
+ double pot = 0.;
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* Distance to the multipole */
const double dx[3] = {gp->x[0] - loc[0], gp->x[1] - loc[1],
gp->x[2] - loc[2]};
+ /* 0th order contributions */
+ pot -= X_000(dx) * lb->F_000;
- /* 0th order interaction */
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ /* 1st order contributions */
a_grav[0] += X_000(dx) * lb->F_100;
a_grav[1] += X_000(dx) * lb->F_010;
a_grav[2] += X_000(dx) * lb->F_001;
+
+ pot -= X_001(dx) * lb->F_001 + X_010(dx) * lb->F_010 + X_100(dx) * lb->F_100;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
- /* 1st order interaction */
+ /* 2nd order contributions */
a_grav[0] +=
X_100(dx) * lb->F_200 + X_010(dx) * lb->F_110 + X_001(dx) * lb->F_101;
a_grav[1] +=
X_100(dx) * lb->F_110 + X_010(dx) * lb->F_020 + X_001(dx) * lb->F_011;
a_grav[2] +=
X_100(dx) * lb->F_101 + X_010(dx) * lb->F_011 + X_001(dx) * lb->F_002;
+
+ pot -= X_002(dx) * lb->F_002 + X_011(dx) * lb->F_011 + X_020(dx) * lb->F_020 +
+ X_101(dx) * lb->F_101 + X_110(dx) * lb->F_110 + X_200(dx) * lb->F_200;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
- /* 2nd order interaction */
+ /* 3rd order contributions */
a_grav[0] +=
X_200(dx) * lb->F_300 + X_020(dx) * lb->F_120 + X_002(dx) * lb->F_102;
a_grav[0] +=
@@ -2299,10 +2307,15 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
X_200(dx) * lb->F_201 + X_020(dx) * lb->F_021 + X_002(dx) * lb->F_003;
a_grav[2] +=
X_110(dx) * lb->F_111 + X_101(dx) * lb->F_102 + X_011(dx) * lb->F_012;
+
+ pot -= X_003(dx) * lb->F_003 + X_012(dx) * lb->F_012 + X_021(dx) * lb->F_021 +
+ X_030(dx) * lb->F_030 + X_102(dx) * lb->F_102 + X_111(dx) * lb->F_111 +
+ X_120(dx) * lb->F_120 + X_201(dx) * lb->F_201 + X_210(dx) * lb->F_210 +
+ X_300(dx) * lb->F_300;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
- /* 3rd order contributions */
+ /* 4th order contributions */
a_grav[0] += X_003(dx) * lb->F_103 + X_012(dx) * lb->F_112 +
X_021(dx) * lb->F_121 + X_030(dx) * lb->F_130 +
X_102(dx) * lb->F_202 + X_111(dx) * lb->F_211 +
@@ -2319,10 +2332,16 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
X_120(dx) * lb->F_121 + X_201(dx) * lb->F_202 +
X_210(dx) * lb->F_211 + X_300(dx) * lb->F_301;
+ pot -= X_004(dx) * lb->F_004 + X_013(dx) * lb->F_013 + X_022(dx) * lb->F_022 +
+ X_031(dx) * lb->F_031 + X_040(dx) * lb->F_040 + X_103(dx) * lb->F_103 +
+ X_112(dx) * lb->F_112 + X_121(dx) * lb->F_121 + X_130(dx) * lb->F_130 +
+ X_202(dx) * lb->F_202 + X_211(dx) * lb->F_211 + X_220(dx) * lb->F_220 +
+ X_301(dx) * lb->F_301 + X_310(dx) * lb->F_310 + X_400(dx) * lb->F_400;
+
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
- /* 4th order contributions */
+ /* 5th order contributions */
a_grav[0] +=
X_004(dx) * lb->F_104 + X_013(dx) * lb->F_113 + X_022(dx) * lb->F_122 +
X_031(dx) * lb->F_131 + X_040(dx) * lb->F_140 + X_103(dx) * lb->F_203 +
@@ -2342,6 +2361,14 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
X_202(dx) * lb->F_203 + X_211(dx) * lb->F_212 + X_220(dx) * lb->F_221 +
X_301(dx) * lb->F_302 + X_310(dx) * lb->F_311 + X_400(dx) * lb->F_401;
+ pot -= X_005(dx) * lb->F_005 + X_014(dx) * lb->F_014 + X_023(dx) * lb->F_023 +
+ X_032(dx) * lb->F_032 + X_041(dx) * lb->F_041 + X_050(dx) * lb->F_050 +
+ X_104(dx) * lb->F_104 + X_113(dx) * lb->F_113 + X_122(dx) * lb->F_122 +
+ X_131(dx) * lb->F_131 + X_140(dx) * lb->F_140 + X_203(dx) * lb->F_203 +
+ X_212(dx) * lb->F_212 + X_221(dx) * lb->F_221 + X_230(dx) * lb->F_230 +
+ X_302(dx) * lb->F_302 + X_311(dx) * lb->F_311 + X_320(dx) * lb->F_320 +
+ X_401(dx) * lb->F_401 + X_410(dx) * lb->F_410 + X_500(dx) * lb->F_500;
+
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
@@ -2351,52 +2378,8 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
gp->a_grav[0] += a_grav[0];
gp->a_grav[1] += a_grav[1];
gp->a_grav[2] += a_grav[2];
+ gp->potential += pot;
}
-
-INLINE static void gravity_M2P(const struct multipole *ma,
- const struct gravity_props *props,
- const double loc[3], struct gpart *gp) {
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
-
- const float eps2 = props->epsilon2;
- const float eps_inv = props->epsilon_inv;
- const float eps_inv3 = props->epsilon_inv3;
-
- /* Distance to the multipole */
- const float dx = gp->x[0] - loc[0];
- const float dy = gp->x[1] - loc[1];
- const float dz = gp->x[2] - loc[2];
- const float r2 = dx * dx + dy * dy + dz * dz;
-
- /* Get the inverse distance */
- const float r_inv = 1.f / sqrtf(r2);
-
- float f, W;
-
- if (r2 >= eps2) {
-
- /* Get Newtonian gravity */
- f = ma->M_000 * r_inv * r_inv * r_inv;
-
- } else {
-
- const float r = r2 * r_inv;
- const float u = r * eps_inv;
-
- kernel_grav_eval(u, &W);
-
- /* Get softened gravity */
- f = ma->M_000 * eps_inv3 * W;
- }
-
- gp->a_grav[0] -= f * dx;
- gp->a_grav[1] -= f * dy;
- gp->a_grav[2] -= f * dz;
-
-#endif
-}
-
/**
* @brief Checks whether a cell-cell interaction can be appromixated by a M-M
* interaction using the distance and cell radius.
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.h b/src/runner_doiact.h
index 2ea4cbf24576e50607379d8659106a57a1f80356..2987d205e5005d4618c087bfe7f18a0b3c12fef5 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -1161,7 +1161,8 @@ void DOPAIR1_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
p->x[1] * runner_shift[sid][1] +
p->x[2] * runner_shift[sid][2];
if (fabsf(d - sort_i[pid].d) - ci->dx_max_sort >
- 1.0e-4 * max(fabsf(d), ci->dx_max_sort_old))
+ 1.0e-4 * max(fabsf(d), ci->dx_max_sort_old) &&
+ fabsf(d - sort_i[pid].d) - ci->dx_max_sort > ci->width[0] * 1.0e-10)
error(
"particle shift diff exceeds dx_max_sort in cell ci. ci->nodeID=%d "
"cj->nodeID=%d d=%e sort_i[pid].d=%e ci->dx_max_sort=%e "
@@ -1174,8 +1175,9 @@ void DOPAIR1_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
const float d = p->x[0] * runner_shift[sid][0] +
p->x[1] * runner_shift[sid][1] +
p->x[2] * runner_shift[sid][2];
- if (fabsf(d - sort_j[pjd].d) - cj->dx_max_sort >
- 1.0e-4 * max(fabsf(d), cj->dx_max_sort_old))
+ if ((fabsf(d - sort_j[pjd].d) - cj->dx_max_sort) >
+ 1.0e-4 * max(fabsf(d), cj->dx_max_sort_old) &&
+ (fabsf(d - sort_j[pjd].d) - cj->dx_max_sort) > cj->width[0] * 1.0e-10)
error(
"particle shift diff exceeds dx_max_sort in cell cj. cj->nodeID=%d "
"ci->nodeID=%d d=%e sort_j[pjd].d=%e cj->dx_max_sort=%e "
@@ -1672,7 +1674,8 @@ void DOPAIR2_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
p->x[1] * runner_shift[sid][1] +
p->x[2] * runner_shift[sid][2];
if (fabsf(d - sort_i[pid].d) - ci->dx_max_sort >
- 1.0e-4 * max(fabsf(d), ci->dx_max_sort_old))
+ 1.0e-4 * max(fabsf(d), ci->dx_max_sort_old) &&
+ fabsf(d - sort_i[pid].d) - ci->dx_max_sort > ci->width[0] * 1.0e-10)
error(
"particle shift diff exceeds dx_max_sort in cell ci. ci->nodeID=%d "
"cj->nodeID=%d d=%e sort_i[pid].d=%e ci->dx_max_sort=%e "
@@ -1686,7 +1689,8 @@ void DOPAIR2_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
p->x[1] * runner_shift[sid][1] +
p->x[2] * runner_shift[sid][2];
if (fabsf(d - sort_j[pjd].d) - cj->dx_max_sort >
- 1.0e-4 * max(fabsf(d), cj->dx_max_sort_old))
+ 1.0e-4 * max(fabsf(d), cj->dx_max_sort_old) &&
+ fabsf(d - sort_j[pjd].d) - cj->dx_max_sort > cj->width[0] * 1.0e-10)
error(
"particle shift diff exceeds dx_max_sort in cell cj. cj->nodeID=%d "
"ci->nodeID=%d d=%e sort_j[pjd].d=%e cj->dx_max_sort=%e "
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 0bc7b851c05e94974ef6337ba3d8df5e5420c9a1..cec450d84449d2b8cbb6790fcf9a596bb24649be 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -21,10 +21,12 @@
#define SWIFT_RUNNER_DOIACT_GRAV_H
/* Includes. */
+#include "active.h"
#include "cell.h"
#include "gravity.h"
#include "inline.h"
#include "part.h"
+#include "timers.h"
/**
* @brief Recursively propagate the multipoles down the tree by applying the
@@ -34,7 +36,8 @@
* @param c The #cell we are working on.
* @param timer Are we timing this ?
*/
-void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
+static INLINE void runner_do_grav_down(struct runner *r, struct cell *c,
+ int timer) {
/* Some constants */
const struct engine *e = r->e;
@@ -125,8 +128,9 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
* @param ci The #cell with field tensor to interact.
* @param cj The #cell with the multipole.
*/
-void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
- struct cell *restrict cj) {
+static INLINE void runner_dopair_grav_mm(const struct runner *r,
+ struct cell *restrict ci,
+ struct cell *restrict cj) {
/* Some constants */
const struct engine *e = r->e;
@@ -204,8 +208,8 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
const float h_inv_i = 1.f / h_i;
const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
- /* Local accumulators for the acceleration */
- float a_x = 0.f, a_y = 0.f, a_z = 0.f;
+ /* Local accumulators for the acceleration and potential */
+ float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
/* Make the compiler understand we are in happy vectorization land */
swift_align_information(float, cj_cache->x, SWIFT_CACHE_ALIGNMENT);
@@ -240,13 +244,15 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
#endif
/* Interact! */
- float f_ij;
- runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij);
+ float f_ij, pot_ij;
+ runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij,
+ &pot_ij);
/* Store it back */
a_x -= f_ij * dx;
a_y -= f_ij * dy;
a_z -= f_ij * dz;
+ pot += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter if it's not a padded gpart */
@@ -258,6 +264,7 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
+ ci_cache->pot[pid] = pot;
}
TIMER_TOC(timer_dopair_grav_pp);
@@ -295,8 +302,8 @@ static INLINE void runner_dopair_grav_pp_truncated(
const float h_inv_i = 1.f / h_i;
const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
- /* Local accumulators for the acceleration */
- float a_x = 0.f, a_y = 0.f, a_z = 0.f;
+ /* Local accumulators for the acceleration and potential */
+ float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
/* Make the compiler understand we are in happy vectorization land */
swift_align_information(float, cj_cache->x, SWIFT_CACHE_ALIGNMENT);
@@ -331,14 +338,15 @@ static INLINE void runner_dopair_grav_pp_truncated(
#endif
/* Interact! */
- float f_ij;
+ float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
- rlr_inv, &f_ij);
+ rlr_inv, &f_ij, &pot_ij);
/* Store it back */
a_x -= f_ij * dx;
a_y -= f_ij * dy;
a_z -= f_ij * dz;
+ pot += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter if it's not a padded gpart */
@@ -350,6 +358,7 @@ static INLINE void runner_dopair_grav_pp_truncated(
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
+ ci_cache->pot[pid] = pot;
}
TIMER_TOC(timer_dopair_grav_pp);
@@ -372,6 +381,7 @@ static INLINE void runner_dopair_grav_pm(
swift_declare_aligned_ptr(float, a_x, ci_cache->a_x, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, a_y, ci_cache->a_y, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, a_z, ci_cache->a_z, SWIFT_CACHE_ALIGNMENT);
+ swift_declare_aligned_ptr(float, pot, ci_cache->pot, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(int, active, ci_cache->active,
SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(int, use_mpole, ci_cache->use_mpole,
@@ -407,14 +417,15 @@ static INLINE void runner_dopair_grav_pm(
const float r2 = dx * dx + dy * dy + dz * dz;
/* Interact! */
- float f_x, f_y, f_z;
- runner_iact_grav_pm(dx, dy, dz, r2, h_i, h_inv_i, multi_j, &f_x, &f_y,
- &f_z);
+ float f_x, f_y, f_z, pot_ij;
+ runner_iact_grav_pm(dx, dy, dz, r2, h_i, h_inv_i, multi_j, &f_x, &f_y, &f_z,
+ &pot_ij);
/* Store it back */
a_x[pid] = f_x;
a_y[pid] = f_y;
a_z[pid] = f_z;
+ pot[pid] = pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter */
@@ -434,7 +445,8 @@ static INLINE void runner_dopair_grav_pm(
* @param ci The first #cell.
* @param cj The other #cell.
*/
-void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
+static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
+ struct cell *cj) {
const struct engine *e = r->e;
@@ -635,7 +647,8 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
*
* @todo Use a local cache for the particles.
*/
-void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
+static INLINE void runner_doself_grav_pp_full(struct runner *r,
+ struct cell *c) {
/* Some constants */
const struct engine *const e = r->e;
@@ -683,7 +696,7 @@ void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
/* Local accumulators for the acceleration */
- float a_x = 0.f, a_y = 0.f, a_z = 0.f;
+ float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
/* Make the compiler understand we are in happy vectorization land */
swift_align_information(float, ci_cache->x, SWIFT_CACHE_ALIGNMENT);
@@ -721,13 +734,15 @@ void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
#endif
/* Interact! */
- float f_ij;
- runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij);
+ float f_ij, pot_ij;
+ runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij,
+ &pot_ij);
/* Store it back */
a_x -= f_ij * dx;
a_y -= f_ij * dy;
a_z -= f_ij * dz;
+ pot += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter if it's not a padded gpart */
@@ -739,6 +754,7 @@ void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
+ ci_cache->pot[pid] = pot;
}
/* Write back to the particles */
@@ -754,7 +770,8 @@ void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
*
* @todo Use a local cache for the particles.
*/
-void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
+static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
+ struct cell *c) {
/* Some constants */
const struct engine *const e = r->e;
@@ -808,8 +825,8 @@ void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
const float h_inv_i = 1.f / h_i;
const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
- /* Local accumulators for the acceleration */
- float a_x = 0.f, a_y = 0.f, a_z = 0.f;
+ /* Local accumulators for the acceleration and potential */
+ float a_x = 0.f, a_y = 0.f, a_z = 0.f, pot = 0.f;
/* Make the compiler understand we are in happy vectorization land */
swift_align_information(float, ci_cache->x, SWIFT_CACHE_ALIGNMENT);
@@ -847,14 +864,15 @@ void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
#endif
/* Interact! */
- float f_ij;
+ float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
- rlr_inv, &f_ij);
+ rlr_inv, &f_ij, &pot_ij);
/* Store it back */
a_x -= f_ij * dx;
a_y -= f_ij * dy;
a_z -= f_ij * dz;
+ pot += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter if it's not a padded gpart */
@@ -866,6 +884,7 @@ void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
+ ci_cache->pot[pid] = pot;
}
/* Write back to the particles */
@@ -879,7 +898,7 @@ void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
* @param r The #runner.
* @param c The #cell.
*/
-void runner_doself_grav_pp(struct runner *r, struct cell *c) {
+static INLINE void runner_doself_grav_pp(struct runner *r, struct cell *c) {
/* Some properties of the space */
const struct engine *e = r->e;
@@ -934,8 +953,8 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
*
* @todo Use a local cache for the particles.
*/
-void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
- int gettimer) {
+static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci,
+ struct cell *cj, int gettimer) {
/* Some constants */
const struct engine *e = r->e;
@@ -1086,7 +1105,8 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
*
* @todo Use a local cache for the particles.
*/
-void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
+static INLINE void runner_doself_grav(struct runner *r, struct cell *c,
+ int gettimer) {
/* Some constants */
const struct engine *e = r->e;
@@ -1137,7 +1157,8 @@ void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
* @param ci The #cell of interest.
* @param timer Are we timing this ?
*/
-void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
+static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
+ int timer) {
/* Some constants */
const struct engine *e = r->e;
diff --git a/src/runner_doiact_vec.c b/src/runner_doiact_vec.c
index aa6e6dc3c15f1846462b900076cf131e726cd883..0daea0cf8a8d868ccb14f7d9b27f04aca5a677c8 100644
--- a/src/runner_doiact_vec.c
+++ b/src/runner_doiact_vec.c
@@ -1111,6 +1111,7 @@ void runner_doself2_force_vec(struct runner *r, struct cell *restrict c) {
#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
const struct engine *e = r->e;
+ const struct cosmology *restrict cosmo = e->cosmology;
const timebin_t max_active_bin = e->max_active_bin;
struct part *restrict parts = c->parts;
const int count = c->count;
@@ -1139,6 +1140,10 @@ void runner_doself2_force_vec(struct runner *r, struct cell *restrict c) {
/* Read the particles from the cell and store them locally in the cache. */
cache_read_force_particles(c, cell_cache);
+ /* 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++) {
@@ -1262,7 +1267,7 @@ void runner_doself2_force_vec(struct runner *r, struct cell *restrict c) {
&cell_cache->vy[pjd], &cell_cache->vz[pjd], &cell_cache->rho[pjd],
&cell_cache->grad_h[pjd], &cell_cache->pOrho2[pjd],
&cell_cache->balsara[pjd], &cell_cache->soundspeed[pjd],
- &cell_cache->m[pjd], v_hi_inv, v_hj_inv, &v_a_hydro_xSum,
+ &cell_cache->m[pjd], v_hi_inv, v_hj_inv, a, H, &v_a_hydro_xSum,
&v_a_hydro_ySum, &v_a_hydro_zSum, &v_h_dtSum, &v_sigSum,
&v_entropy_dtSum, v_doi_mask);
}
@@ -1988,6 +1993,7 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
const struct engine *restrict e = r->e;
+ const struct cosmology *restrict cosmo = e->cosmology;
const timebin_t max_active_bin = e->max_active_bin;
TIMER_TIC;
@@ -2019,6 +2025,10 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
const int active_ci = cell_is_active_hydro(ci, e) && ci_local;
const int active_cj = cell_is_active_hydro(cj, e) && cj_local;
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles have been drifted to the current time */
for (int pid = 0; pid < count_i; pid++)
@@ -2213,7 +2223,7 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
&cj_cache->grad_h[cj_cache_idx], &cj_cache->pOrho2[cj_cache_idx],
&cj_cache->balsara[cj_cache_idx],
&cj_cache->soundspeed[cj_cache_idx], &cj_cache->m[cj_cache_idx],
- v_hi_inv, v_hj_inv, &v_a_hydro_xSum, &v_a_hydro_ySum,
+ v_hi_inv, v_hj_inv, a, H, &v_a_hydro_xSum, &v_a_hydro_ySum,
&v_a_hydro_zSum, &v_h_dtSum, &v_sigSum, &v_entropy_dtSum,
v_doi_mask);
}
@@ -2349,7 +2359,7 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
&ci_cache->grad_h[ci_cache_idx], &ci_cache->pOrho2[ci_cache_idx],
&ci_cache->balsara[ci_cache_idx],
&ci_cache->soundspeed[ci_cache_idx], &ci_cache->m[ci_cache_idx],
- v_hj_inv, v_hi_inv, &v_a_hydro_xSum, &v_a_hydro_ySum,
+ v_hj_inv, v_hi_inv, a, H, &v_a_hydro_xSum, &v_a_hydro_ySum,
&v_a_hydro_zSum, &v_h_dtSum, &v_sigSum, &v_entropy_dtSum,
v_doj_mask);
}
diff --git a/src/scheduler.c b/src/scheduler.c
index 2c0406953581f55a1086ef313da8a389c18a9bd5..756f70c6b02b9e58d0a8a7f1e2b1ad7f11af4ca5 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1435,7 +1435,7 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
switch (t->type) {
case task_type_self:
case task_type_sub_self:
- if (t->subtype == task_subtype_grav)
+ if (t->subtype == task_subtype_grav || t->subtype == task_subtype_external_grav)
qid = t->ci->super_gravity->owner;
else
qid = t->ci->super_hydro->owner;
diff --git a/src/serial_io.c b/src/serial_io.c
index 3e62bd7fda99ed7703a4fe67340d5a03e5dc45b8..aa8d126e4c6f3e19aed7d205c7ced39f7eded6f8 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -197,13 +197,13 @@ void prepareArray(const struct engine* e, hid_t grp, char* fileName,
rank = 2;
shape[0] = N_total;
shape[1] = props.dimension;
- chunk_shape[0] = 1 << 16; /* Just a guess...*/
+ chunk_shape[0] = 1 << 20; /* Just a guess...*/
chunk_shape[1] = props.dimension;
} else {
rank = 1;
shape[0] = N_total;
shape[1] = 0;
- chunk_shape[0] = 1 << 16; /* Just a guess...*/
+ chunk_shape[0] = 1 << 20; /* Just a guess...*/
chunk_shape[1] = 0;
}
@@ -211,7 +211,7 @@ void prepareArray(const struct engine* e, hid_t grp, char* fileName,
if (chunk_shape[0] > N_total) chunk_shape[0] = N_total;
/* Change shape of data space */
- hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, NULL);
+ hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, shape);
if (h_err < 0) {
error("Error while changing data space shape for field '%s'.", props.name);
}
@@ -228,11 +228,15 @@ void prepareArray(const struct engine* e, hid_t grp, char* fileName,
/* Impose data compression */
if (e->snapshotCompression > 0) {
+ h_err = H5Pset_shuffle(h_prop);
+ if (h_err < 0)
+ error("Error while setting shuffling options for field '%s'.",
+ props.name);
+
h_err = H5Pset_deflate(h_prop, e->snapshotCompression);
- if (h_err < 0) {
+ if (h_err < 0)
error("Error while setting compression options for field '%s'.",
props.name);
- }
}
/* Create dataset */
diff --git a/src/single_io.c b/src/single_io.c
index caa39e4a7cf583b210d3e9322abcd53c332cfe70..0fecd724908f7ce2779526792ce8512b9236d58e 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -204,13 +204,13 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
rank = 2;
shape[0] = N;
shape[1] = props.dimension;
- chunk_shape[0] = 1 << 16; /* Just a guess...*/
+ chunk_shape[0] = 1 << 20; /* Just a guess...*/
chunk_shape[1] = props.dimension;
} else {
rank = 1;
shape[0] = N;
shape[1] = 0;
- chunk_shape[0] = 1 << 16; /* Just a guess...*/
+ chunk_shape[0] = 1 << 20; /* Just a guess...*/
chunk_shape[1] = 0;
}
@@ -218,7 +218,7 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
if (chunk_shape[0] > N) chunk_shape[0] = N;
/* Change shape of data space */
- hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, NULL);
+ hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, shape);
if (h_err < 0) {
error("Error while changing data space shape for field '%s'.", props.name);
}
@@ -235,11 +235,15 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
/* Impose data compression */
if (e->snapshotCompression > 0) {
+ h_err = H5Pset_shuffle(h_prop);
+ if (h_err < 0)
+ error("Error while setting shuffling options for field '%s'.",
+ props.name);
+
h_err = H5Pset_deflate(h_prop, e->snapshotCompression);
- if (h_err < 0) {
+ if (h_err < 0)
error("Error while setting compression options for field '%s'.",
props.name);
- }
}
/* Create dataset */
diff --git a/src/statistics.c b/src/statistics.c
index ed3197fbd0538394bb1c40d834c16e174a410068..8d8afae0d1441927e001d7d58a3ee5ff0399f5ed 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -191,7 +191,7 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
stats.E_int += m * u_inter;
stats.E_rad += cooling_get_radiated_energy(xp);
if (gp != NULL) {
- stats.E_pot_self += m * gravity_get_potential(gp) * a_inv;
+ stats.E_pot_self += 0.5f * m * gravity_get_physical_potential(gp, cosmo);
stats.E_pot_ext += m * external_gravity_get_potential_energy(
time, potential, phys_const, gp);
}
@@ -292,7 +292,7 @@ void stats_collect_gpart_mapper(void *map_data, int nr_gparts,
/* Collect energies. */
stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) *
a_inv2; /* 1/2 m a^2 \dot{r}^2 */
- stats.E_pot_self += m * gravity_get_potential(gp) * a_inv;
+ stats.E_pot_self += 0.5f * m * gravity_get_physical_potential(gp, cosmo);
stats.E_pot_ext += m * external_gravity_get_potential_energy(
time, potential, phys_const, gp);
}
diff --git a/src/vector_power.h b/src/vector_power.h
index 6e8377586fdbc072e8f894c0805f43c122f8c54f..ffaa3b3f6cea69626462a9b5f26f94827c6af556 100644
--- a/src/vector_power.h
+++ b/src/vector_power.h
@@ -415,273 +415,273 @@ __attribute__((always_inline)) INLINE static double X_112(const double v[3]) {
/***************************/
/**
-* @brief Compute \f$ \frac{1}{(0,0,5)!}\vec{v}^{(0,0,5)} \f$.
-*
-* Note \f$ \vec{v}^{(0,0,5)} = v_z^5 \f$
-* and \f$ \frac{1}{(0,0,5)!} = 1/(0!*0!*5!) = 1/120 = 8.333333e-03 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(0,0,5)!}\vec{v}^{(0,0,5)} \f$.
+ *
+ * Note \f$ \vec{v}^{(0,0,5)} = v_z^5 \f$
+ * and \f$ \frac{1}{(0,0,5)!} = 1/(0!*0!*5!) = 1/120 = 8.333333e-03 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_005(const double v[3]) {
return 8.333333333333333e-03 * v[2] * v[2] * v[2] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(0,1,4)!}\vec{v}^{(0,1,4)} \f$.
-*
-* Note \f$ \vec{v}^{(0,1,4)} = v_y^1 v_z^4 \f$
-* and \f$ \frac{1}{(0,1,4)!} = 1/(0!*1!*4!) = 1/24 = 4.166667e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(0,1,4)!}\vec{v}^{(0,1,4)} \f$.
+ *
+ * Note \f$ \vec{v}^{(0,1,4)} = v_y^1 v_z^4 \f$
+ * and \f$ \frac{1}{(0,1,4)!} = 1/(0!*1!*4!) = 1/24 = 4.166667e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_014(const double v[3]) {
return 4.166666666666666e-02 * v[1] * v[2] * v[2] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(0,2,3)!}\vec{v}^{(0,2,3)} \f$.
-*
-* Note \f$ \vec{v}^{(0,2,3)} = v_y^2 v_z^3 \f$
-* and \f$ \frac{1}{(0,2,3)!} = 1/(0!*2!*3!) = 1/12 = 8.333333e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(0,2,3)!}\vec{v}^{(0,2,3)} \f$.
+ *
+ * Note \f$ \vec{v}^{(0,2,3)} = v_y^2 v_z^3 \f$
+ * and \f$ \frac{1}{(0,2,3)!} = 1/(0!*2!*3!) = 1/12 = 8.333333e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_023(const double v[3]) {
return 8.333333333333333e-02 * v[1] * v[1] * v[2] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(0,3,2)!}\vec{v}^{(0,3,2)} \f$.
-*
-* Note \f$ \vec{v}^{(0,3,2)} = v_y^3 v_z^2 \f$
-* and \f$ \frac{1}{(0,3,2)!} = 1/(0!*3!*2!) = 1/12 = 8.333333e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(0,3,2)!}\vec{v}^{(0,3,2)} \f$.
+ *
+ * Note \f$ \vec{v}^{(0,3,2)} = v_y^3 v_z^2 \f$
+ * and \f$ \frac{1}{(0,3,2)!} = 1/(0!*3!*2!) = 1/12 = 8.333333e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_032(const double v[3]) {
return 8.333333333333333e-02 * v[1] * v[1] * v[1] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(0,4,1)!}\vec{v}^{(0,4,1)} \f$.
-*
-* Note \f$ \vec{v}^{(0,4,1)} = v_y^4 v_z^1 \f$
-* and \f$ \frac{1}{(0,4,1)!} = 1/(0!*4!*1!) = 1/24 = 4.166667e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(0,4,1)!}\vec{v}^{(0,4,1)} \f$.
+ *
+ * Note \f$ \vec{v}^{(0,4,1)} = v_y^4 v_z^1 \f$
+ * and \f$ \frac{1}{(0,4,1)!} = 1/(0!*4!*1!) = 1/24 = 4.166667e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_041(const double v[3]) {
return 4.166666666666666e-02 * v[1] * v[1] * v[1] * v[1] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(0,5,0)!}\vec{v}^{(0,5,0)} \f$.
-*
-* Note \f$ \vec{v}^{(0,5,0)} = v_y^5 \f$
-* and \f$ \frac{1}{(0,5,0)!} = 1/(0!*5!*0!) = 1/120 = 8.333333e-03 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(0,5,0)!}\vec{v}^{(0,5,0)} \f$.
+ *
+ * Note \f$ \vec{v}^{(0,5,0)} = v_y^5 \f$
+ * and \f$ \frac{1}{(0,5,0)!} = 1/(0!*5!*0!) = 1/120 = 8.333333e-03 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_050(const double v[3]) {
return 8.333333333333333e-03 * v[1] * v[1] * v[1] * v[1] * v[1];
}
/**
-* @brief Compute \f$ \frac{1}{(1,0,4)!}\vec{v}^{(1,0,4)} \f$.
-*
-* Note \f$ \vec{v}^{(1,0,4)} = v_x^1 v_z^4 \f$
-* and \f$ \frac{1}{(1,0,4)!} = 1/(1!*0!*4!) = 1/24 = 4.166667e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(1,0,4)!}\vec{v}^{(1,0,4)} \f$.
+ *
+ * Note \f$ \vec{v}^{(1,0,4)} = v_x^1 v_z^4 \f$
+ * and \f$ \frac{1}{(1,0,4)!} = 1/(1!*0!*4!) = 1/24 = 4.166667e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_104(const double v[3]) {
return 4.166666666666666e-02 * v[0] * v[2] * v[2] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(1,1,3)!}\vec{v}^{(1,1,3)} \f$.
-*
-* Note \f$ \vec{v}^{(1,1,3)} = v_x^1 v_y^1 v_z^3 \f$
-* and \f$ \frac{1}{(1,1,3)!} = 1/(1!*1!*3!) = 1/6 = 1.666667e-01 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(1,1,3)!}\vec{v}^{(1,1,3)} \f$.
+ *
+ * Note \f$ \vec{v}^{(1,1,3)} = v_x^1 v_y^1 v_z^3 \f$
+ * and \f$ \frac{1}{(1,1,3)!} = 1/(1!*1!*3!) = 1/6 = 1.666667e-01 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_113(const double v[3]) {
return 1.666666666666667e-01 * v[0] * v[1] * v[2] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(1,2,2)!}\vec{v}^{(1,2,2)} \f$.
-*
-* Note \f$ \vec{v}^{(1,2,2)} = v_x^1 v_y^2 v_z^2 \f$
-* and \f$ \frac{1}{(1,2,2)!} = 1/(1!*2!*2!) = 1/4 = 2.500000e-01 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(1,2,2)!}\vec{v}^{(1,2,2)} \f$.
+ *
+ * Note \f$ \vec{v}^{(1,2,2)} = v_x^1 v_y^2 v_z^2 \f$
+ * and \f$ \frac{1}{(1,2,2)!} = 1/(1!*2!*2!) = 1/4 = 2.500000e-01 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_122(const double v[3]) {
return 2.500000000000000e-01 * v[0] * v[1] * v[1] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(1,3,1)!}\vec{v}^{(1,3,1)} \f$.
-*
-* Note \f$ \vec{v}^{(1,3,1)} = v_x^1 v_y^3 v_z^1 \f$
-* and \f$ \frac{1}{(1,3,1)!} = 1/(1!*3!*1!) = 1/6 = 1.666667e-01 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(1,3,1)!}\vec{v}^{(1,3,1)} \f$.
+ *
+ * Note \f$ \vec{v}^{(1,3,1)} = v_x^1 v_y^3 v_z^1 \f$
+ * and \f$ \frac{1}{(1,3,1)!} = 1/(1!*3!*1!) = 1/6 = 1.666667e-01 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_131(const double v[3]) {
return 1.666666666666667e-01 * v[0] * v[1] * v[1] * v[1] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(1,4,0)!}\vec{v}^{(1,4,0)} \f$.
-*
-* Note \f$ \vec{v}^{(1,4,0)} = v_x^1 v_y^4 \f$
-* and \f$ \frac{1}{(1,4,0)!} = 1/(1!*4!*0!) = 1/24 = 4.166667e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(1,4,0)!}\vec{v}^{(1,4,0)} \f$.
+ *
+ * Note \f$ \vec{v}^{(1,4,0)} = v_x^1 v_y^4 \f$
+ * and \f$ \frac{1}{(1,4,0)!} = 1/(1!*4!*0!) = 1/24 = 4.166667e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_140(const double v[3]) {
return 4.166666666666666e-02 * v[0] * v[1] * v[1] * v[1] * v[1];
}
/**
-* @brief Compute \f$ \frac{1}{(2,0,3)!}\vec{v}^{(2,0,3)} \f$.
-*
-* Note \f$ \vec{v}^{(2,0,3)} = v_x^2 v_z^3 \f$
-* and \f$ \frac{1}{(2,0,3)!} = 1/(2!*0!*3!) = 1/12 = 8.333333e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(2,0,3)!}\vec{v}^{(2,0,3)} \f$.
+ *
+ * Note \f$ \vec{v}^{(2,0,3)} = v_x^2 v_z^3 \f$
+ * and \f$ \frac{1}{(2,0,3)!} = 1/(2!*0!*3!) = 1/12 = 8.333333e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_203(const double v[3]) {
return 8.333333333333333e-02 * v[0] * v[0] * v[2] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(2,1,2)!}\vec{v}^{(2,1,2)} \f$.
-*
-* Note \f$ \vec{v}^{(2,1,2)} = v_x^2 v_y^1 v_z^2 \f$
-* and \f$ \frac{1}{(2,1,2)!} = 1/(2!*1!*2!) = 1/4 = 2.500000e-01 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(2,1,2)!}\vec{v}^{(2,1,2)} \f$.
+ *
+ * Note \f$ \vec{v}^{(2,1,2)} = v_x^2 v_y^1 v_z^2 \f$
+ * and \f$ \frac{1}{(2,1,2)!} = 1/(2!*1!*2!) = 1/4 = 2.500000e-01 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_212(const double v[3]) {
return 2.500000000000000e-01 * v[0] * v[0] * v[1] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(2,2,1)!}\vec{v}^{(2,2,1)} \f$.
-*
-* Note \f$ \vec{v}^{(2,2,1)} = v_x^2 v_y^2 v_z^1 \f$
-* and \f$ \frac{1}{(2,2,1)!} = 1/(2!*2!*1!) = 1/4 = 2.500000e-01 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(2,2,1)!}\vec{v}^{(2,2,1)} \f$.
+ *
+ * Note \f$ \vec{v}^{(2,2,1)} = v_x^2 v_y^2 v_z^1 \f$
+ * and \f$ \frac{1}{(2,2,1)!} = 1/(2!*2!*1!) = 1/4 = 2.500000e-01 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_221(const double v[3]) {
return 2.500000000000000e-01 * v[0] * v[0] * v[1] * v[1] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(2,3,0)!}\vec{v}^{(2,3,0)} \f$.
-*
-* Note \f$ \vec{v}^{(2,3,0)} = v_x^2 v_y^3 \f$
-* and \f$ \frac{1}{(2,3,0)!} = 1/(2!*3!*0!) = 1/12 = 8.333333e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(2,3,0)!}\vec{v}^{(2,3,0)} \f$.
+ *
+ * Note \f$ \vec{v}^{(2,3,0)} = v_x^2 v_y^3 \f$
+ * and \f$ \frac{1}{(2,3,0)!} = 1/(2!*3!*0!) = 1/12 = 8.333333e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_230(const double v[3]) {
return 8.333333333333333e-02 * v[0] * v[0] * v[1] * v[1] * v[1];
}
/**
-* @brief Compute \f$ \frac{1}{(3,0,2)!}\vec{v}^{(3,0,2)} \f$.
-*
-* Note \f$ \vec{v}^{(3,0,2)} = v_x^3 v_z^2 \f$
-* and \f$ \frac{1}{(3,0,2)!} = 1/(3!*0!*2!) = 1/12 = 8.333333e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(3,0,2)!}\vec{v}^{(3,0,2)} \f$.
+ *
+ * Note \f$ \vec{v}^{(3,0,2)} = v_x^3 v_z^2 \f$
+ * and \f$ \frac{1}{(3,0,2)!} = 1/(3!*0!*2!) = 1/12 = 8.333333e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_302(const double v[3]) {
return 8.333333333333333e-02 * v[0] * v[0] * v[0] * v[2] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(3,1,1)!}\vec{v}^{(3,1,1)} \f$.
-*
-* Note \f$ \vec{v}^{(3,1,1)} = v_x^3 v_y^1 v_z^1 \f$
-* and \f$ \frac{1}{(3,1,1)!} = 1/(3!*1!*1!) = 1/6 = 1.666667e-01 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(3,1,1)!}\vec{v}^{(3,1,1)} \f$.
+ *
+ * Note \f$ \vec{v}^{(3,1,1)} = v_x^3 v_y^1 v_z^1 \f$
+ * and \f$ \frac{1}{(3,1,1)!} = 1/(3!*1!*1!) = 1/6 = 1.666667e-01 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_311(const double v[3]) {
return 1.666666666666667e-01 * v[0] * v[0] * v[0] * v[1] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(3,2,0)!}\vec{v}^{(3,2,0)} \f$.
-*
-* Note \f$ \vec{v}^{(3,2,0)} = v_x^3 v_y^2 \f$
-* and \f$ \frac{1}{(3,2,0)!} = 1/(3!*2!*0!) = 1/12 = 8.333333e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(3,2,0)!}\vec{v}^{(3,2,0)} \f$.
+ *
+ * Note \f$ \vec{v}^{(3,2,0)} = v_x^3 v_y^2 \f$
+ * and \f$ \frac{1}{(3,2,0)!} = 1/(3!*2!*0!) = 1/12 = 8.333333e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_320(const double v[3]) {
return 8.333333333333333e-02 * v[0] * v[0] * v[0] * v[1] * v[1];
}
/**
-* @brief Compute \f$ \frac{1}{(4,0,1)!}\vec{v}^{(4,0,1)} \f$.
-*
-* Note \f$ \vec{v}^{(4,0,1)} = v_x^4 v_z^1 \f$
-* and \f$ \frac{1}{(4,0,1)!} = 1/(4!*0!*1!) = 1/24 = 4.166667e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(4,0,1)!}\vec{v}^{(4,0,1)} \f$.
+ *
+ * Note \f$ \vec{v}^{(4,0,1)} = v_x^4 v_z^1 \f$
+ * and \f$ \frac{1}{(4,0,1)!} = 1/(4!*0!*1!) = 1/24 = 4.166667e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_401(const double v[3]) {
return 4.166666666666666e-02 * v[0] * v[0] * v[0] * v[0] * v[2];
}
/**
-* @brief Compute \f$ \frac{1}{(4,1,0)!}\vec{v}^{(4,1,0)} \f$.
-*
-* Note \f$ \vec{v}^{(4,1,0)} = v_x^4 v_y^1 \f$
-* and \f$ \frac{1}{(4,1,0)!} = 1/(4!*1!*0!) = 1/24 = 4.166667e-02 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(4,1,0)!}\vec{v}^{(4,1,0)} \f$.
+ *
+ * Note \f$ \vec{v}^{(4,1,0)} = v_x^4 v_y^1 \f$
+ * and \f$ \frac{1}{(4,1,0)!} = 1/(4!*1!*0!) = 1/24 = 4.166667e-02 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_410(const double v[3]) {
return 4.166666666666666e-02 * v[0] * v[0] * v[0] * v[0] * v[1];
}
/**
-* @brief Compute \f$ \frac{1}{(5,0,0)!}\vec{v}^{(5,0,0)} \f$.
-*
-* Note \f$ \vec{v}^{(5,0,0)} = v_x^5 \f$
-* and \f$ \frac{1}{(5,0,0)!} = 1/(5!*0!*0!) = 1/120 = 8.333333e-03 \f$
-*
-* @param v vector (\f$ v \f$).
-*/
+ * @brief Compute \f$ \frac{1}{(5,0,0)!}\vec{v}^{(5,0,0)} \f$.
+ *
+ * Note \f$ \vec{v}^{(5,0,0)} = v_x^5 \f$
+ * and \f$ \frac{1}{(5,0,0)!} = 1/(5!*0!*0!) = 1/120 = 8.333333e-03 \f$
+ *
+ * @param v vector (\f$ v \f$).
+ */
__attribute__((always_inline)) INLINE static double X_500(const double v[3]) {
return 8.333333333333333e-03 * v[0] * v[0] * v[0] * v[0] * v[0];
diff --git a/tests/Makefile.am b/tests/Makefile.am
index e830aad453653dd2ffa265627953af1c28d8cf59..8a757e0b87ffae4d6daa1a6f2b67cf9548b868d0 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -15,9 +15,9 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Add the source directory and debug to CFLAGS
-AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS)
+AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS) $(GSL_INCS)
-AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS) $(GRACKLE_LIBS)
+AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS)
# List of programs and scripts to run in the test suite
TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetry \
@@ -26,7 +26,8 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC \
testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
- testPeriodicBC.sh testPeriodicBCPerturbed.sh
+ testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
+ testPotentialPair
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -36,7 +37,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
- testGravityDerivatives
+ testGravityDerivatives testPotentialSelf testPotentialPair
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -100,6 +101,10 @@ testLogger_SOURCES = testLogger.c
testGravityDerivatives_SOURCES = testGravityDerivatives.c
+testPotentialSelf_SOURCES = testPotentialSelf.c
+
+testPotentialPair_SOURCES = testPotentialPair.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 31b7eb431da6fc12e767bc9053a21dc9ef010598..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_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,
+ &(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);
}
diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c
new file mode 100644
index 0000000000000000000000000000000000000000..5fbfb04ecd0b4a8c05c07d622dd77bd38af970ae
--- /dev/null
+++ b/tests/testPotentialPair.c
@@ -0,0 +1,367 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 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/>.
+ *
+ ******************************************************************************/
+#include "../config.h"
+
+/* Some standard headers. */
+#include <fenv.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+
+/* Local headers. */
+#include "runner_doiact_grav.h"
+#include "swift.h"
+
+const int num_tests = 100;
+const double eps = 0.1;
+
+/**
+ * @brief Check that a and b are consistent (up to some relative error)
+ *
+ * @param a First value
+ * @param b Second value
+ * @param s String used to identify this check in messages
+ */
+void check_value(double a, double b, const char *s) {
+ if (fabs(a - b) / fabs(a + b) > 2e-6 && fabs(a - b) > 1.e-6)
+ error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s);
+}
+
+/* Definitions of the potential and force that match
+ exactly the theory document */
+double S(double x) { return exp(x) / (1. + exp(x)); }
+
+double S_prime(double x) { return exp(x) / ((1. + exp(x)) * (1. + exp(x))); }
+
+double potential(double mass, double r, double H, double rlr) {
+
+ const double u = r / H;
+ const double x = r / rlr;
+ double pot;
+ if (u > 1.)
+ pot = -mass / r;
+ else
+ pot = -mass *
+ (-3. * u * u * u * u * u * u * u + 15. * u * u * u * u * u * u -
+ 28. * u * u * u * u * u + 21. * u * u * u * u - 7. * u * u + 3.) /
+ H;
+
+ return pot * (2. - 2. * S(2. * x));
+}
+
+double acceleration(double mass, double r, double H, double rlr) {
+
+ const double u = r / H;
+ const double x = r / rlr;
+ double acc;
+ if (u > 1.)
+ acc = -mass / (r * r * r);
+ else
+ acc = -mass * (21. * u * u * u * u * u - 90. * u * u * u * u +
+ 140. * u * u * u - 84. * u * u + 14.) /
+ (H * H * H);
+
+ return r * acc * (4. * x * S_prime(2 * x) - 2. * S(2. * x) + 2.);
+}
+
+int main() {
+
+ /* Initialize CPU frequency, this also starts time. */
+ unsigned long long cpufreq = 0;
+ clocks_set_cpufreq(cpufreq);
+
+ /* Initialise a few things to get us going */
+ struct engine e;
+ e.max_active_bin = num_time_bins;
+ e.time = 0.1f;
+ e.ti_current = 8;
+ e.time_base = 1e-10;
+
+ struct space s;
+ s.periodic = 0;
+ s.width[0] = 0.2;
+ s.width[1] = 0.2;
+ s.width[2] = 0.2;
+ e.s = &s;
+
+ struct gravity_props props;
+ props.a_smooth = 1.25;
+ props.r_cut_min = 0.;
+ props.theta_crit2 = 0.;
+ e.gravity_properties = &props;
+
+ struct runner r;
+ bzero(&r, sizeof(struct runner));
+ r.e = &e;
+
+ const double rlr = props.a_smooth * s.width[0] * FLT_MAX;
+
+ /* Init the cache for gravity interaction */
+ gravity_cache_init(&r.ci_gravity_cache, num_tests);
+ gravity_cache_init(&r.cj_gravity_cache, num_tests);
+
+ /* Let's create one cell with a massive particle and a bunch of test particles
+ */
+ struct cell ci, cj;
+ bzero(&ci, sizeof(struct cell));
+ bzero(&cj, sizeof(struct cell));
+
+ ci.width[0] = 1.;
+ ci.width[1] = 1.;
+ ci.width[2] = 1.;
+ ci.loc[0] = 0.;
+ ci.loc[1] = 0.;
+ ci.loc[2] = 0.;
+ ci.gcount = 1;
+ ci.ti_old_gpart = 8;
+ ci.ti_old_multipole = 8;
+ ci.ti_gravity_end_min = 8;
+ ci.ti_gravity_end_max = 8;
+
+ cj.width[0] = 1.;
+ cj.width[1] = 1.;
+ cj.width[2] = 1.;
+ cj.loc[0] = 1.;
+ cj.loc[1] = 0.;
+ cj.loc[2] = 0.;
+ cj.gcount = num_tests;
+ cj.ti_old_gpart = 8;
+ cj.ti_old_multipole = 8;
+ cj.ti_gravity_end_min = 8;
+ cj.ti_gravity_end_max = 8;
+
+ /* Allocate multipoles */
+ ci.multipole = malloc(sizeof(struct gravity_tensors));
+ cj.multipole = malloc(sizeof(struct gravity_tensors));
+ bzero(ci.multipole, sizeof(struct gravity_tensors));
+ bzero(cj.multipole, sizeof(struct gravity_tensors));
+
+ /* Set the multipoles */
+ ci.multipole->r_max = 0.1;
+ cj.multipole->r_max = 0.1;
+
+ /* Allocate the particles */
+ if (posix_memalign((void **)&ci.gparts, gpart_align,
+ ci.gcount * sizeof(struct gpart)) != 0)
+ error("Error allocating gparts for cell ci");
+ bzero(ci.gparts, ci.gcount * sizeof(struct gpart));
+
+ if (posix_memalign((void **)&cj.gparts, gpart_align,
+ cj.gcount * sizeof(struct gpart)) != 0)
+ error("Error allocating gparts for cell ci");
+ bzero(cj.gparts, cj.gcount * sizeof(struct gpart));
+
+ /* Create the mass-less test particles */
+ for (int n = 0; n < num_tests; ++n) {
+
+ struct gpart *gp = &cj.gparts[n];
+
+ gp->x[0] = 1. + (n + 1) / ((double)num_tests);
+ gp->x[1] = 0.5;
+ gp->x[2] = 0.5;
+ gp->mass = 0.;
+ gp->epsilon = eps;
+ gp->time_bin = 1;
+ gp->type = swift_type_dark_matter;
+ gp->id_or_neg_offset = n + 1;
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->ti_drift = 8;
+#endif
+ }
+
+ /***********************************************/
+ /* Let's start by testing the P-P interactions */
+ /***********************************************/
+
+ /* Create the massive particle */
+ ci.gparts[0].x[0] = 0.;
+ ci.gparts[0].x[1] = 0.5;
+ ci.gparts[0].x[2] = 0.5;
+ ci.gparts[0].mass = 1.;
+ ci.gparts[0].epsilon = eps;
+ ci.gparts[0].time_bin = 1;
+ ci.gparts[0].type = swift_type_dark_matter;
+ ci.gparts[0].id_or_neg_offset = 1;
+#ifdef SWIFT_DEBUG_CHECKS
+ ci.gparts[0].ti_drift = 8;
+#endif
+
+ /* Now compute the forces */
+ runner_dopair_grav_pp(&r, &ci, &cj);
+
+ /* Verify everything */
+ for (int n = 0; n < num_tests; ++n) {
+ const struct gpart *gp = &cj.gparts[n];
+ const struct gpart *gp2 = &ci.gparts[0];
+
+ double pot_true =
+ potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
+ double acc_true =
+ acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
+
+ message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - gp2->x[0],
+ gp->a_grav[0], acc_true, gp->potential, pot_true);
+
+ check_value(gp->potential, pot_true, "potential");
+ check_value(gp->a_grav[0], acc_true, "acceleration");
+ }
+
+ message("\n\t\t P-P interactions all good\n");
+
+ /* Reset the accelerations */
+ for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
+
+ /**********************************/
+ /* Test the basic PM interactions */
+ /**********************************/
+
+ /* Set an opening angle that allows P-M interactions */
+ props.theta_crit2 = 1.;
+
+ ci.gparts[0].mass = 0.;
+ ci.multipole->CoM[0] = 0.;
+ ci.multipole->CoM[1] = 0.5;
+ ci.multipole->CoM[2] = 0.5;
+
+ bzero(&ci.multipole->m_pole, sizeof(struct multipole));
+ bzero(&cj.multipole->m_pole, sizeof(struct multipole));
+ ci.multipole->m_pole.M_000 = 1.;
+
+ /* Now compute the forces */
+ runner_dopair_grav_pp(&r, &ci, &cj);
+
+ /* Verify everything */
+ for (int n = 0; n < num_tests; ++n) {
+ const struct gpart *gp = &cj.gparts[n];
+ const struct gravity_tensors *mpole = ci.multipole;
+
+ double pot_true = potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0],
+ gp->epsilon, rlr * FLT_MAX);
+ double acc_true =
+ acceleration(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], gp->epsilon,
+ rlr * FLT_MAX);
+
+ message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - mpole->CoM[0],
+ gp->a_grav[0], acc_true, gp->potential, pot_true);
+
+ check_value(gp->potential, pot_true, "potential");
+ check_value(gp->a_grav[0], acc_true, "acceleration");
+ }
+
+ message("\n\t\t basic P-M interactions all good\n");
+
+ /* Reset the accelerations */
+ for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
+
+/***************************************/
+/* Test the high-order PM interactions */
+/***************************************/
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER >= 3
+
+ /* Let's make ci more interesting */
+ free(ci.gparts);
+ ci.gcount = 8;
+ if (posix_memalign((void **)&ci.gparts, gpart_align,
+ ci.gcount * sizeof(struct gpart)) != 0)
+ error("Error allocating gparts for cell ci");
+ bzero(ci.gparts, ci.gcount * sizeof(struct gpart));
+
+ /* Place particles on a simple cube of side-length 0.2 */
+ for (int n = 0; n < 8; ++n) {
+ if (n & 1)
+ ci.gparts[n].x[0] = 0.0 - 0.1;
+ else
+ ci.gparts[n].x[0] = 0.0 + 0.1;
+
+ if (n & 2)
+ ci.gparts[n].x[1] = 0.5 - 0.1;
+ else
+ ci.gparts[n].x[1] = 0.5 + 0.1;
+
+ if (n & 2)
+ ci.gparts[n].x[2] = 0.5 - 0.1;
+ else
+ ci.gparts[n].x[2] = 0.5 + 0.1;
+
+ ci.gparts[n].mass = 1. / 8.;
+
+ ci.gparts[n].epsilon = eps;
+ ci.gparts[n].time_bin = 1;
+ ci.gparts[n].type = swift_type_dark_matter;
+ ci.gparts[n].id_or_neg_offset = 1;
+#ifdef SWIFT_DEBUG_CHECKS
+ ci.gparts[n].ti_drift = 8;
+#endif
+ }
+
+ /* Now let's make a multipole out of it. */
+ gravity_reset(ci.multipole);
+ gravity_P2M(ci.multipole, ci.gparts, ci.gcount);
+
+ // message("CoM=[%e %e %e]", ci.multipole->CoM[0], ci.multipole->CoM[1],
+ // ci.multipole->CoM[2]);
+ gravity_multipole_print(&ci.multipole->m_pole);
+
+ /* Compute the forces */
+ runner_dopair_grav_pp(&r, &ci, &cj);
+
+ /* Verify everything */
+ for (int n = 0; n < num_tests; ++n) {
+ const struct gpart *gp = &cj.gparts[n];
+ const struct gravity_tensors *mpole = ci.multipole;
+
+ double pot_true = 0, acc_true[3] = {0., 0., 0.};
+
+ for (int i = 0; i < 8; ++i) {
+ const struct gpart *gp2 = &ci.gparts[i];
+
+ const double dx[3] = {gp2->x[0] - gp->x[0], gp2->x[1] - gp->x[1],
+ gp2->x[2] - gp->x[2]};
+ const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
+
+ pot_true += potential(gp2->mass, d, gp->epsilon, rlr * FLT_MAX);
+ acc_true[0] -=
+ acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[0] / d;
+ acc_true[1] -=
+ acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[1] / d;
+ acc_true[2] -=
+ acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[2] / d;
+ }
+
+ message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e",
+ gp->x[0] - mpole->CoM[0], gp->a_grav[0], acc_true[0], gp->potential,
+ pot_true, acc_true[1], acc_true[2]);
+
+ // check_value(gp->potential, pot_true, "potential");
+ // check_value(gp->a_grav[0], acc_true[0], "acceleration");
+ }
+
+ message("\n\t\t high-order P-M interactions all good\n");
+
+#endif
+
+ free(ci.multipole);
+ free(cj.multipole);
+ free(ci.gparts);
+ free(cj.gparts);
+ return 0;
+}
diff --git a/tests/testPotentialSelf.c b/tests/testPotentialSelf.c
new file mode 100644
index 0000000000000000000000000000000000000000..33fe574e34bcc49ff8bac6f19f22b55fda19f186
--- /dev/null
+++ b/tests/testPotentialSelf.c
@@ -0,0 +1,187 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 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/>.
+ *
+ ******************************************************************************/
+#include "../config.h"
+
+/* Some standard headers. */
+#include <fenv.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+
+/* Local headers. */
+#include "runner_doiact_grav.h"
+#include "swift.h"
+
+const int num_tests = 100;
+const double eps = 0.02;
+
+/**
+ * @brief Check that a and b are consistent (up to some relative error)
+ *
+ * @param a First value
+ * @param b Second value
+ * @param s String used to identify this check in messages
+ */
+void check_value(double a, double b, const char *s) {
+ if (fabs(a - b) / fabs(a + b) > 2e-6 && fabs(a - b) > 1.e-6)
+ error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s);
+}
+
+/* Definitions of the potential and force that match
+ exactly the theory document */
+double S(double x) { return exp(x) / (1. + exp(x)); }
+
+double S_prime(double x) { return exp(x) / ((1. + exp(x)) * (1. + exp(x))); }
+
+double potential(double mass, double r, double H, double rlr) {
+
+ const double u = r / H;
+ const double x = r / rlr;
+ double pot;
+ if (u > 1.)
+ pot = -mass / r;
+ else
+ pot = -mass *
+ (-3. * u * u * u * u * u * u * u + 15. * u * u * u * u * u * u -
+ 28. * u * u * u * u * u + 21. * u * u * u * u - 7. * u * u + 3.) /
+ H;
+
+ return pot * (2. - 2. * S(2. * x));
+}
+
+double acceleration(double mass, double r, double H, double rlr) {
+
+ const double u = r / H;
+ const double x = r / rlr;
+ double acc;
+ if (u > 1.)
+ acc = -mass / (r * r * r);
+ else
+ acc = -mass * (21. * u * u * u * u * u - 90. * u * u * u * u +
+ 140. * u * u * u - 84. * u * u + 14.) /
+ (H * H * H);
+
+ return r * acc * (4. * x * S_prime(2 * x) - 2. * S(2. * x) + 2.);
+}
+
+int main() {
+
+ /* Initialize CPU frequency, this also starts time. */
+ unsigned long long cpufreq = 0;
+ clocks_set_cpufreq(cpufreq);
+
+ /* Initialise a few things to get us going */
+ struct engine e;
+ e.max_active_bin = num_time_bins;
+ e.time = 0.1f;
+ e.ti_current = 8;
+ e.time_base = 1e-10;
+
+ struct space s;
+ s.width[0] = 0.2;
+ s.width[1] = 0.2;
+ s.width[2] = 0.2;
+ e.s = &s;
+
+ struct gravity_props props;
+ props.a_smooth = 1.25;
+ e.gravity_properties = &props;
+
+ struct runner r;
+ bzero(&r, sizeof(struct runner));
+ r.e = &e;
+
+ const double rlr = props.a_smooth * s.width[0];
+
+ /* Init the cache for gravity interaction */
+ gravity_cache_init(&r.ci_gravity_cache, num_tests * 2);
+
+ /* Let's create one cell with a massive particle and a bunch of test particles
+ */
+ struct cell c;
+ bzero(&c, sizeof(struct cell));
+ c.width[0] = 1.;
+ c.width[1] = 1.;
+ c.width[2] = 1.;
+ c.loc[0] = 0.;
+ c.loc[1] = 0.;
+ c.loc[2] = 0.;
+ c.gcount = 1 + num_tests;
+ c.ti_old_gpart = 8;
+ c.ti_gravity_end_min = 8;
+ c.ti_gravity_end_max = 8;
+
+ if (posix_memalign((void **)&c.gparts, gpart_align,
+ c.gcount * sizeof(struct gpart)) != 0)
+ error("Impossible to allocate memory for the gparts.");
+ bzero(c.gparts, c.gcount * sizeof(struct gpart));
+
+ /* Create the massive particle */
+ c.gparts[0].x[0] = 0.;
+ c.gparts[0].x[1] = 0.5;
+ c.gparts[0].x[2] = 0.5;
+ c.gparts[0].mass = 1.;
+ c.gparts[0].epsilon = eps;
+ c.gparts[0].time_bin = 1;
+ c.gparts[0].type = swift_type_dark_matter;
+ c.gparts[0].id_or_neg_offset = 1;
+#ifdef SWIFT_DEBUG_CHECKS
+ c.gparts[0].ti_drift = 8;
+#endif
+
+ /* Create the mass-less particles */
+ for (int n = 1; n < num_tests + 1; ++n) {
+
+ struct gpart *gp = &c.gparts[n];
+
+ gp->x[0] = n / ((double)num_tests);
+ gp->x[1] = 0.5;
+ gp->x[2] = 0.5;
+ gp->mass = 0.;
+ gp->epsilon = eps;
+ gp->time_bin = 1;
+ gp->type = swift_type_dark_matter;
+ gp->id_or_neg_offset = n + 1;
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->ti_drift = 8;
+#endif
+ }
+
+ /* Now compute the forces */
+ runner_doself_grav_pp_truncated(&r, &c);
+
+ /* Verify everything */
+ for (int n = 1; n < num_tests + 1; ++n) {
+ const struct gpart *gp = &c.gparts[n];
+
+ double pot_true = potential(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
+ double acc_true =
+ acceleration(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
+
+ // message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_true);
+
+ check_value(gp->potential, pot_true, "potential");
+ check_value(gp->a_grav[0], acc_true, "acceleration");
+ }
+
+ free(c.gparts);
+ return 0;
+}
diff --git a/tests/tolerance_27_perturbed_h2.dat b/tests/tolerance_27_perturbed_h2.dat
index 781531d6a0d180d58ab74b9ef12efb927ccc733d..882554741734aeb270427b62580d6077907056ad 100644
--- a/tests/tolerance_27_perturbed_h2.dat
+++ b/tests/tolerance_27_perturbed_h2.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-6 1e-4 5e-4 1.5e-2 1.4e-5 3e-6 3e-6 1e-5
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.5e-6 2.5e-2 1e-5 5.86e-3 4.96e-4 3e-3 4.5e-3 3e-3
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.5e-6 2.5e-2 1e-5 5.86e-3 1.17e-3 3e-3 5.65e-3 3e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e0 1e-6 4e-6 4e-6 4e-6