diff --git a/examples/DiscPatch/HydroStatic/disc-patch-icc.yml b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml
index 9d1fe1b0948594f79c401e50308166313aa39924..6f17cfbb1e0125faf8e47fe4e9e55bfdf4df7b71 100644
--- a/examples/DiscPatch/HydroStatic/disc-patch-icc.yml
+++ b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml
@@ -39,8 +39,8 @@ InitialConditions:
DiscPatchPotential:
surface_density: 10.
scale_height: 100.
- z_disc: 400.
- z_trunc: 300.
- z_max: 350.
+ x_disc: 400.
+ x_trunc: 300.
+ x_max: 350.
timestep_mult: 0.03
growth_time: 5.
diff --git a/examples/DiscPatch/HydroStatic/disc-patch.yml b/examples/DiscPatch/HydroStatic/disc-patch.yml
index 8ed78eaf8fffbb731cab4f583c3b3e8acee2b7be..8816bc17ca526d01b7abcf55bb43287bbb36224a 100644
--- a/examples/DiscPatch/HydroStatic/disc-patch.yml
+++ b/examples/DiscPatch/HydroStatic/disc-patch.yml
@@ -39,7 +39,7 @@ InitialConditions:
DiscPatchPotential:
surface_density: 10.
scale_height: 100.
- z_disc: 400.
- z_trunc: 300.
- z_max: 380.
+ x_disc: 400.
+ x_trunc: 300.
+ x_max: 380.
timestep_mult: 0.03
diff --git a/examples/DiscPatch/HydroStatic/makeIC.py b/examples/DiscPatch/HydroStatic/makeIC.py
index 423dc8a44714a6e710b0c037f8a33d2120fbdcc8..11b482059b494fc9a6b9447fdfe2e7ec543d52ff 100644
--- a/examples/DiscPatch/HydroStatic/makeIC.py
+++ b/examples/DiscPatch/HydroStatic/makeIC.py
@@ -1,22 +1,24 @@
###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
- # Tom Theuns (tom.theuns@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 file is part of SWIFT.
+# Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+# Tom Theuns (tom.theuns@durham.ac.uk)
+# 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+# Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
import h5py
import sys
@@ -37,16 +39,16 @@ import random
# Size of the patch -- side_length
# Parameters of the gas disc
-surface_density = 10.
+surface_density = 10.
scale_height = 100.
gas_gamma = 5./3.
# Parameters of the problem
-z_factor = 2
+x_factor = 2
side_length = 400.
# File
-fileName = "Disc-Patch.hdf5"
+fileName = "Disc-Patch.hdf5"
####################################################################
@@ -64,30 +66,33 @@ unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
unit_velocity_in_cgs = (1e5)
unit_time_in_cgs = unit_length_in_cgs / unit_velocity_in_cgs
-print "UnitMass_in_cgs: %.5e"%unit_mass_in_cgs
+print "UnitMass_in_cgs: %.5e"%unit_mass_in_cgs
print "UnitLength_in_cgs: %.5e"%unit_length_in_cgs
print "UnitVelocity_in_cgs: %.5e"%unit_velocity_in_cgs
print "UnitTime_in_cgs: %.5e"%unit_time_in_cgs
print ""
# Derived units
-const_G = NEWTON_GRAVITY_CGS * unit_mass_in_cgs * unit_time_in_cgs**2 * unit_length_in_cgs**-3
+const_G = NEWTON_GRAVITY_CGS * unit_mass_in_cgs * unit_time_in_cgs**2 * \
+ unit_length_in_cgs**-3
const_mp = PROTON_MASS_IN_CGS * unit_mass_in_cgs**-1
-const_kb = BOLTZMANN_IN_CGS * unit_mass_in_cgs**-1 * unit_length_in_cgs**-2 * unit_time_in_cgs**2
+const_kb = BOLTZMANN_IN_CGS * unit_mass_in_cgs**-1 * unit_length_in_cgs**-2 * \
+ unit_time_in_cgs**2
print "--- Some constants [internal units] ---"
-print "G_Newton: %.5e"%const_G
-print "m_proton: %.5e"%const_mp
-print "k_boltzmann: %.5e"%const_kb
+print "G_Newton: %.5e"%const_G
+print "m_proton: %.5e"%const_mp
+print "k_boltzmann: %.5e"%const_kb
print ""
# derived quantities
-temp = math.pi * const_G * surface_density * scale_height * const_mp / const_kb
-u_therm = const_kb * temp / ((gas_gamma-1) * const_mp)
-v_disp = math.sqrt(2 * u_therm)
-soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma-1.)))
-t_dyn = math.sqrt(scale_height / (const_G * surface_density))
-t_cross = scale_height / soundspeed
+temp = math.pi * const_G * surface_density * scale_height * const_mp / \
+ const_kb
+u_therm = const_kb * temp / ((gas_gamma-1) * const_mp)
+v_disp = math.sqrt(2 * u_therm)
+soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma-1.)))
+t_dyn = math.sqrt(scale_height / (const_G * surface_density))
+t_cross = scale_height / soundspeed
print "--- Properties of the gas [internal units] ---"
print "Gas temperature: %.5e"%temp
@@ -101,18 +106,20 @@ print ""
# Problem properties
boxSize_x = side_length
boxSize_y = boxSize_x
-boxSize_z = boxSize_x * z_factor
+boxSize_z = boxSize_x
+boxSize_x *= x_factor
volume = boxSize_x * boxSize_y * boxSize_z
-M_tot = boxSize_x * boxSize_y * surface_density * math.tanh(boxSize_z / (2. * scale_height))
+M_tot = boxSize_y * boxSize_z * surface_density * \
+ math.tanh(boxSize_x / (2. * scale_height))
density = M_tot / volume
entropy = (gas_gamma - 1.) * u_therm / density**(gas_gamma - 1.)
print "--- Problem properties [internal units] ---"
-print "Box: [%.1f, %.1f, %.1f]"%(boxSize_x, boxSize_y, boxSize_z)
-print "Volume: %.5e"%volume
-print "Total mass: %.5e"%M_tot
-print "Density: %.5e"%density
-print "Entropy: %.5e"%entropy
+print "Box: [%.1f, %.1f, %.1f]"%(boxSize_x, boxSize_y, boxSize_z)
+print "Volume: %.5e"%volume
+print "Total mass: %.5e"%M_tot
+print "Density: %.5e"%density
+print "Entropy: %.5e"%entropy
print ""
####################################################################
@@ -123,34 +130,24 @@ one_glass_pos = infile["/PartType0/Coordinates"][:,:]
one_glass_h = infile["/PartType0/SmoothingLength"][:]
# Rescale to the problem size
-one_glass_pos *= boxSize_x
-one_glass_h *= boxSize_x
-
-#print min(one_glass_p[:,0]), max(one_glass_p[:,0])
-#print min(one_glass_p[:,1]), max(one_glass_p[:,1])
-#print min(one_glass_p[:,2]), max(one_glass_p[:,2])
+one_glass_pos *= side_length
+one_glass_h *= side_length
-# Now create enough copies to fill the volume in z
+# Now create enough copies to fill the volume in x
pos = np.copy(one_glass_pos)
h = np.copy(one_glass_h)
-for i in range(1, z_factor):
-
- one_glass_pos[:,2] += boxSize_x
-
+for i in range(1, x_factor):
+ one_glass_pos[:,0] += side_length
pos = np.append(pos, one_glass_pos, axis=0)
h = np.append(h, one_glass_h, axis=0)
-#print min(pos[:,0]), max(pos[:,0])
-#print min(pos[:,1]), max(pos[:,1])
-#print min(pos[:,2]), max(pos[:,2])
-
# Compute further properties of ICs
numPart = np.size(h)
mass = M_tot / numPart
print "--- Particle properties [internal units] ---"
print "Number part.: ", numPart
-print "Part. mass: %.5e"%mass
+print "Part. mass: %.5e"%mass
print ""
# Create additional arrays
@@ -159,7 +156,6 @@ mass = np.ones(numPart) * mass
vel = np.zeros((numPart, 3))
ids = 1 + np.linspace(0, numPart, numPart, endpoint=False)
-
####################################################################
# Create and write output file
@@ -169,7 +165,7 @@ file = h5py.File(fileName, 'w')
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = unit_length_in_cgs
-grp.attrs["Unit mass in cgs (U_M)"] = unit_mass_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = unit_mass_in_cgs
grp.attrs["Unit time in cgs (U_t)"] = unit_time_in_cgs
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
@@ -200,13 +196,12 @@ ds = grp0.create_dataset('SmoothingLength', (numPart,), 'f', data=h)
ds = grp0.create_dataset('InternalEnergy', (numPart,), 'f', data=u)
ds = grp0.create_dataset('ParticleIDs', (numPart, ), 'L', data=ids)
-
####################################################################
print "--- Runtime parameters (YAML file): ---"
print "DiscPatchPotential:surface_density: ", surface_density
print "DiscPatchPotential:scale_height: ", scale_height
-print "DiscPatchPotential:z_disc: ", boxSize_z / 2.
+print "DiscPatchPotential:x_disc: ", 0.5 * boxSize_x
print ""
print "--- Constant parameters: ---"
diff --git a/examples/DiscPatch/HydroStatic/plotSolution.py b/examples/DiscPatch/HydroStatic/plotSolution.py
index e70f595c645eaadc282e511bbc0d2c510025cc87..681f7d8ab3f2320b5de75e688edcb92efef9d883 100644
--- a/examples/DiscPatch/HydroStatic/plotSolution.py
+++ b/examples/DiscPatch/HydroStatic/plotSolution.py
@@ -1,6 +1,7 @@
################################################################################
# This file is part of SWIFT.
# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+# 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
@@ -34,9 +35,9 @@ import sys
# Parameters
surface_density = 10.
scale_height = 100.
-z_disc = 400.
-z_trunc = 300.
-z_max = 350.
+x_disc = 400.
+x_trunc = 300.
+x_max = 350.
utherm = 20.2678457288
gamma = 5. / 3.
@@ -50,14 +51,14 @@ if len(sys.argv) > 2:
# Get the analytic solution for the density
def get_analytic_density(x):
return 0.5 * surface_density / scale_height / \
- np.cosh( (x - z_disc) / scale_height )**2
+ np.cosh( (x - x_disc) / scale_height )**2
# Get the analytic solution for the (isothermal) pressure
def get_analytic_pressure(x):
return (gamma - 1.) * utherm * get_analytic_density(x)
# Get the data fields to plot from the snapshot file with the given name:
-# snapshot time, z-coord, density, pressure, velocity norm
+# snapshot time, x-coord, density, pressure, velocity norm
def get_data(name):
file = h5py.File(name, "r")
coords = np.array(file["/PartType0/Coordinates"])
@@ -69,7 +70,7 @@ def get_data(name):
vtot = np.sqrt( v[:,0]**2 + v[:,1]**2 + v[:,2]**2 )
- return float(file["/Header"].attrs["Time"]), coords[:,2], rho, P, vtot
+ return float(file["/Header"].attrs["Time"]), coords[:,0], rho, P, vtot
# scan the folder for snapshot files and plot all of them (within the requested
# range)
@@ -80,38 +81,38 @@ for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
print "processing", f, "..."
- zrange = np.linspace(0., 2. * z_disc, 1000)
- time, z, rho, P, v = get_data(f)
+ xrange = np.linspace(0., 2. * x_disc, 1000)
+ time, x, rho, P, v = get_data(f)
fig, ax = pl.subplots(3, 1, sharex = True)
- ax[0].plot(z, rho, "r.")
- ax[0].plot(zrange, get_analytic_density(zrange), "k-")
- ax[0].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
- ax[0].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
- ax[0].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
- ax[0].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
- ax[0].set_ylim(0., 1.2 * get_analytic_density(z_disc))
+ ax[0].plot(x, rho, "r.")
+ ax[0].plot(xrange, get_analytic_density(xrange), "k-")
+ ax[0].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+ ax[0].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+ ax[0].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[0].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[0].set_ylim(0., 1.2 * get_analytic_density(x_disc))
ax[0].set_ylabel("density")
- ax[1].plot(z, v, "r.")
- ax[1].plot(zrange, np.zeros(len(zrange)), "k-")
- ax[1].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
- ax[1].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
- ax[1].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
- ax[1].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
+ ax[1].plot(x, v, "r.")
+ ax[1].plot(xrange, np.zeros(len(xrange)), "k-")
+ ax[1].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+ ax[1].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+ ax[1].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[1].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
ax[1].set_ylim(-0.5, 10.)
ax[1].set_ylabel("velocity norm")
- ax[2].plot(z, P, "r.")
- ax[2].plot(zrange, get_analytic_pressure(zrange), "k-")
- ax[2].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
- ax[2].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
- ax[2].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
- ax[2].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
- ax[2].set_xlim(0., 2. * z_disc)
- ax[2].set_ylim(0., 1.2 * get_analytic_pressure(z_disc))
- ax[2].set_xlabel("z")
+ ax[2].plot(x, P, "r.")
+ ax[2].plot(xrange, get_analytic_pressure(xrange), "k-")
+ ax[2].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+ ax[2].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+ ax[2].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[2].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[2].set_xlim(0., 2. * x_disc)
+ ax[2].set_ylim(0., 1.2 * get_analytic_pressure(x_disc))
+ ax[2].set_xlabel("x")
ax[2].set_ylabel("pressure")
pl.suptitle("t = {0:.2f}".format(time))
diff --git a/examples/DiscPatch/HydroStatic/run.sh b/examples/DiscPatch/HydroStatic/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..e1f47ecad54e7e171d78b7da080d56579e985d1e
--- /dev/null
+++ b/examples/DiscPatch/HydroStatic/run.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_32.hdf5 ]
+then
+ echo "Fetching initial glass file for the disc patch example..."
+ ./getGlass.sh
+fi
+if [ ! -e Disc-Patch.hdf5 ]
+then
+ echo "Generating initial conditions for the disc patch example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../../swift -g -s -t 4 disc-patch-icc.yml 2>&1 | tee output.log
+
+python plotSolution.py
diff --git a/examples/DiscPatch/HydroStatic_1D/disc-patch-icc.yml b/examples/DiscPatch/HydroStatic_1D/disc-patch-icc.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6f17cfbb1e0125faf8e47fe4e9e55bfdf4df7b71
--- /dev/null
+++ b/examples/DiscPatch/HydroStatic_1D/disc-patch-icc.yml
@@ -0,0 +1,46 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.9885e33 # Grams
+ UnitLength_in_cgs: 3.08567758149e18 # Centimeters
+ UnitVelocity_in_cgs: 1e5 # 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: 968. # The end time of the simulation (in internal units).
+ dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 10. # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 12. # Time between statistics output
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: Disc-Patch # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 48. # Time difference between outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 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.
+ max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
+ h_max: 60. # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: Disc-Patch.hdf5 # The file to read
+
+# External potential parameters
+DiscPatchPotential:
+ surface_density: 10.
+ scale_height: 100.
+ x_disc: 400.
+ x_trunc: 300.
+ x_max: 350.
+ timestep_mult: 0.03
+ growth_time: 5.
diff --git a/examples/DiscPatch/HydroStatic_1D/makeIC.py b/examples/DiscPatch/HydroStatic_1D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..1589dfc8c73e5b9bf3c2cad4bcf3029654d9e67e
--- /dev/null
+++ b/examples/DiscPatch/HydroStatic_1D/makeIC.py
@@ -0,0 +1,194 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+# Tom Theuns (tom.theuns@durham.ac.uk)
+# 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+# Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+import h5py
+import sys
+import numpy as np
+import math
+import random
+
+# Generates a disc-patch in hydrostatic equilibrium
+#
+# See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
+#
+#
+# Disc parameters are: surface density -- sigma
+# scale height -- b
+# gas adiabatic index -- gamma
+#
+# Problem parameters are: Ratio height/width of the box -- z_factor
+# Size of the patch -- side_length
+
+# Parameters of the gas disc
+surface_density = 10.
+scale_height = 100.
+gas_gamma = 5./3.
+
+# Parameters of the problem
+x_factor = 2
+side_length = 400.
+numPart = 1000
+
+# File
+fileName = "Disc-Patch.hdf5"
+
+####################################################################
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS = 6.67408e-8
+SOLAR_MASS_IN_CGS = 1.9885e33
+PARSEC_IN_CGS = 3.08567758149e18
+PROTON_MASS_IN_CGS = 1.672621898e-24
+BOLTZMANN_IN_CGS = 1.38064852e-16
+YEAR_IN_CGS = 3.15569252e7
+
+# choice of units
+unit_length_in_cgs = (PARSEC_IN_CGS)
+unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
+unit_velocity_in_cgs = (1e5)
+unit_time_in_cgs = unit_length_in_cgs / unit_velocity_in_cgs
+
+print "UnitMass_in_cgs: %.5e"%unit_mass_in_cgs
+print "UnitLength_in_cgs: %.5e"%unit_length_in_cgs
+print "UnitVelocity_in_cgs: %.5e"%unit_velocity_in_cgs
+print "UnitTime_in_cgs: %.5e"%unit_time_in_cgs
+print ""
+
+# Derived units
+const_G = NEWTON_GRAVITY_CGS * unit_mass_in_cgs * unit_time_in_cgs**2 * \
+ unit_length_in_cgs**-3
+const_mp = PROTON_MASS_IN_CGS * unit_mass_in_cgs**-1
+const_kb = BOLTZMANN_IN_CGS * unit_mass_in_cgs**-1 * unit_length_in_cgs**-2 * \
+ unit_time_in_cgs**2
+
+print "--- Some constants [internal units] ---"
+print "G_Newton: %.5e"%const_G
+print "m_proton: %.5e"%const_mp
+print "k_boltzmann: %.5e"%const_kb
+print ""
+
+# derived quantities
+temp = math.pi * const_G * surface_density * scale_height * const_mp / \
+ const_kb
+u_therm = const_kb * temp / ((gas_gamma-1) * const_mp)
+v_disp = math.sqrt(2 * u_therm)
+soundspeed = math.sqrt(u_therm / (gas_gamma * (gas_gamma-1.)))
+t_dyn = math.sqrt(scale_height / (const_G * surface_density))
+t_cross = scale_height / soundspeed
+
+print "--- Properties of the gas [internal units] ---"
+print "Gas temperature: %.5e"%temp
+print "Gas thermal_energy: %.5e"%u_therm
+print "Dynamical time: %.5e"%t_dyn
+print "Sound crossing time: %.5e"%t_cross
+print "Gas sound speed: %.5e"%soundspeed
+print "Gas 3D vel_disp: %.5e"%v_disp
+print ""
+
+# Problem properties
+boxSize_x = side_length
+boxSize_x *= x_factor
+volume = boxSize_x
+M_tot = surface_density * math.tanh(boxSize_x / (2. * scale_height))
+density = M_tot / volume
+entropy = (gas_gamma - 1.) * u_therm / density**(gas_gamma - 1.)
+
+print "--- Problem properties [internal units] ---"
+print "Box: %.1f"%boxSize_x
+print "Volume: %.5e"%volume
+print "Total mass: %.5e"%M_tot
+print "Density: %.5e"%density
+print "Entropy: %.5e"%entropy
+print ""
+
+####################################################################
+
+# Now create enough copies to fill the volume in x
+pos = np.zeros((numPart, 3))
+h = np.zeros(numPart) + 2. * boxSize_x / numPart
+for i in range(numPart):
+ pos[i, 0] = (i + 0.5) * boxSize_x / numPart
+
+# Compute further properties of ICs
+mass = M_tot / numPart
+
+print "--- Particle properties [internal units] ---"
+print "Number part.: ", numPart
+print "Part. mass: %.5e"%mass
+print ""
+
+# Create additional arrays
+u = np.ones(numPart) * u_therm
+mass = np.ones(numPart) * mass
+vel = np.zeros((numPart, 3))
+ids = 1 + np.linspace(0, numPart, numPart, endpoint=False)
+
+####################################################################
+# Create and write output file
+
+#File
+file = h5py.File(fileName, 'w')
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = unit_length_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = unit_mass_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = unit_time_in_cgs
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [boxSize_x, 1., 1.]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["Dimension"] = 1
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+# write gas particles
+grp0 = file.create_group("/PartType0")
+
+ds = grp0.create_dataset('Coordinates', (numPart, 3), 'f', data=pos)
+ds = grp0.create_dataset('Velocities', (numPart, 3), 'f')
+ds = grp0.create_dataset('Masses', (numPart,), 'f', data=mass)
+ds = grp0.create_dataset('SmoothingLength', (numPart,), 'f', data=h)
+ds = grp0.create_dataset('InternalEnergy', (numPart,), 'f', data=u)
+ds = grp0.create_dataset('ParticleIDs', (numPart, ), 'L', data=ids)
+
+####################################################################
+
+print "--- Runtime parameters (YAML file): ---"
+print "DiscPatchPotential:surface_density: ", surface_density
+print "DiscPatchPotential:scale_height: ", scale_height
+print "DiscPatchPotential:x_disc: ", 0.5 * boxSize_x
+print ""
+
+print "--- Constant parameters: ---"
+print "const_isothermal_internal_energy: %ef"%u_therm
diff --git a/examples/DiscPatch/HydroStatic_1D/plotSolution.py b/examples/DiscPatch/HydroStatic_1D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..681f7d8ab3f2320b5de75e688edcb92efef9d883
--- /dev/null
+++ b/examples/DiscPatch/HydroStatic_1D/plotSolution.py
@@ -0,0 +1,121 @@
+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+# 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/>.
+#
+################################################################################
+
+##
+# This script plots the Disc-Patch_*.hdf5 snapshots.
+# It takes two (optional) parameters: the counter value of the first and last
+# snapshot to plot (default: 0 21).
+##
+
+import numpy as np
+import h5py
+import matplotlib
+matplotlib.use("Agg")
+import pylab as pl
+import glob
+import sys
+
+# Parameters
+surface_density = 10.
+scale_height = 100.
+x_disc = 400.
+x_trunc = 300.
+x_max = 350.
+utherm = 20.2678457288
+gamma = 5. / 3.
+
+start = 0
+stop = 21
+if len(sys.argv) > 1:
+ start = int(sys.argv[1])
+if len(sys.argv) > 2:
+ stop = int(sys.argv[2])
+
+# Get the analytic solution for the density
+def get_analytic_density(x):
+ return 0.5 * surface_density / scale_height / \
+ np.cosh( (x - x_disc) / scale_height )**2
+
+# Get the analytic solution for the (isothermal) pressure
+def get_analytic_pressure(x):
+ return (gamma - 1.) * utherm * get_analytic_density(x)
+
+# Get the data fields to plot from the snapshot file with the given name:
+# snapshot time, x-coord, density, pressure, velocity norm
+def get_data(name):
+ file = h5py.File(name, "r")
+ coords = np.array(file["/PartType0/Coordinates"])
+ rho = np.array(file["/PartType0/Density"])
+ u = np.array(file["/PartType0/InternalEnergy"])
+ v = np.array(file["/PartType0/Velocities"])
+
+ P = (gamma - 1.) * rho * u
+
+ vtot = np.sqrt( v[:,0]**2 + v[:,1]**2 + v[:,2]**2 )
+
+ return float(file["/Header"].attrs["Time"]), coords[:,0], rho, P, vtot
+
+# scan the folder for snapshot files and plot all of them (within the requested
+# range)
+for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
+ num = int(f[-8:-5])
+ if num < start or num > stop:
+ continue
+
+ print "processing", f, "..."
+
+ xrange = np.linspace(0., 2. * x_disc, 1000)
+ time, x, rho, P, v = get_data(f)
+
+ fig, ax = pl.subplots(3, 1, sharex = True)
+
+ ax[0].plot(x, rho, "r.")
+ ax[0].plot(xrange, get_analytic_density(xrange), "k-")
+ ax[0].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+ ax[0].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+ ax[0].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[0].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[0].set_ylim(0., 1.2 * get_analytic_density(x_disc))
+ ax[0].set_ylabel("density")
+
+ ax[1].plot(x, v, "r.")
+ ax[1].plot(xrange, np.zeros(len(xrange)), "k-")
+ ax[1].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+ ax[1].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+ ax[1].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[1].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[1].set_ylim(-0.5, 10.)
+ ax[1].set_ylabel("velocity norm")
+
+ ax[2].plot(x, P, "r.")
+ ax[2].plot(xrange, get_analytic_pressure(xrange), "k-")
+ ax[2].plot([x_disc - x_max, x_disc - x_max], [0, 10], "k--", alpha=0.5)
+ ax[2].plot([x_disc + x_max, x_disc + x_max], [0, 10], "k--", alpha=0.5)
+ ax[2].plot([x_disc - x_trunc, x_disc - x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[2].plot([x_disc + x_trunc, x_disc + x_trunc], [0, 10], "k--", alpha=0.5)
+ ax[2].set_xlim(0., 2. * x_disc)
+ ax[2].set_ylim(0., 1.2 * get_analytic_pressure(x_disc))
+ ax[2].set_xlabel("x")
+ ax[2].set_ylabel("pressure")
+
+ pl.suptitle("t = {0:.2f}".format(time))
+
+ pl.savefig("{name}.png".format(name = f[:-5]))
+ pl.close()
diff --git a/examples/DiscPatch/HydroStatic_1D/run.sh b/examples/DiscPatch/HydroStatic_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..e9d073a6cc7a06ec9ebd9fdb556c44778d32c7f4
--- /dev/null
+++ b/examples/DiscPatch/HydroStatic_1D/run.sh
@@ -0,0 +1,13 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e Disc-Patch.hdf5 ]
+then
+ echo "Generating initial conditions for the disc patch example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../../swift -g -s -t 4 disc-patch-icc.yml 2>&1 | tee output.log
+
+python plotSolution.py
diff --git a/examples/plot_threadpool.py b/examples/plot_threadpool.py
deleted file mode 100755
index 68b038ed406fdac8c4fee3e250fbb397bb51ae21..0000000000000000000000000000000000000000
--- a/examples/plot_threadpool.py
+++ /dev/null
@@ -1,271 +0,0 @@
-#!/usr/bin/env python
-"""
-Usage:
- plot_threadpool.py [options] input.dat output.png
-
-where input.dat is a threadpool info file for a step. Use the '-Y interval'
-flag of the swift command to create these. The output plot will be called
-'output.png'. The --limit option can be used to produce plots with the same
-time span and the --expand option to expand each thread line into '*expand'
-lines, so that adjacent tasks of the same type can be distinguished. Other
-options can be seen using the --help flag.
-
-This file is part of SWIFT.
-Copyright (c) 2015 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- (c) 2017 Peter W. Draper (p.w.draper@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 matplotlib
-matplotlib.use("Agg")
-import matplotlib.collections as collections
-import matplotlib.ticker as plticker
-import pylab as pl
-import sys
-import argparse
-
-# Handle the command line.
-parser = argparse.ArgumentParser(description="Plot threadpool function graphs")
-
-parser.add_argument("input", help="Threadpool data file (-Y output)")
-parser.add_argument("outpng", help="Name for output graphic file (PNG)")
-parser.add_argument("-l", "--limit", dest="limit",
- help="Upper time limit in millisecs (def: depends on data)",
- default=0, type=int)
-parser.add_argument("-e", "--expand", dest="expand",
- help="Thread expansion factor (def: 1)",
- default=1, type=int)
-parser.add_argument("--height", dest="height",
- help="Height of plot in inches (def: 4)",
- default=4., type=float)
-parser.add_argument("--width", dest="width",
- help="Width of plot in inches (def: 16)",
- default=16., type=float)
-parser.add_argument("--nolegend", dest="nolegend",
- help="Whether to show the legend (def: False)",
- default=False, action="store_true")
-parser.add_argument("-v", "--verbose", dest="verbose",
- help="Show colour assignments and other details (def: False)",
- default=False, action="store_true")
-
-args = parser.parse_args()
-infile = args.input
-outpng = args.outpng
-delta_t = args.limit
-expand = args.expand
-
-# Basic plot configuration.
-PLOT_PARAMS = {"axes.labelsize": 10,
- "axes.titlesize": 10,
- "font.size": 12,
- "legend.fontsize": 12,
- "xtick.labelsize": 10,
- "ytick.labelsize": 10,
- "figure.figsize" : (args.width, args.height),
- "figure.subplot.left" : 0.03,
- "figure.subplot.right" : 0.995,
- "figure.subplot.bottom" : 0.09,
- "figure.subplot.top" : 0.99,
- "figure.subplot.wspace" : 0.,
- "figure.subplot.hspace" : 0.,
- "lines.markersize" : 6,
- "lines.linewidth" : 3.
- }
-pl.rcParams.update(PLOT_PARAMS)
-
-# A number of colours for the various types. Recycled when there are
-# more task types than colours...
-colours = ["cyan", "lightgray", "darkblue", "yellow", "tan", "dodgerblue",
- "sienna", "aquamarine", "bisque", "blue", "green", "lightgreen",
- "brown", "purple", "moccasin", "olivedrab", "chartreuse",
- "darksage", "darkgreen", "green", "mediumseagreen",
- "mediumaquamarine", "darkslategrey", "mediumturquoise",
- "black", "cadetblue", "skyblue", "red", "slategray", "gold",
- "slateblue", "blueviolet", "mediumorchid", "firebrick",
- "magenta", "hotpink", "pink", "orange", "lightgreen"]
-maxcolours = len(colours)
-
-# Read header. First two lines.
-with open(infile) as infid:
- head = [next(infid) for x in xrange(2)]
-header = head[1][2:].strip()
-header = eval(header)
-nthread = int(header['num_threads']) + 1
-CPU_CLOCK = float(header['cpufreq']) / 1000.0
-print "Number of threads: ", nthread
-if args.verbose:
- print "CPU frequency:", CPU_CLOCK * 1000.0
-
-# Read input.
-data = pl.genfromtxt(infile, dtype=None, delimiter=" ")
-
-# Mixed types, so need to separate.
-tics = []
-tocs = []
-funcs = []
-threads = []
-chunks = []
-for i in data:
- if i[0] != "#":
- funcs.append(i[0])
- if i[1] < 0:
- threads.append(nthread-1)
- else:
- threads.append(i[1])
- chunks.append(i[2])
- tics.append(i[3])
- tocs.append(i[4])
-tics = pl.array(tics)
-tocs = pl.array(tocs)
-funcs = pl.array(funcs)
-threads = pl.array(threads)
-chunks = pl.array(chunks)
-
-
-# Recover the start and end time
-tic_step = min(tics)
-toc_step = max(tocs)
-
-# Not known.
-
-# Calculate the time range, if not given.
-delta_t = delta_t * CPU_CLOCK
-if delta_t == 0:
- dt = toc_step - tic_step
- if dt > delta_t:
- delta_t = dt
- print "Data range: ", delta_t / CPU_CLOCK, "ms"
-
-# Once more doing the real gather and plots this time.
-start_t = float(tic_step)
-tics -= tic_step
-tocs -= tic_step
-end_t = (toc_step - start_t) / CPU_CLOCK
-
-# Get all "task" names and assign colours.
-TASKTYPES = pl.unique(funcs)
-print TASKTYPES
-
-# Set colours of task/subtype.
-TASKCOLOURS = {}
-ncolours = 0
-for task in TASKTYPES:
- TASKCOLOURS[task] = colours[ncolours]
- ncolours = (ncolours + 1) % maxcolours
-
-# For fiddling with colours...
-if args.verbose:
- print "#Selected colours:"
- for task in sorted(TASKCOLOURS.keys()):
- print "# " + task + ": " + TASKCOLOURS[task]
- for task in sorted(SUBCOLOURS.keys()):
- print "# " + task + ": " + SUBCOLOURS[task]
-
-tasks = {}
-tasks[-1] = []
-for i in range(nthread*expand):
- tasks[i] = []
-
-# Counters for each thread when expanding.
-ecounter = []
-for i in range(nthread):
- ecounter.append(0)
-
-for i in range(len(threads)):
- thread = threads[i]
-
- # Expand to cover extra lines if expanding.
- ethread = thread * expand + (ecounter[thread] % expand)
- ecounter[thread] = ecounter[thread] + 1
- thread = ethread
-
- tasks[thread].append({})
- tasks[thread][-1]["type"] = funcs[i]
- tic = tics[i] / CPU_CLOCK
- toc = tocs[i] / CPU_CLOCK
- tasks[thread][-1]["tic"] = tic
- tasks[thread][-1]["toc"] = toc
- tasks[thread][-1]["colour"] = TASKCOLOURS[funcs[i]]
-
-# Use expanded threads from now on.
-nthread = nthread * expand
-
-typesseen = []
-fig = pl.figure()
-ax = fig.add_subplot(1,1,1)
-ax.set_xlim(-delta_t * 0.01 / CPU_CLOCK, delta_t * 1.01 / CPU_CLOCK)
-ax.set_ylim(0, nthread)
-
-# Fake thread is used to colour the whole range, do that first.
-tictocs = []
-colours = []
-j = 0
-for task in tasks[nthread - expand]:
- tictocs.append((task["tic"], task["toc"] - task["tic"]))
- colours.append(task["colour"])
-ax.broken_barh(tictocs, [0,(nthread-1)], facecolors = colours, linewidth=0, alpha=0.15)
-
-# And we don't plot the fake thread.
-nthread = nthread - expand
-for i in range(nthread):
-
- # Collect ranges and colours into arrays.
- tictocs = []
- colours = []
- j = 0
- for task in tasks[i]:
- tictocs.append((task["tic"], task["toc"] - task["tic"]))
- colours.append(task["colour"])
-
- # Legend support, collections don't add to this.
- qtask = task["type"]
- if qtask not in typesseen:
- pl.plot([], [], color=task["colour"], label=qtask)
- typesseen.append(qtask)
-
- # Now plot.
- ax.broken_barh(tictocs, [i+0.05,0.90], facecolors = colours, linewidth=0)
-
-# Legend and room for it.
-nrow = len(typesseen) / 5
-if not args.nolegend:
- ax.fill_between([0, 0], nthread+0.5, nthread + nrow + 0.5, facecolor="white")
- ax.set_ylim(0, nthread + 0.5)
- ax.legend(loc=1, shadow=True, bbox_to_anchor=(0., 1.05 ,1., 0.2), mode="expand", ncol=5)
- box = ax.get_position()
- ax.set_position([box.x0, box.y0, box.width, box.height*0.8])
-
-# Start and end of time-step
-ax.plot([0, 0], [0, nthread + nrow + 1], 'k--', linewidth=1)
-ax.plot([end_t, end_t], [0, nthread + nrow + 1], 'k--', linewidth=1)
-
-ax.set_xlabel("Wall clock time [ms]", labelpad=0.)
-if expand == 1:
- ax.set_ylabel("Thread ID", labelpad=0 )
-else:
- ax.set_ylabel("Thread ID * " + str(expand), labelpad=0 )
-ax.set_yticks(pl.array(range(nthread)), True)
-
-loc = plticker.MultipleLocator(base=expand)
-ax.yaxis.set_major_locator(loc)
-ax.grid(True, which='major', axis="y", linestyle="-")
-
-pl.show()
-pl.savefig(outpng)
-print "Graphics done, output written to", outpng
-
-sys.exit(0)
diff --git a/examples/process_plot_threadpool b/examples/process_plot_threadpool
deleted file mode 100755
index d7bef48d313d50a2627de69263202fa4782efa74..0000000000000000000000000000000000000000
--- a/examples/process_plot_threadpool
+++ /dev/null
@@ -1,96 +0,0 @@
-#!/bin/bash
-#
-# Usage:
-# process_plot_threadpool nprocess time-range-ms
-#
-# Description:
-# Process all the threadpool info files in the current directory
-# creating function graphs for steps and threads.
-#
-# The input files are created by a run using the "-Y interval" flag and
-# should be named "threadpool_info-step<n>.dat" in the current directory.
-# All located files will be processed using "nprocess" concurrent
-# processes and all plots will have the given time range. An output
-# HTML file "index.html" will be created to view all the plots.
-#
-#
-# This file is part of SWIFT:
-#
-# Copyright (C) 2017 Peter W. Draper (p.w.draper@durham.ac.uk)
-# All Rights Reserved.
-#
-# 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/>.
-
-# Handle command-line
-if test "$2" == ""; then
- echo "Usage: $0 nprocess time-range-ms"
- exit 1
-fi
-NPROCS=$1
-TIMERANGE=$2
-
-# Find all thread info files. Use version sort to get into correct order.
-files=$(ls -v threadpool_info-step*.dat)
-if test $? != 0; then
- echo "Failed to find any threadpool info files"
- exit 1
-fi
-
-# Construct list of names, the step no and names for the graphics.
-list=""
-for f in $files; do
- s=$(echo $f| sed 's,threadpool_info-step\(.*\).dat,\1,')
- list="$list $f $s poolstep${s}r"
-done
-
-# And process them,
-echo "Processing threadpool info files..."
-echo $list | xargs -P $NPROCS -n 3 /bin/bash -c "./plot_threadpool.py --expand 3 --limit $TIMERANGE --width 16 --height 4 \$0 \$2 "
-
-echo "Writing output index.html file"
-# Construct document - serial.
-cat <<EOF > index.html
- <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
-<html>
- <head>
- <title>SWIFT task graphs</title>
- </head>
- <body>
- <h1>SWIFT task graphs</h1>
-EOF
-
-echo $list | xargs -n 3 | while read f s g; do
- cat <<EOF >> index.html
-<h2>Step $s</h2>
-EOF
- cat <<EOF >> index.html
-<a href="poolstep${s}r${i}.html"><img src="poolstep${s}r${i}.png" width=400px/></a>
-EOF
- cat <<EOF > poolstep${s}r${i}.html
- <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
-<html>
-<body>
-<img src="poolstep${s}r${i}.png">
-<pre>
-EOF
-done
-
-cat <<EOF >> index.html
- </body>
-</html>
-EOF
-
-echo "Finished"
-
-exit
diff --git a/src/Makefile.am b/src/Makefile.am
index 2ddcdb0908201c65053d7cc5380a4217277b5c13..3784886278544a1c6b606689bfbae0d6fa52cf19 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -86,6 +86,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
hydro/Gizmo/hydro_slope_limiters_cell.h \
hydro/Gizmo/hydro_slope_limiters_face.h \
hydro/Gizmo/hydro_slope_limiters.h \
+ hydro/Gizmo/hydro_flux_limiters.h \
hydro/Gizmo/hydro_unphysical.h \
hydro/Gizmo/hydro_velocities.h \
hydro/Shadowswift/hydro_debug.h \
diff --git a/src/cell.c b/src/cell.c
index 7edee4d9ae1c02311af295bb6eba5c846a94259c..8386161a15cc110b9abc92e3539776cba0c918ec 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1745,9 +1745,9 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (l == NULL) error("Missing link to send_xv task.");
scheduler_activate(s, l->t);
- /* Drift the cell which will be sent at the level at which it is sent,
- i.e. drift the cell specified in the send task (l->t) itself. */
- cell_activate_drift_part(l->t->ci, s);
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(cj, s);
if (cell_is_active(cj, e)) {
struct link *l = NULL;
@@ -1802,9 +1802,9 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (l == NULL) error("Missing link to send_xv task.");
scheduler_activate(s, l->t);
- /* Drift the cell which will be sent at the level at which it is sent,
- i.e. drift the cell specified in the send task (l->t) itself. */
- cell_activate_drift_part(l->t->ci, s);
+ /* Drift the cell which will be sent; note that not all sent
+ particles will be drifted, only those that are needed. */
+ cell_activate_drift_part(ci, s);
if (cell_is_active(ci, e)) {
diff --git a/src/const.h b/src/const.h
index f5232de0b2b4aba6701c07f1e5f91d5be3d5b6cc..c8060a2be51468a791e192a65a74f1a4d9bc8e30 100644
--- a/src/const.h
+++ b/src/const.h
@@ -49,6 +49,9 @@
#define SLOPE_LIMITER_PER_FACE
#define SLOPE_LIMITER_CELL_WIDE
+/* Types of flux limiter to use (GIZMO_SPH only) */
+#define GIZMO_FLUX_LIMITER
+
/* Options to control the movement of particles for GIZMO_SPH. */
/* This option disables particle movement */
//#define GIZMO_FIX_PARTICLES
diff --git a/src/engine.c b/src/engine.c
index 6bd2fb5ade526e4ff7fd363fa655fd096afcd827..ee48d6565dc62f725b0159dc44e8d9d92d7a4adf 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2402,6 +2402,21 @@ void engine_make_gravityrecursive_tasks(struct engine *e) {
/* } */
}
+void engine_check_sort_tasks(struct engine *e, struct cell *c) {
+
+ /* Find the parent sort task, if any, and copy its flags. */
+ if (c->sorts != NULL) {
+ struct cell *parent = c->parent;
+ while (parent != NULL && parent->sorts == NULL) parent = parent->parent;
+ if (parent != NULL) c->sorts->flags |= parent->sorts->flags;
+ }
+
+ /* Recurse? */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) engine_check_sort_tasks(e, c->progeny[k]);
+}
+
/**
* @brief Fill the #space's task list.
*
@@ -2480,6 +2495,9 @@ void engine_maketasks(struct engine *e) {
for (int k = 0; k < nr_cells; k++)
engine_make_hierarchical_tasks(e, &cells[k]);
+ /* Append hierarchical tasks to each cell. */
+ for (int k = 0; k < nr_cells; k++) engine_check_sort_tasks(e, &cells[k]);
+
/* Run through the tasks and make force tasks for each density task.
Each force task depends on the cell ghosts and unlocks the kick task
of its super-cell. */
@@ -2774,7 +2792,7 @@ int engine_marktasks(struct engine *e) {
/* Run through the tasks and mark as skip or not. */
size_t extra_data[3] = {(size_t)e, rebuild_space, (size_t)&e->sched};
threadpool_map(&e->threadpool, engine_marktasks_mapper, s->tasks, s->nr_tasks,
- sizeof(struct task), 0, extra_data);
+ sizeof(struct task), 10000, extra_data);
rebuild_space = extra_data[1];
if (e->verbose)
@@ -3678,7 +3696,7 @@ void engine_drift_all(struct engine *e) {
#endif
threadpool_map(&e->threadpool, engine_do_drift_all_mapper, e->s->cells_top,
- e->s->nr_cells, sizeof(struct cell), 0, e);
+ e->s->nr_cells, sizeof(struct cell), 1, e);
/* Synchronize particle positions */
space_synchronize_particle_positions(e->s);
@@ -3730,7 +3748,7 @@ void engine_drift_top_multipoles(struct engine *e) {
const ticks tic = getticks();
threadpool_map(&e->threadpool, engine_do_drift_top_multipoles_mapper,
- e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 0, e);
+ e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 10, e);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time. */
@@ -3768,7 +3786,7 @@ void engine_reconstruct_multipoles(struct engine *e) {
const ticks tic = getticks();
threadpool_map(&e->threadpool, engine_do_reconstruct_multipoles_mapper,
- e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 0, e);
+ e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 10, e);
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
diff --git a/src/engine.h b/src/engine.h
index 79ff45a69a39ee7d4d45589df51ce2d53f810fda..c7aaa08b57eb2e61b311deae7a7ccb102f7e3cf8 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -77,7 +77,7 @@ extern const char *engine_policy_names[];
#define engine_queue_scale 1.2
#define engine_maxtaskspercell 96
#define engine_maxproxies 64
-#define engine_tasksreweight 1
+#define engine_tasksreweight 10
#define engine_parts_size_grow 1.05
#define engine_redistribute_alloc_margin 1.2
#define engine_default_energy_file_name "energy"
diff --git a/src/gravity.c b/src/gravity.c
index 49bbaca39b5278009543204a0d9f5e72d69806c4..97b2955b32e1513c3d86d1d1f4da2169130feb77 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -207,7 +207,7 @@ void gravity_exact_force_compute(struct space *s, const struct engine *e) {
data.const_G = e->physical_constants->const_newton_G;
threadpool_map(&s->e->threadpool, gravity_exact_force_compute_mapper,
- s->gparts, s->nr_gparts, sizeof(struct gpart), 0, &data);
+ s->gparts, s->nr_gparts, sizeof(struct gpart), 1000, &data);
message("Computed exact gravity for %d gparts (took %.3f %s). ",
data.counter_global, clocks_from_ticks(getticks() - tic),
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 49e748c71b4544ff6148e6910f66381294608e8f..2c2f54699bb380a491edf61a83ad8a031572c86c 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -49,17 +49,21 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return CFL_condition;
#endif
- if (p->timestepvars.vmax == 0.) {
- /* vmax can be zero in vacuum cells that only have vacuum neighbours */
- /* in this case, the time step should be limited by the maximally
- allowed time step. Since we do not know what that value is here, we set
- the time step to a very large value */
- return FLT_MAX;
- } else {
- const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
- hydro_dimension_inv);
- return 2. * CFL_condition * psize / fabsf(p->timestepvars.vmax);
+ float vrel[3];
+ vrel[0] = p->primitives.v[0] - xp->v_full[0];
+ vrel[1] = p->primitives.v[1] - xp->v_full[1];
+ vrel[2] = p->primitives.v[2] - xp->v_full[2];
+ float vmax =
+ sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
+ sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+ vmax = max(vmax, p->timestepvars.vmax);
+ const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
+ hydro_dimension_inv);
+ float dt = FLT_MAX;
+ if (vmax > 0.) {
+ dt = psize / vmax;
}
+ return CFL_condition * dt;
}
/**
@@ -421,7 +425,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp) {
/* Initialize time step criterion variables */
- p->timestepvars.vmax = 0.0f;
+ p->timestepvars.vmax = 0.;
/* Set the actual velocity of the particle */
hydro_velocities_prepare_force(p, xp);
@@ -638,24 +642,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
a_grav[1] = p->gpart->a_grav[1];
a_grav[2] = p->gpart->a_grav[2];
- /* Store the gravitational acceleration for later use. */
- /* This is used for the prediction step. */
- p->gravity.old_a[0] = a_grav[0];
- p->gravity.old_a[1] = a_grav[1];
- p->gravity.old_a[2] = a_grav[2];
-
/* Make sure the gpart knows the mass has changed. */
p->gpart->mass = p->conserved.mass;
- /* Kick the momentum for half a time step */
- /* Note that this also affects the particle movement, as the velocity for
- the particles is set after this. */
- p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
- p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
- p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
-
#if !defined(EOS_ISOTHERMAL_GAS)
- /* This part still needs to be tested! */
+ /* If the energy needs to be updated, we need to do it before the momentum
+ is updated, as the old value of the momentum enters the equations. */
p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
p->conserved.momentum[1] * a_grav[1] +
p->conserved.momentum[2] * a_grav[2]);
@@ -664,6 +656,13 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
a_grav[1] * p->gravity.mflux[1] +
a_grav[2] * p->gravity.mflux[2]);
#endif
+
+ /* Kick the momentum for half a time step */
+ /* Note that this also affects the particle movement, as the velocity for
+ the particles is set after this. */
+ p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
+ p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
+ p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
}
/* reset fluxes */
diff --git a/src/hydro/Gizmo/hydro_debug.h b/src/hydro/Gizmo/hydro_debug.h
index a05ff9a7d96f04ca3354235540adc31386a2d2e3..17e7f8a08570e355a701f8e165ee8af745fa34ab 100644
--- a/src/hydro/Gizmo/hydro_debug.h
+++ b/src/hydro/Gizmo/hydro_debug.h
@@ -46,7 +46,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"volume=%.3e, "
"matrix_E=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]]}, "
"timestepvars={"
- "vmax=%.3e}, "
+ "vmax=%.3e},"
"density={"
"div_v=%.3e, "
"wcount_dh=%.3e, "
diff --git a/src/hydro/Gizmo/hydro_flux_limiters.h b/src/hydro/Gizmo/hydro_flux_limiters.h
new file mode 100644
index 0000000000000000000000000000000000000000..dc91cf2808e02d903ff97efddc20c164db9c954e
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_flux_limiters.h
@@ -0,0 +1,81 @@
+
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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_HYDRO_FLUX_LIMITERS_H
+#define SWIFT_HYDRO_FLUX_LIMITERS_H
+
+#ifdef GIZMO_FLUX_LIMITER
+
+#define HYDRO_FLUX_LIMITER_IMPLEMENTATION "GIZMO flux limiter"
+
+/**
+ * @brief Limit the flux between two particles.
+ *
+ * @param flux Unlimited flux between the particles.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_flux_limiters_apply(
+ float* flux, struct part* pi, struct part* pj) {
+
+ float flux_limit_factor = 1.;
+ const float timefac = max(pi->force.dt, pj->force.dt);
+ const float areafac = max(pi->geometry.Atot, pj->geometry.Atot);
+ const float totfac = timefac * areafac;
+ if (flux[0] * totfac > pi->conserved.mass) {
+ flux_limit_factor = pi->conserved.mass / (flux[0] * totfac);
+ }
+ if (flux[0] * totfac > pj->conserved.mass) {
+ flux_limit_factor =
+ min(pj->conserved.mass / (flux[0] * totfac), flux_limit_factor);
+ }
+ if (flux[4] * totfac > pi->conserved.energy) {
+ flux_limit_factor =
+ min(pi->conserved.energy / (flux[4] * totfac), flux_limit_factor);
+ }
+ if (flux[4] * totfac > pj->conserved.energy) {
+ flux_limit_factor =
+ min(pj->conserved.energy / (flux[4] * totfac), flux_limit_factor);
+ }
+
+ flux[0] *= flux_limit_factor;
+ flux[1] *= flux_limit_factor;
+ flux[2] *= flux_limit_factor;
+ flux[3] *= flux_limit_factor;
+ flux[4] *= flux_limit_factor;
+}
+
+#else
+
+#define HYDRO_FLUX_LIMITER_IMPLEMENTATION "No flux limiter"
+
+/**
+ * @brief Limit the flux between two particles.
+ *
+ * @param flux Unlimited flux between the particles.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_flux_limiters_apply(
+ float* flux, struct part* pi, struct part* pj) {}
+
+#endif
+
+#endif // SWIFT_HYDRO_FLUX_LIMITERS_H
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index 5ad6d87619a7629a703a8b9c03d089e69ffbdf7d..896128bd45d7964c1f4c8d63564f6fced38db770 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -99,7 +99,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
float xij_j[3];
int k;
float xfac;
- float a_grav_i[3], a_grav_j[3];
/* perform gradient reconstruction in space and time */
/* space */
@@ -141,34 +140,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
pj->primitives.gradients.P[1] * xij_j[1] +
pj->primitives.gradients.P[2] * xij_j[2];
- a_grav_i[0] = pi->gravity.old_a[0];
- a_grav_i[1] = pi->gravity.old_a[1];
- a_grav_i[2] = pi->gravity.old_a[2];
-
- a_grav_i[0] += pi->gravity.grad_a[0][0] * xij_i[0] +
- pi->gravity.grad_a[0][1] * xij_i[1] +
- pi->gravity.grad_a[0][2] * xij_i[2];
- a_grav_i[1] += pi->gravity.grad_a[1][0] * xij_i[0] +
- pi->gravity.grad_a[1][1] * xij_i[1] +
- pi->gravity.grad_a[1][2] * xij_i[2];
- a_grav_i[2] += pi->gravity.grad_a[2][0] * xij_i[0] +
- pi->gravity.grad_a[2][1] * xij_i[1] +
- pi->gravity.grad_a[2][2] * xij_i[2];
-
- a_grav_j[0] = pj->gravity.old_a[0];
- a_grav_j[1] = pj->gravity.old_a[1];
- a_grav_j[2] = pj->gravity.old_a[2];
-
- a_grav_j[0] += pj->gravity.grad_a[0][0] * xij_j[0] +
- pj->gravity.grad_a[0][1] * xij_j[1] +
- pj->gravity.grad_a[0][2] * xij_j[2];
- a_grav_j[1] += pj->gravity.grad_a[1][0] * xij_j[0] +
- pj->gravity.grad_a[1][1] * xij_j[1] +
- pj->gravity.grad_a[1][2] * xij_j[2];
- a_grav_j[2] += pj->gravity.grad_a[2][0] * xij_j[0] +
- pj->gravity.grad_a[2][1] * xij_j[1] +
- pj->gravity.grad_a[2][2] * xij_j[2];
-
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
/* time */
@@ -198,10 +169,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
pi->primitives.gradients.v[1][1] +
pi->primitives.gradients.v[2][2]));
-
- dWi[1] += 0.5 * mindt * a_grav_i[0];
- dWi[2] += 0.5 * mindt * a_grav_i[1];
- dWi[3] += 0.5 * mindt * a_grav_i[2];
}
if (Wj[0] > 0.0f) {
@@ -230,10 +197,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
-
- dWj[1] += 0.5 * mindt * a_grav_j[0];
- dWj[2] += 0.5 * mindt * a_grav_j[1];
- dWj[3] += 0.5 * mindt * a_grav_j[2];
}
Wi[0] += dWi[0];
diff --git a/src/hydro/Gizmo/hydro_gradients_gizmo.h b/src/hydro/Gizmo/hydro_gradients_gizmo.h
index ee3ad6919f81f042ceacc5db8b4e818d63c90266..bc50c10d84cdd6b444887a8bb5fdf7b49a004eb8 100644
--- a/src/hydro/Gizmo/hydro_gradients_gizmo.h
+++ b/src/hydro/Gizmo/hydro_gradients_gizmo.h
@@ -45,18 +45,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
p->primitives.gradients.P[1] = 0.0f;
p->primitives.gradients.P[2] = 0.0f;
- p->gravity.grad_a[0][0] = 0.0f;
- p->gravity.grad_a[0][1] = 0.0f;
- p->gravity.grad_a[0][2] = 0.0f;
-
- p->gravity.grad_a[1][0] = 0.0f;
- p->gravity.grad_a[1][1] = 0.0f;
- p->gravity.grad_a[1][2] = 0.0f;
-
- p->gravity.grad_a[2][0] = 0.0f;
- p->gravity.grad_a[2][1] = 0.0f;
- p->gravity.grad_a[2][2] = 0.0f;
-
hydro_slope_limit_cell_init(p);
}
@@ -157,35 +145,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
(Wi[4] - Wj[4]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->gravity.grad_a[0][0] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[0][1] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[0][2] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->gravity.grad_a[1][0] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[1][1] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[1][2] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->gravity.grad_a[2][0] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[2][1] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[2][2] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
} else {
/* The gradient matrix was not well-behaved, switch to SPH gradients */
@@ -223,27 +182,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
pi->primitives.gradients.P[2] -=
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
- pi->gravity.grad_a[0][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pi->gravity.grad_a[0][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pi->gravity.grad_a[0][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-
- pi->gravity.grad_a[1][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pi->gravity.grad_a[1][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pi->gravity.grad_a[1][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-
- pi->gravity.grad_a[2][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pi->gravity.grad_a[2][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pi->gravity.grad_a[2][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
}
hydro_slope_limit_cell_collect(pi, pj, r);
@@ -306,35 +244,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
(Wi[4] - Wj[4]) * wj *
(Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->gravity.grad_a[0][0] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->gravity.grad_a[0][1] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->gravity.grad_a[0][2] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
- pj->gravity.grad_a[1][0] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->gravity.grad_a[1][1] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->gravity.grad_a[1][2] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
- pj->gravity.grad_a[2][0] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->gravity.grad_a[2][1] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->gravity.grad_a[2][2] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
} else {
/* SPH gradients */
@@ -371,27 +280,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
pj->primitives.gradients.P[2] -=
wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
- pj->gravity.grad_a[0][0] -=
- wj_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pj->gravity.grad_a[0][1] -=
- wj_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pj->gravity.grad_a[0][2] -=
- wj_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-
- pj->gravity.grad_a[1][0] -=
- wj_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pj->gravity.grad_a[1][1] -=
- wj_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pj->gravity.grad_a[1][2] -=
- wj_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-
- pj->gravity.grad_a[2][0] -=
- wj_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pj->gravity.grad_a[2][1] -=
- wj_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pj->gravity.grad_a[2][2] -=
- wj_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
}
hydro_slope_limit_cell_collect(pj, pi, r);
@@ -493,35 +381,6 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
(Wi[4] - Wj[4]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->gravity.grad_a[0][0] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[0][1] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[0][2] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->gravity.grad_a[1][0] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[1][1] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[1][2] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->gravity.grad_a[2][0] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[2][1] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[2][2] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
} else {
/* Gradient matrix is not well-behaved, switch to SPH gradients */
@@ -558,27 +417,6 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
pi->primitives.gradients.P[2] -=
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
- pi->gravity.grad_a[0][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pi->gravity.grad_a[0][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pi->gravity.grad_a[0][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-
- pi->gravity.grad_a[1][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pi->gravity.grad_a[1][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pi->gravity.grad_a[1][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-
- pi->gravity.grad_a[2][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pi->gravity.grad_a[2][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pi->gravity.grad_a[2][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
}
hydro_slope_limit_cell_collect(pi, pj, r);
@@ -618,17 +456,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
p->primitives.gradients.P[1] *= ihdim;
p->primitives.gradients.P[2] *= ihdim;
- p->gravity.grad_a[0][0] *= ihdim;
- p->gravity.grad_a[0][1] *= ihdim;
- p->gravity.grad_a[0][2] *= ihdim;
-
- p->gravity.grad_a[1][0] *= ihdim;
- p->gravity.grad_a[1][1] *= ihdim;
- p->gravity.grad_a[1][2] *= ihdim;
-
- p->gravity.grad_a[2][0] *= ihdim;
- p->gravity.grad_a[2][1] *= ihdim;
- p->gravity.grad_a[2][2] *= ihdim;
} else {
const float ihdimp1 = pow_dimension_plus_one(ih);
@@ -653,18 +480,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
p->primitives.gradients.P[0] *= ihdimp1 * volume;
p->primitives.gradients.P[1] *= ihdimp1 * volume;
p->primitives.gradients.P[2] *= ihdimp1 * volume;
-
- p->gravity.grad_a[0][0] *= ihdimp1 * volume;
- p->gravity.grad_a[0][1] *= ihdimp1 * volume;
- p->gravity.grad_a[0][2] *= ihdimp1 * volume;
-
- p->gravity.grad_a[1][0] *= ihdimp1 * volume;
- p->gravity.grad_a[1][1] *= ihdimp1 * volume;
- p->gravity.grad_a[1][2] *= ihdimp1 * volume;
-
- p->gravity.grad_a[2][0] *= ihdimp1 * volume;
- p->gravity.grad_a[2][1] *= ihdimp1 * volume;
- p->gravity.grad_a[2][2] *= ihdimp1 * volume;
}
hydro_slope_limit_cell(p);
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index c7e86c22da5a8beeece0d7a1fcc9cff560105e6e..0c7c8251b7d1c105dfc0c4b1637724accadaa4ae 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -20,6 +20,7 @@
******************************************************************************/
#include "adiabatic_index.h"
+#include "hydro_flux_limiters.h"
#include "hydro_gradients.h"
#include "riemann.h"
@@ -346,8 +347,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
}
dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
(Wi[3] - Wj[3]) * dx[2];
- if (dvdotdx > 0.) {
- vmax -= dvdotdx / r;
+ dvdotdx = min(dvdotdx, (vi[0] - vj[0]) * dx[0] + (vi[1] - vj[1]) * dx[1] +
+ (vi[2] - vj[2]) * dx[2]);
+ if (dvdotdx < 0.) {
+ /* the magical factor 3 also appears in Gadget2 */
+ vmax -= 3. * dvdotdx / r;
}
pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
if (mode == 1) {
@@ -487,36 +491,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float totflux[5];
riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
- /* Flux limiter */
- float flux_limit_factor = 1.;
- float timefac = max(dti, dtj);
- float areafac = max(pi->geometry.Atot, pj->geometry.Atot);
- if (totflux[0] * areafac * timefac > pi->conserved.mass) {
- flux_limit_factor = pi->conserved.mass / (totflux[0] * areafac * timefac);
- }
- if (totflux[0] * areafac * timefac > pj->conserved.mass) {
- flux_limit_factor =
- min(pj->conserved.mass / (totflux[0] * areafac * timefac),
- flux_limit_factor);
- }
- if (totflux[4] * areafac * timefac > pi->conserved.energy) {
- flux_limit_factor =
- min(pi->conserved.energy / (totflux[4] * areafac * timefac),
- flux_limit_factor);
- }
- if (totflux[4] * areafac * timefac > pj->conserved.energy) {
- flux_limit_factor =
- min(pj->conserved.energy / (totflux[4] * areafac * timefac),
- flux_limit_factor);
- }
- totflux[0] *= flux_limit_factor;
- totflux[1] *= flux_limit_factor;
- totflux[2] *= flux_limit_factor;
- totflux[3] *= flux_limit_factor;
- totflux[4] *= flux_limit_factor;
+ hydro_flux_limiters_apply(totflux, pi, pj);
/* Store mass flux */
- float mflux = mindt * Anorm * totflux[0];
+ float mflux = Anorm * totflux[0];
pi->gravity.mflux[0] += mflux * dx[0];
pi->gravity.mflux[1] += mflux * dx[1];
pi->gravity.mflux[2] += mflux * dx[2];
@@ -554,7 +532,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
if (mode == 1 || pj->force.active == 0) {
/* Store mass flux */
- mflux = mindt * Anorm * totflux[0];
+ mflux = Anorm * totflux[0];
pj->gravity.mflux[0] -= mflux * dx[0];
pj->gravity.mflux[1] -= mflux * dx[1];
pj->gravity.mflux[2] -= mflux * dx[2];
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 3d58be2f47c4e1904aaac5f69d1862f1d453e488..d20f7e2eb1cf50be7690e15a9569d8e9c4605af5 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -18,6 +18,7 @@
******************************************************************************/
#include "adiabatic_index.h"
+#include "hydro_flux_limiters.h"
#include "hydro_gradients.h"
#include "hydro_slope_limiters.h"
#include "io_properties.h"
@@ -127,7 +128,7 @@ float convert_Etot(struct engine* e, 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("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -152,8 +153,6 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
list[9] =
io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
parts, conserved.energy, convert_Etot);
- list[10] = io_make_output_field("GravAcceleration", FLOAT, 3,
- UNIT_CONV_ACCELERATION, parts, gravity.old_a);
}
/**
@@ -171,6 +170,10 @@ void writeSPHflavour(hid_t h_grpsph) {
io_write_attribute_s(h_grpsph, "Piecewise slope limiter model",
HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
+ /* Flux limiter information */
+ io_write_attribute_s(h_grpsph, "Flux limiter model",
+ HYDRO_FLUX_LIMITER_IMPLEMENTATION);
+
/* Riemann solver information */
io_write_attribute_s(h_grpsph, "Riemann solver type",
RIEMANN_SOLVER_IMPLEMENTATION);
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index 6c96004847ae23b46ec3f5182f742e0e84f1118d..47f722c5a2dcce2f3ce603ade3029821d6686067 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -153,10 +153,13 @@ struct part {
} geometry;
- /* Variables used for timestep calculation (currently not used). */
+ /* Variables used for timestep calculation. */
struct {
- /* Maximum fluid velocity among all neighbours. */
+ /* Maximum signal velocity among all the neighbours of the particle. The
+ * signal velocity encodes information about the relative fluid velocities
+ * AND particle velocities of the neighbour and this particle, as well as
+ * the sound speed of both particles. */
float vmax;
} timestepvars;
@@ -201,14 +204,6 @@ struct part {
/* Specific stuff for the gravity-hydro coupling. */
struct {
- /* Previous value of the gravitational acceleration. */
- float old_a[3];
-
- float grad_a[3][3];
-
- /* Previous value of the mass flux vector. */
- float old_mflux[3];
-
/* Current value of the mass flux vector. */
float mflux[3];
diff --git a/src/potential/disc_patch/potential.h b/src/potential/disc_patch/potential.h
index 07f3d7457e7e53debbdda09686ed1166966c1913..ab229d009c692db727e8f2341c3c49813f74f2b8 100644
--- a/src/potential/disc_patch/potential.h
+++ b/src/potential/disc_patch/potential.h
@@ -56,20 +56,20 @@ struct external_potential {
/*! Inverse of disc scale-height (1/b) */
float scale_height_inv;
- /*! Position of the disc along the z-axis */
- float z_disc;
+ /*! Position of the disc along the x-axis */
+ float x_disc;
/*! Position above which the accelerations get truncated */
- float z_trunc;
+ float x_trunc;
/*! Position above which the accelerations are zero */
- float z_max;
+ float x_max;
/*! The truncated transition regime */
- float z_trans;
+ float x_trans;
/*! Inverse of the truncated transition regime */
- float z_trans_inv;
+ float x_trans_inv;
/*! Dynamical time of the system */
float dynamical_time;
@@ -115,36 +115,36 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const float norm = potential->norm;
/* absolute value of height above disc */
- const float dz = fabsf(g->x[2] - potential->z_disc);
+ const float dx = fabsf(g->x[0] - potential->x_disc);
/* vertical acceleration */
- const float z_accel = norm * tanhf(dz * b_inv);
+ const float x_accel = norm * tanhf(dx * b_inv);
float dt = dt_dyn;
/* demand that dt * velocity < fraction of scale height of disc */
- if (g->v_full[2] != 0.f) {
+ if (g->v_full[0] != 0.f) {
- const float dt1 = b / fabsf(g->v_full[2]);
+ const float dt1 = b / fabsf(g->v_full[0]);
dt = min(dt1, dt);
}
/* demand that dt^2 * acceleration < fraction of scale height of disc */
- if (z_accel != 0.f) {
+ if (x_accel != 0.f) {
- const float dt2 = b / fabsf(z_accel);
+ const float dt2 = b / fabsf(x_accel);
if (dt2 < dt * dt) dt = sqrtf(dt2);
}
/* demand that dt^3 * jerk < fraction of scale height of disc */
- if (g->v_full[2] != 0.f) {
+ if (g->v_full[0] != 0.f) {
- const float cosh_dz_inv = 1.f / coshf(dz * b_inv);
- const float cosh_dz_inv2 = cosh_dz_inv * cosh_dz_inv;
- const float dz_accel_over_dt =
- norm * cosh_dz_inv2 * b_inv * fabsf(g->v_full[2]);
+ const float cosh_dx_inv = 1.f / coshf(dx * b_inv);
+ const float cosh_dx_inv2 = cosh_dx_inv * cosh_dx_inv;
+ const float dx_accel_over_dt =
+ norm * cosh_dx_inv2 * b_inv * fabsf(g->v_full[0]);
- const float dt3 = b / fabsf(dz_accel_over_dt);
+ const float dt3 = b / fabsf(dx_accel_over_dt);
if (dt3 < dt * dt * dt) dt = cbrtf(dt3);
}
@@ -152,13 +152,13 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
}
/**
- * @brief Computes the gravitational acceleration along z due to a hydrostatic
+ * @brief Computes the gravitational acceleration along x due to a hydrostatic
* disc
*
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equation 17.
- * We truncate the accelerations beyond z_trunc using a 1-cos(z) function
- * that smoothly brings the accelerations to 0 at z_max.
+ * We truncate the accelerations beyond x_trunc using a 1-cos(x) function
+ * that smoothly brings the accelerations to 0 at x_max.
*
* @param time The current time in internal units.
* @param potential The properties of the potential.
@@ -169,40 +169,40 @@ __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 dz = g->x[2] - potential->z_disc;
- const float abs_dz = fabsf(dz);
+ const float dx = g->x[0] - potential->x_disc;
+ const float abs_dx = fabsf(dx);
const float t_growth = potential->growth_time;
const float t_growth_inv = potential->growth_time_inv;
const float b_inv = potential->scale_height_inv;
- const float z_trunc = potential->z_trunc;
- const float z_max = potential->z_max;
- const float z_trans_inv = potential->z_trans_inv;
+ const float x_trunc = potential->x_trunc;
+ const float x_max = potential->x_max;
+ const float x_trans_inv = potential->x_trans_inv;
const float norm_over_G = potential->norm_over_G;
/* Are we still growing the disc ? */
const float reduction_factor = time < t_growth ? time * t_growth_inv : 1.f;
/* Truncated or not ? */
- float a_z;
- if (abs_dz < z_trunc) {
+ float a_x;
+ if (abs_dx < x_trunc) {
- /* Acc. 2 pi sigma tanh(z/b) */
- a_z = reduction_factor * norm_over_G * tanhf(abs_dz * b_inv);
- } else if (abs_dz < z_max) {
+ /* Acc. 2 pi sigma tanh(x/b) */
+ a_x = reduction_factor * norm_over_G * tanhf(abs_dx * b_inv);
+ } else if (abs_dx < x_max) {
- /* Acc. 2 pi sigma tanh(z/b) [1/2 + 1/2cos((z-zmax)/(pi z_trans))] */
- a_z =
- reduction_factor * norm_over_G * tanhf(abs_dz * b_inv) *
- (0.5f + 0.5f * cosf((float)(M_PI) * (abs_dz - z_trunc) * z_trans_inv));
+ /* Acc. 2 pi sigma tanh(x/b) [1/2 + 1/2cos((x-xmax)/(pi x_trans))] */
+ a_x =
+ reduction_factor * norm_over_G * tanhf(abs_dx * b_inv) *
+ (0.5f + 0.5f * cosf((float)(M_PI) * (abs_dx - x_trunc) * x_trans_inv));
} else {
/* Acc. 0 */
- a_z = 0.f;
+ a_x = 0.f;
}
/* Get the correct sign. Recall G is multipiled in later on */
- if (dz > 0) g->a_grav[2] -= a_z;
- if (dz < 0) g->a_grav[2] += a_z;
+ if (dx > 0) g->a_grav[0] -= a_x;
+ if (dx < 0) g->a_grav[0] += a_x;
}
/**
@@ -211,8 +211,8 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
*
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equation 22.
- * We truncate the accelerations beyond z_trunc using a 1-cos(z) function
- * that smoothly brings the accelerations to 0 at z_max.
+ * We truncate the accelerations beyond x_trunc using a 1-cos(x) function
+ * that smoothly brings the accelerations to 0 at x_max.
*
* @param time The current time.
* @param potential The #external_potential used in the run.
@@ -224,28 +224,28 @@ external_gravity_get_potential_energy(
double time, const struct external_potential* potential,
const struct phys_const* const phys_const, const struct gpart* gp) {
- const float dz = gp->x[2] - potential->z_disc;
- const float abs_dz = fabsf(dz);
+ const float dx = gp->x[0] - potential->x_disc;
+ const float abs_dx = fabsf(dx);
const float t_growth = potential->growth_time;
const float t_growth_inv = potential->growth_time_inv;
const float b = potential->scale_height;
const float b_inv = potential->scale_height_inv;
const float norm = potential->norm;
- const float z_trunc = potential->z_trunc;
- const float z_max = potential->z_max;
+ const float x_trunc = potential->x_trunc;
+ const float x_max = potential->x_max;
/* Are we still growing the disc ? */
const float reduction_factor = time < t_growth ? time * t_growth_inv : 1.f;
/* Truncated or not ? */
float pot;
- if (abs_dz < z_trunc) {
+ if (abs_dx < x_trunc) {
- /* Potential (2 pi G sigma b ln(cosh(z/b)) */
- pot = b * logf(coshf(dz * b_inv));
- } else if (abs_dz < z_max) {
+ /* Potential (2 pi G sigma b ln(cosh(x/b)) */
+ pot = b * logf(coshf(dx * b_inv));
+ } else if (abs_dx < x_max) {
- /* Potential. At z>>b, phi(z) = norm * z / b */
+ /* Potential. At x>>b, phi(x) = norm * x / b */
pot = 0.f;
} else {
@@ -277,14 +277,14 @@ static INLINE void potential_init_backend(
parameter_file, "DiscPatchPotential:surface_density");
potential->scale_height = parser_get_param_double(
parameter_file, "DiscPatchPotential:scale_height");
- potential->z_disc =
- parser_get_param_double(parameter_file, "DiscPatchPotential:z_disc");
- potential->z_trunc = parser_get_opt_param_double(
- parameter_file, "DiscPatchPotential:z_trunc", FLT_MAX);
- potential->z_max = parser_get_opt_param_double(
- parameter_file, "DiscPatchPotential:z_max", FLT_MAX);
- potential->z_disc =
- parser_get_param_double(parameter_file, "DiscPatchPotential:z_disc");
+ potential->x_disc =
+ parser_get_param_double(parameter_file, "DiscPatchPotential:x_disc");
+ potential->x_trunc = parser_get_opt_param_double(
+ parameter_file, "DiscPatchPotential:x_trunc", FLT_MAX);
+ potential->x_max = parser_get_opt_param_double(
+ parameter_file, "DiscPatchPotential:x_max", FLT_MAX);
+ potential->x_disc =
+ parser_get_param_double(parameter_file, "DiscPatchPotential:x_disc");
potential->timestep_mult = parser_get_param_double(
parameter_file, "DiscPatchPotential:timestep_mult");
potential->growth_time = parser_get_opt_param_double(
@@ -299,22 +299,22 @@ static INLINE void potential_init_backend(
potential->growth_time *= potential->dynamical_time;
/* Some cross-checks */
- if (potential->z_trunc > potential->z_max)
- error("Potential truncation z larger than maximal z");
- if (potential->z_trunc < potential->scale_height)
- error("Potential truncation z smaller than scale height");
+ if (potential->x_trunc > potential->x_max)
+ error("Potential truncation x larger than maximal z");
+ if (potential->x_trunc < potential->scale_height)
+ error("Potential truncation x smaller than scale height");
/* Compute derived quantities */
potential->scale_height_inv = 1. / potential->scale_height;
potential->norm =
2. * M_PI * phys_const->const_newton_G * potential->surface_density;
potential->norm_over_G = 2 * M_PI * potential->surface_density;
- potential->z_trans = potential->z_max - potential->z_trunc;
+ potential->x_trans = potential->x_max - potential->x_trunc;
- if (potential->z_trans != 0.f)
- potential->z_trans_inv = 1. / potential->z_trans;
+ if (potential->x_trans != 0.f)
+ potential->x_trans_inv = 1. / potential->x_trans;
else
- potential->z_trans_inv = FLT_MAX;
+ potential->x_trans_inv = FLT_MAX;
if (potential->growth_time != 0.)
potential->growth_time_inv = 1. / potential->growth_time;
@@ -331,14 +331,14 @@ static INLINE void potential_print_backend(
const struct external_potential* potential) {
message(
- "External potential is 'Disk-patch' with Sigma=%f, z_disc=%f, b=%f and "
+ "External potential is 'Disk-patch' with Sigma=%f, x_disc=%f, b=%f and "
"dt_mult=%f.",
- potential->surface_density, potential->z_disc, potential->scale_height,
+ potential->surface_density, potential->x_disc, potential->scale_height,
potential->timestep_mult);
- if (potential->z_max < FLT_MAX)
- message("Potential will be truncated at z_trunc=%f and zeroed at z_max=%f",
- potential->z_trunc, potential->z_max);
+ if (potential->x_max < FLT_MAX)
+ message("Potential will be truncated at x_trunc=%f and zeroed at x_max=%f",
+ potential->x_trunc, potential->x_max);
if (potential->growth_time > 0.)
message("Disc will grow for %f [time_units]. (%f dynamical time)",
diff --git a/src/queue.h b/src/queue.h
index c85cf0cabe30a03d163e2564fdc216c19495761a..951a3e5a056d7ad0c3935f98341a0d93c805e3ad 100644
--- a/src/queue.h
+++ b/src/queue.h
@@ -29,7 +29,7 @@
#define queue_sizeinit 100
#define queue_sizegrow 2
#define queue_search_window 8
-#define queue_incoming_size 10240
+#define queue_incoming_size 1024
#define queue_struct_align 64
/* Counters. */
diff --git a/src/runner.c b/src/runner.c
index 04d230209a27d9c16d5ffbc0fffc9d5519e025d4..c32ba4bf76ae2cd06c28ed1d41b4ebf16abe0885 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1544,11 +1544,6 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
time_bin_min = min(time_bin_min, parts[k].time_bin);
time_bin_max = max(time_bin_max, parts[k].time_bin);
h_max = max(h_max, parts[k].h);
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (parts[k].ti_drift != ti_current)
- error("Received un-drifted particle !");
-#endif
}
/* Convert into a time */
diff --git a/src/scheduler.c b/src/scheduler.c
index 2e645951256d415b072cc964003d4042a842923d..161be418046636c473362f061d4afa1c64ace04a 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -759,7 +759,7 @@ void scheduler_splittasks(struct scheduler *s) {
/* Call the mapper on each current task. */
threadpool_map(s->threadpool, scheduler_splittasks_mapper, s->tasks,
- s->nr_tasks, sizeof(struct task), 0, s);
+ s->nr_tasks, sizeof(struct task), 1000, s);
}
/**
@@ -1174,7 +1174,7 @@ void scheduler_start(struct scheduler *s) {
/* Re-wait the tasks. */
if (s->active_count > 1000) {
threadpool_map(s->threadpool, scheduler_rewait_mapper, s->tid_active,
- s->active_count, sizeof(int), 0, s);
+ s->active_count, sizeof(int), 1000, s);
} else {
scheduler_rewait_mapper(s->tid_active, s->active_count, s);
}
@@ -1250,7 +1250,7 @@ void scheduler_start(struct scheduler *s) {
/* Loop over the tasks and enqueue whoever is ready. */
if (s->active_count > 1000) {
threadpool_map(s->threadpool, scheduler_enqueue_mapper, s->tid_active,
- s->active_count, sizeof(int), 0, s);
+ s->active_count, sizeof(int), 1000, s);
} else {
scheduler_enqueue_mapper(s->tid_active, s->active_count, s);
}
@@ -1363,11 +1363,6 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
} else if (t->subtype == task_subtype_xv ||
t->subtype == task_subtype_rho ||
t->subtype == task_subtype_gradient) {
-#ifdef SWIFT_DEBUG_CHECKS
- for (int k = 0; k < t->ci->count; k++)
- if (t->ci->parts[k].ti_drift != s->space->e->ti_current)
- error("Sending un-drifted particle !");
-#endif
err = MPI_Isend(t->ci->parts, t->ci->count, part_mpi_type,
t->cj->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
// message( "sending %i parts with tag=%i from %i to %i." ,
diff --git a/src/space.c b/src/space.c
index 8ad571be3800fbeebd280dfd09b9ee29158bfdf8..23902a37501c7b13992f5423a3f002d526ba2c27 100644
--- a/src/space.c
+++ b/src/space.c
@@ -378,7 +378,7 @@ void space_regrid(struct space *s, int verbose) {
/* Free the old cells, if they were allocated. */
if (s->cells_top != NULL) {
threadpool_map(&s->e->threadpool, space_rebuild_recycle_mapper,
- s->cells_top, s->nr_cells, sizeof(struct cell), 0, s);
+ s->cells_top, s->nr_cells, sizeof(struct cell), 100, s);
free(s->cells_top);
free(s->multipoles_top);
s->maxdepth = 0;
@@ -491,7 +491,7 @@ void space_regrid(struct space *s, int verbose) {
/* Free the old cells, if they were allocated. */
threadpool_map(&s->e->threadpool, space_rebuild_recycle_mapper,
- s->cells_top, s->nr_cells, sizeof(struct cell), 0, s);
+ s->cells_top, s->nr_cells, sizeof(struct cell), 100, s);
s->maxdepth = 0;
}
@@ -970,7 +970,7 @@ void space_split(struct space *s, struct cell *cells, int nr_cells,
const ticks tic = getticks();
threadpool_map(&s->e->threadpool, space_split_mapper, cells, nr_cells,
- sizeof(struct cell), 0, s);
+ sizeof(struct cell), 1, s);
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -1004,7 +1004,7 @@ void space_sanitize(struct space *s) {
if (s->e->nodeID == 0) message("Cleaning up unreasonable values of h");
threadpool_map(&s->e->threadpool, space_sanitize_mapper, s->cells_top,
- s->nr_cells, sizeof(struct cell), 0, NULL);
+ s->nr_cells, sizeof(struct cell), 1, NULL);
}
/**
@@ -1187,7 +1187,7 @@ void space_parts_get_cell_index(struct space *s, int *ind, struct cell *cells,
data.ind = ind;
threadpool_map(&s->e->threadpool, space_parts_get_cell_index_mapper, s->parts,
- s->nr_parts, sizeof(struct part), 0, &data);
+ s->nr_parts, sizeof(struct part), 1000, &data);
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -1214,7 +1214,7 @@ void space_gparts_get_cell_index(struct space *s, int *gind, struct cell *cells,
data.ind = gind;
threadpool_map(&s->e->threadpool, space_gparts_get_cell_index_mapper,
- s->gparts, s->nr_gparts, sizeof(struct gpart), 0, &data);
+ s->gparts, s->nr_gparts, sizeof(struct gpart), 1000, &data);
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -1241,7 +1241,7 @@ void space_sparts_get_cell_index(struct space *s, int *sind, struct cell *cells,
data.ind = sind;
threadpool_map(&s->e->threadpool, space_sparts_get_cell_index_mapper,
- s->sparts, s->nr_sparts, sizeof(struct spart), 0, &data);
+ s->sparts, s->nr_sparts, sizeof(struct spart), 1000, &data);
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -2501,7 +2501,7 @@ void space_synchronize_particle_positions(struct space *s) {
(s->nr_gparts > 0 && s->nr_sparts > 0))
threadpool_map(&s->e->threadpool,
space_synchronize_particle_positions_mapper, s->gparts,
- s->nr_gparts, sizeof(struct gpart), 0, (void *)s);
+ s->nr_gparts, sizeof(struct gpart), 1000, (void *)s);
}
/**
diff --git a/src/statistics.c b/src/statistics.c
index 5a3f1ff4f9a2232a14817e7e55fd2cff5bdcd80e..57d60bcb1b247c9616c859b7ac8a475acdcd878f 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -271,12 +271,12 @@ void stats_collect(const struct space *s, struct statistics *stats) {
/* Run parallel collection of statistics for parts */
if (s->nr_parts > 0)
threadpool_map(&s->e->threadpool, stats_collect_part_mapper, s->parts,
- s->nr_parts, sizeof(struct part), 0, &extra_data);
+ s->nr_parts, sizeof(struct part), 10000, &extra_data);
/* Run parallel collection of statistics for gparts */
if (s->nr_gparts > 0)
threadpool_map(&s->e->threadpool, stats_collect_gpart_mapper, s->gparts,
- s->nr_gparts, sizeof(struct gpart), 0, &extra_data);
+ s->nr_gparts, sizeof(struct gpart), 10000, &extra_data);
}
/**
diff --git a/src/swift.h b/src/swift.h
index 1d1a7c7d04b3662c524504c292aa7d9eee2c3d09..20397eb24df478cba65a0e35d686b402f1d8ee70 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -57,7 +57,6 @@
#include "sourceterms.h"
#include "space.h"
#include "task.h"
-#include "threadpool.h"
#include "timeline.h"
#include "timers.h"
#include "tools.h"
diff --git a/src/threadpool.c b/src/threadpool.c
index d1b3d0c4269a17022e30c5f5f8da5c8654b79034..3ad02cb369bfcec22689a5aacf97816213eb49f0 100644
--- a/src/threadpool.c
+++ b/src/threadpool.c
@@ -41,7 +41,7 @@
#ifdef SWIFT_DEBUG_THREADPOOL
/**
- * @breif Store a log entry of the given chunk.
+ * @brief Store a log entry of the given chunk.
*/
void threadpool_log(struct threadpool *tp, int tid, size_t chunk_size,
ticks tic, ticks toc) {
diff --git a/tests/test27cells.c b/tests/test27cells.c
index f34d82c587a5a4abb91670b03f02393fa48e57ca..458129edd3b92d78a45102c118fb3e6f44ff604d 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -204,6 +204,10 @@ void zero_particle_fields(struct cell *c) {
void end_calculation(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
hydro_end_density(&c->parts[pid]);
+
+ /* Recover the common "Neighbour number" definition */
+ c->parts[pid].density.wcount *= pow_dimension(c->parts[pid].h);
+ c->parts[pid].density.wcount *= kernel_norm;
}
}
diff --git a/tests/testThreadpool.c b/tests/testThreadpool.c
index 2a9e98c5ca71ba62c5f6a266f4a36f658f61b063..aa65d533a29afbe4e7e8384fb887281822a31e58 100644
--- a/tests/testThreadpool.c
+++ b/tests/testThreadpool.c
@@ -23,7 +23,6 @@
#include <unistd.h>
// Local includes.
-#include "../config.h"
#include "../src/atomic.h"
#include "../src/threadpool.h"
@@ -79,11 +78,6 @@ int main(int argc, char *argv[]) {
threadpool_map(&tp, map_function_first, data, N, sizeof(int), 2, NULL);
}
-/* If logging was enabled, dump the log. */
-#ifdef SWIFT_DEBUG_THREADPOOL
- threadpool_dump_log(&tp, "threadpool_log.txt", 1);
-#endif
-
/* Be clean */
threadpool_clean(&tp);