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/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/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/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)",